+(*)Add option for MEKE to calculate total depth

  Added the new run-time option, MEKE_FIXED_TOTAL_DEPTH, to use the actual total
thickness instead the nominal depth in bathyT in several of the MEKE
calculations.  Also simplified and corrected a minor dimensional inconsistency
in some expressions that effectively set masks when estimating the interface
height diffusivities at tracer points for use in MEKE.  By default, the answers
in the MOM-examples test cases are bitwise identical, but answers could change
in some cases due to the proper thickness weighting in the calculation of
MEKE%Kh_diff.  There are new entries in some MOM_parameter_doc files, so a minor
spelling error correction in a MOM_parameter_doc entry was also included in this
PR.

src/parameterizations/lateral/MOM_MEKE.F90

@@ -82,7 +82,9 @@ module MOM_MEKE
   real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered
                                 !! when computing beta in Rhines scale [nondim]
   real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1].
-
+  logical :: fixed_total_depth  !< If true, use the nominal bathymetric depth as the estimate of
+                        !! the time-varying ocean depth.  Otherwise base the depth on the total
+                        !! ocean mass per unit area.
   logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled.
   logical :: initialize !< If True, invokes a steady state solver to calculate MEKE.
   logical :: debug      !< If true, write out checksums of data for debugging
@@ -283,12 +285,17 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       enddo
     enddo
 
-    !$OMP parallel do default(shared)
-    do j=js-1,je+1 ; do i=is-1,ie+1
-      depth_tot(i,j) = G%bathyT(i,j)
-      !### Try changing this to:
-      ! depth_tot(i,j) = mass(i,j) * I_Rho0
-    enddo ; enddo
+    if (CS%fixed_total_depth) then
+      !$OMP parallel do default(shared)
+      do j=js-1,je+1 ; do i=is-1,ie+1
+        depth_tot(i,j) = G%bathyT(i,j)
+      enddo ; enddo
+    else
+      !$OMP parallel do default(shared)
+      do j=js-1,je+1 ; do i=is-1,ie+1
+        depth_tot(i,j) = mass(i,j) * I_Rho0
+      enddo ; enddo
+    endif
 
     if (CS%initialize) then
       call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot)
@@ -711,8 +718,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m
       if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
         beta_topo_x = 0. ; beta_topo_y = 0.
       else
-        !### Consider different combinations of these estimates of topographic beta, and the use
-        !    of the water column thickness instead of the bathymetric depth.
+        !### Consider different combinations of these estimates of topographic beta.
         beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
                       (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j)  &
                   / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) &
@@ -887,14 +893,13 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, &
       FatH = 0.25* ( ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) ) + &
                      ( G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1) ) )  ! Coriolis parameter at h points
 
-      ! If bathyT is zero, then a division by zero FPE will be raised.  In this
+      ! If depth_tot is zero, then a division by zero FPE will be raised.  In this
       ! case, we apply Adcroft's rule of reciprocals and set the term to zero.
       ! Since zero-bathymetry cells are masked, this should not affect values.
       if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then
         beta_topo_x = 0. ; beta_topo_y = 0.
       else
-        !### Consider different combinations of these estimates of topographic beta, and the use
-        !    of the water column thickness instead of the bathymetric depth.
+        !### Consider different combinations of these estimates of topographic beta.
         beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( &
                       (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j)  &
                  / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) &
@@ -1167,6 +1172,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS)
                  "If positive, is a fixed length contribution to the expression "//&
                  "for mixing length used in MEKE-derived diffusivity.", &
                  units="m", default=0.0, scale=US%m_to_L)
+  call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", CS%fixed_total_depth, &
+                 "If true, use the nominal bathymetric depth as the estimate of the "//&
+                 "time-varying ocean depth.  Otherwise base the depth on the total ocean mass"//&
+                 "per unit area.", default=.true.)
+
   call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", CS%aDeform, &
                  "If positive, is a coefficient weighting the deformation scale "//&
                  "in the expression for mixing length used in MEKE-derived diffusivity.", &

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -145,6 +145,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
     KH_u_CFL      ! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1]
   real, dimension(SZI_(G), SZJB_(G)) :: &
     KH_v_CFL      ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1]
+  real, dimension(SZI_(G), SZJ_(G)) :: &
+    htot          ! The sum of the total layer thicknesses [H ~> m or kg m-2]
   real :: Khth_Loc_u(SZIB_(G), SZJ_(G))
   real :: Khth_Loc_v(SZI_(G), SZJB_(G))
   real :: Khth_Loc(SZIB_(G), SZJB_(G))  ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1]
@@ -479,38 +481,41 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
         ! thicknesses across u and v faces, converted to 0/1 mask
         ! layer average of the interface diffusivities KH_u and KH_v
         do j=js,je ; do I=is-1,ie
-          hu(I,j)       = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect)
-          if (hu(I,j) /= 0.0) hu(I,j) = 1.0
-          !### The same result would be accomplished with the following without a division:
-          ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0
+          ! This expression uses harmonic mean thicknesses:
+          ! hu(I,j)       = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
+          ! This expression is a 0/1 mask based on depths where there are thick layers:
+          hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0
           KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1))
         enddo ; enddo
         do J=js-1,je ; do i=is,ie
-          hv(i,J)       = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
-          if (hv(i,J) /= 0.0) hv(i,J) = 1.0
-          !### The same result would be accomplished with the following without a division:
-          ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0
+          ! This expression uses harmonic mean thicknesses:
+          ! hv(i,J)       = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
+          ! This expression is a 0/1 mask based on depths where there are thick layers:
+          hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0
           KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1))
         enddo ; enddo
-        ! diagnose diffusivity at T-point
-        !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent.
+        ! diagnose diffusivity at T-points
         do j=js,je ; do i=is,ie
-          Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j))  &
-                        +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) &
-                       / (hu(I-1,j)+hu(I,j)+hv(i,J-1)+hv(i,J)+h_neglect)
+          Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j) + hu(I,j)*KH_u_lay(I,j)) + &
+                         (hv(i,J-1)*KH_v_lay(i,J-1) + hv(i,J)*KH_v_lay(i,J))) / &
+                        ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + 1.0e-20)
+          ! Use this denominator instead if hu and hv are actual thicknesses rather than a 0/1 mask:
+          !              ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + h_neglect)
         enddo ; enddo
       enddo
 
       if (CS%Use_KH_in_MEKE) then
         MEKE%Kh_diff(:,:) = 0.0
+        htot(:,:) = 0.0
         do k=1,nz
           do j=js,je ; do i=is,ie
             MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + Kh_t(i,j,k) * h(i,j,k)
+            htot(i,j) = htot(i,j) + h(i,j,k)
           enddo ; enddo
         enddo
 
         do j=js,je ; do i=is,ie
-          MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j))
+          MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0*GV%m_to_H, htot(i,j))
         enddo ; enddo
       endif
 

src/parameterizations/vertical/MOM_ALE_sponge.F90

@@ -218,7 +218,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h,
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
   call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, &
-                 "If true, use the order of arithmetic for horizonal regridding that recovers "//&
+                 "If true, use the order of arithmetic for horizontal regridding that recovers "//&
                  "the answers from the end of 2018.  Otherwise, use rotationally symmetric "//&
                  "forms of the same expressions.", default=default_2018_answers)
   call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, &
@@ -478,7 +478,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
   call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, &
-                 "If true, use the order of arithmetic for horizonal regridding that recovers "//&
+                 "If true, use the order of arithmetic for horizontal regridding that recovers "//&
                  "the answers from the end of 2018 and retain a bug in the 3-dimensional mask "//&
                  "returned in certain cases.  Otherwise, use rotationally symmetric "//&
                  "forms of the same expressions and initialize the mask properly.", &