From 1a4a33b25a37f043760ddd7f3a4b95f7214fe33a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Jun 2016 15:43:01 -0400 Subject: [PATCH 1/5] Added g_Earth to ice_shelf_CS Made g_Earth into an element of the ice_shelf code, to accomodate the removal of g_Earth from the initialization of the ocean_grid_type. Also replaced G%g_Earth with CS%g_Earth throughout MOM_ice_shelf.F90. All answers are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index b4533b3de5..4cfc98198f 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -241,6 +241,7 @@ module MOM_ice_shelf real :: ustar_bg ! A minimum value for ustar under ice shelves, in m s-1. real :: cdrag ! drag coefficient under ice shelves , non-dimensional. + real :: g_Earth ! The gravitational acceleration in m s-2. real :: Cp ! The heat capacity of sea water, in J kg-1 K-1. real :: Rho0 ! A reference ocean density in kg/m3. real :: Cp_ice ! The heat capacity of fresh ice, in J kg-1 K-1. @@ -499,7 +500,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) do j=js,je ! Find the pressure at the ice-ocean interface, averaged only over the ! part of the cell covered by ice shelf. - do i=is,ie ; p_int(i) = G%g_Earth * CS%mass_shelf(i,j) ; enddo + do i=is,ie ; p_int(i) = CS%g_Earth * CS%mass_shelf(i,j) ; enddo ! Calculate insitu densities and expansion coefficients call calculate_density(state%sst(:,j),state%sss(:,j), p_int, & @@ -548,8 +549,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) Sbdry = state%sss(i,j) ; Sb_max_set = .false. ; Sb_min_set = .false. ! Determine the mixed layer buoyancy flux, wB_flux. - dB_dS = (G%g_Earth / Rhoml(i)) * dR0_dS(i) - dB_dT = (G%g_Earth / Rhoml(i)) * dR0_dT(i) + dB_dS = (CS%g_Earth / Rhoml(i)) * dR0_dS(i) + dB_dT = (CS%g_Earth / Rhoml(i)) * dR0_dT(i) ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) do it1 = 1,20 @@ -891,10 +892,10 @@ subroutine add_shelf_flux(G, CS, state, fluxes) if (associated(fluxes%sens)) fluxes%sens(i,j) = -frac_area*CS%t_flux(i,j)*CS%flux_factor if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = frac_area * CS%salt_flux(i,j)*CS%flux_factor - if (associated(fluxes%p_surf)) fluxes%p_surf(i,j) = frac_area * G%g_Earth * CS%mass_shelf(i,j) + if (associated(fluxes%p_surf)) fluxes%p_surf(i,j) = frac_area * CS%g_Earth * CS%mass_shelf(i,j) ! Same for IOB%p if (associated(fluxes%p_surf_full) ) fluxes%p_surf_full(i,j) = & - frac_area * G%g_Earth * CS%mass_shelf(i,j) + frac_area * CS%g_Earth * CS%mass_shelf(i,j) endif enddo ; enddo @@ -1016,10 +1017,10 @@ end subroutine add_shelf_flux ! ! fluxes%salt_flux(i,j) = fluxes%salt_flux(i,j) + frac_area * CS%salt_flux(i,j) ! ! ! Same for IOB%salt_flux. ! fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + & -! frac_area * G%g_Earth * CS%mass_shelf(i,j) +! frac_area * CS%g_Earth * CS%mass_shelf(i,j) ! ! Same for IOB%p ! if (associated(fluxes%p_surf_full)) fluxes%p_surf_full(i,j) = & -! fluxes%p_surf_full(i,j) + frac_area * G%g_Earth * CS%mass_shelf(i,j) +! fluxes%p_surf_full(i,j) + frac_area * CS%g_Earth * CS%mass_shelf(i,j) ! endif ! enddo ; enddo @@ -1165,6 +1166,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti "exchange velocity at the ice-ocean interface.", & units="m s-1", fail_if_missing=.true.) + call get_param(param_file, mod, "G_EARTH", CS%g_Earth, & + "The gravitational acceleration of the Earth.", & + units="m s-2", default = 9.80) call get_param(param_file, mod, "C_P", CS%Cp, & "The heat capacity of sea water.", units="J kg-1 K-1", & fail_if_missing=.true.) @@ -1569,10 +1573,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti !if (.not. solo_ice_sheet) then if (G%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = CS%area_shelf_h(i,j) / G%areaT(i,j) if (associated(fluxes%p_surf)) & - fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + fluxes%frac_shelf_h(i,j) * (G%g_Earth * CS%mass_shelf(i,j)) + fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + & + fluxes%frac_shelf_h(i,j) * (CS%g_Earth * CS%mass_shelf(i,j)) if (associated(fluxes%p_surf_full)) & fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + & - fluxes%frac_shelf_h(i,j) * (G%g_Earth * CS%mass_shelf(i,j)) + fluxes%frac_shelf_h(i,j) * (CS%g_Earth * CS%mass_shelf(i,j)) !endif enddo ; enddo @@ -1806,9 +1811,9 @@ subroutine initialize_shelf_mass(G, param_file, CS, new_sim) CS%area_shelf_h(i,j) = 0.0 enddo ; enddo - case ("USER") - call USER_initialize_shelf_mass(CS%mass_shelf, CS%area_shelf_h, CS%h_shelf, CS%hmask, G, CS%user_CS, param_file, new_sim_2) + call USER_initialize_shelf_mass(CS%mass_shelf, CS%area_shelf_h, & + CS%h_shelf, CS%hmask, G, CS%user_CS, param_file, new_sim_2) case default ; call MOM_error(FATAL,"initialize_ice_shelf: "// & "Unrecognized ice shelf setup "//trim(config)) @@ -1820,7 +1825,8 @@ subroutine update_shelf_mass(CS, Time) type(ice_shelf_CS), pointer :: CS type(time_type), intent(in) :: Time - call USER_update_shelf_mass(CS%mass_shelf, CS%area_shelf_h, CS%h_shelf, CS%hmask, CS%grid, CS%user_CS, Time, .true.) + call USER_update_shelf_mass(CS%mass_shelf, CS%area_shelf_h, CS%h_shelf, & + CS%hmask, CS%grid, CS%user_CS, Time, .true.) end subroutine update_shelf_mass From f926dae091bfe7542ea1e8c616df96d385811c86 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Jun 2016 15:49:45 -0400 Subject: [PATCH 2/5] +Initialize ocean_grid from dyn_horgrid Completed the initialization of fixed elements of MOM6 using a shareable dynamic horizontal grid. Changes include: - Moved BATHYMETRY_AT_VEL, GLOBAL_INDEXING and FIRST_DIRECTION into MOM.F90 - Eliminated g_Earth from the dyn_horgrid. - Added optional arguments to MOM_grid_init and create_dyn_horgrid to set the parameters read in by initialize_MOM. - Allocation and deallocation of Dopen_u, etc., in dyn_horgrid is conditional upon BATHYMETRY_AT_VEL. All solutions are bitwise identical, although several interfaces and all MOM_parameter_doc files have changed. --- src/core/MOM.F90 | 110 +++++++++++++++++++------------ src/core/MOM_dyn_horgrid.F90 | 29 +++++--- src/core/MOM_grid.F90 | 56 +++++++--------- src/core/MOM_transcribe_grid.F90 | 2 - src/core/MOM_verticalGrid.F90 | 2 +- 5 files changed, 115 insertions(+), 84 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 87df795436..ba9e2941ec 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -85,8 +85,9 @@ module MOM use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_EOS, only : EOS_init use MOM_error_checking, only : check_redundant -use MOM_grid, only : MOM_grid_init, ocean_grid_type, MOM_grid_end -use MOM_hor_index, only : hor_index_type +use MOM_grid, only : ocean_grid_type, set_first_direction +use MOM_grid, only : MOM_grid_init, MOM_grid_end +use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_hor_visc, only : horizontal_viscosity, hor_visc_init use MOM_interface_heights, only : find_eta use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init @@ -1365,10 +1366,11 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) !! model is not being started from a restart file ! local - type(ocean_grid_type), pointer :: G ! pointer to a structure with metrics and related + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to a structure with metrics and related + type(hor_index_type) :: HI ! A hor_index_type for array extents type(verticalGrid_type), pointer :: GV => NULL() type(dyn_horgrid_type), pointer :: dG => NULL() - type(diag_ctrl), pointer :: diag + type(diag_ctrl), pointer :: diag character(len=4), parameter :: vers_num = 'v2.0' @@ -1384,7 +1386,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) real, allocatable, dimension(:,:) :: eta ! free surface height (m) or bottom press (Pa) type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() - integer :: nkml, nkbl, verbosity, write_geom real :: default_val ! default value for a parameter logical :: write_geom_files ! If true, write out the grid geometry files. logical :: new_sim @@ -1394,6 +1395,15 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) logical :: save_IC ! If true, save the initial conditions. logical :: do_unit_tests ! If true, call unit tests. logical :: test_grid_copy = .false. + logical :: global_indexing ! If true use global horizontal index values instead + ! of having the data domain on each processor start at 1. + logical :: bathy_at_vel ! If true, also define bathymetric fields at the + ! the velocity points. + integer :: first_direction ! An integer that indicates which direction is to be + ! updated first in directionally split parts of the + ! calculation. This can be altered during the course + ! of the run via calls to set_first_direction. + integer :: nkml, nkbl, verbosity, write_geom type(time_type) :: Start_time type(ocean_internal_state) :: MOM_internal_state @@ -1486,6 +1496,11 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "If true, do thickness diffusion before dynamics.\n"//& "This is only used if THICKNESSDIFFUSE is true.", & default=.false.) + call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, & + "If true, there are separate values for the basin depths \n"//& + "at velocity points. Otherwise the effects of topography \n"//& + "are entirely determined from thickness points.", & + default=.false.) call get_param(param_file, "MOM", "DEBUG", CS%debug, & "If true, write out verbose debugging data.", default=.false.) @@ -1581,6 +1596,21 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) default=2) endif + call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, & + "If true, use a global lateral indexing convention, so \n"//& + "that corresponding points on different processors have \n"//& + "the same index. This does not work with static memory.", & + default=.false., layoutParam=.true.) +#ifdef STATIC_MEMORY_ + if (global_indexing) call MOM_error(FATAL, "initialize_MOM: "//& + "GLOBAL_INDEXING can not be true with STATIC_MEMORY.") +#endif + call get_param(param_file, "MOM", "FIRST_DIRECTION", first_direction, & + "An integer that indicates which direction goes first \n"//& + "in parts of the code that use directionally split \n"//& + "updates, with even numbers (or 0) used for x- first \n"//& + "and odd numbers used for y-first.", default=0) + call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", & CS%check_bad_surface_vals, & "If true, check the surface state for ridiculous values.", & @@ -1631,6 +1661,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.") if (CS%adiabatic .and. CS%bulkmixedlayer) call MOM_error(FATAL, & "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.") + if (CS%bulkmixedlayer .and. .not.use_EOS) call MOM_error(FATAL, & + "initialize_MOM: A bulk mixed layer can only be used with T & S as "//& + "state variables. Add USE_EOS = True to MOM_input.") call callTree_waypoint("MOM parameters read (initialize_MOM)") @@ -1651,20 +1684,18 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("domains initialized (initialize_MOM)") call MOM_checksums_init(param_file) - call diag_mediator_infrastructure_init() call MOM_io_init(param_file) - call MOM_grid_init(G, param_file) - call create_dyn_horgrid(dG, G%HI) - dG%first_direction = G%first_direction - dG%bathymetry_at_vel = G%bathymetry_at_vel + call hor_index_init(G%Domain, HI, param_file, & + local_indexing=.not.global_indexing) + + call create_dyn_horgrid(dG, HI, bathymetry_at_vel=bathy_at_vel) call clone_MOM_domain(G%Domain, dG%Domain) call verticalGridInit( param_file, CS%GV ) GV => CS%GV - dG%g_Earth = GV%g_Earth ; G%g_Earth = GV%g_Earth - +! dG%g_Earth = GV%g_Earth ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. if (CS%debug .or. dG%symmetric) & @@ -1677,10 +1708,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call tracer_registry_init(param_file, CS%tracer_Reg) - ! Copy a common variable from the vertical grid to the horizontal grid. - ! Consider removing this later? -! G%ke = GV%ke - is = dG%isc ; ie = dG%iec ; js = dG%jsc ; je = dG%jec ; nz = GV%ke isd = dG%isd ; ied = dG%ied ; jsd = dG%jsd ; jed = dG%jed IsdB = dG%IsdB ; IedB = dG%IedB ; JsdB = dG%JsdB ; JedB = dG%JedB @@ -1701,8 +1728,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) CS%vd_S = var_desc(name="S",units="PPT",longname="Salinity",& cmor_field_name="so",cmor_units="ppt", & conversion=0.001) - call register_tracer(CS%tv%T, CS%vd_T, param_file, G%HI, GV, CS%tracer_Reg, CS%vd_T) - call register_tracer(CS%tv%S, CS%vd_S, param_file, G%HI, GV, CS%tracer_Reg, CS%vd_S) + call register_tracer(CS%tv%T, CS%vd_T, param_file, dG%HI, GV, CS%tracer_Reg, CS%vd_T) + call register_tracer(CS%tv%S, CS%vd_S, param_file, dG%HI, GV, CS%tracer_Reg, CS%vd_S) endif if (CS%use_frazil) then allocate(CS%tv%frazil(isd:ied,jsd:jed)) ; CS%tv%frazil(:,:) = 0.0 @@ -1712,11 +1739,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) endif if (CS%bulkmixedlayer) then - if (.not.use_EOS) call MOM_error(FATAL, & - "initialize_MOM: A bulk mixed layer can only be used with T & S as "//& - "state variables. Add #define USE_EOS.") - GV%nkml = nkml - GV%nk_rho_varies = nkml + nkbl + GV%nkml = nkml ; GV%nk_rho_varies = nkml + nkbl allocate(CS%tv%Hml(isd:ied,jsd:jed)) ; CS%tv%Hml(:,:) = 0.0 else GV%nkml = 0 ; GV%nk_rho_varies = 0 @@ -1769,30 +1792,30 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call set_restart_fields(GV, param_file, CS) if (CS%split) then if (CS%legacy_split) then - call register_restarts_dyn_legacy_split(G%HI, GV, param_file, & + call register_restarts_dyn_legacy_split(dG%HI, GV, param_file, & CS%dyn_legacy_split_CSp, CS%restart_CSp, CS%uh, CS%vh) else - call register_restarts_dyn_split_RK2(G%HI, GV, param_file, & + call register_restarts_dyn_split_RK2(dG%HI, GV, param_file, & CS%dyn_split_RK2_CSp, CS%restart_CSp, CS%uh, CS%vh) endif else if (CS%use_RK2) then - call register_restarts_dyn_unsplit_RK2(G%HI, GV, param_file, & + call register_restarts_dyn_unsplit_RK2(dG%HI, GV, param_file, & CS%dyn_unsplit_RK2_CSp, CS%restart_CSp) else - call register_restarts_dyn_unsplit(G%HI, GV, param_file, & + call register_restarts_dyn_unsplit(dG%HI, GV, param_file, & CS%dyn_unsplit_CSp, CS%restart_CSp) endif endif ! This subroutine calls user-specified tracer registration routines. ! Additional calls can be added to MOM_tracer_flow_control.F90. - call call_tracer_register(G%HI, GV, param_file, CS%tracer_flow_CSp, & + call call_tracer_register(dG%HI, GV, param_file, CS%tracer_flow_CSp, & CS%tracer_Reg, CS%restart_CSp) - call MEKE_alloc_register_restart(G%HI, param_file, CS%MEKE, CS%restart_CSp) - call set_visc_register_restarts(G%HI, GV, param_file, CS%visc, CS%restart_CSp) - call mixedlayer_restrat_register_restarts(G%HI, param_file, CS%mixedlayer_restrat_CSp, CS%restart_CSp) + call MEKE_alloc_register_restart(dG%HI, param_file, CS%MEKE, CS%restart_CSp) + call set_visc_register_restarts(dG%HI, GV, param_file, CS%visc, CS%restart_CSp) + call mixedlayer_restrat_register_restarts(dG%HI, param_file, CS%mixedlayer_restrat_CSp, CS%restart_CSp) ! Initialize fields call callTree_waypoint("restart registration complete (initialize_MOM)") @@ -1809,16 +1832,22 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("returned from ALE_init() (initialize_MOM)") endif - ! Shift from using the temporary dynamic grid type to using the final (potentially - ! static) ocean grid type. - ! call clone_MOM_domain(dG%Domain, CS%G%Domain) - ! call MOM_grid_init(CS%G, param_file) + ! Shift from using the temporary dynamic grid type to using the final + ! (potentially static) ocean-specific grid type. + ! The next line would be needed if G%Domain had not already been init'd above: + ! call clone_MOM_domain(dG%Domain, G%Domain) + call MOM_grid_init(G, param_file, HI, bathymetry_at_vel=bathy_at_vel) + call copy_dyngrid_to_MOM_grid(dG, G) + call destroy_dyn_horgrid(dG) - call copy_dyngrid_to_MOM_grid(dg, G) + ! Set a few remaining fields that are specific to the ocean grid type. + call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. if (CS%debug .or. G%symmetric) & call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) - G%ke = GV%ke + ! Copy common variables from the vertical grid to the horizontal grid. + ! Consider removing this later? + G%ke = GV%ke ; G%g_Earth = GV%g_Earth call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, & dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & @@ -1827,7 +1856,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") ! From this point, there may be pointers being set, so the final grid type - ! that will persist through the run has to be used. + ! that will persist throughout the run has to be used. if (test_grid_copy) then ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G. @@ -1839,17 +1868,14 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call copy_MOM_grid_to_dyngrid(G, dg) call copy_dyngrid_to_MOM_grid(dg, CS%G) - ! Copy a common variable from the vertical grid to the horizontal grid. - ! Consider removing this later? - CS%G%ke = GV%ke call destroy_dyn_horgrid(dG) call MOM_grid_end(G) ; deallocate(G) G => CS%G - if (CS%debug .or. CS%G%symmetric) & call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) + G%ke = GV%ke ; G%g_Earth = GV%g_Earth endif diff --git a/src/core/MOM_dyn_horgrid.F90 b/src/core/MOM_dyn_horgrid.F90 index 498c2bfa23..2d42e7f13a 100644 --- a/src/core/MOM_dyn_horgrid.F90 +++ b/src/core/MOM_dyn_horgrid.F90 @@ -125,7 +125,7 @@ module MOM_dyn_horgrid CoriolisBu ! The Coriolis parameter at corner points, in s-1. real, allocatable, dimension(:,:) :: & dF_dx, dF_dy ! Derivatives of f (Coriolis parameter) at h-points, in s-1 m-1. - real :: g_Earth ! The gravitational acceleration in m s-2. +! real :: g_Earth ! The gravitational acceleration in m s-2. ! These variables are global sums that are useful for 1-d diagnostics real :: areaT_global ! Global sum of h-cell area in m2 @@ -146,9 +146,13 @@ module MOM_dyn_horgrid !--------------------------------------------------------------------- !> Allocate memory used by the dyn_horgrid_type and related structures. -subroutine create_dyn_horgrid(G, HI) +subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel) type(dyn_horgrid_type), pointer :: G !< A pointer to the dynamic horizontal grid type type(hor_index_type), intent(in) :: HI !< A hor_index_type for array extents + logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are + !! separate values for the basin depths at velocity + !! points. Otherwise the effects of topography are + !! entirely determined from thickness points. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg ! This subroutine allocates the lateral elements of the dyn_horgrid_type that @@ -174,6 +178,9 @@ subroutine create_dyn_horgrid(G, HI) G%isd_global = G%isd + HI%idg_offset ; G%jsd_global = G%jsd + HI%jdg_offset G%symmetric = HI%symmetric + G%bathymetry_at_vel = .false. + if (present(bathymetry_at_vel)) G%bathymetry_at_vel = bathymetry_at_vel + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg @@ -232,10 +239,12 @@ subroutine create_dyn_horgrid(G, HI) allocate(G%sin_rot(isd:ied,jsd:jed)) ; G%sin_rot(:,:) = 0.0 allocate(G%cos_rot(isd:ied,jsd:jed)) ; G%cos_rot(:,:) = 1.0 - allocate(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = 0.0 - allocate(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = 0.0 - allocate(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = 0.0 - allocate(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = 0.0 + if (G%bathymetry_at_vel) then + allocate(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = 0.0 + allocate(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = 0.0 + allocate(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = 0.0 + allocate(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = 0.0 + endif allocate(G%gridLonT(isg:ieg)) ; G%gridLonT(:) = 0.0 allocate(G%gridLonB(G%IsgB:G%IegB)) ; G%gridLonB(:) = 0.0 @@ -333,8 +342,10 @@ subroutine destroy_dyn_horgrid(G) deallocate(G%dF_dx) ; deallocate(G%dF_dy) deallocate(G%sin_rot) ; deallocate(G%cos_rot) - deallocate(G%Dblock_u) ; deallocate(G%Dopen_u) - deallocate(G%Dblock_v) ; deallocate(G%Dopen_v) + if (allocated(G%Dblock_u)) deallocate(G%Dblock_u) + if (allocated(G%Dopen_u)) deallocate(G%Dopen_u) + if (allocated(G%Dblock_v)) deallocate(G%Dblock_v) + if (allocated(G%Dopen_v)) deallocate(G%Dopen_v) deallocate(G%gridLonT) ; deallocate(G%gridLatT) deallocate(G%gridLonB) ; deallocate(G%gridLatB) @@ -342,6 +353,8 @@ subroutine destroy_dyn_horgrid(G) ! Do not deallocate G%Domain%mpp_domain or deallocate(G%Domain) because ! these are pointers to types used elsewhere. G%Domain => NULL() + ! deallocate(G%Domain%mpp_domain) + ! deallocate(G%Domain) deallocate(G) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index de44bd8f7f..cae09ffcbd 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -151,13 +151,18 @@ module MOM_grid !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> MOM_grid_init initializes the ocean grid array sizes and grid memory. -subroutine MOM_grid_init(G, param_file, HI) +subroutine MOM_grid_init(G, param_file, HI, global_indexing, bathymetry_at_vel) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(hor_index_type), optional, intent(in) :: HI !< A hor_index_type for array extents -! Arguments: G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. + type(hor_index_type), & + optional, intent(in) :: HI !< A hor_index_type for array extents + logical, optional, intent(in) :: global_indexing !< If true use global index + !! values instead of having the data domain on each + !! processor start at 1. + logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are + !! separate values for the ocean bottom depths at + !! velocity points. Otherwise the effects of topography + !! are entirely determined from thickness points. ! This include declares and sets the variable "version". #include "version_variable.h" @@ -165,8 +170,9 @@ subroutine MOM_grid_init(G, param_file, HI) integer :: IsdB, IedB, JsdB, JedB integer :: ied_max, jed_max integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j - logical :: global_indexing ! If true use global index values instead of having + logical :: local_indexing ! If false use global index values instead of having ! the data domain on each processor start at 1. + integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend character(len=40) :: mod_nm = "MOM_grid" ! This module's name. @@ -174,29 +180,7 @@ subroutine MOM_grid_init(G, param_file, HI) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mod_nm, version, & "Parameters providing information about the lateral grid.") -! call get_param(param_file, "MOM", "G_EARTH", G%g_Earth, & -! "The gravitational acceleration of the Earth.", & -! units="m s-2", default = 9.80) - call get_param(param_file, mod_nm, "GLOBAL_INDEXING", global_indexing, & - "If true, use a global lateral indexing convention, so \n"//& - "that corresponding points on different processors have \n"//& - "the same index. This does not work with static memory.", & - default=.false., layoutParam=.true.) -#ifdef STATIC_MEMORY_ - if (global_indexing) call MOM_error(FATAL, "MOM_grid_init : "//& - "GLOBAL_INDEXING can not be true with STATIC_MEMORY.") -#endif - call get_param(param_file, mod_nm, "FIRST_DIRECTION", G%first_direction, & - "An integer that indicates which direction goes first \n"//& - "in parts of the code that use directionally split \n"//& - "updates, with even numbers (or 0) used for x- first \n"//& - "and odd numbers used for y-first.", default=0) - - call get_param(param_file, mod_nm, "BATHYMETRY_AT_VEL", G%bathymetry_at_vel, & - "If true, there are separate values for the basin depths \n"//& - "at velocity points. Otherwise the effects of of \n"//& - "topography are entirely determined from thickness points.", & - default=.false.) + call get_param(param_file, mod_nm, "NIBLOCK", niblock, "The number of blocks "// & "in the x-direction on each processor (for openmp).", default=1, & @@ -220,8 +204,10 @@ subroutine MOM_grid_init(G, param_file, HI) G%isd_global = G%isd + HI%idg_offset ; G%jsd_global = G%jsd + HI%jdg_offset G%symmetric = HI%symmetric else + local_indexing = .true. + if (present(global_indexing)) local_indexing = .not.global_indexing call hor_index_init(G%Domain, G%HI, param_file, & - local_indexing=.not.global_indexing) + local_indexing=local_indexing) ! get_domain_extent ensures that domains start at 1 for compatibility between ! static and dynamically allocated arrays, unless global_indexing is true. @@ -229,7 +215,7 @@ subroutine MOM_grid_init(G, param_file, HI) G%isd, G%ied, G%jsd, G%jed, & G%isg, G%ieg, G%jsg, G%jeg, & G%idg_offset, G%jdg_offset, G%symmetric, & - local_indexing=.not.global_indexing) + local_indexing=local_indexing) G%isd_global = G%isd+G%idg_offset ; G%jsd_global = G%jsd+G%jdg_offset endif @@ -254,6 +240,9 @@ subroutine MOM_grid_init(G, param_file, HI) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + G%bathymetry_at_vel = .false. + if (present(bathymetry_at_vel)) G%bathymetry_at_vel = bathymetry_at_vel if (G%bathymetry_at_vel) then ALLOC_(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = 0.0 ALLOC_(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = 0.0 @@ -535,6 +524,11 @@ subroutine MOM_grid_end(G) DEALLOC_(G%dF_dx) ; DEALLOC_(G%dF_dy) DEALLOC_(G%sin_rot) ; DEALLOC_(G%cos_rot) + if (G%bathymetry_at_vel) then + DEALLOC_(G%Dblock_u) ; DEALLOC_(G%Dopen_u) + DEALLOC_(G%Dblock_v) ; DEALLOC_(G%Dopen_v) + endif + deallocate(G%gridLonT) ; deallocate(G%gridLatT) deallocate(G%gridLonB) ; deallocate(G%gridLatB) diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index def4c1f606..b21b96e58f 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -139,7 +139,6 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG) oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon oG%Rad_Earth = dG%Rad_Earth ; oG%max_depth = dG%max_depth - oG%g_Earth = dG%g_Earth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(oG%areaT, oG%Domain) @@ -290,7 +289,6 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG) dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth - dG%g_Earth = oG%g_Earth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(dG%areaT, dG%Domain) diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index b13c443743..51034fac3a 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -100,7 +100,7 @@ subroutine verticalGridInit( param_file, GV ) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mod, version, & "Parameters providing information about the vertical grid.") - call get_param(param_file, "MOM", "G_EARTH", GV%g_Earth, & + call get_param(param_file, mod, "G_EARTH", GV%g_Earth, & "The gravitational acceleration of the Earth.", & units="m s-2", default = 9.80) call get_param(param_file, mod, "RHO_0", GV%Rho0, & From 984cc750c431e9c1c43c595ba9d3fcf61cc2577f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Jun 2016 16:41:49 -0400 Subject: [PATCH 3/5] Moved MOM_dyn_horgrid.F90 into src/framework Moved MOM_dyn_horgrid.F90 into the MOM6/src/framework directory to facilitate sharing this file with other components. All answers are bitwise identical. --- src/{core => framework}/MOM_dyn_horgrid.F90 | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{core => framework}/MOM_dyn_horgrid.F90 (100%) diff --git a/src/core/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 similarity index 100% rename from src/core/MOM_dyn_horgrid.F90 rename to src/framework/MOM_dyn_horgrid.F90 From 60f68089ee4c9b0de968851ea73961287ca77176 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Jun 2016 16:42:29 -0400 Subject: [PATCH 4/5] Deallocate G%Domain in destroy_dyn_horgrid Added calls to deallocate G%Domain and G%Domain%mpp_domain in destroy_dyn_horgrid. This eliminates a potential memory leak. It had not been done before because there was some concern that these pointers might have been set to point to an element of another type, but this is now all being done via a copy, so they can (and should) be deallocated. All answers are bitwise identical. --- src/framework/MOM_dyn_horgrid.F90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 2d42e7f13a..8a57ac6c90 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -125,7 +125,6 @@ module MOM_dyn_horgrid CoriolisBu ! The Coriolis parameter at corner points, in s-1. real, allocatable, dimension(:,:) :: & dF_dx, dF_dy ! Derivatives of f (Coriolis parameter) at h-points, in s-1 m-1. -! real :: g_Earth ! The gravitational acceleration in m s-2. ! These variables are global sums that are useful for 1-d diagnostics real :: areaT_global ! Global sum of h-cell area in m2 @@ -350,11 +349,8 @@ subroutine destroy_dyn_horgrid(G) deallocate(G%gridLonT) ; deallocate(G%gridLatT) deallocate(G%gridLonB) ; deallocate(G%gridLatB) - ! Do not deallocate G%Domain%mpp_domain or deallocate(G%Domain) because - ! these are pointers to types used elsewhere. - G%Domain => NULL() - ! deallocate(G%Domain%mpp_domain) - ! deallocate(G%Domain) + deallocate(G%Domain%mpp_domain) + deallocate(G%Domain) deallocate(G) From 14200dad29cff986a0cebcff1d231094b3340b13 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Jun 2016 17:57:15 -0400 Subject: [PATCH 5/5] Eliminate MOM_memory.h from fixed_init code MOM_memory.h sets various macros, but these should no longer be used in the strictly dynamic MOM_fixed_initialization and MOM_grid_initialize code, and hence the lines to #include have been eliminated. A array size few declarations were changed to avoid using these macros, and one routine was dOxyGenized. All answers and interfaces are bitwise identical. --- .../MOM_fixed_initialization.F90 | 23 ++++++++++--------- src/initialization/MOM_grid_initialize.F90 | 3 --- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 786cb19b16..252c4778c1 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -13,7 +13,6 @@ module MOM_fixed_initialization use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_file_parser, only : log_version -! use MOM_grid, only : ocean_grid_type use MOM_io, only : close_file, create_file, fieldtype, file_exists use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field, var_desc @@ -33,8 +32,6 @@ module MOM_fixed_initialization implicit none ; private -#include - public MOM_initialize_fixed, MOM_initialize_rotation, MOM_initialize_topography character(len=40) :: mod = "MOM_fixed_initialization" ! This module's name. @@ -271,10 +268,13 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) end subroutine MOM_initialize_topography ! ----------------------------------------------------------------------------- + +!> Return the global maximum ocean bottom depth in m. function diagnoseMaximumDepth(D,G) - type(dyn_horgrid_type), intent(in) :: G - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: D - real :: diagnoseMaximumDepth + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(in) :: D !< Ocean bottom depth in m + real :: diagnoseMaximumDepth !< The global maximum ocean bottom depth in m ! Local variables integer :: i,j diagnoseMaximumDepth=D(G%isc,G%jsc) @@ -877,7 +877,7 @@ end subroutine reset_face_lengths_file ! ----------------------------------------------------------------------------- subroutine reset_face_lengths_list(G, param_file) type(dyn_horgrid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file + type(param_file_type), intent(in) :: param_file ! This subroutine sets the open face lengths at selected points to restrict ! passages to their observed widths. @@ -1212,10 +1212,11 @@ subroutine write_ocean_geometry_file(G, param_file, directory) integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB logical :: multiple_files - real :: out_h(SZI_(G),SZJ_(G)) - real :: out_u(SZIB_(G),SZJ_(G)) - real :: out_v(SZI_(G),SZJB_(G)) - real :: out_q(SZIB_(G),SZJB_(G)) + real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q + real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u + real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 46f4dc9c18..61a6f5ed66 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -76,7 +76,6 @@ module MOM_grid_initialize use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -! use MOM_grid, only : ocean_grid_type, set_derived_metrics use MOM_io, only : read_data, slasher, file_exists use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE @@ -84,8 +83,6 @@ module MOM_grid_initialize implicit none ; private -#include - public set_grid_metrics, initialize_masks, Adcroft_reciprocal type, public :: GPS ; private