Changed constraint on Leith viscosity.

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -221,16 +221,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
     FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction (W/m2)
     Leith_Kh_h, & ! Leith Laplacian viscosity at h-points (m2 s-1)
     Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points (m4 s-1)
-    pl            ! Planetary number (nondim)
-
+    pl_h          ! Planetary number (nondim)
+ 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain (s-1)
     sh_xy,  &     ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms
     str_xy, &     ! str_xy is the cross term in the stress tensor (H m2 s-2)
     bhstr_xy, &   ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2)
     vort_xy, & ! Vertical vorticity (dv/dx - du/dy) (s-1)
     Leith_Kh_q, & ! Leith Laplacian viscosity at q-points (m2 s-1)
-    Leith_Ah_q    ! Leith bi-harmonic viscosity at q-points (m4 s-1)
+    Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points (m4 s-1)
+    pl_q          ! Planetary number (nondim)
 
   real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
     Ah_q, &   ! biharmonic viscosity at corner points (m4/s)
@@ -272,7 +273,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
   real :: RoScl     ! The scaling function for MEKE source term
   real :: FatH      ! abs(f) at h-point for MEKE source term (s-1)
   real :: local_strain ! Local variable for interpolating computed strain rates (s-1).
-  real :: beta, u_scale, epsilon, grid_sp_h2, pl_u, pl_v, mod_Leith_pl
+  real :: beta, u_scale, epsilon, grid_sp_h2, grid_sp_q2
   real :: DY_dxBu, DX_dyBu
   logical :: rescale_Kh, legacy_bound
   logical :: find_FrictWork
@@ -578,28 +579,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
           beta = sqrt( G%dF_dx(i,j)**2 + G%dF_dy(i,j)**2 )
           u_scale = sqrt((0.5*(u(I,j,k)+u(I-1,j,k)))**2 + (0.5*(v(i,J,k)+v(i,J-1,k)))**2)
           grid_sp_h2 = (2.0*CS%DX2h(i,j)*CS%DY2h(i,j)) / (CS%DX2h(i,j) + CS%DY2h(i,j))
-          pl(i,j) = beta * grid_sp_h2 / (u_scale + epsilon)
+          pl_h(i,j) = beta * grid_sp_h2 / (u_scale + epsilon)
         enddo; enddo
 
-        mod_Leith_pl = 1.0
+        do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+          beta = sqrt( (0.25*(G%dF_dx(i,j)+G%dF_dx(i+1,j)+G%dF_dx(i,j+1)+G%dF_dx(i+1,j+1))**2) + &
+                       (0.25*(G%dF_dy(i,j)+G%dF_dy(i+1,j)+G%dF_dy(i,j+1)+G%dF_dy(i+1,j+1))**2) )
+          u_scale = sqrt((0.5*(u(I,j,k)+u(I,j+1,k)))**2 + (0.5*(v(i,J,k)+v(i,J+1,k)))**2)
+          grid_sp_q2 = (2.0*CS%DX2q(I,J)*CS%DY2q(I,J)) / (CS%DX2q(I,J) + CS%DY2q(I,J))         
+          pl_q(I,J) = beta * grid_sp_q2 / (u_scale + epsilon)
+        enddo ; enddo
+
         do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
-          pl_v = 0.5 * (pl(i,j) + pl(i,j+1))
-          if (pl_v > 1.0) then
-            vort_xy_dx(i,J) = 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1))
-            mod_Leith_pl = 0.0
-          else
             vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1))
-          endif
         enddo ; enddo
         do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
-          pl_u = 0.5 * (pl(i,j) + pl(i+1,j))
-          if (pl_u > 1.0) then
-            vort_xy_dy(I,j) = 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j))
-            mod_Leith_pl = 0.0
-          else
             vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j))
-          endif
         enddo ; enddo
+
       endif
 
       mod_Leith = 0.; if (CS%modified_Leith) mod_Leith = 1.0
@@ -616,11 +613,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
                 (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1))))
       endif
       if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
-        vert_vort_mag = sqrt( &
+        if (pl_h(i,j) > 1) then
+          vert_vort_mag = sqrt( G%dF_dx(i,j)**2 + G%dF_dy(i,j)**2 )
+        else
+          vert_vort_mag = sqrt( &
           0.5*((vort_xy_dx(i,J-1)*vort_xy_dx(i,J-1) + vort_xy_dx(i,J)*vort_xy_dx(i,J)) + &
               (vort_xy_dy(I-1,j)*vort_xy_dy(I-1,j) + vort_xy_dy(I,j)*vort_xy_dy(I,j))) + &
           mod_Leith*mod_Leith_pl*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I-1,j)* &
           div_xx_dx(I-1,j)) + (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i,J-1)*div_xx_dy(i,J-1))))
+        endif
       endif
       if (CS%better_bound_Ah .or. CS%better_bound_Kh) then
         hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / &
@@ -749,7 +750,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
                   (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j))))
       endif
       if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
-        vert_vort_mag = sqrt( &
+        if (pl_q(I,J) > 1) then 
+          vert_vort_mag = sqrt( &
+                (0.25*(G%dF_dx(i,j)+G%dF_dx(i+1,j)+G%dF_dx(i,j+1)+G%dF_dx(i+1,j+1))**2) + &
+                (0.25*(G%dF_dy(i,j)+G%dF_dy(i+1,j)+G%dF_dy(i,j+1)+G%dF_dy(i+1,j+1))**2) )
+        else
+          vert_vort_mag = sqrt( &
           0.5*((vort_xy_dx(i,J)*vort_xy_dx(i,J) + vort_xy_dx(i+1,J)*vort_xy_dx(i+1,J)) + &
               (vort_xy_dy(I,j)*vort_xy_dy(I,j) + vort_xy_dy(I,j+1)*vort_xy_dy(I,j+1))) + &
           mod_Leith*mod_Leith_pl*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)* &