diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fe44a68e0b..d53e788a04 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -86,8 +86,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 @@ -1366,10 +1367,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' @@ -1385,7 +1387,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 @@ -1395,6 +1396,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 @@ -1487,6 +1497,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.) @@ -1582,6 +1597,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.", & @@ -1632,6 +1662,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)") @@ -1652,20 +1685,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) & @@ -1678,10 +1709,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 @@ -1702,8 +1729,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 @@ -1713,11 +1740,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 @@ -1770,30 +1793,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)") @@ -1810,16 +1833,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, & @@ -1828,7 +1857,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. @@ -1840,17 +1869,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_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, & diff --git a/src/core/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 similarity index 94% rename from src/core/MOM_dyn_horgrid.F90 rename to src/framework/MOM_dyn_horgrid.F90 index 498c2bfa23..8a57ac6c90 100644 --- a/src/core/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 @@ -146,9 +145,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 +177,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 +238,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,15 +341,16 @@ 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) - ! 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/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 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