From d4770e88fbbafe270f45f67b8a1336cc8679dcaa Mon Sep 17 00:00:00 2001
From: Brandon Reichl <breichl@users.noreply.github.com>
Date: Wed, 5 Feb 2025 18:03:37 -0500
Subject: [PATCH] Correct bug in kappa shear viscosity with vertex shear
 option. (#824)

* Correct bug in kappa shear viscosity with vertex shear option.

- Viscosities at vertices along the coast were incorrectly zero'd out.  This commit removes that mask so the non-zero shear driven viscosities can be interpolated from in the model.  This bug caused diffusivities to be very large in channels and potentially near coastlines.
- A bugfix flag is added with an option to use the current behavior for legacy purposes.

* Fix missing paranthesis in previous commit (VS_viscosity_bug)

* Update logging of vertex shear viscosity bug parameter
---
 .../vertical/MOM_kappa_shear.F90              | 21 +++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90
index 05f312225f..8efde6352f 100644
--- a/src/parameterizations/vertical/MOM_kappa_shear.F90
+++ b/src/parameterizations/vertical/MOM_kappa_shear.F90
@@ -97,6 +97,8 @@ module MOM_kappa_shear
                              !! time average TKE when there is mass in all layers.  Otherwise always
                              !! report the time-averaged TKE, as is currently done when there
                              !! are some massless layers.
+  logical :: VS_viscosity_bug !< If true, use a bug in the calculation of the viscosity that sets
+                             !! it to zero for all vertices that are on a coastline.
   logical :: restrictive_tolerance_check !< If false, uses the less restrictive tolerance check to
                              !! determine if a timestep is acceptable for the KS_it outer iteration
                              !! loop, as the code was originally written.  True uses the more
@@ -607,10 +609,17 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_
       enddo
     endif ; enddo ! i-loop
 
-    do K=1,nz+1 ; do I=IsB,IeB
-      tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K)
-      kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_2d(I,K,J2) ) * CS%Prandtl_turb
-    enddo ; enddo
+    if (CS%VS_viscosity_bug) then
+      do K=1,nz+1 ; do I=IsB,IeB
+        tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K)
+        kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_2d(I,K,J2) ) * CS%Prandtl_turb
+      enddo; enddo
+    else
+      do K=1,nz+1 ; do I=IsB,IeB
+        tke_io(I,J,K) = tke_2d(I,K)
+        kv_io(I,J,K) = kappa_2d(I,K,J2) * CS%Prandtl_turb
+      enddo; enddo
+    endif
     if (J>=G%jsc) then ; do K=1,nz+1 ; do i=G%isc,G%iec
       ! Set the diffusivities in tracer columns from the values at vertices.
       kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * &
@@ -1873,6 +1882,10 @@ function kappa_shear_init(Time, G, GV, US, param_file, diag, CS)
                  "If true, do the calculations of the shear-driven mixing "//&
                  "at the cell vertices (i.e., the vorticity points).", &
                  default=.false., do_not_log=just_read)
+  call get_param(param_file, mdl, "VERTEX_SHEAR_VISCOSITY_BUG", CS%VS_viscosity_bug, &
+                 "If true, use a bug in vertex shear that zeros out viscosities at "//&
+                 "vertices on coastlines.", &
+                 default=.true., do_not_log=just_read.or.(.not.CS%KS_at_vertex))
   call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, &
                  "The critical Richardson number for shear mixing.", &
                  units="nondim", default=0.25, do_not_log=just_read)