Options to enforce KHTR_MIN and KHTH_MIN in the whole water column an…

…d fix naming inconsistency (#294)

* Rename KhTh_use_ebt_struct to KhTr_use_ebt_struct

The parameter KHTR_USE_EBT_STRUCT was introduced with an incorrect
variable name within the local control structure. The variable was
named KhTh_use_ebt_struct, indicating thickness diffusivity, instead
of the correct KhTr_use_ebt_struct, which indicates tracer diffusivity.
This commit rectifies the naming inconsistency.

* Option to apply KHTR_MIN in the whole water column

This commit introduces a new parameter, FULL_DEPTH_KHTR_MIN, which
enforces a user-specified minimum diffusivity (KHTR_MIN) to be
applied throughout the entire water column, instead of only at
the surface. This option is available only when
KHTR_USE_EBT_STRUCT=True and KHTR_MIN > 0.

* Option to apply KHTH_MIN in the whole water column

This commit introduces a new parameter, FULL_DEPTH_KHTH_MIN, which
enforces a user-specified minimum diffusivity (KHTH_MIN) to be
applied throughout the entire water column instead of only at
the surface. This option is available only when
KHTH_USE_EBT_STRUCT=True and KHTH_MIN > 0.

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -49,6 +49,9 @@ module MOM_thickness_diffuse
   real    :: kappa_smooth        !< Vertical diffusivity used to interpolate more sensible values
                                  !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
   logical :: thickness_diffuse   !< If true, interfaces heights are diffused.
+  logical :: full_depth_khth_min !< If true, KHTH_MIN is enforced throughout the whole water column.
+                                 !! Otherwise, KHTH_MIN is only enforced at the surface. This parameter
+                                 !! is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0.
   logical :: use_FGNV_streamfn   !< If true, use the streamfunction formulation of
                                  !! Ferrari et al., 2010, which effectively emphasizes
                                  !! graver vertical modes by smoothing in the vertical.
@@ -301,10 +304,18 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
   enddo ; enddo
 
   if (khth_use_ebt_struct) then
-    !$OMP do
-    do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
-      KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
-    enddo ; enddo ; enddo
+    if (CS%full_depth_khth_min) then
+      !$OMP do
+      do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
+        KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
+        KH_u(I,j,K) = max(KH_u(I,j,K), CS%Khth_Min)
+      enddo ; enddo ; enddo
+    else
+      !$OMP do
+      do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
+        KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
+      enddo ; enddo ; enddo
+    endif
   else
     !$OMP do
     do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
@@ -397,10 +408,18 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
   endif
 
   if (khth_use_ebt_struct) then
-    !$OMP do
-    do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
-      KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
-    enddo ; enddo ; enddo
+    if (CS%full_depth_khth_min) then
+      !$OMP do
+      do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
+        KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
+        KH_v(i,J,K) = max(KH_v(i,J,K), CS%Khth_Min)
+      enddo ; enddo ; enddo
+    else
+      !$OMP do
+      do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
+        KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
+      enddo ; enddo ; enddo
+    endif
   else
     !$OMP do
     do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
@@ -2169,7 +2188,11 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
                        ! rotation [nondim].
   real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
                         ! temperature variance [nondim]
-  integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags.
+  logical :: khth_use_ebt_struct ! If true, uses the equivalent barotropic structure
+                                 ! as the vertical structure of thickness diffusivity.
+                                 ! Used to determine if FULL_DEPTH_KHTH_MIN should be
+                                 ! available.
+  integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
   integer :: i, j
 
   CS%initialized = .true.
@@ -2215,6 +2238,17 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
   call get_param(param_file, mdl, "KHTH_MIN", CS%KHTH_Min, &
                  "The minimum horizontal thickness diffusivity.", &
                  default=0.0, units="m2 s-1", scale=US%m_to_L**2*US%T_to_s)
+  call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", khth_use_ebt_struct, &
+                 "If true, uses the equivalent barotropic structure "//&
+                 "as the vertical structure of thickness diffusivity.",&
+                 default=.false., do_not_log=.true.)
+  if (khth_use_ebt_struct .and. CS%KHTH_Min>0.0) then
+    call get_param(param_file, mdl, "FULL_DEPTH_KHTH_MIN", CS%full_depth_khth_min, &
+                   "If true, KHTH_MIN is enforced throughout the whole water column. "//&
+                   "Otherwise, KHTH_MIN is only enforced at the surface. This parameter "//&
+                   "is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0.",      &
+                   default=.false.)
+  endif
   call get_param(param_file, mdl, "KHTH_MAX", CS%KHTH_Max, &
                  "The maximum horizontal thickness diffusivity.", &
                  default=0.0, units="m2 s-1", scale=US%m_to_L**2*US%T_to_s)

src/tracer/MOM_tracer_hor_diff.F90

@@ -52,8 +52,11 @@ module MOM_tracer_hor_diff
   real    :: max_diff_CFL         !< If positive, locally limit the along-isopycnal
                                   !! tracer diffusivity to keep the diffusive CFL
                                   !! locally at or below this value [nondim].
-  logical :: KhTh_use_ebt_struct  !< If true, uses the equivalent barotropic structure
+  logical :: KhTr_use_ebt_struct  !< If true, uses the equivalent barotropic structure
                                   !! as the vertical structure of tracer diffusivity.
+  logical :: full_depth_khtr_min  !< If true, KHTR_MIN is enforced throughout the whole water column.
+                                  !! Otherwise, KHTR_MIN is only enforced at the surface. This parameter
+                                  !! is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.
   logical :: Diffuse_ML_interior  !< If true, diffuse along isopycnals between
                                   !! the mixed layer and the interior.
   logical :: check_diffusive_CFL  !< If true, automatically iterate the diffusion
@@ -422,21 +425,40 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
         enddo
       enddo
     enddo
-    if (CS%KhTh_use_ebt_struct) then
-      do K=2,nz+1
-        do J=js-1,je
-          do i=is,ie
-            Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
+    if (CS%KhTr_use_ebt_struct) then
+      if (CS%full_depth_khtr_min) then
+        do K=2,nz+1
+          do J=js-1,je
+            do i=is,ie
+              Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
+              Coef_y(i,J,K) = max(Coef_y(i,J,K), CS%KhTr_min)
+            enddo
           enddo
         enddo
-      enddo
-      do k=2,nz+1
-        do j=js,je
-          do I=is-1,ie
-            Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
+        do k=2,nz+1
+          do j=js,je
+            do I=is-1,ie
+              Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
+              Coef_x(I,j,K) = max(Coef_x(I,j,K), CS%KhTr_min)
+            enddo
           enddo
         enddo
-      enddo
+      else
+        do K=2,nz+1
+          do J=js-1,je
+            do i=is,ie
+              Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
+            enddo
+          enddo
+        enddo
+        do k=2,nz+1
+          do j=js,je
+            do I=is-1,ie
+              Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
+            enddo
+          enddo
+        enddo
+      endif
     endif
 
     do itt=1,num_itts
@@ -478,7 +500,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
         enddo
       enddo
     enddo
-    if (CS%KhTh_use_ebt_struct) then
+    if (CS%KhTr_use_ebt_struct) then
       do K=2,nz+1
         do J=js-1,je
           do i=is,ie
@@ -605,7 +627,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
     do j=js,je ; do I=is-1,ie
       Kh_u(I,j,:) = G%mask2dCu(I,j)*Kh_u(I,j,1)
     enddo ; enddo
-    if (CS%KhTh_use_ebt_struct) then
+    if (CS%KhTr_use_ebt_struct) then
       do K=2,nz+1
         do j=js,je
           do I=is-1,ie
@@ -621,7 +643,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
     do J=js-1,je ; do i=is,ie
       Kh_v(i,J,:) = G%mask2dCv(i,J)*Kh_v(i,J,1)
     enddo ; enddo
-    if (CS%KhTh_use_ebt_struct) then
+    if (CS%KhTr_use_ebt_struct) then
       do K=2,nz+1
         do J=js-1,je
           do i=is,ie
@@ -647,7 +669,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_
                          (G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37)
       Kh_h(i,j,:) = normalize*G%mask2dT(i,j)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + &
                                              (Kh_v(i,J-1,1)+Kh_v(i,J,1)))
-      if (CS%KhTh_use_ebt_struct) then
+      if (CS%KhTr_use_ebt_struct) then
         do K=2,nz+1
           Kh_h(i,j,K) = normalize*G%mask2dT(i,j)*VarMix%ebt_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + &
                                                                             (Kh_v(i,J-1,1)+Kh_v(i,J,1)))
@@ -1630,7 +1652,7 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic
   call get_param(param_file, mdl, "KHTR", CS%KhTr, &
                  "The background along-isopycnal tracer diffusivity.", &
                  units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s)
-  call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, &
+  call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTr_use_ebt_struct, &
                  "If true, uses the equivalent barotropic structure "//&
                  "as the vertical structure of the tracer diffusivity.",&
                  default=.false.)
@@ -1642,6 +1664,13 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic
   call get_param(param_file, mdl, "KHTR_MIN", CS%KhTr_Min, &
                  "The minimum along-isopycnal tracer diffusivity.", &
                  units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s)
+  if (CS%KhTr_use_ebt_struct .and. CS%KhTr_Min > 0.0) then
+    call get_param(param_file, mdl, "FULL_DEPTH_KHTR_MIN", CS%full_depth_khtr_min, &
+                   "If true, KHTR_MIN is enforced throughout the whole water column. "//&
+                   "Otherwise, KHTR_MIN is only enforced at the surface. This parameter "//&
+                   "is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.",      &
+                   default=.false.)
+  endif
   call get_param(param_file, mdl, "KHTR_MAX", CS%KhTr_Max, &
                  "The maximum along-isopycnal tracer diffusivity.", &
                  units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s)