Bounds vort_xy_* in Leith.

NOAA-GFDL · Oct 17, 2018 · 92e2875 · 92e2875
1 parent d02fd47
commit 92e2875
Showing 1 changed file with 29 additions and 7 deletions.
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -220,7 +220,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
     bhstr_xx,&    ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2)
     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)
+    Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points (m4 s-1)
+    pl            ! Planetary number (nondim)
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain (s-1)
@@ -271,6 +272,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 :: DY_dxBu, DX_dyBu
   logical :: rescale_Kh, legacy_bound
   logical :: find_FrictWork
@@ -572,11 +574,31 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
 
       ! Add in beta for the Leith viscosity
       if (CS%use_beta_in_Leith) then
+        do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+          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)
+        enddo; enddo
+
+        mod_Leith_pl = 1.0
         do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
-          vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+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
-          vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j) )
+          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
 
@@ -597,8 +619,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
         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*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))))
+          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
       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)) / &
@@ -730,8 +752,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
         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*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)*div_xx_dx(I,j+1)) + &
-              (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J))))
+          mod_Leith*mod_Leith_pl*0.5*((div_xx_dx(I,j)*div_xx_dx(I,j) + div_xx_dx(I,j+1)* &
+          div_xx_dx(I,j+1)) + (div_xx_dy(i,J)*div_xx_dy(i,J) + div_xx_dy(i+1,J)*div_xx_dy(i+1,J))))
       endif
       h2uq = 4.0 * h_u(I,j) * h_u(I,j+1)
       h2vq = 4.0 * h_v(i,J) * h_v(i+1,J)