From e2d90161adb5707b327975c686a1db1e722f27f9 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Thu, 23 Jun 2016 17:04:32 -0400
Subject: [PATCH 1/5] +Expanded use of GV in framework code

  Increased the use of the vertical grid type, and eliminated the use of the
horizontal grid type to specify the number of layers (i.e., the use of G%ke and
SZK_(G)) throughout the MOM6 framework code.  Also, renamed the variable "sum"
as "asum" to avoid name-space confusion with the F90 intrisic function sum() in
two routines in MOM_spatial_means.F90 and revised global_volume_mean to avoid
the unnecessary use of 3-d arrays.  All answers are bitwise identical, although
an argument was added to global_layer_mean, global_volume_mean,
diag_mediator_init, and diag_masks_set.
---
 src/core/MOM.F90                    |  4 +-
 src/diagnostics/MOM_diagnostics.F90 |  8 +--
 src/framework/MOM_checksums.F90     | 10 ++--
 src/framework/MOM_diag_mediator.F90 | 59 ++++++++++--------
 src/framework/MOM_restart.F90       |  2 +-
 src/framework/MOM_spatial_means.F90 | 92 +++++++++++++++--------------
 6 files changed, 95 insertions(+), 80 deletions(-)

diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 3afd5de7de..54b46ef5b5 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -1884,12 +1884,12 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
 
    diag    => CS%diag
   ! Initialize the diag mediator.
-  call diag_mediator_init(G, param_file, diag, doc_file_dir=dirs%output_directory)
+  call diag_mediator_init(G, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
 
   ! Initialize the diagnostics mask arrays.
   ! This step has to be done after call to MOM_initialize_state
   ! and before MOM_diagnostics_init
-  call diag_masks_set(G, CS%missing, diag)
+  call diag_masks_set(G, GV%ke, CS%missing, diag)
 
   ! Set up a pointers h within diag mediator control structure,
   ! this needs to occur _after_ CS%h has been allocated.
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 5dda3dce55..86e6812969 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -326,7 +326,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
 
   ! volume mean potential temperature 
   if (CS%id_thetaoga>0) then
-    thetaoga = global_volume_mean(tv%T, h, G)
+    thetaoga = global_volume_mean(tv%T, h, G, GV)
     call post_data(CS%id_thetaoga, thetaoga, CS%diag)
   endif
 
@@ -341,7 +341,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
 
   ! volume mean salinity 
   if (CS%id_soga>0) then
-    soga = global_volume_mean(tv%S, h, G)
+    soga = global_volume_mean(tv%S, h, G, GV)
     call post_data(CS%id_soga, soga, CS%diag)
   endif
 
@@ -356,13 +356,13 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
 
   ! layer mean potential temperature 
   if (CS%id_temp_layer_ave>0) then
-    temp_layer_ave = global_layer_mean(tv%T, h, G)
+    temp_layer_ave = global_layer_mean(tv%T, h, G, GV)
     call post_data_1d_k(CS%id_temp_layer_ave, temp_layer_ave, CS%diag)
   endif
 
   ! layer mean salinity 
   if (CS%id_salt_layer_ave>0) then
-    salt_layer_ave = global_layer_mean(tv%S, h, G)
+    salt_layer_ave = global_layer_mean(tv%S, h, G, GV)
     call post_data_1d_k(CS%id_salt_layer_ave, salt_layer_ave, CS%diag)
   endif
 
diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90
index c30f4a38d2..30a0c1398d 100644
--- a/src/framework/MOM_checksums.F90
+++ b/src/framework/MOM_checksums.F90
@@ -1156,10 +1156,11 @@ function totalStuff(G, hThick, stuff)
   real, dimension(G%isd:,G%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights
   real, dimension(G%isd:,G%jsd:,:), intent(in) :: stuff  !< The array of stuff to be summed 
   real                                         :: totalStuff
-  integer :: i, j, k
+  integer :: i, j, k, nz
 
+  nz = size(hThick,3)
   totalStuff = 0.
-  do k = 1, G%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
+  do k = 1, nz ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
     totalStuff = totalStuff + hThick(i,j,k) * stuff(i,j,k) * G%areaT(i,j)
   enddo ; enddo ; enddo
   call sum_across_PEs(totalStuff)
@@ -1185,10 +1186,11 @@ subroutine totalTandS(G, hThick, temperature, salinity, mesg)
 
   logical, save :: firstCall = .true.
   real :: thisH, thisT, thisS, delH, delT, delS
-  integer :: i, j, k
+  integer :: i, j, k, nz
 
+  nz = size(hThick,3)
   thisH = 0.
-  do k = 1, G%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
+  do k = 1, nz ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
     thisH = thisH + hThick(i,j,k) * G%areaT(i,j)
   enddo ; enddo ; enddo
   call sum_across_PEs(thisH)
diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90
index 233877f152..dba5253cba 100644
--- a/src/framework/MOM_diag_mediator.F90
+++ b/src/framework/MOM_diag_mediator.F90
@@ -204,7 +204,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical)
   integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_zzl, id_zzi
   integer :: k, nz
   integer :: nzi(4)
-  real :: zlev(G%ke), zinter(G%ke+1)
+  real :: zlev(GV%ke), zinter(GV%ke+1)
   logical :: set_vert
   character(len=200) :: inputdir, string, filename, varname, dimname
 
@@ -234,7 +234,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical)
               'h point nominal latitude', Domain2=G%Domain%mpp_domain)
 
   if (set_vert) then
-    nz = G%ke
+    nz = GV%ke
     zinter(1:nz+1) = GV%sInterface(1:nz+1)
     zlev(1:nz) = GV%sLayer(1:nz)
     id_zl = diag_axis_init('zl', zlev, trim(GV%zAxisUnits), 'z', &
@@ -262,7 +262,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical)
   if (len_trim(string) > 0) then
     if (trim(string) == 'UNIFORM') then
       ! initialise a uniform coordinate with depth
-      nzi(1) = G%ke + 1
+      nzi(1) = GV%ke + 1
       allocate(diag_cs%zi_remap(nzi(1)))
       allocate(diag_cs%zl_remap(nzi(1) - 1))
 
@@ -828,14 +828,14 @@ end subroutine remap_diag_to_z
 !! height changes.
 subroutine diag_update_target_grids(diag_cs)
   type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure
+
   ! Local variables
-  real, dimension(size(diag_cs%h, 3)) :: h_src
   type(ocean_grid_type), pointer :: G
   real :: depth
   integer :: nz_src, nz_dest
   integer :: i, j, k
   logical :: force, h_changed
-  real, dimension(diag_cs%G%ke) :: h_src_1d
+  real, dimension(size(diag_cs%h, 3)) :: h_src_1d
   real, dimension(diag_cs%nz_remap) :: h_zout_1d
   real, dimension(diag_cs%nz_remap+1) :: z_out_1d
 
@@ -1654,11 +1654,16 @@ subroutine diag_mediator_infrastructure_init(err_msg)
   call diag_manager_init(err_msg=err_msg)
 end subroutine diag_mediator_infrastructure_init
 
-subroutine diag_mediator_init(G, param_file, diag_cs, doc_file_dir)
-  type(ocean_grid_type), target, intent(inout) :: G
-  type(param_file_type),      intent(in)    :: param_file
-  type(diag_ctrl),            intent(inout) :: diag_cs
-  character(len=*), optional, intent(in)    :: doc_file_dir
+!> diag_mediator_init initializes the MOM diag_mediator and opens the available
+!! diagnostics file, if appropriate.
+subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir)
+  type(ocean_grid_type), target, intent(inout) :: G  !< The ocean grid type.
+  integer,                    intent(in)    :: nz    !< The number of layers in the model's native grid.
+  type(param_file_type),      intent(in)    :: param_file !< Parameter file structure
+  type(diag_ctrl),            intent(inout) :: diag_cs !< A pointer to a type with many variables
+                                                     !! used for diagnostics
+  character(len=*), optional, intent(in)    :: doc_file_dir !< A directory in which to create the
+                                                     !! file 
 
   ! This subroutine initializes the diag_mediator and the diag_manager.
   ! The grid type should have its dimensions set by this point, but it
@@ -1693,7 +1698,7 @@ subroutine diag_mediator_init(G, param_file, diag_cs, doc_file_dir)
   diag_cs%do_z_remapping_on_T = .false.
   diag_cs%remapping_initialized = .false.
 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
-  allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,G%ke))
+  allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz))
   diag_cs%h_old(:,:,:) = 0.0
 #endif
 
@@ -1755,11 +1760,13 @@ subroutine diag_set_thickness_ptr(h, diag_cs)
 
 end subroutine
 
-subroutine diag_masks_set(G, missing_value, diag_cs)
-! Setup the 2d masks for diagnostics
-  type(ocean_grid_type), target, intent(in) :: G
-  real,                          intent(in) :: missing_value
-  type(diag_ctrl),          pointer    :: diag_cs
+!> diag_masks_set sets up the 2d and 3d masks for diagnostics
+subroutine diag_masks_set(G, nz, missing_value, diag_cs)
+  type(ocean_grid_type), target, intent(in) :: G  !< The ocean grid type.
+  integer,                       intent(in) :: nz !< The number of layers in the model's native grid.
+  real,                          intent(in) :: missing_value !< A value to use for masked points.
+  type(diag_ctrl),               pointer    :: diag_cs !< A pointer to a type with many variables
+                                                  !! used for diagnostics
   ! Local variables
   integer :: k
 
@@ -1767,21 +1774,21 @@ subroutine diag_masks_set(G, missing_value, diag_cs)
   diag_cs%mask2dBu=> G%mask2dBu
   diag_cs%mask2dCu=> G%mask2dCu
   diag_cs%mask2dCv=> G%mask2dCv
-  allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:G%ke))
-  allocate(diag_cs%mask3dBuL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:G%ke))
-  allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:G%ke))
-  allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:G%ke))
-  do k = 1,G%ke
+  allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:nz))
+  allocate(diag_cs%mask3dBuL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz))
+  allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:nz))
+  allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:nz))
+  do k=1,nz
     diag_cs%mask3dTL(:,:,k)  = diag_cs%mask2dT (:,:)
     diag_cs%mask3dBuL(:,:,k) = diag_cs%mask2dBu(:,:)
     diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:)
     diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:)
   enddo
-  allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:G%ke+1))
-  allocate(diag_cs%mask3dBui(G%IsdB:G%IedB,G%JsdB:G%JedB,1:G%ke+1))
-  allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:G%ke+1))
-  allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:G%ke+1))
-  do k = 1,G%ke+1
+  allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:nz+1))
+  allocate(diag_cs%mask3dBui(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz+1))
+  allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:nz+1))
+  allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:nz+1))
+  do k=1,nz+1
     diag_cs%mask3dTi(:,:,k)  = diag_cs%mask2dT (:,:)
     diag_cs%mask3dBui(:,:,k) = diag_cs%mask2dBu(:,:)
     diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:)
diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90
index 57a8be9451..a564bd9cfa 100644
--- a/src/framework/MOM_restart.F90
+++ b/src/framework/MOM_restart.F90
@@ -759,7 +759,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV)
 
   num_files = 0
   next_var = 0
-  nz = G%ke
+  nz = 1 ; if (present(GV)) nz = GV%ke
 
   call get_time(time,seconds,days)
   restart_time = real(days) + real(seconds)/86400.0
diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90
index 0bce719151..54f58b7269 100644
--- a/src/framework/MOM_spatial_means.F90
+++ b/src/framework/MOM_spatial_means.F90
@@ -27,6 +27,7 @@ module MOM_spatial_means
 use MOM_error_handler, only : MOM_error, NOTE, WARNING, FATAL, is_root_pe
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_grid, only : ocean_grid_type
+use MOM_verticalGrid, only : verticalGrid_type
 
 implicit none ; private
 
@@ -74,21 +75,22 @@ function global_area_integral(var,G)
 
 end function global_area_integral
 
-function global_layer_mean(var,h,G)
-  type(ocean_grid_type),                       intent(in)  :: G
-  real, dimension(SZI_(G), SZJ_(G), SZK_(G)),  intent(in)  :: var
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)  :: h
-  real, dimension(SZK_(G))                   :: global_layer_mean
+function global_layer_mean(var, h, G, GV)
+  type(ocean_grid_type),                     intent(in)  :: G
+  type(verticalGrid_type),                   intent(in)  :: GV
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: var
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h
+  real, dimension(SZK_(GV))                   :: global_layer_mean
 
-  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: tmpForSumming, weight
-  real, dimension(SZK_(G)) :: scalarij, weightij
-  real, dimension(SZK_(G)) :: global_temp_scalar, global_weight_scalar
+  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: tmpForSumming, weight
+  real, dimension(SZK_(GV)) :: scalarij, weightij
+  real, dimension(SZK_(GV)) :: global_temp_scalar, global_weight_scalar
   integer :: i, j, k, is, ie, js, je, nz
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.
 
-  do k=1, nz ; do j=js,je ; do i=is, ie
+  do k=1,nz ; do j=js,je ; do i=is,ie
     weight(i,j,k)  =  h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j))
     tmpForSumming(i,j,k) =  var(i,j,k) * weight(i,j,k)
   enddo ; enddo ; enddo
@@ -102,23 +104,27 @@ function global_layer_mean(var,h,G)
 
 end function global_layer_mean
 
-function global_volume_mean(var,h,G)
-  type(ocean_grid_type),                       intent(in)  :: G
-  real, dimension(SZI_(G), SZJ_(G), SZK_(G)),  intent(in)  :: var
-  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)  :: h
+function global_volume_mean(var, h, G, GV)
+  type(ocean_grid_type),                     intent(in)  :: G
+  type(verticalGrid_type),                   intent(in)  :: GV
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: var
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h
   real :: global_volume_mean
 
-  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: tmpForSumming, weight
+  real :: weight_here
+  real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming, sum_weight
   integer :: i, j, k, is, ie, js, je, nz
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
-  tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.
+  tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0.
 
-  do k=1, nz ; do j=js,je ; do i=is, ie
-    weight(i,j,k)  =  h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j))
-    tmpForSumming(i,j,k) =  var(i,j,k) * weight(i,j,k)
+  do k=1,nz ; do j=js,je ; do i=is,ie
+    weight_here  =  h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j))
+    tmpForSumming(i,j) = tmpForSumming(i,j) + var(i,j,k) * weight_here
+    sum_weight(i,j) = sum_weight(i,j) + weight_here
   enddo ; enddo ; enddo
-  global_volume_mean = (reproducing_sum(sum(tmpForSumming,3))) / (reproducing_sum(sum(weight,3)))
+  global_volume_mean = (reproducing_sum(tmpForSumming)) / &
+                       (reproducing_sum(sum_weight))
 
 end function global_volume_mean
 
@@ -137,7 +143,7 @@ subroutine global_i_mean(array, i_mean, G, mask)
 !  (in)      G - The ocean's grid structure.
 !  (in)      mask - An array used for weighting the i-mean.
 
-  type(EFP_type), allocatable, dimension(:) :: sum, mask_sum
+  type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
   real :: mask_sum_r
   integer :: is, ie, js, je, idg_off, jdg_off
   integer :: i, j
@@ -147,23 +153,23 @@ subroutine global_i_mean(array, i_mean, G, mask)
 
   call reset_EFP_overflow_error()
 
-  allocate(sum(G%jsg:G%jeg))
+  allocate(asum(G%jsg:G%jeg))
   if (present(mask)) then
     allocate(mask_sum(G%jsg:G%jeg))
 
     do j=G%jsg,G%jeg
-      sum(j) = real_to_EFP(0.0) ; mask_sum(j) = real_to_EFP(0.0)
+      asum(j) = real_to_EFP(0.0) ; mask_sum(j) = real_to_EFP(0.0)
     enddo
 
     do i=is,ie ; do j=js,je
-      sum(j+jdg_off) = sum(j+jdg_off) + real_to_EFP(array(i,j)*mask(i,j))
+      asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(array(i,j)*mask(i,j))
       mask_sum(j+jdg_off) = mask_sum(j+jdg_off) + real_to_EFP(mask(i,j))
     enddo ; enddo
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_i_mean overflow error occurred before sums across PEs.")
 
-    call EFP_list_sum_across_PEs(sum(G%jsg:G%jeg), G%jeg-G%jsg+1)
+    call EFP_list_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1)
     call EFP_list_sum_across_PEs(mask_sum(G%jsg:G%jeg), G%jeg-G%jsg+1)
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
@@ -172,32 +178,32 @@ subroutine global_i_mean(array, i_mean, G, mask)
     do j=js,je
       mask_sum_r = EFP_to_real(mask_sum(j+jdg_off))
       if (mask_sum_r == 0.0 ) then ; i_mean(j) = 0.0 ; else
-        i_mean(j) = EFP_to_real(sum(j+jdg_off)) / mask_sum_r
+        i_mean(j) = EFP_to_real(asum(j+jdg_off)) / mask_sum_r
       endif
     enddo
 
     deallocate(mask_sum)
   else
-    do j=G%jsg,G%jeg ; sum(j) = real_to_EFP(0.0) ; enddo
+    do j=G%jsg,G%jeg ; asum(j) = real_to_EFP(0.0) ; enddo
 
     do i=is,ie ; do j=js,je
-      sum(j+jdg_off) = sum(j+jdg_off) + real_to_EFP(array(i,j))
+      asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(array(i,j))
     enddo ; enddo
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_i_mean overflow error occurred before sum across PEs.")
 
-    call EFP_list_sum_across_PEs(sum(G%jsg:G%jeg), G%jeg-G%jsg+1)
+    call EFP_list_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1)
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_i_mean overflow error occurred during sum across PEs.")
 
     do j=js,je
-      i_mean(j) = EFP_to_real(sum(j+jdg_off)) / real(G%ieg-G%isg+1)
+      i_mean(j) = EFP_to_real(asum(j+jdg_off)) / real(G%ieg-G%isg+1)
     enddo
   endif
 
-  deallocate(sum)
+  deallocate(asum)
 
 end subroutine global_i_mean
 
@@ -215,7 +221,7 @@ subroutine global_j_mean(array, j_mean, G, mask)
 !  (in)      G - The ocean's grid structure.
 !  (in)      mask - An array used for weighting the j-mean.
 
-  type(EFP_type), allocatable, dimension(:) :: sum, mask_sum
+  type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
   real :: mask_sum_r
   integer :: is, ie, js, je, idg_off, jdg_off
   integer :: i, j
@@ -225,23 +231,23 @@ subroutine global_j_mean(array, j_mean, G, mask)
 
   call reset_EFP_overflow_error()
 
-  allocate(sum(G%isg:G%ieg))
+  allocate(asum(G%isg:G%ieg))
   if (present(mask)) then
     allocate (mask_sum(G%isg:G%ieg))
 
     do i=G%isg,G%ieg
-      sum(i) = real_to_EFP(0.0) ; mask_sum(i) = real_to_EFP(0.0)
+      asum(i) = real_to_EFP(0.0) ; mask_sum(i) = real_to_EFP(0.0)
     enddo
 
     do i=is,ie ; do j=js,je
-      sum(i+idg_off) = sum(i+idg_off) + real_to_EFP(array(i,j)*mask(i,j))
+      asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(array(i,j)*mask(i,j))
       mask_sum(i+idg_off) = mask_sum(i+idg_off) + real_to_EFP(mask(i,j))
     enddo ; enddo
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_j_mean overflow error occurred before sums across PEs.")
 
-    call EFP_list_sum_across_PEs(sum(G%isg:G%ieg), G%ieg-G%isg+1)
+    call EFP_list_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1)
     call EFP_list_sum_across_PEs(mask_sum(G%isg:G%ieg), G%ieg-G%isg+1)
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
@@ -250,32 +256,32 @@ subroutine global_j_mean(array, j_mean, G, mask)
     do i=is,ie
       mask_sum_r = EFP_to_real(mask_sum(i+idg_off))
       if (mask_sum_r == 0.0 ) then ; j_mean(i) = 0.0 ; else
-        j_mean(i) = EFP_to_real(sum(i+idg_off)) / mask_sum_r
+        j_mean(i) = EFP_to_real(asum(i+idg_off)) / mask_sum_r
       endif
     enddo
 
     deallocate(mask_sum)
   else
-    do i=G%isg,G%ieg ; sum(i) = real_to_EFP(0.0) ; enddo
+    do i=G%isg,G%ieg ; asum(i) = real_to_EFP(0.0) ; enddo
 
     do i=is,ie ; do j=js,je
-      sum(i+idg_off) = sum(i+idg_off) + real_to_EFP(array(i,j))
+      asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(array(i,j))
     enddo ; enddo
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_j_mean overflow error occurred before sum across PEs.")
 
-    call EFP_list_sum_across_PEs(sum(G%isg:G%ieg), G%ieg-G%isg+1)
+    call EFP_list_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1)
 
     if (query_EFP_overflow_error()) call MOM_error(FATAL, &
       "global_j_mean overflow error occurred during sum across PEs.")
 
     do i=is,ie
-      j_mean(i) = EFP_to_real(sum(i+idg_off)) / real(G%jeg-G%jsg+1)
+      j_mean(i) = EFP_to_real(asum(i+idg_off)) / real(G%jeg-G%jsg+1)
     enddo
   endif
 
-  deallocate(sum)
+  deallocate(asum)
 
 end subroutine global_j_mean
 

From 7a4803e438fc8f341f894a7923585da05e86a8f6 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Thu, 23 Jun 2016 18:20:59 -0400
Subject: [PATCH 2/5] +Revised create_file to use ocean_grid or dyn_grid

  Modified create_file and reopen_file to permit the use of either an ocean grid
type variable or dyn_horgrid_type variable.  Both have similar information, and
either can be provided as an optional argument.  All results are bitwise
identical, but because the ocean_grid_type is now an optional argument to
reopen_file, the order of its arguments has changed.
---
 src/diagnostics/MOM_sum_output.F90 |   3 +-
 src/framework/MOM_io.F90           | 101 +++++++++++++++++++++--------
 2 files changed, 76 insertions(+), 28 deletions(-)

diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index d8101ca9c3..d0aa3cd65c 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -592,7 +592,8 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp)
     energypath_nc = trim(CS%energyfile) // ".nc"
     if (day > CS%Start_time) then
       call reopen_file(CS%fileenergy_nc, trim(energypath_nc), vars, &
-                       num_nc_fields, G, CS%fields, SINGLE_FILE, CS%timeunit, GV=GV)
+                       num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, &
+                       G=G, GV=GV)
     else
       call create_file(CS%fileenergy_nc, trim(energypath_nc), vars, &
                        num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, &
diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90
index 7a1741d7b3..91cf9137bf 100644
--- a/src/framework/MOM_io.F90
+++ b/src/framework/MOM_io.F90
@@ -7,8 +7,9 @@ module MOM_io
 use MOM_error_handler,    only : MOM_error, NOTE, FATAL, WARNING
 use MOM_domains,          only : MOM_domain_type
 use MOM_file_parser,      only : log_version, param_file_type
-use MOM_string_functions, only : lowercase
 use MOM_grid,             only : ocean_grid_type
+use MOM_dyn_horgrid,      only : dyn_horgrid_type
+use MOM_string_functions, only : lowercase
 use MOM_verticalGrid,     only : verticalGrid_type
 
 use ensemble_manager_mod, only : get_ensemble_id
@@ -76,7 +77,7 @@ module MOM_io
 !> Routine creates a new NetCDF file.  It also sets up
 !! structures that describe this file and variables that will
 !! later be written to this file. Type for describing a variable, typically a tracer
-subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, GV)
+subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV)
   integer,               intent(out)   :: unit       !< unit id of an open file or -1 on a
                                                      !! nonwriting PE with single file output
   character(len=*),      intent(in)    :: filename   !< full path to the file to create
@@ -86,8 +87,11 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
   integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE
   real, optional,        intent(in)    :: timeunit   !< length, in seconds, of the units for time. The
                                                      !! default value is 86400.0, for 1 day.
-  type(ocean_grid_type),   optional, intent(in) :: G  !< ocean grid structure structure, which is
-                                                     !! required if the new file uses any
+  type(ocean_grid_type),   optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
+                                                     !! is required if the new file uses any
+                                                     !! horizontal grid axes.
+  type(dyn_horgrid_type),  optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG
+                                                     !! is required if the new file uses any
                                                      !! horizontal grid axes.
   type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is
                                                      !! required if the new file uses any
@@ -95,15 +99,21 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
 
   logical        :: use_lath, use_lonh, use_latq, use_lonq, use_time
   logical        :: use_layer, use_int, use_periodic
-  logical        :: one_file
+  logical        :: one_file, domain_set
   type(axistype) :: axis_lath, axis_latq, axis_lonh, axis_lonq
   type(axistype) :: axis_layer, axis_int, axis_time, axis_periodic
   type(axistype) :: axes(4)
+  type(MOM_domain_type), pointer :: Domain => NULL()
   type(domain1d) :: x_domain, y_domain
   integer        :: numaxes, pack, thread, k
+  integer        :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB
   integer        :: var_periods, num_periods=0
   real, dimension(:), allocatable :: period_val
-  character(len=40) :: time_units
+  real, pointer, dimension(:) :: &
+    gridLatT => NULL(), & ! The latitude or longitude of T or B points for
+    gridLatB => NULL(), & ! the purpose of labeling the output axes.
+    gridLonT => NULL(), gridLonB => NULL()
+  character(len=40) :: time_units, x_axis_units, y_axis_units
   character(len=8)  :: t_grid, t_grid_read
 
   use_lath  = .false. ; use_lonh     = .false.
@@ -114,15 +124,32 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
   thread = SINGLE_FILE
   if (PRESENT(threading)) thread = threading
 
-  one_file = .true.
+  domain_set = .false.
   if (present(G)) then
-    one_file = ((thread == SINGLE_FILE) .or. .not.G%Domain%use_io_layout)
+    domain_set = .true. ; Domain => G%Domain
+    gridLatT => G%gridLatT ; gridLatB => G%gridLatB
+    gridLonT => G%gridLonT ; gridLonB => G%gridLonB
+    x_axis_units = G%x_axis_units ; y_axis_units = G%y_axis_units
+    isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg
+    IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB
+  elseif (present(dG)) then
+    domain_set = .true. ; Domain => dG%Domain
+    gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB
+    gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB
+    x_axis_units = dG%x_axis_units ; y_axis_units = dG%y_axis_units
+    isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg
+    IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB
+  endif
+
+  one_file = .true.
+  if (domain_set) then
+    one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout)
   endif
 
   if (one_file) then
     call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, threading=thread)
   else
-    call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, domain=G%Domain%mpp_domain)
+    call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, domain=Domain%mpp_domain)
   endif
 
 ! Define the coordinates.
@@ -183,10 +210,10 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
   end do
 
   if ((use_lath .or. use_lonh .or. use_latq .or. use_lonq)) then
-    if (.not.present(G)) call MOM_error(FATAL, "create_file: "//&
-      "An ocean_grid_type is required to create a file with a vertical coordinate.")
+    if (.not.domain_set) call MOM_error(FATAL, "create_file: "//&
+      "An ocean_grid_type or dyn_horgrid_type is required to create a file with a horizontal coordinate.")
 
-    call mpp_get_domain_components(G%Domain%mpp_domain, x_domain, y_domain)
+    call mpp_get_domain_components(Domain%mpp_domain, x_domain, y_domain)
   endif
   if ((use_layer .or. use_int) .and. .not.present(GV)) call MOM_error(FATAL, &
     "create_file: A vertical grid type is required to create a file with a vertical coordinate.")
@@ -194,20 +221,20 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit
 ! Specify all optional arguments to mpp_write_meta: name, units, longname, cartesian, calendar, sense, domain, data, min)
 ! Otherwise if optional arguments are added to mpp_write_meta the compiler may (and in case of GNU is) get confused and crash.
   if (use_lath) &
-    call mpp_write_meta(unit, axis_lath, name="lath", units=G%y_axis_units, longname="Latitude", &
-                   cartesian='Y', domain = y_domain, data=G%gridLatT(G%jsg:G%jeg))
+    call mpp_write_meta(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", &
+                   cartesian='Y', domain = y_domain, data=gridLatT(jsg:jeg))
 
   if (use_lonh) &
-    call mpp_write_meta(unit, axis_lonh, name="lonh", units=G%x_axis_units, longname="Longitude", &
-                   cartesian='X', domain = x_domain, data=G%gridLonT(G%isg:G%ieg))
+    call mpp_write_meta(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", &
+                   cartesian='X', domain = x_domain, data=gridLonT(isg:ieg))
 
   if (use_latq) &
-    call mpp_write_meta(unit, axis_latq, name="latq", units=G%y_axis_units, longname="Latitude", &
-                   cartesian='Y', domain = y_domain, data=G%gridLatB(G%JsgB:G%JegB))
+    call mpp_write_meta(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", &
+                   cartesian='Y', domain = y_domain, data=gridLatB(JsgB:JegB))
 
   if (use_lonq) &
-    call mpp_write_meta(unit, axis_lonq, name="lonq", units=G%x_axis_units, longname="Longitude", &
-                   cartesian='X', domain = x_domain, data=G%gridLonB(G%IsgB:G%IegB))
+    call mpp_write_meta(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", &
+                   cartesian='X', domain = x_domain, data=gridLonB(IsgB:IegB))
 
   if (use_layer) &
     call mpp_write_meta(unit, axis_layer, name="Layer", units=trim(GV%zAxisUnits), &
@@ -305,24 +332,30 @@ end subroutine create_file
 !! does not find the file, a new file is created.  It also sets up
 !! structures that describe this file and the variables that will
 !! later be written to this file.
-subroutine reopen_file(unit, filename, vars, novars, G, fields, threading, timeunit, GV)
+subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV)
   integer,               intent(out)   :: unit       !< unit id of an open file or -1 on a
                                                      !! nonwriting PE with single file output
   character(len=*),      intent(in)    :: filename   !< full path to the file to create
   type(vardesc),         intent(in)    :: vars(:)    !< structures describing fields written to filename
   integer,               intent(in)    :: novars     !< number of fields written to filename
-  type(ocean_grid_type), intent(in)    :: G          !< ocean grid structure
   type(fieldtype),       intent(inout) :: fields(:)  !< array of fieldtypes for each variable
   integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE
   real, optional,        intent(in)    :: timeunit   !< length, in seconds, of the units for time. The
                                                      !! default value is 86400.0, for 1 day.
+  type(ocean_grid_type),   optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
+                                                     !! is required if a new file uses any
+                                                     !! horizontal grid axes.
+  type(dyn_horgrid_type),  optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG
+                                                     !! is required if a new file uses any
+                                                     !! horizontal grid axes.
   type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is
-                                                     !! required if the new file uses any
+                                                     !! required if a new file uses any
                                                      !! vertical grid axes.
 
+  type(MOM_domain_type), pointer :: Domain => NULL()
   character(len=200) :: check_name, mesg
   integer :: length, ndim, nvar, natt, ntime, thread
-  logical :: exists
+  logical :: exists, one_file, domain_set
 
   thread = SINGLE_FILE
   if (PRESENT(threading)) thread = threading
@@ -335,12 +368,26 @@ subroutine reopen_file(unit, filename, vars, novars, G, fields, threading, timeu
   inquire(file=check_name,EXIST=exists)
 
   if (.not.exists) then
-    call create_file(unit, filename, vars, novars, fields, threading, timeunit, G=G, GV=GV)
+    call create_file(unit, filename, vars, novars, fields, threading, timeunit, &
+                     G=G, dG=dG, GV=GV)
   else
-    if ((thread == SINGLE_FILE) .or. .not.G%Domain%use_io_layout) then
+
+    domain_set = .false.
+    if (present(G)) then
+      domain_set = .true. ; Domain => G%Domain
+    elseif (present(dG)) then
+      domain_set = .true. ; Domain => dG%Domain
+    endif
+
+    one_file = .true.
+    if (domain_set) then
+      one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout)
+    endif
+
+    if (one_file) then
       call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, threading=thread)
     else
-      call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, domain=G%Domain%mpp_domain)
+      call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, domain=Domain%mpp_domain)
     endif
     if (unit < 0) return
 

From 7ea0255abf47f6cf5187a24a8a8dc7c20fb913ea Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Fri, 24 Jun 2016 17:30:38 -0400
Subject: [PATCH 3/5] +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                              |  50 +++++----
 src/core/MOM_transcribe_grid.F90              |   4 +-
 src/ice_shelf/MOM_ice_shelf.F90               |  44 +++++---
 .../MOM_fixed_initialization.F90              | 106 +++++++++---------
 src/initialization/MOM_grid_initialize.F90    |  37 +++---
 src/user/DOME2d_initialization.F90            |  12 +-
 src/user/DOME_initialization.F90              |  13 ++-
 src/user/ISOMIP_initialization.F90            |  14 +--
 src/user/Phillips_initialization.F90          |  19 ++--
 src/user/benchmark_initialization.F90         |  13 ++-
 src/user/seamount_initialization.F90          |  16 +--
 src/user/sloshing_initialization.F90          |  16 +--
 src/user/user_initialization.F90              |  16 +--
 13 files changed, 199 insertions(+), 161 deletions(-)

diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 54b46ef5b5..f888599b6e 100644
--- a/src/core/MOM.F90
+++ b/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,8 +1663,8 @@ 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)
 
@@ -1665,11 +1672,11 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
 
   ! 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)
diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90
index a3d6e75375..def4c1f606 100644
--- a/src/core/MOM_transcribe_grid.F90
+++ b/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)
diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90
index ef2b752cd9..14cde0c674 100644
--- a/src/ice_shelf/MOM_ice_shelf.F90
+++ b/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
diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90
index 5ebe25a844..b7e5aa3008 100644
--- a/src/initialization/MOM_fixed_initialization.F90
+++ b/src/initialization/MOM_fixed_initialization.F90
@@ -8,11 +8,12 @@ module MOM_fixed_initialization
 use MOM_coms, only : max_across_PEs
 use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
 use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
 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, isPointInCell
+! 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
@@ -44,7 +45,7 @@ module MOM_fixed_initialization
 !> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's
 !!   horizontal grid, bathymetry, and the Coriolis parameter.
 subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir)
-  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
+  type(dyn_horgrid_type),  intent(inout) :: G    !< The ocean's grid structure.
   type(param_file_type),   intent(in)    :: PF   !< A structure indicating the open file
                                                  !! to parse for model parameter values.
   logical,                 intent(in)    :: write_geom !< If true, write grid geometry files.
@@ -147,15 +148,11 @@ subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir)
 end subroutine MOM_initialize_fixed
 ! -----------------------------------------------------------------------------
 
-
+!> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter.
 subroutine MOM_initialize_rotation(f, G, PF)
-  type(ocean_grid_type),                        intent(in)  :: G
-  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f
-  type(param_file_type),                        intent(in)  :: PF
-! Arguments: f  - the Coriolis parameter in s-1. Intent out.
-!  (in)      G  - The ocean's grid structure.
-!  (in)      PF - A structure indicating the open file to parse for
-!                         model parameter values.
+  type(dyn_horgrid_type),                       intent(in)  :: G  !< The dynamic horizontal grid type
+  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f  !< The Coriolis parameter in s-1
+  type(param_file_type),                        intent(in)  :: PF !< Parameter file structure
 
 !   This subroutine makes the appropriate call to set up the Coriolis parameter.
 ! This is a separate subroutine so that it can be made public and shared with
@@ -185,9 +182,11 @@ end subroutine MOM_initialize_rotation
 
 !> Calculates the components of grad f (Coriolis parameter)
 subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G)
-  type(ocean_grid_type),            intent(inout) :: G !< Grid type
-  real, dimension(SZI_(G),SZJ_(G)), intent(out)   :: dF_dx !< x-component of grad f
-  real, dimension(SZI_(G),SZJ_(G)), intent(out)   :: dF_dy !< y-component of grad f
+  type(dyn_horgrid_type),             intent(inout) :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out)   :: dF_dx !< x-component of grad f
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out)   :: dF_dy !< y-component of grad f
   ! Local variables
   integer :: i,j
   real :: f1, f2
@@ -202,15 +201,13 @@ subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G)
   call pass_vector(dF_dx, dF_dy, G%Domain, stagger=AGRID)
 end subroutine MOM_calculate_grad_Coriolis
 
+!> MOM_initialize_topography makes the appropriate call to set up the bathymetry.
 subroutine MOM_initialize_topography(D, max_depth, G, PF)
-  type(ocean_grid_type),            intent(in)  :: G
-  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D
-  real,                             intent(out) :: max_depth
-  type(param_file_type),            intent(in)  :: PF
-! Arguments: D  - the bottom depth in m. Intent out.
-!  (in)      G  - The ocean's grid structure.
-!  (in)      PF - A structure indicating the open file to parse for
-!                         model parameter values.
+  type(dyn_horgrid_type),           intent(in)  :: G  !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                    intent(out) :: D  !< Ocean bottom depth in m
+  type(param_file_type),            intent(in)  :: PF !< Parameter file structure
+  real,                             intent(out) :: max_depth !< Maximum depth of model in m
 
 !  This subroutine makes the appropriate call to set up the bottom depth.
 !  This is a separate subroutine so that it can be made public and shared with
@@ -249,13 +246,13 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF)
     case ("bowl");      call initialize_topography_named(D, G, PF, config, max_depth)
     case ("halfpipe");  call initialize_topography_named(D, G, PF, config, max_depth)
     case ("DOME");      call DOME_initialize_topography(D, G, PF, max_depth)
-    case ("ISOMIP");      call ISOMIP_initialize_topography(D, G, PF, max_depth)
+    case ("ISOMIP");    call ISOMIP_initialize_topography(D, G, PF, max_depth)
     case ("benchmark"); call benchmark_initialize_topography(D, G, PF, max_depth)
     case ("DOME2D");    call DOME2d_initialize_topography(D, G, PF, max_depth)
     case ("sloshing");  call sloshing_initialize_topography(D, G, PF, max_depth)
     case ("seamount");  call seamount_initialize_topography(D, G, PF, max_depth)
-    case ("Phillips");  call Phillips_initialize_topography(D, G, PF)
-    case ("USER");      call user_initialize_topography(D, G, PF)
+    case ("Phillips");  call Phillips_initialize_topography(D, G, PF, max_depth)
+    case ("USER");      call user_initialize_topography(D, G, PF, max_depth)
     case default ;      call MOM_error(FATAL,"MOM_initialize_topography: "// &
       "Unrecognized topography setup '"//trim(config)//"'")
   end select
@@ -275,7 +272,7 @@ end subroutine MOM_initialize_topography
 
 ! -----------------------------------------------------------------------------
 function diagnoseMaximumDepth(D,G)
-  type(ocean_grid_type),            intent(in) :: G
+  type(dyn_horgrid_type),            intent(in) :: G
   real, dimension(SZI_(G),SZJ_(G)), intent(in) :: D
   real :: diagnoseMaximumDepth
   ! Local variables
@@ -292,8 +289,9 @@ end function diagnoseMaximumDepth
 
 !> Read gridded depths from file
 subroutine initialize_topography_from_file(D, G, param_file)
-  type(ocean_grid_type),            intent(in)  :: G !< Grid structure
-  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Bottom depth (positive, in m)
+  type(dyn_horgrid_type),           intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                    intent(out) :: D !< Ocean bottom depth in m
   type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
   ! Local variables
   character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
@@ -333,9 +331,11 @@ end subroutine initialize_topography_from_file
 
 !> Applies a list of topography overrides read from a netcdf file
 subroutine apply_topography_edits_from_file(D, G, param_file)
-  type(ocean_grid_type),            intent(in)    :: G !< Grid structure
-  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: D !< Bottom depth (positive, in m)
+  type(dyn_horgrid_type),           intent(in)    :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                    intent(inout) :: D !< Ocean bottom depth in m
   type(param_file_type),            intent(in)    :: param_file !< Parameter file structure
+
   ! Local variables
   character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
   character(len=40)  :: mod = "apply_topography_edits_from_file" ! This subroutine's name.
@@ -441,12 +441,16 @@ subroutine apply_topography_edits_from_file(D, G, param_file)
 end subroutine apply_topography_edits_from_file
 
 ! -----------------------------------------------------------------------------
+!> initialized the bathymetry based on one of several named idealized configurations
 subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth)
-  type(ocean_grid_type),            intent(in)  :: G
-  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D
-  type(param_file_type),            intent(in)  :: param_file
-  character(len=*),                 intent(in)  :: topog_config
-  real,                             intent(in)  :: max_depth
+  type(dyn_horgrid_type),           intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                    intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
+  character(len=*),                 intent(in)  :: topog_config !< The name of an idealized
+                                                              !! topographic configuration
+  real,                             intent(in)  :: max_depth  !< Maximum depth of model in m
+
 ! Arguments: D          - the bottom depth in m. Intent out.
 !  (in)      G          - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
@@ -553,11 +557,13 @@ end subroutine initialize_topography_named
 ! -----------------------------------------------------------------------------
 
 ! -----------------------------------------------------------------------------
+!> limit_topography ensures that  min_depth < D(x,y) < max_depth
 subroutine limit_topography(D, G, param_file, max_depth)
-  type(ocean_grid_type),            intent(in)    :: G
-  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: D
-  type(param_file_type),            intent(in)    :: param_file
-  real,                             intent(in)    :: max_depth
+  type(dyn_horgrid_type), intent(in)    :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                          intent(inout) :: D !< Ocean bottom depth in m
+  type(param_file_type),  intent(in)    :: param_file !< Parameter file structure
+  real,                   intent(in)    :: max_depth  !< Maximum depth of model in m
 ! Arguments: D          - the bottom depth in m. Intent in/out.
 !  (in)      G          - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
@@ -602,7 +608,7 @@ end subroutine limit_topography
 
 ! -----------------------------------------------------------------------------
 subroutine set_rotation_planetary(f, G, param_file)
-  type(ocean_grid_type),                        intent(in)  :: G
+  type(dyn_horgrid_type),                        intent(in)  :: G
   real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f
   type(param_file_type),                        intent(in)  :: param_file
 ! Arguments: f          - Coriolis parameter (vertical component) in s^-1
@@ -631,7 +637,7 @@ end subroutine set_rotation_planetary
 
 ! -----------------------------------------------------------------------------
 subroutine set_rotation_beta_plane(f, G, param_file)
-  type(ocean_grid_type),                        intent(in)  :: G
+  type(dyn_horgrid_type),                        intent(in)  :: G
   real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f
   type(param_file_type),                        intent(in)  :: param_file
 ! Arguments: f          - Coriolis parameter (vertical component) in s^-1
@@ -677,7 +683,7 @@ end subroutine set_rotation_beta_plane
 
 ! -----------------------------------------------------------------------------
 subroutine reset_face_lengths_named(G, param_file, name)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   type(param_file_type), intent(in)    :: param_file
   character(len=*),      intent(in)    :: name
 !   This subroutine sets the open face lengths at selected points to restrict
@@ -802,7 +808,7 @@ end subroutine reset_face_lengths_named
 
 ! -----------------------------------------------------------------------------
 subroutine reset_face_lengths_file(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   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.
@@ -870,7 +876,7 @@ end subroutine reset_face_lengths_file
 
 ! -----------------------------------------------------------------------------
 subroutine reset_face_lengths_list(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   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.
@@ -1129,7 +1135,7 @@ end subroutine read_face_length_list
 
 ! -----------------------------------------------------------------------------
 subroutine set_velocity_depth_max(G)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   ! This subroutine sets the 4 bottom depths at velocity points to be the
   ! maximum of the adjacent depths.
   integer :: i, j
@@ -1147,7 +1153,7 @@ end subroutine set_velocity_depth_max
 
 ! -----------------------------------------------------------------------------
 subroutine compute_global_grid_integrals(G)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   ! Subroutine to pre-compute global integrals of grid quantities for
   ! later use in reporting diagnostics
   integer :: i,j
@@ -1168,7 +1174,7 @@ end subroutine compute_global_grid_integrals
 
 ! -----------------------------------------------------------------------------
 subroutine set_velocity_depth_min(G)
-  type(ocean_grid_type), intent(inout) :: G
+  type(dyn_horgrid_type), intent(inout) :: G
   ! This subroutine sets the 4 bottom depths at velocity points to be the
   ! minimum of the adjacent depths.
   integer :: i, j
@@ -1186,9 +1192,9 @@ end subroutine set_velocity_depth_min
 
 ! -----------------------------------------------------------------------------
 subroutine write_ocean_geometry_file(G, param_file, directory)
-  type(ocean_grid_type), intent(inout) :: G
-  type(param_file_type), intent(in)    :: param_file
-  character(len=*),      intent(in)    :: directory
+  type(dyn_horgrid_type), intent(inout) :: G
+  type(param_file_type),  intent(in)    :: param_file
+  character(len=*),       intent(in)    :: directory
 !   This subroutine writes out a file containing all of the ocean geometry
 ! and grid data uses by the MOM ocean model.
 ! Arguments: G - The ocean's grid structure.  Effectively intent in.
@@ -1269,7 +1275,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory)
   if (multiple_files) file_threading = MULTIPLE
 
   call create_file(unit, trim(filepath), vars, nFlds_used, fields, &
-                   file_threading, G=G)
+                   file_threading, dG=G)
 
   do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLatBu(I,J); enddo; enddo
   call write_field(unit, fields(1), G%Domain%mpp_domain, out_q)
diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90
index 03332babf1..46f4dc9c18 100644
--- a/src/initialization/MOM_grid_initialize.F90
+++ b/src/initialization/MOM_grid_initialize.F90
@@ -72,10 +72,11 @@ module MOM_grid_initialize
 use MOM_domains, only : To_North, To_South, To_East, To_West
 use MOM_domains, only : MOM_define_domain, MOM_define_IO_domain
 use MOM_domains, only : MOM_domain_type
+use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid
 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_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
 
@@ -107,7 +108,7 @@ module MOM_grid_initialize
 !!   grid.  The bathymetry, land-sea mask and any restricted channel widths are
 !!   not know yet, so these are set later.
 subroutine set_grid_metrics(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G           !< The horizontal grid structure
+  type(dyn_horgrid_type), intent(inout) :: G          !< The dynamic horizontal grid type
   type(param_file_type), intent(in)    :: param_file  !< Parameter file structure
 ! Arguments:
 !  (inout)   G - The ocean's grid structure.
@@ -155,7 +156,8 @@ subroutine set_grid_metrics(G, param_file)
 
 ! Calculate derived metrics (i.e. reciprocals and products)
   call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90")
-  call set_derived_metrics(G)
+!  call set_derived_metrics(G)
+  call set_derived_dyn_horgrid(G)
   call callTree_leave("set_derived_metrics()")
 
   if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics',G)
@@ -169,7 +171,7 @@ end subroutine set_grid_metrics
 !!   debugging.
 subroutine grid_metrics_chksum(parent, G)
   character(len=*),      intent(in) :: parent  !< A string identifying the caller
-  type(ocean_grid_type), intent(in) :: G       !< The horizontal grid structure
+  type(dyn_horgrid_type), intent(in) :: G      !< The dynamic horizontal grid type
 
   real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH
   real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ
@@ -271,9 +273,10 @@ end subroutine grid_metrics_chksum
 
 ! ------------------------------------------------------------------------------
 
-subroutine set_grid_metrics_from_mosaic(G,param_file)
-  type(ocean_grid_type), intent(inout) :: G
-  type(param_file_type), intent(in)    :: param_file
+!>  et_grid_metrics_from_mosaic sets the grid metrics from a mosaic file.
+subroutine set_grid_metrics_from_mosaic(G, param_file)
+  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
+  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure
 !   This subroutine sets the grid metrics from a mosaic file.
 !  
 ! Arguments:
@@ -511,8 +514,9 @@ end subroutine set_grid_metrics_from_mosaic
 ! ------------------------------------------------------------------------------
 
 subroutine set_grid_metrics_cartesian(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G
-  type(param_file_type), intent(in)    :: param_file
+  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
+  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure
+
 ! Arguments:
 !  (inout)   G - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
@@ -646,8 +650,9 @@ end subroutine set_grid_metrics_cartesian
 ! ------------------------------------------------------------------------------
 
 subroutine set_grid_metrics_spherical(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G
-  type(param_file_type), intent(in)    :: param_file
+  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
+  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure
+
 ! Arguments:
 !  (inout)   G - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
@@ -785,8 +790,9 @@ end subroutine set_grid_metrics_spherical
 ! ------------------------------------------------------------------------------
 
 subroutine set_grid_metrics_mercator(G, param_file)
-  type(ocean_grid_type), intent(inout) :: G
-  type(param_file_type), intent(in)    :: param_file
+  type(dyn_horgrid_type), intent(inout) :: G           !< The dynamic horizontal grid type
+  type(param_file_type), intent(in)     :: param_file  !< Parameter file structure
+
 ! Arguments:
 !  (inout)   G - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
@@ -1327,8 +1333,9 @@ end function Adcroft_reciprocal
 !> initialize_masks initializes the grid masks and any metrics that come
 !!    with masks already applied.
 subroutine initialize_masks(G, PF)
-  type(ocean_grid_type), intent(inout) :: G  !< The ocean's grid structure
-  type(param_file_type), intent(in)    :: PF !< The param_file handle
+  type(dyn_horgrid_type), intent(inout) :: G   !< The dynamic horizontal grid type
+  type(param_file_type), intent(in)     :: PF  !< Parameter file structure
+
 ! Arguments:
 !  (inout)   G - The ocean's grid structure.
 !  (in)      PF - A structure indicating the open file to parse for
diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90
index 0e96df8f28..36c3996bc4 100644
--- a/src/user/DOME2d_initialization.F90
+++ b/src/user/DOME2d_initialization.F90
@@ -1,6 +1,7 @@
 module DOME2d_initialization
 
 use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
@@ -33,10 +34,11 @@ module DOME2d_initialization
 !> Initialize topography with a shelf and slope in a 2D domain
 subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth )
   ! Arguments
-  type(ocean_grid_type),            intent(in)  :: G !< Ocean grid structure
-  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Ocean bottom depth
-  type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
-  real,                             intent(in)  :: max_depth !< Maximum depth of model
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
   ! Local variables
   integer :: i, j
   real    :: x, bay_depth, l1, l2
@@ -82,8 +84,8 @@ end subroutine DOME2d_initialize_topography
 !> Initialize thicknesses according to coordinate mode
 subroutine DOME2d_initialize_thickness ( h, G, GV, param_file )
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure
-  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Layer thicknesses
   type(verticalGrid_type),                   intent(in)  :: GV !< Vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< Layer thicknesses
   type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
   ! Local variables
   real :: e0(SZK_(G))     ! The resting interface heights, in m, usually !
diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90
index 829912a338..2a8c22c04f 100644
--- a/src/user/DOME_initialization.F90
+++ b/src/user/DOME_initialization.F90
@@ -20,6 +20,7 @@ module DOME_initialization
 !***********************************************************************
 
 use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
@@ -43,11 +44,11 @@ module DOME_initialization
 ! -----------------------------------------------------------------------------
 !> This subroutine sets up the DOME topography
 subroutine DOME_initialize_topography(D, G, param_file, max_depth)
-  type(ocean_grid_type), intent(in)           :: G !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m.
-  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
-                                                  !! to parse for model parameter values.
-  real,                  intent(in)               :: max_depth !< Maximum depth.
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
 
   real :: min_depth ! The minimum and maximum depths in m.
 ! This include declares and sets the variable "version".
@@ -90,7 +91,7 @@ end subroutine DOME_initialize_topography
 subroutine DOME_initialize_thickness(h, G, GV, param_file)
   type(ocean_grid_type),   intent(in) :: G  !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being initialized.
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being initialized.
   type(param_file_type),   intent(in) :: param_file !< A structure indicating the open file
                                                     !! to parse for model parameter values.
 
diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90
index d38d9db420..c47671a2cc 100644
--- a/src/user/ISOMIP_initialization.F90
+++ b/src/user/ISOMIP_initialization.F90
@@ -20,6 +20,7 @@ module ISOMIP_initialization
 !***********************************************************************
 
 use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe, WARNING
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
@@ -58,12 +59,11 @@ module ISOMIP_initialization
 
 !> Initialization of topography
 subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth)
-  type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D        !< The bottom depth in m.
-  type(param_file_type), intent(in)           :: param_file !< A structure indicating the
-                                                            !! open file to parse for model
-                                                            !! parameter values.
-  real,                  intent(in)           :: max_depth  !< Maximum depth.
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
 
 ! This subroutine sets up the ISOMIP topography
   real :: min_depth ! The minimum and maximum depths in m.
@@ -146,7 +146,7 @@ end subroutine ISOMIP_initialize_topography
 subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv )
   type(ocean_grid_type), intent(in) :: G                !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV             !< The ocean's vertical grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being
                                                         !! initialized.
   type(param_file_type), intent(in) :: param_file       !< A structure indicating the
                                                         !! open file to parse for model
diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90
index 492fd2e8ee..16970cb1a0 100644
--- a/src/user/Phillips_initialization.F90
+++ b/src/user/Phillips_initialization.F90
@@ -20,6 +20,7 @@ module Phillips_initialization
 !***********************************************************************
 
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
 use MOM_grid, only : ocean_grid_type
@@ -51,9 +52,9 @@ module Phillips_initialization
 !> Initialize thickness field.
 subroutine Phillips_initialize_thickness(h, G, GV, param_file)
   type(ocean_grid_type),   intent(in) :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is
-                                                    !! being initialized.
   type(verticalGrid_type), intent(in) :: GV         !< The ocean's vertical grid structure.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(out) :: h         !< The thickness that is being initialized.
   type(param_file_type),   intent(in) :: param_file !< A structure indicating the
                                                     !! open file to parse for model
                                                     !! parameter values.
@@ -275,12 +276,12 @@ function sech(x)
 end function sech
 
 !> Initialize topography.
-subroutine Phillips_initialize_topography(D, G, param_file)
-  type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D        !< The bottom depth in m.
-  type(param_file_type), intent(in)           :: param_file !< A structure indicating the
-                                                            !! open file to parse for model
-                                                            !! parameter values.
+subroutine Phillips_initialize_topography(D, G, param_file, max_depth)
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
 
   real :: PI, Htop, Wtop, Ltop, offset, dist, &
           x1, x2, x3, x4, y1, y2
@@ -316,7 +317,7 @@ subroutine Phillips_initialize_topography(D, G, param_file)
        D(i,j) = 2.0/3.0*Htop*sin(PI*(G%geoLonT(i,j)-x3)/(x4-x3))**2 &
                     *sin(PI*(G%geoLatT(i,j)-y1)/(y2-y1))**2
      end if
-     D(i,j)=G%max_depth-D(i,j)
+     D(i,j)=max_depth-D(i,j)
   enddo; enddo
 
 end subroutine Phillips_initialize_topography
diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90
index d003773f08..e7953f86d1 100644
--- a/src/user/benchmark_initialization.F90
+++ b/src/user/benchmark_initialization.F90
@@ -20,6 +20,7 @@ module benchmark_initialization
 !***********************************************************************
 
 use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
@@ -43,11 +44,11 @@ module benchmark_initialization
 ! -----------------------------------------------------------------------------
 !> This subroutine sets up the benchmark test case topography.
 subroutine benchmark_initialize_topography(D, G, param_file, max_depth)
-  type(ocean_grid_type),           intent(in) :: G          !< The ocean's grid structure.
-  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D        !< the bottom depth in m.
-  type(param_file_type),           intent(in) :: param_file !< A structure indicating the open
-                                                   !! file to parse for model parameter values.
-  real,                            intent(in) :: max_depth  !< The Maximum depth.
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
 
 ! This subroutine sets up the benchmark test case topography
   real :: min_depth            ! The minimum and maximum depths in m.
@@ -94,7 +95,7 @@ end subroutine benchmark_initialize_topography
 subroutine benchmark_initialize_thickness(h, G, GV, param_file, eqn_of_state, P_ref)
   type(ocean_grid_type),   intent(in) :: G                    !< The ocean's grid structure.
   type(verticalGrid_type), intent(in) :: GV                   !< The ocean's vertical grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being
                                                               !! initialized.
   type(param_file_type),   intent(in) :: param_file           !< A structure indicating the open
                                                               !! file to parse for model
diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90
index 7cb9ba4612..79768bca5e 100644
--- a/src/user/seamount_initialization.F90
+++ b/src/user/seamount_initialization.F90
@@ -20,6 +20,7 @@ module seamount_initialization
 !***********************************************************************
 
 use MOM_domains, only : sum_across_PEs
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
 use MOM_file_parser, only : get_param, param_file_type
 use MOM_get_input, only : directories
@@ -57,12 +58,11 @@ module seamount_initialization
 !> Initialization of topography.
 subroutine seamount_initialize_topography ( D, G, param_file, max_depth )
   ! Arguments 
-  type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D        !< The bottom depth in m.
-  type(param_file_type), intent(in)           :: param_file !< A structure indicating the
-                                                            !! open file to parse for model
-                                                            !! parameter values.
-  real,                  intent(in)           :: max_depth  !< Maximum depth.
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
   
   ! Local variables 
   integer   :: i, j
@@ -99,9 +99,9 @@ end subroutine seamount_initialize_topography
 !! This subroutine initializes the layer thicknesses to be uniform.
 subroutine seamount_initialize_thickness ( h, G, GV, param_file )
   type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being
-                                                            !! initialized.
   type(verticalGrid_type), intent(in)         :: GV         !< The ocean's vertical grid structure.
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being
+                                                            !! initialized.
   type(param_file_type), intent(in)           :: param_file !< A structure indicating the
                                                             !! open file to parse for model
                                                             !! parameter values.
diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90
index 791c1ddcd6..a24c8f845a 100644
--- a/src/user/sloshing_initialization.F90
+++ b/src/user/sloshing_initialization.F90
@@ -20,6 +20,7 @@ module sloshing_initialization
 !***********************************************************************
 
 use MOM_domains, only : sum_across_PEs
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
 use MOM_file_parser, only : get_param, param_file_type
 use MOM_get_input, only : directories
@@ -53,12 +54,11 @@ module sloshing_initialization
 
 !> Initialization of topography.
 subroutine sloshing_initialize_topography ( D, G, param_file, max_depth )
-  type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D        !< The bottom depth in m.
-  type(param_file_type), intent(in)           :: param_file !< A structure indicating the
-                                                            !! open file to parse for model
-                                                            !! parameter values.
-  real,                  intent(in)           :: max_depth  !< Maximum depth.
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
   
   ! Local variables 
   integer   :: i, j
@@ -84,9 +84,9 @@ end subroutine sloshing_initialize_topography
 !! left and minimum value on the right. This sets off a regular sloshing motion.
 subroutine sloshing_initialize_thickness ( h, G, GV, param_file )
   type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being
-                                                            !! initialized.
   type(verticalGrid_type), intent(in)         :: GV         !< The ocean's vertical grid structure.
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being
+                                                            !! initialized.
   type(param_file_type), intent(in)           :: param_file !< A structure indicating the
                                                             !! open file to parse for model
                                                             !! parameter values.
diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90
index 3a01978e67..a96b820c65 100644
--- a/src/user/user_initialization.F90
+++ b/src/user/user_initialization.F90
@@ -20,6 +20,7 @@ module user_initialization
 !***********************************************************************
 
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_get_input, only : directories
 use MOM_grid, only : ocean_grid_type
@@ -69,12 +70,13 @@ subroutine USER_set_coord(Rlay, g_prime, GV, param_file, eqn_of_state)
 end subroutine USER_set_coord
 
 !> Initialize topography.
-subroutine USER_initialize_topography(D, G, param_file)
-  type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D        !< The bottom depth in m.
-  type(param_file_type), intent(in)           :: param_file !< A structure indicating the
-                                                            !! open file to parse for model
-                                                            !! parameter values.
+subroutine USER_initialize_topography(D, G, param_file, max_depth)
+  type(dyn_horgrid_type),             intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                      intent(out) :: D !< Ocean bottom depth in m
+  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
+  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m
+
   call MOM_error(FATAL, &
    "USER_initialization.F90, USER_initialize_topography: " // &
    "Unmodified user routine called - you must edit the routine to use it")
@@ -88,7 +90,7 @@ end subroutine USER_initialize_topography
 !> initialize thicknesses.
 subroutine USER_initialize_thickness(h, G, param_file, T)
   type(ocean_grid_type), intent(in)           :: G          !< The ocean's grid structure.
-  real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being
+  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h !< The thicknesses being
                                                             !! initialized.
   type(param_file_type), intent(in)           :: param_file !< A structure indicating the
                                                             !! open file to parse for model

From c88a8520cde68508d9b5b72c1484374e169c49a6 Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Fri, 24 Jun 2016 18:11:56 -0400
Subject: [PATCH 4/5] Eliminated the use of SZD... macros in EOS code

  The various SZDI_, etc., macros always resolve to the same thing, so in the
interest of code clarity, these macros have been eliminated throughout the MOM6
equation_of_state code.  All answers are bitwise identical.
---
 src/equation_of_state/MOM_EOS.F90        | 151 +++++++++++++----------
 src/equation_of_state/MOM_EOS_Wright.F90 |  44 ++++---
 src/equation_of_state/MOM_EOS_linear.F90 |  42 ++++---
 3 files changed, 140 insertions(+), 97 deletions(-)

diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90
index e05d868c88..20f72c8364 100644
--- a/src/equation_of_state/MOM_EOS.F90
+++ b/src/equation_of_state/MOM_EOS.F90
@@ -311,34 +311,34 @@ end subroutine calculate_2_densities
 subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, &
                                dza, intp_dza, intx_dza, inty_dza, halo_size)
   !> The horizontal index structure
-    type(hor_index_type),                  intent(in)  :: HI
+    type(hor_index_type),                          intent(in)  :: HI
   !> Potential temperature referenced to the surface (degC)
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: T
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  intent(in)  :: T
   !> Salinity (PSU)
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: S
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  intent(in)  :: S
   !> Pressure at the top of the layer in Pa.
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: p_t
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  intent(in)  :: p_t
   !> Pressure at the bottom of the layer in Pa.
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: p_b
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  intent(in)  :: p_b
   !> A mean specific volume that is subtracted out to reduce the magnitude of
   !! each of the integrals, m3 kg-1. The calculation is mathematically identical
   !! with different values of alpha_ref, but this reduces the effects of roundoff.
-    real,                            intent(in)  :: alpha_ref
+    real,                                          intent(in)  :: alpha_ref
   !> Equation of state structure
-    type(EOS_type),                  pointer     :: EOS
+    type(EOS_type),                                pointer     :: EOS
   !> The change in the geopotential anomaly across the layer, in m2 s-2.
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(out) :: dza
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  intent(out) :: dza
   !> The integral in pressure through the layer of the geopotential anomaly
   !! relative to the anomaly at the bottom of the layer, in Pa m2 s-2.
-    real, dimension(SZDI_(HI),SZDJ_(HI)),  optional, intent(out) :: intp_dza
+    real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  optional, intent(out) :: intp_dza
   !> The integral in x of the difference between the geopotential anomaly at the
   !! top and bottom of the layer divided by the x grid spacing, in m2 s-2.
-    real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza
+    real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), optional, intent(out) :: intx_dza
   !> The integral in y of the difference between the geopotential anomaly at the
   !! top and bottom of the layer divided by the y grid spacing, in m2 s-2.
-    real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza
+    real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), optional, intent(out) :: inty_dza
   !> The width of halo points on which to calculate dza.
-    integer,                         optional, intent(in)  :: halo_size
+    integer,                             optional, intent(in)  :: halo_size
 
   if (.not.associated(EOS)) call MOM_error(FATAL, &
     "int_specific_vol_dp called with an unassociated EOS_type EOS.")
@@ -374,13 +374,13 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, &
   !> Ocean horizontal index structures for the output arrays
     type(hor_index_type),            intent(in)  :: HIO
   !> Potential temperature referenced to the surface (degC)
-    real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T
+    real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in)  :: T
   !> Salinity (PSU)
-    real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: S
+    real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in)  :: S
   !> Height at the top of the layer in m.
-    real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: z_t
+    real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in)  :: z_t
   !> Height at the bottom of the layer in m.
-    real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: z_b
+    real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in)  :: z_b
   !> A mean density, in kg m-3, that is subtracted out to reduce the magnitude
   !! of each of the integrals. (The pressure is calculated as p~=-z*rho_0*G_e.)
     real,                            intent(in)  :: rho_ref
@@ -392,16 +392,16 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, &
   !> Equation of state structure
     type(EOS_type),                  pointer     :: EOS
   !> The change in the pressure anomaly across the layer, in Pa.
-    real, dimension(SZDI_(HIO),SZDJ_(HIO)),            intent(out) :: dpa
+    real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), intent(out) :: dpa
   !> The integral through the thickness of the layer of the pressure anomaly
   !! relative to the anomaly at the top of the layer, in Pa m.
-    real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
+    real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed),  optional, intent(out) :: intz_dpa
   !> The integral in x of the difference between the pressure anomaly at the
   !! top and bottom of the layer divided by the x grid spacing, in Pa.
-    real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa
+    real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), optional, intent(out) :: intx_dpa
   !> The integral in y of the difference between the pressure anomaly at the
   !! top and bottom of the layer divided by the y grid spacing, in Pa.
-    real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
+    real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), optional, intent(out) :: inty_dpa
 
   if (.not.associated(EOS)) call MOM_error(FATAL, &
     "int_density_dz called with an unassociated EOS_type EOS.")
@@ -561,14 +561,19 @@ end subroutine EOS_use_linear
 
 subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
                                   EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
-  type(hor_index_type),                    intent(in)  :: HII, HIO
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T, S, z_t, z_b
-  real,                                    intent(in)  :: rho_ref, rho_0, G_e
-  type(EOS_type),                          pointer     :: EOS !< Equation of state structure
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  intent(out) :: dpa
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
-  real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
+  type(hor_index_type), intent(in)  :: HII, HIO
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
+                        intent(in)  :: T, S, z_t, z_b
+  real,                 intent(in)  :: rho_ref, rho_0, G_e
+  type(EOS_type),       pointer     :: EOS !< Equation of state structure
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                        intent(out) :: dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intz_dpa
+  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intx_dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
+              optional, intent(out) :: inty_dpa
 !   This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.  The one
@@ -897,16 +902,21 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, &
                                        rho_0, G_e, H_subroundoff, bathyT, HII, HIO, EOS, dpa, &
                                        intz_dpa, intx_dpa, inty_dpa, &
                                        useMassWghtInterp)
-  type(hor_index_type),                    intent(in)  :: HII, HIO
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T_t, T_b, S_t, S_b, z_t, z_b
-  real,                                    intent(in)  :: rho_ref, rho_0, G_e, H_subroundoff
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: bathyT
-  type(EOS_type), pointer                              :: EOS !< Equation of state structure
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  intent(out) :: dpa
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
-  real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
-  logical,                                 optional, intent(in)  :: useMassWghtInterp
+  type(hor_index_type),   intent(in)  :: HII, HIO
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
+                          intent(in)  :: T_t, T_b, S_t, S_b, z_t, z_b
+  real,                   intent(in)  :: rho_ref, rho_0, G_e, H_subroundoff
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed),  intent(in)  :: bathyT
+  type(EOS_type),         pointer     :: EOS !< Equation of state structure
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                          intent(out) :: dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                optional, intent(out) :: intz_dpa
+  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
+                optional, intent(out) :: intx_dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
+                optional, intent(out) :: inty_dpa
+  logical,      optional, intent(in)  :: useMassWghtInterp
 ! This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.  The one
@@ -1294,16 +1304,19 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, &
                                        z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
                                        EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
   
-  type(hor_index_type),                    intent(in)  :: HII, HIO
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T, T_t, T_b, S, S_t, S_b, &
-                                                        z_t, z_b
+  type(hor_index_type), intent(in)  :: HII, HIO
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
+                        intent(in)  :: T, T_t, T_b, S, S_t, S_b, z_t, z_b
   real,                                    intent(in)  :: rho_ref, rho_0, G_e
-  type(EOS_type), pointer                            :: EOS !< Equation of state structure
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  intent(out) :: dpa
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
-
-  real, dimension(SZDIB_(HIO),SZDJB_(HIO)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
+  type(EOS_type),       pointer     :: EOS !< Equation of state structure
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                        intent(out) :: dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intz_dpa
+  real, dimension(HIO%IsdB:HIO%IedB,HIO%JsdB:HIO%JedB), &
+              optional, intent(out) :: intx_dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
+              optional, intent(out) :: inty_dpa
 ! This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.  The one
@@ -1543,14 +1556,19 @@ end subroutine int_density_dz_generic_ppm
 subroutine int_density_dz_generic_plm_analytic (T_t, T_b, S_t, S_b, z_t, &
             z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, intz_dpa, intx_dpa, inty_dpa)
   
-  type(hor_index_type),                  intent(in)  :: HI
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: T_t, T_b, S_t, S_b, z_t, z_b
-  real,                            intent(in)  :: rho_ref, rho_0, G_e
-  type(EOS_type), pointer                      :: EOS !< Equation of state structure
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(out) :: dpa
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  optional, intent(out) :: intz_dpa
-  real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dpa
+  type(hor_index_type), intent(in)  :: HI
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(in)  :: T_t, T_b, S_t, S_b, z_t, z_b
+  real,                 intent(in)  :: rho_ref, rho_0, G_e
+  type(EOS_type),       pointer     :: EOS !< Equation of state structure
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(out) :: dpa
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(out) :: intz_dpa
+  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
+              optional, intent(out) :: intx_dpa
+  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
+              optional, intent(out) :: inty_dpa
 ! This subroutine calculates (by numerical quadrature) integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.  The one
@@ -1949,15 +1967,20 @@ end subroutine evaluate_shape_quadratic
 
 subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, &
                                    dza, intp_dza, intx_dza, inty_dza, halo_size)
-  type(hor_index_type),                  intent(in)  :: HI
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: T, S, p_t, p_b
-  real,                                  intent(in)  :: alpha_ref
-  type(EOS_type),                        pointer     :: EOS !< Equation of state structure
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(out) :: dza
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  optional, intent(out) :: intp_dza
-  real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza
-  real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza
-  integer,                               optional, intent(in)  :: halo_size
+  type(hor_index_type), intent(in)  :: HI
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(in)  :: T, S, p_t, p_b
+  real,                 intent(in)  :: alpha_ref
+  type(EOS_type),       pointer     :: EOS !< Equation of state structure
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(out) :: dza
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(out) :: intp_dza
+  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
+              optional, intent(out) :: intx_dza
+  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
+              optional, intent(out) :: inty_dza
+  integer,    optional, intent(in)  :: halo_size
 !   This subroutine calculates analytical and nearly-analytical integrals in
 ! pressure across layers of geopotential anomalies, which are required for
 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90
index 97e0b90e8b..778dff90ae 100644
--- a/src/equation_of_state/MOM_EOS_Wright.F90
+++ b/src/equation_of_state/MOM_EOS_Wright.F90
@@ -259,13 +259,18 @@ end subroutine calculate_2_densities_wright
 
 subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, &
                                  dpa, intz_dpa, intx_dpa, inty_dpa)
-  type(hor_index_type),                    intent(in)  :: HII, HIO
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T, S, z_t, z_b
-  real,                                    intent(in)  :: rho_ref, rho_0, G_e
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  intent(out) :: dpa
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
-  real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
+  type(hor_index_type), intent(in)  :: HII, HIO
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
+                        intent(in)  :: T, S, z_t, z_b
+  real,                 intent(in)  :: rho_ref, rho_0, G_e
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                        intent(out) :: dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intz_dpa
+  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intx_dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), &
+              optional, intent(out) :: inty_dpa
 !   This subroutine calculates analytical and nearly-analytical integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.
@@ -293,7 +298,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO,
 !                       pressure anomaly at the top and bottom of the layer
 !                       divided by the y grid spacing, in Pa.
 
-  real, dimension(SZDI_(HII),SZDJ_(HII)) :: al0_2d, p0_2d, lambda_2d
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed) :: al0_2d, p0_2d, lambda_2d
   real :: al0, p0, lambda
   real :: eps, eps2, rho_anom, rem
   real :: w_left, w_right, intz(5)
@@ -393,14 +398,19 @@ end subroutine int_density_dz_wright
 
 subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, &
                                   intp_dza, intx_dza, inty_dza, halo_size)
-  type(hor_index_type),                  intent(in)  :: HI
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: T, S, p_t, p_b
-  real,                                  intent(in)  :: alpha_ref
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(out) :: dza
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  optional, intent(out) :: intp_dza
-  real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza
-  real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza
-  integer,                               optional, intent(in)  :: halo_size
+  type(hor_index_type), intent(in)  :: HI
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(in)  :: T, S, p_t, p_b
+  real,                 intent(in)  :: alpha_ref
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(out) :: dza
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(out) :: intp_dza
+  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
+              optional, intent(out) :: intx_dza
+  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
+              optional, intent(out) :: inty_dza
+  integer,    optional, intent(in)  :: halo_size
 !   This subroutine calculates analytical and nearly-analytical integrals in
 ! pressure across layers of geopotential anomalies, which are required for
 ! calculating the finite-volume form pressure accelerations in a non-Boussinesq
@@ -431,7 +441,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, &
 !                       divided by the y grid spacing, in m2 s-2.
 !  (in,opt)  halo_size - The width of halo points on which to calculate dza.
 
-  real, dimension(SZDI_(HI),SZDJ_(HI)) :: al0_2d, p0_2d, lambda_2d
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d
   real :: al0, p0, lambda
   real :: alpha_anom, dp, p_ave
   real :: rem, eps, eps2
diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90
index 090e12f5ac..96e2468ae5 100644
--- a/src/equation_of_state/MOM_EOS_linear.F90
+++ b/src/equation_of_state/MOM_EOS_linear.F90
@@ -213,14 +213,19 @@ end subroutine calculate_2_densities_linear
 
 subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HII, HIO, &
                  Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa)
-  type(hor_index_type),                    intent(in)  :: HII, HIO
-  real, dimension(SZDI_(HII),SZDJ_(HII)),  intent(in)  :: T, S, z_t, z_b
-  real,                                    intent(in)  :: rho_ref, rho_0_pres, G_e
-  real,                                    intent(in)  :: Rho_T0_S0, dRho_dT, dRho_dS
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  intent(out) :: dpa
-  real, dimension(SZDI_(HIO),SZDJ_(HIO)),  optional, intent(out) :: intz_dpa
-  real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa
-  real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa
+  type(hor_index_type), intent(in)  :: HII, HIO
+  real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), &
+                        intent(in)  :: T, S, z_t, z_b
+  real,                 intent(in)  :: rho_ref, rho_0_pres, G_e
+  real,                 intent(in)  :: Rho_T0_S0, dRho_dT, dRho_dS
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+                        intent(out) :: dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), &
+              optional, intent(out) :: intz_dpa
+  real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed),  &
+              optional, intent(out) :: intx_dpa
+  real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB),  &
+              optional, intent(out) :: inty_dpa
 !   This subroutine calculates analytical and nearly-analytical integrals of
 ! pressure anomalies across layers, which are required for calculating the
 ! finite-volume form pressure accelerations in a Boussinesq model.
@@ -292,14 +297,19 @@ end subroutine int_density_dz_linear
 
 subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, &
                dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size)
-  type(hor_index_type),                  intent(in)  :: HI
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(in)  :: T, S, p_t, p_b
-  real,             intent(in)  :: alpha_ref
-  real,             intent(in)  :: Rho_T0_S0, dRho_dT, dRho_dS
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  intent(out) :: dza
-  real, dimension(SZDI_(HI),SZDJ_(HI)),  optional, intent(out) :: intp_dza
-  real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza
-  real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza
+  type(hor_index_type), intent(in)  :: HI
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed),  &
+                        intent(in)  :: T, S, p_t, p_b
+  real,                 intent(in)  :: alpha_ref
+  real,                 intent(in)  :: Rho_T0_S0, dRho_dT, dRho_dS
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+                        intent(out) :: dza
+  real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), &
+              optional, intent(out) :: intp_dza
+  real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), &
+              optional, intent(out) :: intx_dza
+  real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), &
+              optional, intent(out) :: inty_dza
   integer,                         optional, intent(in)  :: halo_size
 !   This subroutine calculates analytical and nearly-analytical integrals in
 ! pressure across layers of geopotential anomalies, which are required for

From 5d7504ea91228edf9845921a262870ad095a1c8b Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Mon, 27 Jun 2016 19:00:35 -0400
Subject: [PATCH 5/5] +Replaced ML_USE_OMEGA with ML_OMEGA_FRAC

  Replaced the logical test to toggle between using the local value of f or
omega in the rotational scaling for the decay of turbulence.  ML_USE_OMEGA=true
is replicated by ML_OMEGA_FRAC=1.0, and ML_USE_OMEGA=false (the default) is the
same as ML_OMEGA_FRAC=0.0 (the default), but now intermediate values can be
used as well.  Due to clever design, all test cases are bitwise identical, but
run-time parameter names and MOM_parameter_doc files have changed.
---
 .../vertical/MOM_bulk_mixed_layer.F90         | 30 +++++++++++----
 .../vertical/MOM_energetic_PBL.F90            | 27 +++++++++----
 .../vertical/MOM_set_diffusivity.F90          | 25 +++++++++---
 .../vertical/MOM_set_viscosity.F90            | 38 ++++++++++++++-----
 4 files changed, 90 insertions(+), 30 deletions(-)

diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
index 9a813d95f6..815aebec9a 100644
--- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
+++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
@@ -118,9 +118,9 @@ module MOM_bulk_mixed_layer
   integer :: ML_presort_nz_conv_adj ! If ML_resort is true, do convective
                              ! adjustment on this many layers (starting from the
                              ! top) before sorting the remaining layers.
-  logical :: use_omega       !   If true, use the absolute rotation rate instead
-                             ! of the vertical component of rotation when
-                             ! setting the decay scale for turbulence.
+  real    :: omega_frac      !   When setting the decay scale for turbulence, use
+                             ! this fraction of the absolute rotation rate blended
+                             ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   logical :: correct_absorption ! If true, the depth at which penetrating
                              ! shortwave radiation is absorbed is corrected by
                              ! moving some of the heating upward in the water
@@ -1368,7 +1368,7 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA,
   is = G%isc ; ie = G%iec ; nz = G%ke
   diag_wt = dt * Idt_diag
 
-  if (CS%use_omega) absf = 2.0*CS%omega
+  if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega
   do i=is,ie
     U_Star = fluxes%ustar(i,j)
     if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
@@ -1378,9 +1378,12 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA,
     endif
 
     if (U_Star < CS%ustar_min) U_Star = CS%ustar_min
-    if (.not.CS%use_omega) &
+    if (CS%omega_frac < 1.0) then
       absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
                    (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J))))
+      if (CS%omega_frac > 0.0) &
+        absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
+    endif
     absf_Ustar = absf / U_Star
     Idecay_len_TKE(i) = (absf_Ustar * CS%TKE_decay) * GV%H_to_m
 
@@ -3355,8 +3358,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS)
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_mixed_layer"  ! This module's name.
+  real :: omega_frac_dflt
   integer :: isd, ied, jsd, jed
-  logical :: use_temperature
+  logical :: use_temperature, use_omega
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
   if (associated(CS)) then
@@ -3446,10 +3450,20 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS)
   call get_param(param_file, mod, "OMEGA",CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)
-  call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+  call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+  omega_frac_dflt = 0.0
+  if (use_omega) then
+    call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+    omega_frac_dflt = 1.0
+  endif
+  call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                 "When setting the decay scale for turbulence, use this \n"//&
+                 "fraction of the absolute rotation rate blended with the \n"//&
+                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                 units="nondim", default=omega_frac_dflt)
   call get_param(param_file, mod, "ML_RESORT", CS%ML_resort, &
                  "If true, resort the topmost layers by potential density \n"//&
                  "before the mixed layer calculations.", default=.false.)
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index eb1d15fced..f56aca2e39 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -99,9 +99,9 @@ module MOM_energetic_PBL
                              ! problems, in m s-1.  If the value is small enough,
                              ! this should not affect the solution.
   real    :: omega           !   The Earth's rotation rate, in s-1.
-  logical :: use_omega       !   If true, use the absolute rotation rate instead
-                             ! of the vertical component of rotation when
-                             ! setting the decay scale for turbulence.
+  real    :: omega_frac      !   When setting the decay scale for turbulence, use
+                             ! this fraction of the absolute rotation rate blended
+                             ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   real    :: wstar_ustar_coef ! A ratio relating the efficiency with which
                              ! convectively released energy is converted to a
                              ! turbulent velocity, relative to mechanically
@@ -443,10 +443,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, &
       endif
       if (U_Star < CS%ustar_min) U_Star = CS%ustar_min
 
-      if (CS%use_omega) then ; absf(i) = 2.0*CS%omega
+      if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega
       else
         absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
                      (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J))))
+        if (CS%omega_frac > 0.0) &
+          absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
       endif
 
       mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3))
@@ -1117,8 +1119,9 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS)
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_energetic_PBL"  ! This module's name.
+  real :: omega_frac_dflt
   integer :: isd, ied, jsd, jed
-  logical :: use_temperature
+  logical :: use_temperature, use_omega
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
   if (associated(CS)) then
@@ -1157,10 +1160,20 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS)
   call get_param(param_file, mod, "OMEGA",CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)
-  call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+  call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+  omega_frac_dflt = 0.0
+  if (use_omega) then
+    call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+    omega_frac_dflt = 1.0
+  endif
+  call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                 "When setting the decay scale for turbulence, use this \n"//&
+                 "fraction of the absolute rotation rate blended with the \n"//&
+                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                 units="nondim", default=omega_frac_dflt)
   call get_param(param_file, mod, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, &
                  "A ratio relating the efficiency with which convectively\n"//&
                  "released energy is converted to a turbulent velocity,\n"//&
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index cc8059355a..171bdfbf12 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -255,6 +255,9 @@ module MOM_set_diffusivity
   logical :: ML_use_omega     ! If true, use absolute rotation rate instead
                               ! of the vertical component of rotation when
                               ! setting the decay scale for mixed layer turbulence.
+  real    :: ML_omega_frac    !   When setting the decay scale for turbulence, use
+                              ! this fraction of the absolute rotation rate blended
+                              ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2.
   logical :: user_change_diff ! If true, call user-defined code to change diffusivity.
   logical :: useKappaShear    ! If true, use the kappa_shear module to find the
                               ! shear-driven diapycnal diffusivity.
@@ -1791,11 +1794,13 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int)
   do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + GV%H_to_m*h(i,j,k) ; enddo ; enddo
 
   do i=is,ie ; if (do_i(i)) then
-    if (CS%ML_use_omega) then
+    if (CS%ML_omega_frac >= 1.0) then
       f_sq = 4.0*Omega2
     else
       f_sq = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
                    (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
+      if (CS%ML_omega_frac > 0.0) &
+        f_sq = CS%ML_omega_frac*4.0*Omega2 + (1.0-CS%ML_omega_frac)*f_sq
     endif
 
     ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2
@@ -2509,13 +2514,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp
 
   real :: decay_length, utide, zbot, hamp
   type(vardesc) :: vd
-  logical :: read_tideamp
+  logical :: read_tideamp, ML_use_omega
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_set_diffusivity"  ! This module's name.
   character(len=20)  :: tmpstr
   character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file
   real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data
+  real :: omega_frac_dflt
   integer :: i, j, is, ie, js, je
 
   if (associated(CS)) then
@@ -2587,11 +2593,20 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp
     call get_param(param_file, mod, "TKE_DECAY", CS%TKE_decay, &
                  "The ratio of the natural Ekman depth to the TKE decay scale.", &
                  units="nondim", default=2.5)
-    call get_param(param_file, mod, "ML_USE_OMEGA", CS%ML_use_omega, &
+    call get_param(param_file, mod, "ML_USE_OMEGA", ML_use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence in the mixed layer. \n"//&
-                 "This is only used if ML_RADIATION is true.", default=.false.)
+                 "scale for turbulence.", default=.false., do_not_log=.true.)
+    omega_frac_dflt = 0.0
+    if (ML_use_omega) then
+      call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+      omega_frac_dflt = 1.0
+    endif
+    call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%ML_omega_frac, &
+                   "When setting the decay scale for turbulence, use this \n"//&
+                   "fraction of the absolute rotation rate blended with the \n"//&
+                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                   units="nondim", default=omega_frac_dflt)
   endif
 
   call get_param(param_file, mod, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 84fb3d2aa8..2d7813586f 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -116,9 +116,9 @@ module MOM_set_visc
                             ! this should not affect the solution.
   real    :: TKE_decay      ! The ratio of the natural Ekman depth to the TKE
                             ! decay scale, nondimensional.
-  logical :: use_omega      !   If true, use the absolute rotation rate instead
-                            ! of the vertical component of rotation when
-                            ! setting the decay scale for turbulence.
+  real    :: omega_frac     !   When setting the decay scale for turbulence, use
+                            ! this fraction of the absolute rotation rate blended
+                            ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
   logical :: debug          ! If true, write verbose checksums for debugging purposes.
   type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
@@ -1064,8 +1064,11 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS)
           vhtot(I) = 0.25 * dt_Rho0 * ((fluxes%tauy(i,J) + fluxes%tauy(i+1,J-1)) + &
                                        (fluxes%tauy(i,J-1) + fluxes%tauy(i+1,J)))
 
-          if (CS%use_omega) then ; absf = 2.0*CS%omega
-          else ; absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1))) ; endif
+          if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
+            absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1)))
+            if (CS%omega_frac > 0.0) &
+              absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)
+          endif
           U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i+1,j)))
           Idecay_len_TKE(I) = ((absf / U_Star) * CS%TKE_decay) * H_to_m
         endif
@@ -1304,8 +1307,12 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS)
           uhtot(i) = 0.25 * dt_Rho0 * ((fluxes%taux(I,j) + fluxes%tauy(I-1,j+1)) + &
                                        (fluxes%taux(I-1,j) + fluxes%tauy(I,j+1)))
 
-         if (CS%use_omega) then ; absf = 2.0*CS%omega
-         else ; absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ; endif
+         if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
+           absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
+           if (CS%omega_frac > 0.0) &
+             absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2)        
+         endif
+ 
          U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i,j+1)))
          Idecay_len_TKE(i) = ((absf / U_Star) * CS%TKE_decay) * H_to_m
 
@@ -1617,8 +1624,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS)
 !                 for this module
   real    :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt
   real    :: Kv_background
+  real    :: omega_frac_dflt
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
-  logical :: use_kappa_shear, adiabatic, differential_diffusion
+  logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mod = "MOM_set_visc"  ! This module's name.
@@ -1695,10 +1703,20 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS)
                  "mixed layer viscosity.  By default, \n"//&
                  "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
                  default=TKE_decay_dflt)
-    call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, &
+    call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, &
                  "If true, use the absolute rotation rate instead of the \n"//&
                  "vertical component of rotation when setting the decay \n"//&
-                 "scale for turbulence.", default=.false.)
+                   "scale for turbulence.", default=.false., do_not_log=.true.)
+    omega_frac_dflt = 0.0
+    if (use_omega) then
+      call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
+      omega_frac_dflt = 1.0
+    endif
+    call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, &
+                   "When setting the decay scale for turbulence, use this \n"//&
+                   "fraction of the absolute rotation rate blended with the \n"//&
+                   "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
+                   units="nondim", default=omega_frac_dflt)
     call get_param(param_file, mod, "OMEGA", CS%omega, &
                  "The rotation rate of the earth.", units="s-1", &
                  default=7.2921e-5)