+Use dyn_horgrid_type in fixed initialization

  Made extensive changes to replace the ocean_grid_type with a dyn_horgrid_type
during the fixed initialization of MOM6.  This is a big step toward being able
to share MOM6 fixed initialization software (of masks, topography, and the
Coriolis parameter) with other components that have their own layouts.  This
changes the type of a number of arguments, and adds some arguments in a few
cases of ..._initialize_topography files (for standardization).  Also added
dOxyGen comments for some subroutines.  All answers are bitwise identical.

src/core/MOM.F90

@@ -1443,8 +1443,15 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   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 clone_MOM_domain(G%Domain, dG%Domain)
+
   call verticalGridInit( param_file, CS%GV )
   GV => CS%GV
+  dG%g_Earth = GV%g_Earth
 
   ! Read relevant parameters and write them to the model log.
   call log_version(param_file, "MOM", version, "")
@@ -1656,20 +1663,20 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
     "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
 
   ! 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.)
+  if (CS%debug .or. dG%symmetric) &
+    call clone_MOM_domain(dG%Domain, dG%Domain_aux, symmetric=.false.)
 
   call MOM_timing_init(CS)
 
   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
+!  G%ke = GV%ke
 
-  is   = G%isc   ; ie   = G%iec  ; js   = G%jsc  ; je   = G%jec ; nz = GV%ke
-  isd  = G%isd   ; ied  = G%ied  ; jsd  = G%jsd  ; jed  = G%jed
-  IsdB = G%IsdB  ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
+  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
 
   ! Allocate and initialize space for primary MOM variables.
   ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz))   ; CS%u(:,:,:) = 0.0
@@ -1784,17 +1791,28 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   call callTree_waypoint("restart registration complete (initialize_MOM)")
 
   call cpu_clock_begin(id_clock_MOM_init)
-  call MOM_initialize_fixed(G, param_file, write_geom_files, dirs%output_directory)
+  call MOM_initialize_fixed(dG, param_file, write_geom_files, dirs%output_directory)
   call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
   call MOM_initialize_coord(GV, param_file, write_geom_files, &
-                            dirs%output_directory, CS%tv, G%max_depth)
+                            dirs%output_directory, CS%tv, dG%max_depth)
   call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
 
   if (CS%use_ALE_algorithm) then
-    call ALE_init(param_file, GV, G%max_depth, CS%ALE_CSp)
+    call ALE_init(param_file, GV, dG%max_depth, CS%ALE_CSp)
     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)
+
+  call copy_dyngrid_to_MOM_grid(dg, G)
+  ! 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
+
   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, &
                             CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in)
@@ -1804,20 +1822,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   ! From this point, there may be pointers being set, so the final grid type
   ! that will persist through the run has to be used.
 
-  !   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)
-  ! call copy_dyngrid_to_MOM_grid(dg, CS%G)
-  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
-  ! if (CS%debug .or. CS%G%symmetric) &
-  !  call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.)
-
-  ! ! Copy a common variable from the vertical grid to the horizontal grid.
-  ! ! Consider removing this later?
-  ! CS%G%ke = GV%ke
-  ! G => CS%G
-
   if (test_grid_copy) then
     !  Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
     call create_dyn_horgrid(dG, G%HI)

src/core/MOM_transcribe_grid.F90

@@ -139,6 +139,7 @@ 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)
@@ -288,7 +289,8 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG)
   dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global
   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%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)

src/ice_shelf/MOM_ice_shelf.F90

@@ -99,7 +99,9 @@ module MOM_ice_shelf
 use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
 use MOM_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid
 use MOM_diag_mediator, only : diag_ctrl, time_type, enable_averaging, disable_averaging
-use MOM_domains, only : MOM_domains_init, pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE
+use MOM_domains, only : MOM_domains_init, clone_MOM_domain
+use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE
+use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
 use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
 use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
 use MOM_grid, only : MOM_grid_init, ocean_grid_type
@@ -113,6 +115,7 @@ module MOM_ice_shelf
 use MOM_restart, only : register_restart_field, query_initialized, save_restart
 use MOM_restart, only : restart_init, restore_state, MOM_restart_CS
 use MOM_time_manager, only : time_type, set_time, time_type_to_real
+use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid
 !use MOM_variables, only : forcing, surface
 use MOM_variables, only : surface
 use MOM_forcing_type, only : forcing, allocate_forcing_type
@@ -152,6 +155,7 @@ module MOM_ice_shelf
 type, public :: ice_shelf_CS ; private
   type(MOM_restart_CS), pointer :: restart_CSp => NULL()
   type(ocean_grid_type) :: grid !< Grid for the ice-shelf model
+!  type(dyn_horgrid_type), pointer :: dG  !< Dynamic grid for the ice-shelf model
   type(ocean_grid_type), pointer :: ocn_grid => NULL() !< A pointer to the ocean model grid
   ! The rest is private
   real ::   flux_factor = 1.0
@@ -685,7 +689,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS)
   fluxes%iceshelf_melt = CS%lprec  * (86400.0*365.0/CS%density_ice)
 
   if (CS%DEBUG) then
-   call hchksum (CS%h_shelf, "melt rate", G, haloshift=0)
+   call hchksum (CS%h_shelf, "melt rate", G%HI, haloshift=0)
   endif
 
   if (CS%shelf_mass_is_dynamic) then
@@ -830,10 +834,10 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
 
   if (CS%debug) then
     if (associated(state%taux_shelf)) then
-      call uchksum(state%taux_shelf, "taux_shelf", G, haloshift=0)
+      call uchksum(state%taux_shelf, "taux_shelf", G%HI, haloshift=0)
     endif
     if (associated(state%tauy_shelf)) then
-      call vchksum(state%tauy_shelf, "tauy_shelf", G, haloshift=0)
+      call vchksum(state%tauy_shelf, "tauy_shelf", G%HI, haloshift=0)
     endif
   endif
 
@@ -896,7 +900,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes)
     endif
   enddo ; enddo
   if (CS%debug) then
-    call hchksum(fluxes%ustar_shelf, "ustar_shelf", G, haloshift=0)
+    call hchksum(fluxes%ustar_shelf, "ustar_shelf", G%HI, haloshift=0)
   endif
 
   ! If the shelf mass is changing, the fluxes%rigidity_ice_[uv] needs to be
@@ -973,10 +977,10 @@ end subroutine add_shelf_flux
 
 !   if (CS%debug) then
 !     if (associated(state%taux_shelf)) then
-!       call uchksum(state%taux_shelf, "taux_shelf", G, haloshift=0)
+!       call uchksum(state%taux_shelf, "taux_shelf", G%HI, haloshift=0)
 !     endif
 !     if (associated(state%tauy_shelf)) then
-!       call vchksum(state%tauy_shelf, "tauy_shelf", G, haloshift=0)
+!       call vchksum(state%tauy_shelf, "tauy_shelf", G%HI, haloshift=0)
 !     endif
 !   endif
 
@@ -1021,7 +1025,7 @@ end subroutine add_shelf_flux
 !   enddo ; enddo
 
 !   if (CS%debug) then
-!     call hchksum(fluxes%ustar_shelf, "ustar_shelf", G, haloshift=0)
+!     call hchksum(fluxes%ustar_shelf, "ustar_shelf", G%HI, haloshift=0)
 !   endif
 
 !   ! If the shelf mass is changing, the fluxes%rigidity_ice_[uv] needs to be
@@ -1057,6 +1061,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti
   type(ocean_grid_type), pointer :: G, OG ! Convenience pointers
   type(directories)  :: dirs
   type(vardesc) :: vd
+  type(dyn_horgrid_type), pointer :: dG => NULL()
   real :: cdrag, drag_bg_vel
   logical :: new_sim, save_IC, var_force
 ! This include declares and sets the variable "version".
@@ -1089,7 +1094,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti
   call MOM_domains_init(CS%grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_)
 ! call diag_mediator_init(CS%grid,param_file,CS%diag) ! this needs to be fixed - will probably break when not using coupled driver 0
   call MOM_grid_init(CS%grid, param_file)
-  call set_grid_metrics(CS%grid, param_file)
+
+  call create_dyn_horgrid(dG, CS%grid%HI)
+  call clone_MOM_domain(CS%grid%Domain, dG%Domain)
+
+  call set_grid_metrics(dG, param_file)
 ! call set_diag_mediator_grid(CS%grid, CS%diag)
 
   ! The ocean grid is possibly different
@@ -1396,9 +1405,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti
   endif
 
 ! Set up the bottom depth, G%D either analytically or from file
-  call MOM_initialize_topography(G%bathyT, G%max_depth, G, param_file)
+  call MOM_initialize_topography(G%bathyT, G%max_depth, dG, param_file)
 ! Set up the Coriolis parameter, G%f, usually analytically.
-  call MOM_initialize_rotation(G%CoriolisBu, G, param_file)
+  call MOM_initialize_rotation(G%CoriolisBu, dG, param_file)
+  call copy_dyngrid_to_MOM_grid(dG, CS%grid)
+
+  call destroy_dyn_horgrid(dG)
 
   ! Set up the restarts.
   call restart_init(param_file, CS%restart_CSp, "Shelf.res")
@@ -2031,8 +2043,8 @@ subroutine ice_shelf_advect(CS, time_step, melt_rate, Time)
   call pass_var(CS%hmask, G%domain)
 
   if (CS%DEBUG) then
-    call hchksum (CS%h_shelf, "h after front", G, haloshift=3)
-    call hchksum (CS%h_shelf, "shelf area after front", G, haloshift=3)
+    call hchksum (CS%h_shelf, "h after front", G%HI, haloshift=3)
+    call hchksum (CS%h_shelf, "shelf area after front", G%HI, haloshift=3)
   endif
 
 
@@ -2257,8 +2269,8 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time)
 
 
     if (CS%DEBUG) then
-      call qchksum (u, "u shelf", G, haloshift=2)
-      call qchksum (v, "v shelf", G, haloshift=2)
+      call qchksum (u, "u shelf", G%HI, haloshift=2)
+      call qchksum (v, "v shelf", G%HI, haloshift=2)
     endif
 
     if (is_root_pe()) print *,"linear solve done",iters," iterations"
@@ -5920,7 +5932,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time)
   call pass_var(CS%tmask, G%domain)
 
   if (CS%DEBUG) then
-    call hchksum (CS%t_shelf, "temp after front", G, haloshift=3)
+    call hchksum (CS%t_shelf, "temp after front", G%HI, haloshift=3)
   endif
 
 end subroutine ice_shelf_temp