Preparing lateral_coeffs and thickness_diffuse for QG Leith calculati…

…ons.

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -587,6 +587,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
 
       else
 
+        do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+          div_xx_dx(I,j) = 0.0
+        enddo ; enddo
+        do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2
+          div_xx_dy(i,J) = 0.0
+        enddo ; enddo
         do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
           grad_div_mag_h(i,j) = 0.0
         enddo ; enddo
@@ -616,7 +622,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS,
 
       if (CS%use_QG_Leith_visc) then
         call calc_QG_Leith_viscosity(VarMix, G, GV, h, k, vort_xy_dx, vort_xy_dy, &
-                                     grad_div_mag_h, grad_div_mag_q)
+                                     div_xx_dx, div_xx_dy)
       endif
 
       do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -91,6 +91,12 @@ module MOM_lateral_mixing_coeffs
     slope_y => NULL(), &  !< Meridional isopycnal slope (non-dimensional)
     ebt_struct => NULL()  !< Vertical structure function to scale diffusivities with (non-dim)
 
+  real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
+    Laplac3_const_u       !< Laplacian  metric-dependent constants (nondim)
+
+  real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: &
+    Laplac3_const_v       !< Laplacian  metric-dependent constants (nondim)
+
   ! Parameters
   integer :: VarMix_Ktop  !< Top layer to start downward integrals
   real :: Visbeck_L_scale !< Fixed length scale in Visbeck formula
@@ -111,23 +117,7 @@ module MOM_lateral_mixing_coeffs
 
   ! Leith parameters
   logical :: use_QG_Leith_GM    !< If true, uses the QG Leith viscosity as the GM coefficient
-!  logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity
-!  logical :: Leith_Kh          !< If true, enables the Leith scheme
-!  logical :: modified_Leith    !< if true, include the divergence contribution to Leith viscosity
-!  logical :: Leith_Ah          !< If true, enables the bi-harmonic Leith scheme
-!  logical :: no_slip           !< If true, no slip boundary conditions are used.
-!                               !! Otherwise free slip boundary conditions are assumed.
-!                               !! The implementation of the free slip boundary
-!                               !! conditions on a C-grid is much cleaner than the
-!                               !! no slip boundary conditions. The use of free slip
-!                               !! b.c.s is strongly encouraged. The no slip b.c.s
-!                               !! are not implemented with the biharmonic viscosity.
-!  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
-!    Laplac3_const_xx, & ! Laplacian  metric-dependent constants (nondim)
-!    biharm5_const_xx    ! Biharmonic metric-dependent constants (nondim)
-!  real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
-!    Laplac3_const_xy, & ! Laplacian  metric-dependent constants (nondim)
-!    biharm5_const_xy    ! Biharmonic metric-dependent constants (nondim)
+  logical :: use_beta_in_QG_Leith ! If true, includes the beta term in the QG Leith GM coefficient
 
   ! Diagnostics
   !>@{
@@ -743,7 +733,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, CS, e, calculate_slopes)
 end subroutine calc_slope_functions_using_just_e
 
 !> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
-subroutine calc_QG_Leith_viscosity(CS, G, GV, h, k, vort_xy_dx, vort_xy_dy)
+subroutine calc_QG_Leith_viscosity(CS, G, GV, h, k, vort_xy_dx, vort_xy_dy, div_xx_dx, div_xx_dy)
   type(VarMix_CS),                           pointer     :: CS !< Variable mixing coefficients
   type(ocean_grid_type),                     intent(in)  :: G !< Ocean grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< The ocean's vertical grid structure.
@@ -753,6 +743,8 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, h, k, vort_xy_dx, vort_xy_dy)
   integer,                                   intent(in)  :: k !< Layer for which to calculate vorticity magnitude
   real, dimension(SZI_(G),SZJB_(G)),         intent(out) :: vort_xy_dx !< x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 s-1)
   real, dimension(SZIB_(G),SZJ_(G)),         intent(out) :: vort_xy_dy !< y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 s-1)
+  real, dimension(SZIB_(G),SZJ_(G)),         intent(out) :: div_xx_dx  ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 s-1)
+  real, dimension(SZI_(G),SZJB_(G)),         intent(out) :: div_xx_dy  ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 s-1)
 !  real, dimension(SZI_(G),SZJ_(G)),          intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity at h-points (m2 s-1)
 !  real, dimension(SZIB_(G),SZJB_(G)),        intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity at q-points (m2 s-1)
 !  real, dimension(SZI_(G),SZJ_(G)),          intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity at h-points (m4 s-1)
@@ -809,17 +801,54 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, h, k, vort_xy_dx, vort_xy_dy)
       vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * &
             ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) &
             + ( h_at_u(I-1,j) * dslopex_dz(I-1,j) + h_at_u(I,j+1) * dslopex_dz(I,j+1) ) ) / &
-              ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) )
+              ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff)
     enddo ; enddo
     do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
       f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
       vort_xy_dy(I,j) = vort_xy_dx(I,j) - f * &
             ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) &
             + ( h_at_v(i,J-1) * dslopey_dz(i,J-1) + h_at_v(i+1,J) * dslopey_dz(i+1,J) ) ) / &
-              ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) )
+              ( ( h_at_v(i,J) + h_at_v(i+1,J-1) ) + ( h_at_v(i,J-1) + h_at_v(i+1,J) ) + GV%H_subroundoff)
     enddo ; enddo
 
 
+    if (CS%use_QG_Leith_GM) then
+      if (CS%use_beta_in_QG_Leith) then
+        do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+          beta_u(I,j) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + &
+                        (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2) )
+        enddo ; enddo
+        do J=js-2,Jeq+1 ; do i=Isq-1,Ieq+2  
+          beta_v(i,J) = sqrt( (0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + &
+                        (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2) )
+        enddo ; enddo
+      endif
+
+      do j=js-1,Jeq+1 ; do I=is-2,Ieq
+        grad_vort_mag_u(I,j) = SQRT(vort_xy_dy(I,j)**2  + (0.25*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J) &
+                             + vort_xy_dx(i,J-1) + vort_xy_dx(i+1,J-1)))**2)
+        grad_div_mag_u(I,j) = SQRT(div_xx_dx(I,j)**2  + (0.25*(div_xx_dy(i,J) + div_xx_dy(i+1,J) &
+                             + div_xx_dy(i,J-1) + div_xx_dy(i+1,J-1)))**2) 
+        if (CS%use_beta_in_QG_Leith) then
+          vert_vort_mag = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), beta_u(I,j)**3)
+        else
+          vert_vort_mag = grad_vort_mag_u(I,j) + grad_div_mag_u(I,j)
+        endif
+      enddo ; enddo
+
+      do J=js-2,Jeq ; do i=is-1,Ieq+1
+        grad_vort_mag_v(i,J) = SQRT(vort_xy_dx(i,J)**2  + (0.25*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j) &
+                             + vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j+1)))**2)
+        grad_div_mag_v(i,J) = SQRT(div_xx_dy(i,J)**2  + (0.25*(div_xx_dx(I,j) + div_xx_dx(I-1,j) &
+                             + div_xx_dx(I,j+1) + div_xx_dx(I-1,j+1)))**2)  
+        if (CS%use_beta_in_QG_Leith) then
+          vert_vort_mag = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), beta_v(i,J)**3)
+        else
+          vert_vort_mag = grad_vort_mag_v(i,J) + grad_div_mag_v(i,J)
+        endif
+      enddo ; enddo
+    endif
+
 end subroutine calc_QG_Leith_viscosity
 
 subroutine VarMix_init(Time, G, param_file, diag, CS)
@@ -837,8 +866,8 @@ subroutine VarMix_init(Time, G, param_file, diag, CS)
   logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
   real :: MLE_front_length
   real :: Leith_Lap_const      ! The non-dimensional coefficient in the Leith viscosity
-  real :: Leith_bi_const       ! The non-dimensional coefficient in the bi-harmonic Leith viscosity
-  real :: DX2, DY2, grid_sp_2, grid_sp_3 ! Intermediate quantities for Leith metrics
+  real :: grid_sp_u2, grid_sp_u3
+  real :: grid_sp_v2, grid_sp_v3 ! Intermediate quantities for Leith metrics
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
   character(len=40)  :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
@@ -1122,83 +1151,39 @@ subroutine VarMix_init(Time, G, param_file, diag, CS)
   endif
 
   ! Leith parameters
-!  call get_param(param_file, mdl, "USE_QG_LEITH", CS%use_QG_Leith, &
-!               "If true, use the QG Leith nonlinear eddy viscosity.", &
-!               default=.false.)
-!  call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, &
-!               "If true, use a Leith nonlinear eddy viscosity.", &
-!               default=CS%use_QG_Leith)
-!  if (CS%use_QG_Leith .and. .not. CS%Leith_Kh) call MOM_error(FATAL, &
-!           "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
-!           "LEITH_KH must be True when USE_QG_LEITH=True.")
-!  call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, &
-!               "If true, include the beta term in the QG Leith nonlinear eddy viscosity.", &
-!               default=CS%use_QG_Leith)
-!  call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, &
-!               "If true, use a biharmonic Leith nonlinear eddy \n"//&
-!               "viscosity.", default=.false.)
-!  call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, &
-!               "If true, add a term to Leith viscosity which is \n"//&
-!               "proportional to the gradient of divergence.", &
-!               default=.false.)
-!  call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
-!               "The nondimensional Laplacian Leith constant, \n"//&
-!               "often set to 1.0", units="nondim", default=0.0, &
-!                fail_if_missing = CS%Leith_Kh)
-!  call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, &
-!               "The nondimensional biharmonic Leith constant, \n"//&
-!               "typical values are thus far undetermined.", units="nondim", default=0.0, &
-!               fail_if_missing = CS%Leith_Ah)
-!  if (CS%Leith_Kh .or. CS%Leith_Ah) then
-!    in_use = .true.
-!    call get_param(param_file, mdl, "NOSLIP", CS%no_slip, &
-!                   "If true, no slip boundary conditions are used; otherwise \n"//&
-!                   "free slip boundary conditions are assumed. The \n"//&
-!                   "implementation of the free slip BCs on a C-grid is much \n"//&
-!                   "cleaner than the no slip BCs. The use of free slip BCs \n"//&
-!                   "is strongly encouraged, and no slip BCs are not used with \n"//&
-!                   "the biharmonic viscosity.", default=.false.)
+  call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%use_QG_Leith_GM, &
+               "If true, use the QG Leith viscosity as the GM coefficient.", &
+               default=.false.)
+
   if (CS%Use_QG_Leith_GM) then
+    call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
+               "The nondimensional Laplacian Leith constant, \n"//&
+               "often set to 1.0", units="nondim", default=0.0)
+
+    call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_QG_Leith, &
+               "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
+               default=.true.)
+
+    allocate(CS%Laplac3_const_u(IsdB:IedB,jsd:jed)) ; CS%Laplac3_const_u(:,:) = 0.0
+    allocate(CS%Laplac3_const_v(isd:ied,JsdB:JedB)) ; CS%Laplac3_const_v(:,:) = 0.0
+
+    do j=Jsq,Jeq+1 ; do I=is-1,Ieq
+      ! Static factors in the Leith schemes
+      grid_sp_u2 = G%dyCu(I,j)*G%dxCu(I,j)
+      grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
+      CS%Laplac3_const_u(I,j) = Leith_Lap_const * grid_sp_v3
+    enddo ; enddo
+    do j=js-1,Jeq ; do I=Isq,Ieq+1
+      ! Static factors in the Leith schemes
+      grid_sp_v2 = G%dyCv(i,J)*G%dxCu(i,J)
+      grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
+      CS%Laplac3_const_v(i,J) = Leith_Lap_const * grid_sp_v3
+    enddo ; enddo
+
     if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
            "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
            "USE_STORED_SLOPES must be True when using QG Leith.")
   endif
-!  if (CS%Leith_Kh) then
-!    allocate(CS%Laplac3_const_xx(isd:ied,jsd:jed)) ; CS%Laplac3_const_xx(:,:) = 0.0
-!    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-!      DX2 = G%dxT(i,j)*G%dxT(i,j)
-!      DY2 = G%dyT(i,j)*G%dyT(i,j)
-!      grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2)
-!      grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
-!      CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_3
-!    enddo ; enddo
-!    allocate(CS%Laplac3_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%Laplac3_const_xy(:,:) = 0.0
-!    do J=js-1,Jeq ; do I=is-1,Ieq
-!      DX2 = G%dxBu(I,J)*G%dxBu(I,J)
-!      DY2 = G%dyBu(I,J)*G%dyBu(I,J)
-!      grid_sp_2 = (2.0*DX2*DY2) / (DX2 + DY2)
-!      grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
-!      CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_3
-!    enddo ; enddo
-!  endif
-!  if (CS%Leith_Ah) then
-!    allocate(CS%biharm5_const_xx(isd:ied,jsd:jed)) ; CS%biharm5_const_xx(:,:) = 0.0
-!    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-!      DX2 = G%dxT(i,j)*G%dxT(i,j)
-!      DY2 = G%dyT(i,j)*G%dyT(i,j)
-!      grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2)
-!      grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
-!      CS%biharm5_const_xx(i,j) = Leith_bi_const * (grid_sp_2 * grid_sp_3)
-!    enddo ; enddo
-!    allocate(CS%biharm5_const_xy(IsdB:IedB,JsdB:JedB)) ; CS%biharm5_const_xy(:,:) = 0.0
-!    do J=js-1,Jeq ; do I=is-1,Ieq
-!      DX2 = G%dxBu(I,J)*G%dxBu(I,J)
-!      DY2 = G%dyBu(I,J)*G%dyBu(I,J)
-!      grid_sp_2 = (2.0*DX2*DY2) / (DX2+DY2)
-!      grid_sp_3 = grid_sp_2*sqrt(grid_sp_2)
-!      CS%biharm5_const_xy(I,J) = Leith_bi_const * (grid_sp_2 * grid_sp_3)
-!    enddo ; enddo
-!  endif
 
   ! If nothing is being stored in this class then deallocate
   if (in_use) then

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -52,6 +52,7 @@ module MOM_thickness_diffuse
                                  !! longer than DT, or 0 (the default) to use DT.
   integer :: nkml                !< number of layers within mixed layer
   logical :: debug               !< write verbose checksums for debugging purposes
+  logical :: QG_Leith_GM         !< If true, uses the QG Leith viscosity as the GM coefficient
   type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics
   real, pointer :: GMwork(:,:)       => NULL()  !< Work by thickness diffusivity (W m-2)
   real, pointer :: diagSlopeX(:,:,:) => NULL()  !< Diagnostic: zonal neutral slope (nondim)
@@ -1728,6 +1729,9 @@ subroutine thickness_diffuse_init(Time, G, GV, param_file, diag, CDp, CS)
                  "marginally unstable value in a pure layered model, but \n"//&
                  "much smaller numbers (e.g. 0.1) seem to work better for \n"//&
                  "ALE-based models.", units = "nondimensional", default=0.8)
+  call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%QG_Leith_GM, &
+               "If true, use the QG Leith viscosity as the GM coefficient.", &
+               default=.false.)
   if (CS%max_Khth_CFL < 0.0) CS%max_Khth_CFL = 0.0
   call get_param(param_file, mdl, "DETANGLE_INTERFACES", CS%detangle_interfaces, &
                  "If defined add 3-d structured enhanced interface height \n"//&