MOM_hor_visc.F90

!> Calculates horizontal viscosity and viscous stresses
module mom_hor_visc

! This file is part of MOM6. See LICENSE.md for the license.

use mom_checksums,             only : hchksum, bchksum
use mom_coms,                  only : min_across_pes
use mom_diag_mediator,         only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator,         only : diag_ctrl, time_type
use mom_domains,               only : pass_var, corner, pass_vector, agrid, bgrid_ne
use mom_error_handler,         only : mom_error, fatal, warning, is_root_pe
use mom_file_parser,           only : get_param, log_version, param_file_type
use mom_grid,                  only : ocean_grid_type
use mom_lateral_mixing_coeffs, only : varmix_cs, calc_qg_leith_viscosity
use mom_barotropic,            only : barotropic_cs, barotropic_get_tav
use mom_thickness_diffuse,     only : thickness_diffuse_cs, thickness_diffuse_get_kh
use mom_io,                    only : mom_read_data, slasher
use mom_meke_types,            only : meke_type
use mom_open_boundary,         only : ocean_obc_type, obc_direction_e, obc_direction_w
use mom_open_boundary,         only : obc_direction_n, obc_direction_s, obc_none
use mom_unit_scaling,          only : unit_scale_type
use mom_verticalgrid,          only : verticalgrid_type
use mom_variables,             only : accel_diag_ptrs

implicit none ; private

#include <MOM_memory.h>

public horizontal_viscosity, hor_visc_init, hor_visc_end

!> Control structure for horizontal viscosity
type, public :: hor_visc_cs ; private
  logical :: laplacian       !< Use a Laplacian horizontal viscosity if true.
  logical :: biharmonic      !< Use a biharmonic horizontal viscosity if true.
  logical :: debug           !< If true, write verbose checksums for debugging purposes.
  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.
  logical :: bound_kh        !< If true, the Laplacian coefficient is locally
                             !! limited to guarantee stability.
  logical :: better_bound_kh !< If true, use a more careful bounding of the
                             !! Laplacian viscosity to guarantee stability.
  logical :: bound_ah        !< If true, the biharmonic coefficient is locally
                             !! limited to guarantee stability.
  logical :: better_bound_ah !< If true, use a more careful bounding of the
                             !! biharmonic viscosity to guarantee stability.
  real    :: re_ah           !! If nonzero, the biharmonic coefficient is scaled
                             !< so that the biharmonic Reynolds number is equal to this.
  real    :: bound_coef      !< The nondimensional coefficient of the ratio of
                             !! the viscosity bounds to the theoretical maximum
                             !! for stability without considering other terms [nondim].
                             !! The default is 0.8.
  logical :: smagorinsky_kh  !< If true, use Smagorinsky nonlinear eddy
                             !! viscosity. KH is the background value.
  logical :: smagorinsky_ah  !< If true, use a biharmonic form of Smagorinsky
                             !! nonlinear eddy viscosity. AH is the background.
  logical :: leith_kh        !< If true, use 2D Leith nonlinear eddy
                             !! viscosity. KH is the background value.
  logical :: modified_leith  !< If true, use extra component of Leith viscosity
                             !! to damp divergent flow. To use, still set Leith_Kh=.TRUE.
  logical :: use_beta_in_leith !< If true, includes the beta term in the Leith viscosity
  logical :: leith_ah        !< If true, use a biharmonic form of 2D Leith
                             !! nonlinear eddy viscosity. AH is the background.
  logical :: use_qg_leith_visc    !< If true, use QG Leith nonlinear eddy viscosity.
                             !! KH is the background value.
  logical :: bound_coriolis  !< If true & SMAGORINSKY_AH is used, the biharmonic
                             !! viscosity is modified to include a term that
                             !! scales quadratically with the velocity shears.
  logical :: use_kh_bg_2d    !< Read 2d background viscosity from a file.
  real    :: kh_bg_min       !< The minimum value allowed for Laplacian horizontal
                             !! viscosity [L2 T-1 ~> m2 s-1]. The default is 0.0.
  logical :: use_land_mask   !< Use the land mask for the computation of thicknesses
                             !! at velocity locations. This eliminates the dependence on
                             !! arbitrary values over land or outside of the domain.
                             !! Default is False to maintain answers with legacy experiments
                             !! but should be changed to True for new experiments.
  logical :: anisotropic     !< If true, allow anisotropic component to the viscosity.
  logical :: add_les_viscosity!< If true, adds the viscosity from Smagorinsky and Leith to
                             !! the background viscosity instead of taking the maximum.
  real    :: kh_aniso        !< The anisotropic viscosity [L2 T-1 ~> m2 s-1].
  logical :: dynamic_aniso   !< If true, the anisotropic viscosity is recomputed as a function
                             !! of state. This is set depending on ANISOTROPIC_MODE.
  logical :: res_scale_meke  !< If true, the viscosity contribution from MEKE is scaled by
                             !! the resolution function.
  logical :: use_gme         !< If true, use GME backscatter scheme.
  logical :: answers_2018    !< If true, use the order of arithmetic and expressions that recover the
                             !! answers from the end of 2018.  Otherwise, use updated and more robust
                             !! forms of the same expressions.
  real    :: gme_h0          !< The strength of GME tapers quadratically to zero when the bathymetric
                             !! depth is shallower than GME_H0 [Z ~> m]
  real    :: gme_efficiency  !< The nondimensional prefactor multiplying the GME coefficient [nondim]
  real    :: gme_limiter     !< The absolute maximum value the GME coefficient is allowed to take [L2 T-1 ~> m2 s-1].
  real    :: min_grid_kh     !< Minimum horizontal Laplacian viscosity used to
                             !! limit the grid Reynolds number [L2 T-1 ~> m2 s-1]
  real    :: min_grid_ah     !< Minimun horizontal biharmonic viscosity used to
                             !! limit grid Reynolds number [L4 T-1 ~> m4 s-1]

  real allocable_, dimension(NIMEM_,NJMEM_) :: kh_bg_xx
                      !< The background Laplacian viscosity at h points [L2 T-1 ~> m2 s-1].
                      !! The actual viscosity may be the larger of this
                      !! viscosity and the Smagorinsky and Leith viscosities.
  real allocable_, dimension(NIMEM_,NJMEM_) :: kh_bg_2d
                      !< The background Laplacian viscosity at h points [L2 T-1 ~> m2 s-1].
                      !! The actual viscosity may be the larger of this
                      !! viscosity and the Smagorinsky and Leith viscosities.
  real allocable_, dimension(NIMEM_,NJMEM_) :: ah_bg_xx
                      !< The background biharmonic viscosity at h points [L4 T-1 ~> m4 s-1].
                      !! The actual viscosity may be the larger of this
                      !! viscosity and the Smagorinsky and Leith viscosities.
  real allocable_, dimension(NIMEM_,NJMEM_) :: reduction_xx
                      !< The amount by which stresses through h points are reduced
                      !! due to partial barriers [nondim].
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    kh_max_xx,      & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
    ah_max_xx,      & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
    n1n2_h,         & !< Factor n1*n2 in the anisotropic direction tensor at h-points
    n1n1_m_n2n2_h,  & !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points
    grid_sp_h2,     & !< Harmonic mean of the squares of the grid [L2 ~> m2]
    grid_sp_h3        !< Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: kh_bg_xy
                      !< The background Laplacian viscosity at q points [L2 T-1 ~> m2 s-1].
                      !! The actual viscosity may be the larger of this
                      !! viscosity and the Smagorinsky and Leith viscosities.
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: ah_bg_xy
                      !< The background biharmonic viscosity at q points [L4 T-1 ~> m4 s-1].
                      !! The actual viscosity may be the larger of this
                      !! viscosity and the Smagorinsky and Leith viscosities.
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: reduction_xy
                      !< The amount by which stresses through q points are reduced
                      !! due to partial barriers [nondim].
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    kh_max_xy,      & !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
    ah_max_xy,      & !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
    n1n2_q,         & !< Factor n1*n2 in the anisotropic direction tensor at q-points
    n1n1_m_n2n2_q     !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points

  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    dx2h,   & !< Pre-calculated dx^2 at h points [L2 ~> m2]
    dy2h,   & !< Pre-calculated dy^2 at h points [L2 ~> m2]
    dx_dyt, & !< Pre-calculated dx/dy at h points [nondim]
    dy_dxt    !< Pre-calculated dy/dx at h points [nondim]
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    dx2q,    & !< Pre-calculated dx^2 at q points [L2 ~> m2]
    dy2q,    & !< Pre-calculated dy^2 at q points [L2 ~> m2]
    dx_dybu, & !< Pre-calculated dx/dy at q points [nondim]
    dy_dxbu    !< Pre-calculated dy/dx at q points [nondim]
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    idx2dycu, & !< 1/(dx^2 dy) at u points [L-3 ~> m-3]
    idxdy2u     !< 1/(dx dy^2) at u points [L-3 ~> m-3]
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    idx2dycv, & !< 1/(dx^2 dy) at v points [L-3 ~> m-3]
    idxdy2v     !< 1/(dx dy^2) at v points [L-3 ~> m-3]

  ! The following variables are precalculated time-invariant combinations of
  ! parameters and metric terms.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    laplac2_const_xx, & !< Laplacian  metric-dependent constants [L2 ~> m2]
    biharm6_const_xx, & !< Biharmonic metric-dependent constants [L6 ~> m6]
    laplac3_const_xx, & !< Laplacian  metric-dependent constants [L3 ~> m3]
    biharm_const_xx,  & !< Biharmonic metric-dependent constants [L4 ~> m4]
    biharm_const2_xx, & !< Biharmonic metric-dependent constants [T L4 ~> s m4]
    re_ah_const_xx      !< Biharmonic metric-dependent constants [L3 ~> m3]

  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    laplac2_const_xy, & !< Laplacian  metric-dependent constants [L2 ~> m2]
    biharm6_const_xy, & !< Biharmonic metric-dependent constants [L6 ~> m6]
    laplac3_const_xy, & !< Laplacian  metric-dependent constants [L3 ~> m3]
    biharm_const_xy,  & !< Biharmonic metric-dependent constants [L4 ~> m4]
    biharm_const2_xy, & !< Biharmonic metric-dependent constants [T L4 ~> s m4]
    re_ah_const_xy      !< Biharmonic metric-dependent constants [L3 ~> m3]

  type(diag_ctrl), pointer :: diag => null() !< structure to regulate diagnostics

  ! real, pointer :: hf_diffu(:,:,:)   => NULL() ! Zonal hor. visc. accel. x fract. thickness [L T-2 ~> m s-2].
  ! real, pointer :: hf_diffv(:,:,:)   => NULL() ! Merdional hor. visc. accel. x fract. thickness [L T-2 ~> m s-2].
  ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
  ! The code is retained for degugging purposes in the future.

  !>@{
  !! Diagnostic id
  integer :: id_grid_re_ah = -1, id_grid_re_kh   = -1
  integer :: id_diffu     = -1, id_diffv         = -1
  ! integer :: id_hf_diffu  = -1, id_hf_diffv      = -1
  integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
  integer :: id_ah_h      = -1, id_ah_q          = -1
  integer :: id_kh_h      = -1, id_kh_q          = -1
  integer :: id_gme_coeff_h = -1, id_gme_coeff_q = -1
  integer :: id_vort_xy_q = -1, id_div_xx_h      = -1
  integer :: id_sh_xy_q = -1,    id_sh_xx_h      = -1
  integer :: id_frictwork = -1, id_frictworkintz = -1
  integer :: id_frictwork_gme = -1
  !>@}


end type hor_visc_cs

contains

!> Calculates the acceleration due to the horizontal viscosity.
!!
!! A combination of biharmonic and Laplacian forms can be used. The coefficient
!! may either be a constant or a shear-dependent form. The biharmonic is
!! determined by twice taking the divergence of an appropriately defined stress
!! tensor. The Laplacian is determined by doing so once.
!!
!! To work, the following fields must be set outside of the usual
!! is:ie range before this subroutine is called:
!!   u[is-2:ie+2,js-2:je+2]
!!   v[is-2:ie+2,js-2:je+2]
!!   h[is-1:ie+1,js-1:je+1]
subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
                                CS, OBC, BT, TD, ADp)
  type(ocean_grid_type),         intent(in)  :: G      !< The ocean's grid structure.
  type(verticalgrid_type),       intent(in)  :: GV     !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                 intent(in)  :: u      !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                 intent(in)  :: v      !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                 intent(inout) :: h    !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                 intent(out) :: diffu  !< Zonal acceleration due to convergence of
                                                       !! along-coordinate stress tensor [L T-2 ~> m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                 intent(out) :: diffv  !< Meridional acceleration due to convergence
                                                       !! of along-coordinate stress tensor [L T-2 ~> m s-2].
  type(meke_type),               pointer     :: MEKE   !< Pointer to a structure containing fields
                                                       !! related to Mesoscale Eddy Kinetic Energy.
  type(varmix_cs),               pointer     :: VarMix !< Pointer to a structure with fields that
                                                       !! specify the spatially variable viscosities
  type(unit_scale_type),         intent(in)  :: US     !< A dimensional unit scaling type
  type(hor_visc_cs),             pointer     :: CS     !< Control structure returned by a previous
                                                       !! call to hor_visc_init.
  type(ocean_obc_type), optional, pointer    :: OBC    !< Pointer to an open boundary condition type
  type(barotropic_cs),  optional, pointer    :: BT     !< Pointer to a structure containing
                                                       !! barotropic velocities.
  type(thickness_diffuse_cs), optional, pointer :: TD  !< Pointer to a structure containing
                                                       !! thickness diffusivities.
  type(accel_diag_ptrs), optional, pointer :: ADp      !< Acceleration diagnostic pointers

  ! Local variables
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    Del2u, &      ! The u-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]
    h_u, &        ! Thickness interpolated to u points [H ~> m or kg m-2].
    vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
    div_xx_dx, &  ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
    ubtav         ! zonal barotropic vel. ave. over baroclinic time-step [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G)) :: &
    Del2v, &      ! The v-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]
    h_v, &        ! Thickness interpolated to v points [H ~> m or kg m-2].
    vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
    div_xx_dy, &  ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
    vbtav         ! meridional barotropic vel. ave. over baroclinic time-step [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G)) :: &
    dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [T-1 ~> s-1]
    div_xx, &     ! Estimate of horizontal divergence at h-points [T-1 ~> s-1]
    sh_xx, &      ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
    sh_xx_bt, &   ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
    str_xx,&      ! str_xx is the diagonal term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2]
    str_xx_GME,&  ! smoothed diagonal term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2]
    bhstr_xx, &   ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]
    FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2]
    ! Leith_Kh_h, & ! Leith Laplacian viscosity at h-points [L2 T-1 ~> m2 s-1]
    ! Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points [L4 T-1 ~> m4 s-1]
    grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]
    grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]
    del2vort_h, & ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1]
    grad_div_mag_h, &     ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1]
    dudx, dvdy, &    ! components in the horizontal tension [T-1 ~> s-1]
    grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2]
    grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2]
    grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2]
    boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]

  real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2]
  real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2]

  real, dimension(SZIB_(G),SZJB_(G)) :: &
    dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
    dDel2vdx, dDel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1]
    dvdx_bt, dudy_bt,   & ! components in the barotropic shearing strain [T-1 ~> s-1]
    sh_xy,  &     ! horizontal shearing strain (du/dy + dv/dx) including metric terms [T-1 ~> s-1]
    sh_xy_bt, &   ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [T-1 ~> s-1]
    str_xy, &     ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2]
    str_xy_GME, & ! smoothed cross term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2]
    bhstr_xy, &   ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]
    vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1]
    Leith_Kh_q, & ! Leith Laplacian viscosity at q-points [L2 T-1 ~> m2 s-1]
    Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points [L4 T-1 ~> m4 s-1]
    grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]
    grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]
    Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1]
    grad_div_mag_q, &  ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1]
    grad_vel_mag_q, &  ! Magnitude of the velocity gradient tensor squared at q-points [T-2 ~> s-2]
    hq, &  ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2]
           ! This form guarantees that hq/hu < 4.
    grad_vel_mag_bt_q, &  ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2]
    boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim]

  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
    Ah_q, &      ! biharmonic viscosity at corner points [L4 T-1 ~> m4 s-1]
    Kh_q, &      ! Laplacian viscosity at corner points [L2 T-1 ~> m2 s-1]
    vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1]
    sh_xy_q,   & ! horizontal shearing strain at corner points [T-1 ~> s-1]
    GME_coeff_q, &  !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]
    max_diss_rate_q ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3]

  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: &
    KH_u_GME  !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: &
    KH_v_GME  !< interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    Ah_h, &          ! biharmonic viscosity at thickness points [L4 T-1 ~> m4 s-1]
    Kh_h, &          ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1]
    max_diss_rate_h, & ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3]
    FrictWork, &     ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2]
    FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2]
    div_xx_h,      & ! horizontal divergence [T-1 ~> s-1]
    sh_xx_h          ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    grid_Re_Kh, &    !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim]
    grid_Re_Ah, &    !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim]
    GME_coeff_h      !< GME coeff. at h-points [L2 T-1 ~> m2 s-1]
  real :: Ah         ! biharmonic viscosity [L4 T-1 ~> m4 s-1]
  real :: Kh         ! Laplacian  viscosity [L2 T-1 ~> m2 s-1]
  real :: AhSm       ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1]
  real :: AhLth      ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1]
  real :: mod_Leith  ! nondimensional coefficient for divergence part of modified Leith
                     ! viscosity. Here set equal to nondimensional Laplacian Leith constant.
                     ! This is set equal to zero if modified Leith is not used.
  real :: Shear_mag  ! magnitude of the shear [T-1 ~> s-1]
  real :: vert_vort_mag ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1]
  real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4].
  real :: hu, hv     ! Thicknesses interpolated by arithmetic means to corner
                     ! points; these are first interpolated to u or v velocity
                     ! points where masks are applied [H ~> m or kg m-2].
  real :: h_neglect  ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2]
  real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6]
  real :: hrat_min   ! minimum thicknesses at the 4 neighboring
                     ! velocity points divided by the thickness at the stress
                     ! point (h or q point) [nondim]
  real :: visc_bound_rem ! fraction of overall viscous bounds that
                         ! remain to be applied [nondim]
  real :: Kh_scale  ! A factor between 0 and 1 by which the horizontal
                    ! Laplacian viscosity is rescaled [nondim]
  real :: RoScl     ! The scaling function for MEKE source term [nondim]
  real :: FatH      ! abs(f) at h-point for MEKE source term [T-1 ~> s-1]
  real :: local_strain ! Local variable for interpolating computed strain rates [T-1 ~> s-1].
  real :: meke_res_fn ! A copy of the resolution scaling factor if being applied to MEKE. Otherwise =1.
  real :: GME_coeff ! The GME (negative) viscosity coefficient [L2 T-1 ~> m2 s-1]
  real :: GME_coeff_limiter ! Maximum permitted value of the GME coefficient [L2 T-1 ~> m2 s-1]
  real :: FWfrac    ! Fraction of maximum theoretical energy transfer to use when scaling GME coefficient [nondim]
  real :: DY_dxBu   ! Ratio of meridional over zonal grid spacing at vertices [nondim]
  real :: DX_dyBu   ! Ratio of zonal over meridiononal grid spacing at vertices [nondim]
  real :: DY_dxCv   ! Ratio of meridional over zonal grid spacing at faces [nondim]
  real :: DX_dyCu   ! Ratio of zonal over meridional grid spacing at faces [nondim]
  real :: Sh_F_pow  ! The ratio of shear over the absolute value of f raised to some power and rescaled [nondim]
  real :: backscat_subround ! The ratio of f over Shear_mag that is so small that the backscatter
                    ! calculation gives the same value as if f were 0 [nondim].
  real :: H0_GME    ! Depth used to scale down GME coefficient in shallow areas [Z ~> m]
  real :: KE        ! Local kinetic energy [L2 T-2 ~> m2 s-2]

  logical :: rescale_Kh, legacy_bound
  logical :: find_FrictWork
  logical :: apply_OBC = .false.
  logical :: use_MEKE_Ku
  logical :: use_MEKE_Au
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: i, j, k, n
  real :: inv_PI3, inv_PI2, inv_PI6
  is  = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  h_neglect  = gv%H_subroundoff
  h_neglect3 = h_neglect**3
  inv_pi3 = 1.0/((4.0*atan(1.0))**3)
  inv_pi2 = 1.0/((4.0*atan(1.0))**2)
  inv_pi6 = inv_pi3 * inv_pi3

  ah_h(:,:,:) = 0.0
  kh_h(:,:,:) = 0.0

  if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then
    apply_obc = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally
    apply_obc = .true.
  endif ; endif ; endif

  if (.not.associated(cs)) call mom_error(fatal, &
         "MOM_hor_visc: Module must be initialized before it is used.")
  if (.not.(cs%Laplacian .or. cs%biharmonic)) return

  find_frictwork = (cs%id_FrictWork > 0)
  if (cs%id_FrictWorkIntz > 0) find_frictwork = .true.
  if (associated(meke)) then
    if (associated(meke%mom_src)) find_frictwork = .true.
    backscat_subround = 0.0
    if (find_frictwork .and. associated(meke%mom_src) .and. (meke%backscatter_Ro_c > 0.0) .and. &
        (meke%backscatter_Ro_Pow /= 0.0)) &
      backscat_subround = (1.0e-16/meke%backscatter_Ro_c)**(1.0/meke%backscatter_Ro_Pow)
  endif

  rescale_kh = .false.
  if (associated(varmix)) then
    rescale_kh = varmix%Resoln_scaled_Kh
    if (rescale_kh .and. &
    (.not.associated(varmix%Res_fn_h) .or. .not.associated(varmix%Res_fn_q))) &
      call mom_error(fatal, "MOM_hor_visc: VarMix%Res_fn_h and " //&
        "VarMix%Res_fn_q both need to be associated with Resoln_scaled_Kh.")
  endif
  legacy_bound = (cs%Smagorinsky_Kh .or. cs%Leith_Kh) .and. &
                 (cs%bound_Kh .and. .not.cs%better_bound_Kh)

  ! Toggle whether to use a Laplacian viscosity derived from MEKE
  use_meke_ku = associated(meke%Ku)
  use_meke_au = associated(meke%Au)

  if (cs%use_GME) then
    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
      boundary_mask_h(i,j) = (g%mask2dCu(i,j) * g%mask2dCv(i,j) * g%mask2dCu(i-1,j) * g%mask2dCv(i,j-1))
    enddo ; enddo

    do j=js-2,jeq+1 ; do i=is-2,ieq+1
      boundary_mask_q(i,j) = (g%mask2dCv(i,j) * g%mask2dCv(i+1,j) * g%mask2dCu(i,j) * g%mask2dCu(i,j-1))
    enddo; enddo

    ! initialize diag. array with zeros
    gme_coeff_h(:,:,:) = 0.0
    gme_coeff_q(:,:,:) = 0.0
    str_xx_gme(:,:) = 0.0
    str_xy_gme(:,:) = 0.0

    ! Get barotropic velocities and their gradients
    call barotropic_get_tav(bt, ubtav, vbtav, g, us)
    call pass_vector(ubtav, vbtav, g%Domain)

    do j=js-1,je+2 ; do i=is-1,ie+2
      dudx_bt(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * ubtav(i,j) - &
                                     g%IdyCu(i-1,j) * ubtav(i-1,j))
      dvdy_bt(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * vbtav(i,j) - &
                                     g%IdxCv(i,j-1) * vbtav(i,j-1))
    enddo; enddo

    do j=jsq,jeq+1 ; do i=isq,ieq+1
      sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
    enddo ; enddo

    ! Components for the barotropic shearing strain
    do j=js-2,jeq+1 ; do i=is-2,ieq+1
      dvdx_bt(i,j) = cs%DY_dxBu(i,j)*(vbtav(i+1,j)*g%IdyCv(i+1,j) &
                                    - vbtav(i,j)*g%IdyCv(i,j))
      dudy_bt(i,j) = cs%DX_dyBu(i,j)*(ubtav(i,j+1)*g%IdxCu(i,j+1) &
                                    - ubtav(i,j)*g%IdxCu(i,j))
    enddo ; enddo

    call pass_vector(dudx_bt, dvdy_bt, g%Domain, stagger=bgrid_ne)
    call pass_vector(dvdx_bt, dudy_bt, g%Domain, stagger=agrid)

    if (cs%no_slip) then
      do j=js-1,jeq ; do i=is-1,ieq
        sh_xy_bt(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
      enddo ; enddo
    else
      do j=js-1,jeq ; do i=is-1,ieq
        sh_xy_bt(i,j) = g%mask2dBu(i,j) * ( dvdx_bt(i,j) + dudy_bt(i,j) )
      enddo ; enddo
    endif

    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
      grad_vel_mag_bt_h(i,j) = boundary_mask_h(i,j) * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
            (0.25*((dvdx_bt(i,j)+dvdx_bt(i-1,j-1))+(dvdx_bt(i,j-1)+dvdx_bt(i-1,j))))**2 + &
            (0.25*((dudy_bt(i,j)+dudy_bt(i-1,j-1))+(dudy_bt(i,j-1)+dudy_bt(i-1,j))))**2)
    enddo ; enddo

    do j=js-2,jeq+1 ; do i=is-2,ieq+1
      grad_vel_mag_bt_q(i,j) = boundary_mask_q(i,j) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &
            (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + &
            (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)
    enddo ; enddo

  endif ! use_GME

  !$OMP parallel do default(none) &
  !$OMP shared( &
  !$OMP   CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, &
  !$OMP   is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
  !$OMP   apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
  !$OMP   use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, &
  !$OMP   backscat_subround, GME_coeff_limiter, &
  !$OMP   h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, &
  !$OMP   diffu, diffv, max_diss_rate_h, max_diss_rate_q, &
  !$OMP   Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, &
  !$OMP   div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, &
  !$OMP   TD, KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah &
  !$OMP ) &
  !$OMP private( &
  !$OMP   i, j, k, n, &
  !$OMP   dudx, dudy, dvdx, dvdy, sh_xx, sh_xy, h_u, h_v, &
  !$OMP   Del2u, Del2v, DY_dxBu, DX_dyBu, sh_xx_bt, sh_xy_bt, &
  !$OMP   str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, &
  !$OMP   vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, &
  !$OMP   grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, &
  !$OMP   grad_vort_mag_h_2d, grad_vort_mag_q_2d, &
  !$OMP   grad_vel_mag_h, grad_vel_mag_q, &
  !$OMP   grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, &
  !$OMP   meke_res_fn, Shear_mag, vert_vort_mag, hrat_min, visc_bound_rem, &
  !$OMP   Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, &
  !$OMP   dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, &
  !$OMP   h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff &
  !$OMP )
  do k=1,nz

    ! The following are the forms of the horizontal tension and horizontal
    ! shearing strain advocated by Smagorinsky (1993) and discussed in
    ! Griffies and Hallberg (2000).

    ! Calculate horizontal tension
    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
      dudx(i,j) = cs%DY_dxT(i,j)*(g%IdyCu(i,j) * u(i,j,k) - &
                                  g%IdyCu(i-1,j) * u(i-1,j,k))
      dvdy(i,j) = cs%DX_dyT(i,j)*(g%IdxCv(i,j) * v(i,j,k) - &
                                  g%IdxCv(i,j-1) * v(i,j-1,k))
      sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
    enddo ; enddo

    ! Components for the shearing strain
    do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2
      dvdx(i,j) = cs%DY_dxBu(i,j)*(v(i+1,j,k)*g%IdyCv(i+1,j) - v(i,j,k)*g%IdyCv(i,j))
      dudy(i,j) = cs%DX_dyBu(i,j)*(u(i,j+1,k)*g%IdxCu(i,j+1) - u(i,j,k)*g%IdxCu(i,j))
    enddo ; enddo

    ! Interpolate the thicknesses to velocity points.
    ! The extra wide halos are to accommodate the cross-corner-point projections
    ! in OBCs, which are not ordinarily be necessary, and might not be necessary
    ! even with OBCs if the accelerations are zeroed at OBC points, in which
    ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH
    if (cs%use_land_mask) then
      do j=js-2,je+2 ; do i=isq-1,ieq+1
        h_u(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k))
      enddo ; enddo
      do j=jsq-1,jeq+1 ; do i=is-2,ie+2
        h_v(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k))
      enddo ; enddo
    else
      do j=js-2,je+2 ; do i=isq-1,ieq+1
        h_u(i,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
      enddo ; enddo
      do j=jsq-1,jeq+1 ; do i=is-2,ie+2
        h_v(i,j) = 0.5 * (h(i,j,k) + h(i,j+1,k))
      enddo ; enddo
    endif

    ! Adjust contributions to shearing strain and interpolated values of
    ! thicknesses on open boundaries.
    if (apply_obc) then ; do n=1,obc%number_of_segments
      j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
      if (obc%zero_strain .or. obc%freeslip_strain .or. obc%computed_strain) then
        if (obc%segment(n)%is_N_or_S .and. (j >= js-2) .and. (j <= jeq+1)) then
          do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
            if (obc%zero_strain) then
              dvdx(i,j) = 0. ; dudy(i,j) = 0.
            elseif (obc%freeslip_strain) then
              dudy(i,j) = 0.
            elseif (obc%computed_strain) then
              if (obc%segment(n)%direction == obc_direction_n) then
                dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &
                            (obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%IdxCu(i,j)
              else
                dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &
                            (u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdxCu(i,j+1)
              endif
            elseif (obc%specified_strain) then
              if (obc%segment(n)%direction == obc_direction_n) then
                dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j)*g%dxBu(i,j)
              else
                dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j+1)*g%dxBu(i,j)
              endif
            endif
          enddo
        elseif (obc%segment(n)%is_E_or_W .and. (i >= is-2) .and. (i <= ieq+1)) then
          do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
            if (obc%zero_strain) then
              dvdx(i,j) = 0. ; dudy(i,j) = 0.
            elseif (obc%freeslip_strain) then
              dvdx(i,j) = 0.
            elseif (obc%computed_strain) then
              if (obc%segment(n)%direction == obc_direction_e) then
                dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &
                            (obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%IdyCv(i,j)
              else
                dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &
                            (v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdyCv(i+1,j)
              endif
            elseif (obc%specified_strain) then
              if (obc%segment(n)%direction == obc_direction_e) then
                dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i,j)*g%dxBu(i,j)
              else
                dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i+1,j)*g%dxBu(i,j)
              endif
            endif
          enddo
        endif
      endif


      if (obc%segment(n)%direction == obc_direction_n) then
        ! There are extra wide halos here to accommodate the cross-corner-point
        ! OBC projections, but they might not be necessary if the accelerations
        ! are always zeroed out at OBC points, in which case the i-loop below
        ! becomes do i=is-1,ie+1. -RWH
        if ((j >= jsq-1) .and. (j <= jeq+1)) then
          do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
            h_v(i,j) = h(i,j,k)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_s) then
        if ((j >= jsq-1) .and. (j <= jeq+1)) then
          do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)
            h_v(i,j) = h(i,j+1,k)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_e) then
        if ((i >= isq-1) .and. (i <= ieq+1)) then
          do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
            h_u(i,j) = h(i,j,k)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_w) then
        if ((i >= isq-1) .and. (i <= ieq+1)) then
          do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)
            h_u(i,j) = h(i+1,j,k)
          enddo
        endif
      endif
    enddo ; endif
    ! Now project thicknesses across corner points on OBCs.
    if (apply_obc) then ; do n=1,obc%number_of_segments
      j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
      if (obc%segment(n)%direction == obc_direction_n) then
        if ((j >= js-2) .and. (j <= je)) then
          do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
            h_u(i,j+1) = h_u(i,j)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_s) then
        if ((j >= js-1) .and. (j <= je+1)) then
          do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+1,obc%segment(n)%HI%ied)
            h_u(i,j) = h_u(i,j+1)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_e) then
        if ((i >= is-2) .and. (i <= ie)) then
          do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
            h_v(i+1,j) = h_v(i,j)
          enddo
        endif
      elseif (obc%segment(n)%direction == obc_direction_w) then
        if ((i >= is-1) .and. (i <= ie+1)) then
          do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)
            h_v(i,j) = h_v(i+1,j)
          enddo
        endif
      endif
    enddo ; endif

    ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).
    ! dudy and dvdx include modifications at OBCs from above.
    if (cs%no_slip) then
      do j=js-2,jeq+1 ; do i=is-2,ieq+1
        sh_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) + dudy(i,j) )
      enddo ; enddo
    else
      do j=js-2,jeq+1 ; do i=is-2,ieq+1
        sh_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) + dudy(i,j) )
      enddo ; enddo
    endif

    !  Evaluate Del2u = x.Div(Grad u) and Del2v = y.Div( Grad u)
    if (cs%biharmonic) then
      do j=js-1,jeq+1 ; do i=isq-1,ieq+1
        del2u(i,j) = cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*sh_xx(i+1,j) - cs%dy2h(i,j)*sh_xx(i,j)) + &
                     cs%Idx2dyCu(i,j)*(cs%dx2q(i,j)*sh_xy(i,j) - cs%dx2q(i,j-1)*sh_xy(i,j-1))
      enddo ; enddo
      do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
        del2v(i,j) = cs%Idxdy2v(i,j)*(cs%dy2q(i,j)*sh_xy(i,j) - cs%dy2q(i-1,j)*sh_xy(i-1,j)) - &
                     cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*sh_xx(i,j+1) - cs%dx2h(i,j)*sh_xx(i,j))
      enddo ; enddo
      if (apply_obc) then; if (obc%zero_biharmonic) then
        do n=1,obc%number_of_segments
          i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
          if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
            do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
              del2v(i,j) = 0.
            enddo
          elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
            do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
              del2u(i,j) = 0.
            enddo
          endif
        enddo
      endif; endif
    endif

    ! Vorticity
    if (cs%no_slip) then
      do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2
        vort_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) - dudy(i,j) )
      enddo ; enddo
    else
      do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2
        vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )
      enddo ; enddo
    endif

    ! Divergence
    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
      div_xx(i,j) = dudx(i,j) + dvdy(i,j)
    enddo ; enddo

    if ((cs%Leith_Kh) .or. (cs%Leith_Ah)) then

      ! Vorticity gradient
      do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
        dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)
        vort_xy_dx(i,j) = dy_dxbu * (vort_xy(i,j) * g%IdyCu(i,j) - vort_xy(i-1,j) * g%IdyCu(i-1,j))
      enddo ; enddo

      do j=jsq-1,jeq+2 ; do i=isq-1,ieq+1
        dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)
        vort_xy_dy(i,j) = dx_dybu * (vort_xy(i,j) * g%IdxCv(i,j) - vort_xy(i,j-1) * g%IdxCv(i,j-1))
      enddo ; enddo

      ! Laplacian of vorticity
      do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
        dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)
        dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)

        del2vort_q(i,j) = dy_dxbu * (vort_xy_dx(i+1,j) * g%IdyCv(i+1,j) - vort_xy_dx(i,j) * g%IdyCv(i,j)) + &
                        dx_dybu * (vort_xy_dy(i,j+1) * g%IdyCu(i,j+1) - vort_xy_dy(i,j) * g%IdyCu(i,j))
      enddo ; enddo
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        del2vort_h(i,j) = 0.25*(del2vort_q(i,j) + del2vort_q(i-1,j) + del2vort_q(i,j-1) + del2vort_q(i-1,j-1))
      enddo ; enddo

      if (cs%modified_Leith) then

        ! Divergence gradient
        do j=jsq-1,jeq+2 ; do i=isq-1,ieq+1
          div_xx_dx(i,j) = g%IdxCu(i,j)*(div_xx(i+1,j) - div_xx(i,j))
        enddo ; enddo
        do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
          div_xx_dy(i,j) = g%IdyCv(i,j)*(div_xx(i,j+1) - div_xx(i,j))
        enddo ; enddo

        ! Magnitude of divergence gradient
        do j=jsq,jeq+1 ; do i=isq,ieq+1
          grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i-1,j)))**2 + &
          (0.5 * (div_xx_dy(i,j) + div_xx_dy(i,j-1)))**2)
        enddo ; enddo
        !do J=js-1,Jeq ; do I=is-1,Ieq
        do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
          grad_div_mag_q(i,j) =sqrt((0.5*(div_xx_dx(i,j) + div_xx_dx(i,j+1)))**2 + &
          (0.5 * (div_xx_dy(i,j) + div_xx_dy(i+1,j)))**2)
        enddo ; enddo

      else

        do j=jsq-1,jeq+2 ; do i=is-2,ieq+1
          div_xx_dx(i,j) = 0.0
        enddo ; enddo
        do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
          div_xx_dy(i,j) = 0.0
        enddo ; enddo
        do j=jsq,jeq+1 ; do i=isq,ieq+1
          grad_div_mag_h(i,j) = 0.0
        enddo ; enddo
        do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
          grad_div_mag_q(i,j) = 0.0
        enddo ; enddo

      endif ! CS%modified_Leith

      ! Add in beta for the Leith viscosity
      if (cs%use_beta_in_Leith) then
        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))
        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))
        enddo ; enddo
      endif ! CS%use_beta_in_Leith

      if (cs%use_QG_Leith_visc) then

        do j=jsq,jeq+1 ; do i=isq,ieq+1
          grad_vort_mag_h_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
                                         (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
        enddo ; enddo
        do j=js-1,jeq ; do i=is-1,ieq
          grad_vort_mag_q_2d(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
                                         (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
        enddo ; enddo

        ! This accumulates terms, some of which are in VarMix, so rescaling can not be done here.
        call calc_qg_leith_viscosity(varmix, g, gv, us, h, k, div_xx_dx, div_xx_dy, &
                                     vort_xy_dx, vort_xy_dy)

      endif

      do j=jsq,jeq+1 ; do i=isq,ieq+1
        grad_vort_mag_h(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2 + &
                                    (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2 )
      enddo ; enddo
      do j=js-1,jeq ; do i=is-1,ieq
        grad_vort_mag_q(i,j) = sqrt((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2 + &
                                    (0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2 )
      enddo ; enddo

    endif ! CS%Leith_Kh

    meke_res_fn = 1.

    do j=jsq,jeq+1 ; do i=isq,ieq+1
      if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then
        shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
          0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
                (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
        if (cs%use_QG_Leith_visc) then
          vert_vort_mag = min(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3.*grad_vort_mag_h_2d(i,j))
        else
          vert_vort_mag = (grad_vort_mag_h(i,j) + grad_div_mag_h(i,j))
        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)) / &
                            (h(i,j,k) + h_neglect) )
        visc_bound_rem = 1.0
      endif

      if (cs%Laplacian) then
        ! Determine the Laplacian viscosity at h points, using the
        ! largest value from several parameterizations.
        kh = cs%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity
        if (cs%add_LES_viscosity) then
          if (cs%Smagorinsky_Kh) kh = kh + cs%Laplac2_const_xx(i,j) * shear_mag
          if (cs%Leith_Kh) kh = kh + cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3
        else
          if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xx(i,j) * shear_mag )
          if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3)
        endif
        ! All viscosity contributions above are subject to resolution scaling
        if (rescale_kh) kh = varmix%Res_fn_h(i,j) * kh
        if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_h(i,j)
        ! Older method of bounding for stability
        if (legacy_bound) kh = min(kh, cs%Kh_Max_xx(i,j))
        kh = max( kh, cs%Kh_bg_min ) ! Place a floor on the viscosity, if desired.
        if (use_meke_ku) &
          kh = kh + meke%Ku(i,j) * meke_res_fn ! *Add* the MEKE contribution (might be negative)
        if (cs%anisotropic) kh = kh + cs%Kh_aniso * ( 1. - cs%n1n2_h(i,j)**2 ) ! *Add* the tension component
                                                                               ! of anisotropic viscosity

        ! Newer method of bounding for stability
        if (cs%better_bound_Kh) then
          if (kh >= hrat_min*cs%Kh_Max_xx(i,j)) then
            visc_bound_rem = 0.0
            kh = hrat_min*cs%Kh_Max_xx(i,j)
          else
            visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xx(i,j))
          endif
        endif

        if ((cs%id_Kh_h>0) .or. find_frictwork .or. cs%debug) kh_h(i,j,k) = kh

        if (cs%id_grid_Re_Kh>0) then
          ke = 0.125*((u(i,j,k)+u(i-1,j,k))**2 + (v(i,j,k)+v(i,j-1,k))**2)
          grid_re_kh(i,j,k) = (sqrt(ke) * sqrt(cs%grid_sp_h2(i,j))) &
              / max(kh, cs%min_grid_Kh)
        endif

        if (cs%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j)
        if (cs%id_sh_xx_h>0) sh_xx_h(i,j,k) = sh_xx(i,j)

        str_xx(i,j) = -kh * sh_xx(i,j)
      else   ! not Laplacian
        str_xx(i,j) = 0.0
      endif ! Laplacian

      if (cs%anisotropic) then
        ! Shearing-strain averaged to h-points
        local_strain = 0.25 * ( (sh_xy(i,j) + sh_xy(i-1,j-1)) + (sh_xy(i-1,j) + sh_xy(i,j-1)) )
        ! *Add* the shear-strain contribution to the xx-component of stress
        str_xx(i,j) = str_xx(i,j) - cs%Kh_aniso * cs%n1n2_h(i,j) * cs%n1n1_m_n2n2_h(i,j) * local_strain
      endif

      if (cs%biharmonic) then
        ! Determine the biharmonic viscosity at h points, using the
        ! largest value from several parameterizations.
        ahsm = 0.0; ahlth = 0.0
        if ((cs%Smagorinsky_Ah) .or. (cs%Leith_Ah)) then
          if (cs%Smagorinsky_Ah) then
            if (cs%bound_Coriolis) then
              ahsm = shear_mag * (cs%Biharm_const_xx(i,j) + &
                                  cs%Biharm_const2_xx(i,j)*shear_mag)
            else
              ahsm = cs%Biharm_const_xx(i,j) * shear_mag
            endif
          endif
          if (cs%Leith_Ah) ahlth = cs%Biharm6_const_xx(i,j) * abs(del2vort_h(i,j)) * inv_pi6
          ah = max(max(cs%Ah_bg_xx(i,j), ahsm), ahlth)
          if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
            ah = min(ah, cs%Ah_Max_xx(i,j))
        else
          ah = cs%Ah_bg_xx(i,j)
        endif ! Smagorinsky_Ah or Leith_Ah

        if (use_meke_au) ah = ah + meke%Au(i,j) ! *Add* the MEKE contribution

        if (cs%Re_Ah > 0.0) then
          ke = 0.125*((u(i,j,k)+u(i-1,j,k))**2 + (v(i,j,k)+v(i,j-1,k))**2)
          ah = sqrt(ke) * cs%Re_Ah_const_xx(i,j)
        endif

        if (cs%better_bound_Ah) then
          ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xx(i,j))
        endif

        if ((cs%id_Ah_h>0) .or. find_frictwork .or. cs%debug) ah_h(i,j,k) = ah

        if (cs%id_grid_Re_Ah>0) then
          ke = 0.125*((u(i,j,k)+u(i-1,j,k))**2 + (v(i,j,k)+v(i,j-1,k))**2)
          grid_re_ah(i,j,k) = (sqrt(ke) * cs%grid_sp_h3(i,j)) &
              / max(ah, cs%min_grid_Ah)
        endif

        str_xx(i,j) = str_xx(i,j) + ah * &
          (cs%DY_dxT(i,j) * (g%IdyCu(i,j)*del2u(i,j) - g%IdyCu(i-1,j)*del2u(i-1,j)) - &
           cs%DX_dyT(i,j) * (g%IdxCv(i,j)*del2v(i,j) - g%IdxCv(i,j-1)*del2v(i,j-1)))

        ! Keep a copy of the biharmonic contribution for backscatter parameterization
        bhstr_xx(i,j) = ah * &
          (cs%DY_dxT(i,j) * (g%IdyCu(i,j)*del2u(i,j) - g%IdyCu(i-1,j)*del2u(i-1,j)) - &
           cs%DX_dyT(i,j) * (g%IdxCv(i,j)*del2v(i,j) - g%IdxCv(i,j-1)*del2v(i,j-1)))
        bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))

      endif  ! biharmonic

    enddo ; enddo

    if (cs%biharmonic) then
      ! Gradient of Laplacian, for use in bi-harmonic term
      do j=js-1,jeq ; do i=is-1,ieq
        ddel2vdx(i,j) = cs%DY_dxBu(i,j)*(del2v(i+1,j)*g%IdyCv(i+1,j) - del2v(i,j)*g%IdyCv(i,j))
        ddel2udy(i,j) = cs%DX_dyBu(i,j)*(del2u(i,j+1)*g%IdxCu(i,j+1) - del2u(i,j)*g%IdxCu(i,j))
      enddo ; enddo
      ! Adjust contributions to shearing strain on open boundaries.
      if (apply_obc) then ; if (obc%zero_strain .or. obc%freeslip_strain) then
        do n=1,obc%number_of_segments
          j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB
          if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= jeq)) then
            do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
              if (obc%zero_strain) then
                ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.
              elseif (obc%freeslip_strain) then
                ddel2udy(i,j) = 0.
              endif
            enddo
          elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ieq)) then
            do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
              if (obc%zero_strain) then
                ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.
              elseif (obc%freeslip_strain) then
                ddel2vdx(i,j) = 0.
              endif
            enddo
          endif
        enddo
      endif ; endif
    endif

    meke_res_fn = 1.

    do j=js-1,jeq ; do i=is-1,ieq
      if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then
        shear_mag = sqrt(sh_xy(i,j)*sh_xy(i,j) + &
            0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + &
                  (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
        if (cs%use_QG_Leith_visc) then
          vert_vort_mag = min(grad_vort_mag_q(i,j) + grad_div_mag_q(i,j), 3.*grad_vort_mag_q_2d(i,j))
        else
          vert_vort_mag = (grad_vort_mag_q(i,j) + grad_div_mag_q(i,j))
        endif
      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)
      hq(i,j) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
              ((h_u(i,j) + h_u(i,j+1)) + (h_v(i,j) + h_v(i+1,j))))

      if (cs%better_bound_Ah .or. cs%better_bound_Kh) then
        hrat_min = min(1.0, min(h_u(i,j), h_u(i,j+1), h_v(i,j), h_v(i+1,j)) / &
                            (hq(i,j) + h_neglect) )
        visc_bound_rem = 1.0
      endif

      if (cs%no_slip .and. (g%mask2dBu(i,j) < 0.5)) then
        if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &
            (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0) then
          ! This is a coastal vorticity point, so modify hq and hrat_min.

          hu = g%mask2dCu(i,j) * h_u(i,j) + g%mask2dCu(i,j+1) * h_u(i,j+1)
          hv = g%mask2dCv(i,j) * h_v(i,j) + g%mask2dCv(i+1,j) * h_v(i+1,j)
          if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) * &
              (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) == 0.0) then
            ! Only one of hu and hv is nonzero, so just add them.
            hq(i,j) = hu + hv
            hrat_min = 1.0
          else
            ! Both hu and hv are nonzero, so take the harmonic mean.
            hq(i,j) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
            hrat_min = min(1.0, min(hu, hv) / (hq(i,j) + h_neglect) )
          endif
        endif
      endif

      if (cs%Laplacian) then
        ! Determine the Laplacian viscosity at q points, using the
        ! largest value from several parameterizations.
        kh = cs%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity
        if (cs%add_LES_viscosity) then
          if (cs%Smagorinsky_Kh) kh = kh + cs%Laplac2_const_xx(i,j) * shear_mag
          if (cs%Leith_Kh) kh = kh + cs%Laplac3_const_xx(i,j) * vert_vort_mag*inv_pi3
        else
          if (cs%Smagorinsky_Kh) kh = max( kh, cs%Laplac2_const_xy(i,j) * shear_mag )
          if (cs%Leith_Kh) kh = max( kh, cs%Laplac3_const_xy(i,j) * vert_vort_mag*inv_pi3)
        endif
        ! All viscosity contributions above are subject to resolution scaling
        if (rescale_kh) kh = varmix%Res_fn_q(i,j) * kh
        if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_q(i,j)
        ! Older method of bounding for stability
        if (legacy_bound) kh = min(kh, cs%Kh_Max_xy(i,j))
        kh = max( kh, cs%Kh_bg_min ) ! Place a floor on the viscosity, if desired.
        if (use_meke_ku) then ! *Add* the MEKE contribution (might be negative)
          kh = kh + 0.25*( (meke%Ku(i,j) + meke%Ku(i+1,j+1)) + &
                           (meke%Ku(i+1,j) + meke%Ku(i,j+1)) ) * meke_res_fn
        endif
        ! Older method of bounding for stability
        if (cs%anisotropic) kh = kh + cs%Kh_aniso * cs%n1n2_q(i,j)**2 ! *Add* the shear component
                                                                      ! of anisotropic viscosity

        ! Newer method of bounding for stability
        if (cs%better_bound_Kh) then
          if (kh >= hrat_min*cs%Kh_Max_xy(i,j)) then
            visc_bound_rem = 0.0
            kh = hrat_min*cs%Kh_Max_xy(i,j)
          elseif (cs%Kh_Max_xy(i,j)>0.) then
            visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xy(i,j))
          endif
        endif

        if (cs%id_Kh_q>0 .or. cs%debug) kh_q(i,j,k) = kh
        if (cs%id_vort_xy_q>0) vort_xy_q(i,j,k) = vort_xy(i,j)
        if (cs%id_sh_xy_q>0) sh_xy_q(i,j,k) = sh_xy(i,j)

        str_xy(i,j) = -kh * sh_xy(i,j)
      else   ! not Laplacian
        str_xy(i,j) = 0.0
      endif ! Laplacian

      if (cs%anisotropic) then
        ! Horizontal-tension averaged to q-points
        local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )
        ! *Add* the tension contribution to the xy-component of stress
        str_xy(i,j) = str_xy(i,j) - cs%Kh_aniso * cs%n1n2_q(i,j) * cs%n1n1_m_n2n2_q(i,j) * local_strain
      endif

      if (cs%biharmonic) then
      ! Determine the biharmonic viscosity at q points, using the
      ! largest value from several parameterizations.
        ahsm = 0.0 ; ahlth = 0.0
        if (cs%Smagorinsky_Ah .or. cs%Leith_Ah) then
          if (cs%Smagorinsky_Ah) then
            if (cs%bound_Coriolis) then
              ahsm = shear_mag * (cs%Biharm_const_xy(i,j) + &
                                  cs%Biharm_const2_xy(i,j)*shear_mag)
            else
              ahsm = cs%Biharm_const_xy(i,j) * shear_mag
            endif
          endif
          if (cs%Leith_Ah) ahlth = cs%Biharm6_const_xy(i,j) * abs(del2vort_q(i,j)) * inv_pi6
          ah = max(max(cs%Ah_bg_xy(i,j), ahsm), ahlth)
          if (cs%bound_Ah .and. .not.cs%better_bound_Ah) &
            ah = min(ah, cs%Ah_Max_xy(i,j))
        else
          ah = cs%Ah_bg_xy(i,j)
        endif ! Smagorinsky_Ah or Leith_Ah

        if (use_meke_au) then ! *Add* the MEKE contribution
          ah = ah + 0.25*( (meke%Au(i,j) + meke%Au(i+1,j+1)) +  &
                           (meke%Au(i+1,j) + meke%Au(i,j+1)) )
        endif

        if (cs%Re_Ah > 0.0) then
          ke = 0.125*((u(i,j,k)+u(i,j+1,k))**2 + (v(i,j,k)+v(i+1,j,k))**2)
          ah = sqrt(ke) * cs%Re_Ah_const_xy(i,j)
        endif

        if (cs%better_bound_Ah) then
          ah = min(ah, visc_bound_rem*hrat_min*cs%Ah_Max_xy(i,j))
        endif

        if (cs%id_Ah_q>0 .or. cs%debug) ah_q(i,j,k) = ah

        str_xy(i,j) = str_xy(i,j) + ah * ( ddel2vdx(i,j) + ddel2udy(i,j) )

        ! Keep a copy of the biharmonic contribution for backscatter parameterization
        bhstr_xy(i,j) = ah * ( ddel2vdx(i,j) + ddel2udy(i,j) ) * &
                        (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))

      endif  ! biharmonic

    enddo ; enddo

    if (cs%use_GME) then
      call thickness_diffuse_get_kh(td, kh_u_gme, kh_v_gme, g)
      call pass_vector(kh_u_gme, kh_v_gme, g%Domain)

      do j=jsq,jeq+1 ; do i=isq,ieq+1
        if (grad_vel_mag_bt_h(i,j)>0) then
          gme_coeff = cs%GME_efficiency * (min(g%bathyT(i,j)/cs%GME_h0,1.0)**2) * &
            (0.25*(kh_u_gme(i,j,k)+kh_u_gme(i-1,j,k)+kh_v_gme(i,j,k)+kh_v_gme(i,j-1,k)))
        else
          gme_coeff = 0.0
        endif

        ! apply mask
        gme_coeff = gme_coeff * boundary_mask_h(i,j)
        gme_coeff = min(gme_coeff, cs%GME_limiter)

        if ((cs%id_GME_coeff_h>0) .or. find_frictwork) gme_coeff_h(i,j,k) = gme_coeff
        str_xx_gme(i,j) = gme_coeff * sh_xx_bt(i,j)

      enddo ; enddo

      do j=js-1,jeq ; do i=is-1,ieq
        if (grad_vel_mag_bt_q(i,j)>0) then
          gme_coeff = cs%GME_efficiency * (min(g%bathyT(i,j)/cs%GME_h0,1.0)**2) * &
              (0.25*(kh_u_gme(i,j,k)+kh_u_gme(i,j+1,k)+kh_v_gme(i,j,k)+kh_v_gme(i+1,j,k)))
        else
          gme_coeff = 0.0
        endif

        ! apply mask
        gme_coeff = gme_coeff * boundary_mask_q(i,j)
        gme_coeff = min(gme_coeff, cs%GME_limiter)

        if (cs%id_GME_coeff_q>0) gme_coeff_q(i,j,k) = gme_coeff
        str_xy_gme(i,j) = gme_coeff * sh_xy_bt(i,j)

      enddo ; enddo

      ! Applying GME diagonal term.  This is linear and the arguments can be rescaled.
      call smooth_gme(cs,g,gme_flux_h=str_xx_gme)
      call smooth_gme(cs,g,gme_flux_q=str_xy_gme)

      do j=jsq,jeq+1 ; do i=isq,ieq+1
        str_xx(i,j) = (str_xx(i,j) + str_xx_gme(i,j)) * (h(i,j,k) * cs%reduction_xx(i,j))
      enddo ; enddo

      do j=js-1,jeq ; do i=is-1,ieq
        ! GME is applied below
        if (cs%no_slip) then
          str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * cs%reduction_xy(i,j))
        else
          str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
        endif
      enddo ; enddo

      if (associated(meke%GME_snk)) then
        do j=js,je ; do i=is,ie
          frictwork_gme(i,j,k) = gme_coeff_h(i,j,k) * h(i,j,k) * gv%H_to_kg_m2 * grad_vel_mag_bt_h(i,j)
        enddo ; enddo
      endif

    else ! use_GME
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
      enddo ; enddo

      do j=js-1,jeq ; do i=is-1,ieq
        if (cs%no_slip) then
          str_xy(i,j) = str_xy(i,j) * (hq(i,j) * cs%reduction_xy(i,j))
        else
          str_xy(i,j) = str_xy(i,j) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
        endif
      enddo ; enddo

    endif ! use_GME

    ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
    do j=js,je ; do i=isq,ieq
      diffu(i,j,k) = ((g%IdyCu(i,j)*(cs%dy2h(i,j) *str_xx(i,j) - &
                                     cs%dy2h(i+1,j)*str_xx(i+1,j)) + &
                       g%IdxCu(i,j)*(cs%dx2q(i,j-1)*str_xy(i,j-1) - &
                                     cs%dx2q(i,j) *str_xy(i,j))) * &
                     g%IareaCu(i,j)) / (h_u(i,j) + h_neglect)

    enddo ; enddo
    if (apply_obc) then
      ! This is not the right boundary condition. If all the masking of tendencies are done
      ! correctly later then eliminating this block should not change answers.
      do n=1,obc%number_of_segments
        if (obc%segment(n)%is_E_or_W) then
          i = obc%segment(n)%HI%IsdB
          do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
            diffu(i,j,k) = 0.
          enddo
        endif
      enddo
    endif

    ! Evaluate 1/h y.Div(h Grad u) or the biharmonic equivalent.
    do j=jsq,jeq ; do i=is,ie
      diffv(i,j,k) = ((g%IdyCv(i,j)*(cs%dy2q(i-1,j)*str_xy(i-1,j) - &
                                    cs%dy2q(i,j) *str_xy(i,j)) - &
                       g%IdxCv(i,j)*(cs%dx2h(i,j) *str_xx(i,j) - &
                                    cs%dx2h(i,j+1)*str_xx(i,j+1))) * &
                     g%IareaCv(i,j)) / (h_v(i,j) + h_neglect)
    enddo ; enddo
    if (apply_obc) then
      ! This is not the right boundary condition. If all the masking of tendencies are done
      ! correctly later then eliminating this block should not change answers.
      do n=1,obc%number_of_segments
        if (obc%segment(n)%is_N_or_S) then
          j = obc%segment(n)%HI%JsdB
          do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
            diffv(i,j,k) = 0.
          enddo
        endif
      enddo
    endif

    if (find_frictwork) then ; do j=js,je ; do i=is,ie
      ! Diagnose   str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
      ! This is the old formulation that includes energy diffusion
      frictwork(i,j,k) = gv%H_to_RZ * ( &
              (str_xx(i,j)*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j)     &
              -str_xx(i,j)*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j))    &
       +0.25*((str_xy(i,j)*(                                     &
                   (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j)            &
                  +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) )          &
              +str_xy(i-1,j-1)*(                                 &
                   (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1)    &
                  +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) )) &
             +(str_xy(i-1,j)*(                                   &
                   (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j)      &
                  +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) )        &
              +str_xy(i,j-1)*(                                   &
                   (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1)          &
                  +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
    enddo ; enddo ; endif

    ! Make a similar calculation as for FrictWork above but accumulating into
    ! the vertically integrated MEKE source term, and adjusting for any
    ! energy loss seen as a reduction in the (biharmonic) frictional source term.
    if (find_frictwork .and. associated(meke)) then ; if (associated(meke%mom_src)) then
      if (k==1) then
        do j=js,je ; do i=is,ie
          meke%mom_src(i,j) = 0.
        enddo ; enddo
        if (associated(meke%GME_snk)) then
          do j=js,je ; do i=is,ie
            meke%GME_snk(i,j) = 0.
          enddo ; enddo
        endif
      endif
      if (meke%backscatter_Ro_c /= 0.) then
        do j=js,je ; do i=is,ie
          fath = 0.25*( (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
                        (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))) )
          shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + &
            0.25*((sh_xy(i-1,j-1)*sh_xy(i-1,j-1) + sh_xy(i,j)*sh_xy(i,j)) + &
                  (sh_xy(i-1,j)*sh_xy(i-1,j) + sh_xy(i,j-1)*sh_xy(i,j-1))))
          if (cs%answers_2018) then
            fath = (us%s_to_T*fath)**meke%backscatter_Ro_pow ! f^n
            ! Note the hard-coded dimensional constant in the following line that can not
            ! be rescaled for dimensional consistency.
            shear_mag = ( ( (us%s_to_T*shear_mag)**meke%backscatter_Ro_pow ) + 1.e-30 ) &
                        * meke%backscatter_Ro_c ! c * D^n
            ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n)
            ! RoScl = 1 - g(Ro)
            roscl = shear_mag / ( fath + shear_mag ) ! = 1 - f^n/(f^n+c*D^n)
          else
            if (fath <= backscat_subround*shear_mag) then
              roscl = 1.0
            else
              sh_f_pow = meke%backscatter_Ro_c * (shear_mag / fath)**meke%backscatter_Ro_pow
              roscl = sh_f_pow / (1.0 + sh_f_pow) ! = 1 - f^n/(f^n+c*D^n)
            endif
          endif


          meke%mom_src(i,j) = meke%mom_src(i,j) + gv%H_to_RZ * ( &
                ((str_xx(i,j)-roscl*bhstr_xx(i,j))*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j)  &
                -(str_xx(i,j)-roscl*bhstr_xx(i,j))*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j)) &
         +0.25*(((str_xy(i,j)-roscl*bhstr_xy(i,j))*(                                  &
                     (u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j)                               &
                    +(v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j) )                             &
                +(str_xy(i-1,j-1)-roscl*bhstr_xy(i-1,j-1))*(                          &
                     (u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1)                       &
                    +(v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1) ))                    &
                +((str_xy(i-1,j)-roscl*bhstr_xy(i-1,j))*(                             &
                     (u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j)                         &
                    +(v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j) )                           &
                +(str_xy(i,j-1)-roscl*bhstr_xy(i,j-1))*(                              &
                     (u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1)                             &
                    +(v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1) )) ) )
        enddo ; enddo
      endif ! MEKE%backscatter_Ro_c

      do j=js,je ; do i=is,ie
        meke%mom_src(i,j) = meke%mom_src(i,j) + frictwork(i,j,k)
      enddo ; enddo

      if (cs%use_GME .and. associated(meke)) then ; if (associated(meke%GME_snk)) then
        do j=js,je ; do i=is,ie
          meke%GME_snk(i,j) = meke%GME_snk(i,j) + frictwork_gme(i,j,k)
        enddo ; enddo
      endif ; endif

    endif ; endif ! find_FrictWork and associated(mom_src)

  enddo ! end of k loop

  ! Offer fields for diagnostic averaging.
  if (cs%id_diffu>0)     call post_data(cs%id_diffu, diffu, cs%diag)
  if (cs%id_diffv>0)     call post_data(cs%id_diffv, diffv, cs%diag)
  if (cs%id_FrictWork>0) call post_data(cs%id_FrictWork, frictwork, cs%diag)
  if (cs%id_FrictWork_GME>0) call post_data(cs%id_FrictWork_GME, frictwork_gme, cs%diag)
  if (cs%id_Ah_h>0)      call post_data(cs%id_Ah_h, ah_h, cs%diag)
  if (cs%id_grid_Re_Ah>0) call post_data(cs%id_grid_Re_Ah, grid_re_ah, cs%diag)
  if (cs%id_div_xx_h>0)  call post_data(cs%id_div_xx_h, div_xx_h, cs%diag)
  if (cs%id_vort_xy_q>0) call post_data(cs%id_vort_xy_q, vort_xy_q, cs%diag)
  if (cs%id_sh_xx_h>0)  call post_data(cs%id_sh_xx_h, sh_xx_h, cs%diag)
  if (cs%id_sh_xy_q>0) call post_data(cs%id_sh_xy_q, sh_xy_q, cs%diag)
  if (cs%id_Ah_q>0)      call post_data(cs%id_Ah_q, ah_q, cs%diag)
  if (cs%id_Kh_h>0)      call post_data(cs%id_Kh_h, kh_h, cs%diag)
  if (cs%id_grid_Re_Kh>0) call post_data(cs%id_grid_Re_Kh, grid_re_kh, cs%diag)
  if (cs%id_Kh_q>0)      call post_data(cs%id_Kh_q, kh_q, cs%diag)
  if (cs%id_GME_coeff_h > 0)  call post_data(cs%id_GME_coeff_h, gme_coeff_h, cs%diag)
  if (cs%id_GME_coeff_q > 0)  call post_data(cs%id_GME_coeff_q, gme_coeff_q, cs%diag)

  if (cs%debug) then
    if (cs%Laplacian) then
      call hchksum(kh_h, "Kh_h", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
      call bchksum(kh_q, "Kh_q", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
    endif
    if (cs%biharmonic) call hchksum(ah_h, "Ah_h", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
    if (cs%biharmonic) call bchksum(ah_q, "Ah_q", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
  endif

  if (cs%id_FrictWorkIntz > 0) then
    do j=js,je
      do i=is,ie ; frictworkintz(i,j) = frictwork(i,j,1) ; enddo
      do k=2,nz ; do i=is,ie
        frictworkintz(i,j) = frictworkintz(i,j) + frictwork(i,j,k)
      enddo ; enddo
    enddo
    call post_data(cs%id_FrictWorkIntz, frictworkintz, cs%diag)
  endif

  ! Diagnostics for terms multiplied by fractional thicknesses

  ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
  ! The code is retained for degugging purposes in the future.
  !if (present(ADp) .and. (CS%id_hf_diffu > 0)) then
  !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
  !    CS%hf_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hfrac_u(I,j,k)
  !  enddo ; enddo ; enddo
  !  call post_data(CS%id_hf_diffu, CS%hf_diffu, CS%diag)
  !endif
  !if (present(ADp) .and. (CS%id_hf_diffv > 0)) then
  !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
  !    CS%hf_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hfrac_v(i,J,k)
  !  enddo ; enddo ; enddo
  !  call post_data(CS%id_hf_diffv, CS%hf_diffv, CS%diag)
  !endif
  if (present(adp) .and. (cs%id_hf_diffu_2d > 0)) then
    allocate(hf_diffu_2d(g%IsdB:g%IedB,g%jsd:g%jed))
    hf_diffu_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=isq,ieq
      hf_diffu_2d(i,j) = hf_diffu_2d(i,j) + diffu(i,j,k) * adp%diag_hfrac_u(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_diffu_2d, hf_diffu_2d, cs%diag)
    deallocate(hf_diffu_2d)
  endif
  if (present(adp) .and. (cs%id_hf_diffv_2d > 0)) then
    allocate(hf_diffv_2d(g%isd:g%ied,g%JsdB:g%JedB))
    hf_diffv_2d(:,:) = 0.0
    do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      hf_diffv_2d(i,j) = hf_diffv_2d(i,j) + diffv(i,j,k) * adp%diag_hfrac_v(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_diffv_2d, hf_diffv_2d, cs%diag)
    deallocate(hf_diffv_2d)
  endif

end subroutine horizontal_viscosity

!> Allocates space for and calculates static variables used by horizontal_viscosity().
!! hor_visc_init calculates and stores the values of a number of metric functions that
!! are used in horizontal_viscosity().
subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp)
  type(time_type),         intent(in)    :: Time !< Current model time.
  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.
  type(hor_visc_cs), pointer             :: CS   !< Pointer to the control structure for this module
  type(meke_type), pointer               :: MEKE !< MEKE data
  type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers
  ! Local variables
  real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v
  real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v
                ! u0v is the Laplacian sensitivities to the v velocities
                ! at u points [L-2 ~> m-2], with u0u, v0u, and v0v defined similarly.
  real :: grid_sp_h2       ! Harmonic mean of the squares of the grid [L2 ~> m2]
  real :: grid_sp_h3       ! Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
  real :: grid_sp_q2       ! spacings at h and q points [L2 ~> m2]
  real :: grid_sp_q3       ! spacings at h and q points^(3/2) [L3 ~> m3]
  real :: min_grid_sp_h2   ! Minimum value of grid_sp_h2 [L2 ~> m2]
  real :: min_grid_sp_h4   ! Minimum value of grid_sp_h2**2 [L4 ~> m4]
  real :: Kh_Limit         ! A coefficient [T-1 ~> s-1] used, along with the
                           ! grid spacing, to limit Laplacian viscosity.
  real :: fmax             ! maximum absolute value of f at the four
                           ! vorticity points around a thickness point [T-1 ~> s-1]
  real :: BoundCorConst    ! A constant used when using viscosity to bound the Coriolis accelerations
                           ! [T2 L-2 ~> s2 m-2]
  real :: Ah_Limit         ! coefficient [T-1 ~> s-1] used, along with the
                           ! grid spacing, to limit biharmonic viscosity
  real :: Kh               ! Lapacian horizontal viscosity [L2 T-1 ~> m2 s-1]
  real :: Ah               ! biharmonic horizontal viscosity [L4 T-1 ~> m4 s-1]
  real :: Kh_vel_scale     ! this speed [L T-1 ~> m s-1] times grid spacing gives Lap visc
  real :: Ah_vel_scale     ! this speed [L T-1 ~> m s-1] times grid spacing cubed gives bih visc
  real :: Ah_time_scale    ! damping time-scale for biharmonic visc [T ~> s]
  real :: Smag_Lap_const   ! nondimensional Laplacian Smagorinsky constant
  real :: Smag_bi_const    ! nondimensional biharmonic Smagorinsky constant
  real :: Leith_Lap_const  ! nondimensional Laplacian Leith constant
  real :: Leith_bi_const   ! nondimensional biharmonic Leith constant
  real :: dt               ! The dynamics time step [T ~> s]
  real :: Idt              ! The inverse of dt [T-1 ~> s-1]
  real :: denom            ! work variable; the denominator of a fraction
  real :: maxvel           ! largest permitted velocity components [m s-1]
  real :: bound_Cor_vel    ! grid-scale velocity variations at which value
                           ! the quadratically varying biharmonic viscosity
                           ! balances Coriolis acceleration [L T-1 ~> m s-1]
  real :: Kh_sin_lat       ! Amplitude of latitudinally dependent viscosity [L2 T-1 ~> m2 s-1]
  real :: Kh_pwr_of_sine   ! Power used to raise sin(lat) when using Kh_sin_lat
  logical :: bound_Cor_def ! parameter setting of BOUND_CORIOLIS
  logical :: get_all       ! If true, read and log all parameters, regardless of
                           ! whether they are used, to enable spell-checking of
                           ! valid parameters.
  logical :: split         ! If true, use the split time stepping scheme.
                           ! If false and USE_GME = True, issue a FATAL error.
  logical :: default_2018_answers
  character(len=64) :: inputdir, filename
  real    :: deg2rad       ! Converts degrees to radians
  real    :: slat_fn       ! sin(lat)**Kh_pwr_of_sine
  real    :: aniso_grid_dir(2) ! Vector (n1,n2) for anisotropic direction
  integer :: aniso_mode    ! Selects the mode for setting the anisotropic direction
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  integer :: i, j
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_hor_visc"  ! module name
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  if (associated(cs)) then
    call mom_error(warning, "hor_visc_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)
  cs%diag => diag
  ! Read parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  !   It is not clear whether these initialization lines are needed for the
  ! cases where the corresponding parameters are not read.
  cs%bound_Kh = .false. ; cs%better_bound_Kh = .false. ; cs%Smagorinsky_Kh = .false. ; cs%Leith_Kh = .false.
  cs%bound_Ah = .false. ; cs%better_bound_Ah = .false. ; cs%Smagorinsky_Ah = .false. ; cs%Leith_Ah = .false.
  cs%use_QG_Leith_visc = .false.
  cs%bound_Coriolis = .false.
  cs%Modified_Leith = .false.
  cs%anisotropic = .false.
  cs%dynamic_aniso = .false.
  kh = 0.0 ; ah = 0.0
  !   If GET_ALL_PARAMS is true, all parameters are read in all cases to enable
  ! parameter spelling checks.
  call get_param(param_file, mdl, "GET_ALL_PARAMS", get_all, default=.false.)
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "HOR_VISC_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
  call get_param(param_file, mdl, "LAPLACIAN", cs%Laplacian, &
                 "If true, use a Laplacian horizontal viscosity.", &
                 default=.false.)
  if (cs%Laplacian .or. get_all) then
    call get_param(param_file, mdl, "KH", kh,                      &
                 "The background Laplacian horizontal viscosity.", &
                 units = "m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
    call get_param(param_file, mdl, "KH_BG_MIN", cs%Kh_bg_min, &
                 "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
                 units = "m2 s-1",  default=0.0, scale=us%m_to_L**2*us%T_to_s)
    call get_param(param_file, mdl, "KH_VEL_SCALE", kh_vel_scale, &
                 "The velocity scale which is multiplied by the grid "//&
                 "spacing to calculate the Laplacian viscosity. "//&
                 "The final viscosity is the largest of this scaled "//&
                 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
                 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
    call get_param(param_file, mdl, "KH_SIN_LAT", kh_sin_lat, &
                 "The amplitude of a latitudinally-dependent background "//&
                 "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &
                 units = "m2 s-1",  default=0.0, scale=us%m_to_L**2*us%T_to_s)
    if (kh_sin_lat>0. .or. get_all) &
      call get_param(param_file, mdl, "KH_PWR_OF_SINE", kh_pwr_of_sine, &
                 "The power used to raise SIN(LAT) when using a latitudinally "//&
                 "dependent background viscosity.", &
                 units = "nondim",  default=4.0)
    call get_param(param_file, mdl, "SMAGORINSKY_KH", cs%Smagorinsky_Kh, &
                 "If true, use a Smagorinsky nonlinear eddy viscosity.", &
                 default=.false.)
    if (cs%Smagorinsky_Kh .or. get_all) &
      call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_lap_const, &
                 "The nondimensional Laplacian Smagorinsky constant, "//&
                 "often 0.15.", units="nondim", default=0.0, &
                  fail_if_missing = cs%Smagorinsky_Kh)
    call get_param(param_file, mdl, "LEITH_KH", cs%Leith_Kh, &
                 "If true, use a Leith nonlinear eddy viscosity.", &
                 default=.false.)
    call get_param(param_file, mdl, "MODIFIED_LEITH", cs%Modified_Leith, &
                 "If true, add a term to Leith viscosity which is "//&
                 "proportional to the gradient of divergence.", &
                 default=.false.)
    call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", cs%res_scale_MEKE, &
                 "If true, the viscosity contribution from MEKE is scaled by "//&
                 "the resolution function.", default=.false.)
    if (cs%Leith_Kh .or. get_all) then
      call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
                 "The nondimensional Laplacian Leith constant, "//&
                 "often set to 1.0", units="nondim", default=0.0, &
                  fail_if_missing = cs%Leith_Kh)
      call get_param(param_file, mdl, "USE_QG_LEITH_VISC", cs%use_QG_Leith_visc, &
                 "If true, use QG Leith nonlinear eddy viscosity.", &
                 default=.false.)
      if (cs%use_QG_Leith_visc .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_VISC=True.")
    endif
    if (cs%Leith_Kh .or. cs%Leith_Ah .or. get_all) then
      call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_Leith, &
                 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
                 default=cs%Leith_Kh)
      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.)
    endif
    call get_param(param_file, mdl, "BOUND_KH", cs%bound_Kh, &
                 "If true, the Laplacian coefficient is locally limited "//&
                 "to be stable.", default=.true.)
    call get_param(param_file, mdl, "BETTER_BOUND_KH", cs%better_bound_Kh, &
                 "If true, the Laplacian coefficient is locally limited "//&
                 "to be stable with a better bounding than just BOUND_KH.", &
                 default=cs%bound_Kh)
    call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", cs%anisotropic, &
                 "If true, allow anistropic viscosity in the Laplacian "//&
                 "horizontal viscosity.", default=.false.)
    call get_param(param_file, mdl, "ADD_LES_VISCOSITY", cs%add_LES_viscosity, &
                 "If true, adds the viscosity from Smagorinsky and Leith to the "//&
                 "background viscosity instead of taking the maximum.", default=.false.)
  endif
  if (cs%anisotropic .or. get_all) then
    call get_param(param_file, mdl, "KH_ANISO", cs%Kh_aniso, &
                 "The background Laplacian anisotropic horizontal viscosity.", &
                 units = "m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
    call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, &
                 "Selects the mode for setting the direction of anistropy.\n"//&
                 "\t 0 - Points along the grid i-direction.\n"//&
                 "\t 1 - Points towards East.\n"//&
                 "\t 2 - Points along the flow direction, U/|U|.", &
                 default=0)
    select case (aniso_mode)
      case (0)
        call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
                 "The vector pointing in the direction of anistropy for "//&
                 "horizont viscosity. n1,n2 are the i,j components relative "//&
                 "to the grid.", units = "nondim", fail_if_missing=.true.)
      case (1)
        call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
                 "The vector pointing in the direction of anistropy for "//&
                 "horizont viscosity. n1,n2 are the i,j components relative "//&
                 "to the spherical coordinates.", units = "nondim", fail_if_missing=.true.)
    end select
  endif
  call get_param(param_file, mdl, "BIHARMONIC", cs%biharmonic, &
                 "If true, use a biharmonic horizontal viscosity. "//&
                 "BIHARMONIC may be used with LAPLACIAN.", &
                 default=.true.)
  if (cs%biharmonic .or. get_all) then
    call get_param(param_file, mdl, "AH", ah, &
                 "The background biharmonic horizontal viscosity.", &
                 units = "m4 s-1", default=0.0, scale=us%m_to_L**4*us%T_to_s)
    call get_param(param_file, mdl, "AH_VEL_SCALE", ah_vel_scale, &
                 "The velocity scale which is multiplied by the cube of "//&
                 "the grid spacing to calculate the biharmonic viscosity. "//&
                 "The final viscosity is the largest of this scaled "//&
                 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
                 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
    call get_param(param_file, mdl, "AH_TIME_SCALE", ah_time_scale, &
                 "A time scale whose inverse is multiplied by the fourth "//&
                 "power of the grid spacing to calculate biharmonic viscosity. "//&
                 "The final viscosity is the largest of all viscosity "//&
                 "formulations in use. 0.0 means that it's not used.", &
                 units="s", default=0.0, scale=us%s_to_T)
    call get_param(param_file, mdl, "SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
                 "If true, use a biharmonic Smagorinsky nonlinear eddy "//&
                 "viscosity.", default=.false.)
    call get_param(param_file, mdl, "LEITH_AH", cs%Leith_Ah, &
                 "If true, use a biharmonic Leith nonlinear eddy "//&
                 "viscosity.", default=.false.)
    call get_param(param_file, mdl, "BOUND_AH", cs%bound_Ah, &
                 "If true, the biharmonic coefficient is locally limited "//&
                 "to be stable.", default=.true.)
    call get_param(param_file, mdl, "BETTER_BOUND_AH", cs%better_bound_Ah, &
                 "If true, the biharmonic coefficient is locally limited "//&
                 "to be stable with a better bounding than just BOUND_AH.", &
                 default=cs%bound_Ah)
    call get_param(param_file, mdl, "RE_AH", cs%Re_Ah, &
                 "If nonzero, the biharmonic coefficient is scaled "//&
                 "so that the biharmonic Reynolds number is equal to this.", &
                 units="nondim", default=0.0)

    if (cs%Smagorinsky_Ah .or. get_all) then
      call get_param(param_file, mdl, "SMAG_BI_CONST",smag_bi_const, &
                 "The nondimensional biharmonic Smagorinsky constant, "//&
                 "typically 0.015 - 0.06.", units="nondim", default=0.0, &
                 fail_if_missing = cs%Smagorinsky_Ah)
      call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_cor_def, default=.false.)
      call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &
                 "If true use a viscosity that increases with the square "//&
                 "of the velocity shears, so that the resulting viscous "//&
                 "drag is of comparable magnitude to the Coriolis terms "//&
                 "when the velocity differences between adjacent grid "//&
                 "points is 0.5*BOUND_CORIOLIS_VEL.  The default is the "//&
                 "value of BOUND_CORIOLIS (or false).", default=bound_cor_def)
      if (cs%bound_Coriolis .or. get_all) then
        call get_param(param_file, mdl, "MAXVEL", maxvel, default=3.0e8)
        bound_cor_vel = maxvel
        call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_cor_vel, &
                 "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&
                 "the biharmonic drag to have comparable magnitude to the "//&
                 "Coriolis acceleration.  The default is set by MAXVEL.", &
                 units="m s-1", default=maxvel, scale=us%m_s_to_L_T)
      endif
    endif
    if (cs%Leith_Ah .or. get_all) &
      call get_param(param_file, mdl, "LEITH_BI_CONST", leith_bi_const, &
                 "The nondimensional biharmonic Leith constant, "//&
                 "typical values are thus far undetermined.", units="nondim", default=0.0, &
                 fail_if_missing = cs%Leith_Ah)
  endif
  call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", cs%use_land_mask, &
                 "If true, use Use the land mask for the computation of thicknesses "//&
                 "at velocity locations. This eliminates the dependence on arbitrary "//&
                 "values over land or outside of the domain.", default=.true.)
  if (cs%better_bound_Ah .or. cs%better_bound_Kh .or. get_all) &
    call get_param(param_file, mdl, "HORVISC_BOUND_COEF", cs%bound_coef, &
                 "The nondimensional coefficient of the ratio of the "//&
                 "viscosity bounds to the theoretical maximum for "//&
                 "stability without considering other terms.", units="nondim", &
                 default=0.8)
  call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
                 "If true, no slip boundary conditions are used; otherwise "//&
                 "free slip boundary conditions are assumed. The "//&
                 "implementation of the free slip BCs on a C-grid is much "//&
                 "cleaner than the no slip BCs. The use of free slip BCs "//&
                 "is strongly encouraged, and no slip BCs are not used with "//&
                 "the biharmonic viscosity.", default=.false.)
  call get_param(param_file, mdl, "USE_KH_BG_2D", cs%use_Kh_bg_2d, &
                 "If true, read a file containing 2-d background harmonic "//&
                 "viscosities. The final viscosity is the maximum of the other "//&
                 "terms and this background value.", default=.false.)

  call get_param(param_file, mdl, "USE_GME", cs%use_GME, &
                 "If true, use the GM+E backscatter scheme in association \n"//&
                 "with the Gent and McWilliams parameterization.", default=.false.)
  if (cs%use_GME) then
    call get_param(param_file, mdl, "SPLIT", split, &
                 "Use the split time stepping if true.", default=.true., &
                  do_not_log=.true.)
    if (.not. split) call mom_error(fatal,"ERROR: Currently, USE_GME = True "// &
                                           "cannot be used with SPLIT=False.")
    call get_param(param_file, mdl, "GME_H0", cs%GME_h0, &
                 "The strength of GME tapers quadratically to zero when the bathymetric "//&
                 "depth is shallower than GME_H0.", units="m", scale=us%m_to_Z, &
                 default=1000.0)
    call get_param(param_file, mdl, "GME_EFFICIENCY", cs%GME_efficiency, &
                 "The nondimensional prefactor multiplying the GME coefficient.", &
                 units="nondim", default=1.0)
    call get_param(param_file, mdl, "GME_LIMITER", cs%GME_limiter, &
                 "The absolute maximum value the GME coefficient is allowed to take.", &
                 units="m2 s-1", scale=us%m_to_L**2*us%T_to_s, default=1.0e7)
  endif
  if (cs%Laplacian .or. cs%biharmonic) then
    call get_param(param_file, mdl, "DT", dt, &
                 "The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T, &
                 fail_if_missing=.true.)
    idt = 1.0 / dt
  endif
  if (cs%no_slip .and. cs%biharmonic) &
    call mom_error(fatal,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
                         "at the same time in MOM.")
  if (.not.(cs%Laplacian .or. cs%biharmonic)) then
    ! Only issue inviscid warning if not in single column mode (usually 2x2 domain)
    if ( max(g%domain%niglobal, g%domain%njglobal)>2 ) call mom_error(warning, &
      "hor_visc_init:  It is usually a very bad idea not to use either "//&
      "LAPLACIAN or BIHARMONIC viscosity.")
    return ! We are not using either Laplacian or Bi-harmonic lateral viscosity
  endif
  deg2rad = atan(1.0) / 45.
  alloc_(cs%dx2h(isd:ied,jsd:jed))        ; cs%dx2h(:,:)    = 0.0
  alloc_(cs%dy2h(isd:ied,jsd:jed))        ; cs%dy2h(:,:)    = 0.0
  alloc_(cs%dx2q(isdb:iedb,jsdb:jedb))    ; cs%dx2q(:,:)    = 0.0
  alloc_(cs%dy2q(isdb:iedb,jsdb:jedb))    ; cs%dy2q(:,:)    = 0.0
  alloc_(cs%dx_dyT(isd:ied,jsd:jed))      ; cs%dx_dyT(:,:)  = 0.0
  alloc_(cs%dy_dxT(isd:ied,jsd:jed))      ; cs%dy_dxT(:,:)  = 0.0
  alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0
  alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0
  if (cs%Laplacian) then
    alloc_(cs%grid_sp_h2(isd:ied,jsd:jed))   ; cs%grid_sp_h2(:,:) = 0.0
    alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed))     ; cs%Kh_bg_xx(:,:) = 0.0
    alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0
    if (cs%bound_Kh .or. cs%better_bound_Kh) then
      alloc_(cs%Kh_Max_xx(isd:ied,jsd:jed)) ; cs%Kh_Max_xx(:,:) = 0.0
      alloc_(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_Max_xy(:,:) = 0.0
    endif
    if (cs%Smagorinsky_Kh) then
      alloc_(cs%Laplac2_const_xx(isd:ied,jsd:jed))     ; cs%Laplac2_const_xx(:,:) = 0.0
      alloc_(cs%Laplac2_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac2_const_xy(:,:) = 0.0
    endif
    if (cs%Leith_Kh) then
      alloc_(cs%Laplac3_const_xx(isd:ied,jsd:jed)) ; cs%Laplac3_const_xx(:,:) = 0.0
      alloc_(cs%Laplac3_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac3_const_xy(:,:) = 0.0
    endif
  endif
  alloc_(cs%reduction_xx(isd:ied,jsd:jed))     ; cs%reduction_xx(:,:) = 0.0
  alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0
  if (cs%anisotropic) then
    alloc_(cs%n1n2_h(isd:ied,jsd:jed)) ; cs%n1n2_h(:,:) = 0.0
    alloc_(cs%n1n1_m_n2n2_h(isd:ied,jsd:jed)) ; cs%n1n1_m_n2n2_h(:,:) = 0.0
    alloc_(cs%n1n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n2_q(:,:) = 0.0
    alloc_(cs%n1n1_m_n2n2_q(isdb:iedb,jsdb:jedb)) ; cs%n1n1_m_n2n2_q(:,:) = 0.0
    select case (aniso_mode)
      case (0)
        call align_aniso_tensor_to_grid(cs, aniso_grid_dir(1), aniso_grid_dir(2))
      case (1)
      ! call align_aniso_tensor_to_grid(CS, aniso_grid_dir(1), aniso_grid_dir(2))
      case (2)
        cs%dynamic_aniso = .true.
      case default
        call mom_error(fatal, "MOM_hor_visc: "//&
             "Runtime parameter ANISOTROPIC_MODE is out of range.")
    end select
  endif
  if (cs%use_Kh_bg_2d) then
    alloc_(cs%Kh_bg_2d(isd:ied,jsd:jed))     ; cs%Kh_bg_2d(:,:) = 0.0
    call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, &
                 'The filename containing a 2d map of "Kh".', &
                 default='KH_background_2d.nc')
    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
    inputdir = slasher(inputdir)
    call mom_read_data(trim(inputdir)//trim(filename), 'Kh', cs%Kh_bg_2d, &
                       g%domain, timelevel=1, scale=us%m_to_L**2*us%T_to_s)
    call pass_var(cs%Kh_bg_2d, g%domain)
  endif
  if (cs%biharmonic) then
    alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0
    alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0
    alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed))  ; cs%Idxdy2u(:,:)  = 0.0
    alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb))  ; cs%Idxdy2v(:,:)  = 0.0
    alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed))     ; cs%Ah_bg_xx(:,:) = 0.0
    alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0
    alloc_(cs%grid_sp_h3(isd:ied,jsd:jed))   ; cs%grid_sp_h3(:,:) = 0.0
    if (cs%bound_Ah .or. cs%better_bound_Ah) then
      alloc_(cs%Ah_Max_xx(isd:ied,jsd:jed))     ; cs%Ah_Max_xx(:,:) = 0.0
      alloc_(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_Max_xy(:,:) = 0.0
    endif
    if (cs%Smagorinsky_Ah) then
      alloc_(cs%Biharm_const_xx(isd:ied,jsd:jed))     ; cs%Biharm_const_xx(:,:) = 0.0
      alloc_(cs%Biharm_const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const_xy(:,:) = 0.0
      if (cs%bound_Coriolis) then
        alloc_(cs%Biharm_const2_xx(isd:ied,jsd:jed))     ; cs%Biharm_const2_xx(:,:) = 0.0
        alloc_(cs%Biharm_const2_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm_const2_xy(:,:) = 0.0
      endif
    endif
    if (cs%Leith_Ah) then
        alloc_(cs%biharm6_const_xx(isd:ied,jsd:jed)) ; cs%biharm6_const_xx(:,:) = 0.0
        alloc_(cs%biharm6_const_xy(isdb:iedb,jsdb:jedb)) ; cs%biharm6_const_xy(:,:) = 0.0
    endif
    if (cs%Re_Ah > 0.0) then
      alloc_(cs%Re_Ah_const_xx(isd:ied,jsd:jed)); cs%Re_Ah_const_xx(:,:) = 0.0
      alloc_(cs%Re_Ah_const_xy(isdb:iedb,jsdb:jedb)); cs%Re_Ah_const_xy(:,:) = 0.0
    endif
  endif
  do j=js-2,jeq+1 ; do i=is-2,ieq+1
    cs%dx2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%dy2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)
    cs%DX_dyBu(i,j) = g%dxBu(i,j)*g%IdyBu(i,j) ; cs%DY_dxBu(i,j) = g%dyBu(i,j)*g%IdxBu(i,j)
  enddo ; enddo
  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
    cs%dx2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%dy2h(i,j) = g%dyT(i,j)*g%dyT(i,j)
    cs%DX_dyT(i,j) = g%dxT(i,j)*g%IdyT(i,j) ; cs%DY_dxT(i,j) = g%dyT(i,j)*g%IdxT(i,j)
  enddo ; enddo
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    cs%reduction_xx(i,j) = 1.0
    if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
        (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &
      cs%reduction_xx(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
    if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &
        (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &
      cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / (g%dyCu(i-1,j))
    if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
        (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &
      cs%reduction_xx(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
    if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &
        (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &
      cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / (g%dxCv(i,j-1))
  enddo ; enddo
  do j=js-1,jeq ; do i=is-1,ieq
    cs%reduction_xy(i,j) = 1.0
    if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &
        (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &
      cs%reduction_xy(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))
    if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &
        (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &
      cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / (g%dyCu(i,j+1))
    if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &
        (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &
      cs%reduction_xy(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))
    if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &
        (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &
      cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / (g%dxCv(i+1,j))
  enddo ; enddo
  if (cs%Laplacian) then
   ! The 0.3 below was 0.4 in MOM1.10.  The change in hq requires
   ! this to be less than 1/3, rather than 1/2 as before.
    if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)
    ! Calculate and store the background viscosity at h-points

    min_grid_sp_h2 = huge(1.)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      ! Static factors in the Smagorinsky and Leith schemes
      grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j) + cs%dy2h(i,j))
      cs%grid_sp_h2(i,j) = grid_sp_h2
      grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
      if (cs%Smagorinsky_Kh) cs%Laplac2_const_xx(i,j) = smag_lap_const * grid_sp_h2
      if (cs%Leith_Kh)       cs%Laplac3_const_xx(i,j) = leith_lap_const * grid_sp_h3
      ! Maximum of constant background and MICOM viscosity
      cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))
      ! Use the larger of the above and values read from a file
      if (cs%use_Kh_bg_2d) cs%Kh_bg_xx(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xx(i,j))
      ! Use the larger of the above and a function of sin(latitude)
      if (kh_sin_lat>0.) then
        slat_fn = abs( sin( deg2rad * g%geoLatT(i,j) ) ) ** kh_pwr_of_sine
        cs%Kh_bg_xx(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xx(i,j))
      endif
      if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
        ! Limit the background viscosity to be numerically stable
        cs%Kh_Max_xx(i,j) = kh_limit * grid_sp_h2
        cs%Kh_bg_xx(i,j) = min(cs%Kh_bg_xx(i,j), cs%Kh_Max_xx(i,j))
      endif
      min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2)
    enddo ; enddo
    call min_across_pes(min_grid_sp_h2)

    ! Calculate and store the background viscosity at q-points
    do j=js-1,jeq ; do i=is-1,ieq
      ! Static factors in the Smagorinsky and Leith schemes
      grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j) + cs%dy2q(i,j))
      grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
      if (cs%Smagorinsky_Kh) cs%Laplac2_const_xy(i,j) = smag_lap_const * grid_sp_q2
      if (cs%Leith_Kh)       cs%Laplac3_const_xy(i,j) = leith_lap_const * grid_sp_q3
      ! Maximum of constant background and MICOM viscosity
      cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))
      ! Use the larger of the above and values read from a file
      if (cs%use_Kh_bg_2d) then
        cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_xy(i,j), &
            0.25*((cs%Kh_bg_2d(i,j) + cs%Kh_bg_2d(i+1,j+1)) + &
                  (cs%Kh_bg_2d(i+1,j) + cs%Kh_bg_2d(i,j+1))) )
      endif

      ! Use the larger of the above and a function of sin(latitude)
      if (kh_sin_lat>0.) then
        slat_fn = abs( sin( deg2rad * g%geoLatBu(i,j) ) ) ** kh_pwr_of_sine
        cs%Kh_bg_xy(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xy(i,j))
      endif
      if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
        ! Limit the background viscosity to be numerically stable
        cs%Kh_Max_xy(i,j) = kh_limit * grid_sp_q2
        cs%Kh_bg_xy(i,j) = min(cs%Kh_bg_xy(i,j), cs%Kh_Max_xy(i,j))
      endif
    enddo ; enddo
  endif
  if (cs%biharmonic) then
    do j=js-1,jeq+1 ; do i=isq-1,ieq+1
      cs%Idx2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)
      cs%Idxdy2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))
    enddo ; enddo
    do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
      cs%Idx2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)
      cs%Idxdy2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))
    enddo ; enddo
    cs%Ah_bg_xy(:,:) = 0.0
   ! The 0.3 below was 0.4 in MOM1.10.  The change in hq requires
   ! this to be less than 1/3, rather than 1/2 as before.
    if (cs%better_bound_Ah .or. cs%bound_Ah) ah_limit = 0.3 / (dt*64.0)
    if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &
      boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))

    min_grid_sp_h4 = huge(1.)
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j)+cs%dy2h(i,j))
      grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
      cs%grid_sp_h3(i,j) = grid_sp_h3
      if (cs%Smagorinsky_Ah) then
        cs%Biharm_const_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)
        if (cs%bound_Coriolis) then
          fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &
                     abs(g%CoriolisBu(i-1,j)),   abs(g%CoriolisBu(i,j)))
          cs%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
                                     (fmax * boundcorconst)
        endif
      endif
      if (cs%Leith_Ah) then
         cs%biharm6_const_xx(i,j) = leith_bi_const * (grid_sp_h3 * grid_sp_h3)
      endif
      cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
      if (cs%Re_Ah > 0.0) cs%Re_Ah_const_xx(i,j) = grid_sp_h3 / cs%Re_Ah
      if (ah_time_scale > 0.) cs%Ah_bg_xx(i,j) = &
            max(cs%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / ah_time_scale)
      if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
        cs%Ah_Max_xx(i,j) = ah_limit * (grid_sp_h2 * grid_sp_h2)
        cs%Ah_bg_xx(i,j) = min(cs%Ah_bg_xx(i,j), cs%Ah_Max_xx(i,j))
      endif
      min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4)
    enddo ; enddo
    call min_across_pes(min_grid_sp_h4)

    do j=js-1,jeq ; do i=is-1,ieq
      grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j)+cs%dy2q(i,j))
      grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
      if (cs%Smagorinsky_Ah) then
        cs%Biharm_const_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)
        if (cs%bound_Coriolis) then
          cs%Biharm_const2_xy(i,j) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
                                     (abs(g%CoriolisBu(i,j)) * boundcorconst)
        endif
      endif
      if (cs%Leith_Ah) then
         cs%biharm6_const_xy(i,j) = leith_bi_const * (grid_sp_q3 * grid_sp_q3)
      endif
      cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
      if (cs%Re_Ah > 0.0) cs%Re_Ah_const_xy(i,j) = grid_sp_q3 / cs%Re_Ah
      if (ah_time_scale > 0.) cs%Ah_bg_xy(i,j) = &
           max(cs%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / ah_time_scale)
      if (cs%bound_Ah .and. .not.cs%better_bound_Ah) then
        cs%Ah_Max_xy(i,j) = ah_limit * (grid_sp_q2 * grid_sp_q2)
        cs%Ah_bg_xy(i,j) = min(cs%Ah_bg_xy(i,j), cs%Ah_Max_xy(i,j))
      endif
    enddo ; enddo
  endif
  ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.
  if (cs%Laplacian .and. cs%better_bound_Kh) then
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      denom = max( &
         (cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &
          max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &
         (cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &
          max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
      cs%Kh_Max_xx(i,j) = 0.0
      if (denom > 0.0) &
        cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom
    enddo ; enddo
    do j=js-1,jeq ; do i=is-1,ieq
      denom = max( &
         (cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &
          max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &
         (cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &
          max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
      cs%Kh_Max_xy(i,j) = 0.0
      if (denom > 0.0) &
        cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom
    enddo ; enddo
    if (cs%debug) then
      call hchksum(cs%Kh_Max_xx, "Kh_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
      call bchksum(cs%Kh_Max_xy, "Kh_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
    endif
  endif
  ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but
  ! empirically work for CS%bound_coef <~ 1.0
  if (cs%biharmonic .and. cs%better_bound_Ah) then
    do j=js-1,jeq+1 ; do i=isq-1,ieq+1
      u0u(i,j) = (cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*cs%DY_dxT(i+1,j)*(g%IdyCu(i+1,j) + g%IdyCu(i,j))   + &
                                   cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) + &
                 cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) + &
                                   cs%dx2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1)) ) )
      u0v(i,j) = (cs%Idxdy2u(i,j)*(cs%dy2h(i+1,j)*cs%DX_dyT(i+1,j)*(g%IdxCv(i+1,j) + g%IdxCv(i+1,j-1)) + &
                                   cs%dy2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) )   + &
                 cs%Idx2dyCu(i,j)*(cs%dx2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j))   + &
                                   cs%dx2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1)) ) )
    enddo ; enddo
    do j=jsq-1,jeq+1 ; do i=is-1,ieq+1
      v0u(i,j) = (cs%Idxdy2v(i,j)*(cs%dy2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j))       + &
                                   cs%dy2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j)) ) + &
                 cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*cs%DY_dxT(i,j+1)*(g%IdyCu(i,j+1) + g%IdyCu(i-1,j+1))   + &
                                   cs%dx2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) ) )
      v0v(i,j) = (cs%Idxdy2v(i,j)*(cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j))   + &
                                   cs%dy2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j)) ) + &
                 cs%Idx2dyCv(i,j)*(cs%dx2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j))   + &
                                   cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) ) )
    enddo ; enddo
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      denom = max( &
         (cs%dy2h(i,j) * &
          (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0u(i,j) + g%IdyCu(i-1,j)*u0u(i-1,j))  + &
           cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0u(i,j) + g%IdxCv(i,j-1)*v0u(i,j-1))) * &
          max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ),   &
         (cs%dx2h(i,j) * &
          (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0v(i,j) + g%IdyCu(i-1,j)*u0v(i-1,j))  + &
           cs%DX_dyT(i,j)*(g%IdxCv(i,j)*v0v(i,j) + g%IdxCv(i,j-1)*v0v(i,j-1))) * &
          max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )
      cs%Ah_Max_xx(i,j) = 0.0
      if (denom > 0.0) &
        cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom
    enddo ; enddo
    do j=js-1,jeq ; do i=is-1,ieq
      denom = max( &
         (cs%dx2q(i,j) * &
          (cs%DX_dyBu(i,j)*(u0u(i,j+1)*g%IdxCu(i,j+1) + u0u(i,j)*g%IdxCu(i,j))  + &
           cs%DY_dxBu(i,j)*(v0u(i+1,j)*g%IdyCv(i+1,j) + v0u(i,j)*g%IdyCv(i,j))) * &
          max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ),    &
         (cs%dy2q(i,j) * &
          (cs%DX_dyBu(i,j)*(u0v(i,j+1)*g%IdxCu(i,j+1) + u0v(i,j)*g%IdxCu(i,j))  + &
           cs%DY_dxBu(i,j)*(v0v(i+1,j)*g%IdyCv(i+1,j) + v0v(i,j)*g%IdyCv(i,j))) * &
          max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )
      cs%Ah_Max_xy(i,j) = 0.0
      if (denom > 0.0) &
        cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom
    enddo ; enddo
    if (cs%debug) then
      call hchksum(cs%Ah_Max_xx, "Ah_Max_xx", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
      call bchksum(cs%Ah_Max_xy, "Ah_Max_xy", g%HI, haloshift=0, scale=us%L_to_m**4*us%s_to_T)
    endif
  endif
  ! Register fields for output from this module.
  cs%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, time, &
      'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, time, &
      'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

  !CS%id_hf_diffu = register_diag_field('ocean_model', 'hf_diffu', diag%axesCuL, Time, &
  !    'Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', 'm s-2', &
  !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if ((CS%id_hf_diffu > 0) .and. (present(ADp))) then
  !  call safe_alloc_ptr(CS%hf_diffu,G%IsdB,G%IedB,G%jsd,G%jed,G%ke)
  !  call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,G%ke)
  !endif

  !CS%id_hf_diffv = register_diag_field('ocean_model', 'hf_diffv', diag%axesCvL, Time, &
  !    'Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', 'm s-2', &
  !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if ((CS%id_hf_diffv > 0) .and. (present(ADp))) then
  !  call safe_alloc_ptr(CS%hf_diffv,G%isd,G%ied,G%JsdB,G%JedB,G%ke)
  !  call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,G%ke)
  !endif

  cs%id_hf_diffu_2d = register_diag_field('ocean_model', 'hf_diffu_2d', diag%axesCu1, time, &
      'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  if ((cs%id_hf_diffu_2d > 0) .and. (present(adp))) then
    call safe_alloc_ptr(adp%diag_hfrac_u,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
  endif

  cs%id_hf_diffv_2d = register_diag_field('ocean_model', 'hf_diffv_2d', diag%axesCv1, time, &
      'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  if ((cs%id_hf_diffv_2d > 0) .and. (present(adp))) then
    call safe_alloc_ptr(adp%diag_hfrac_v,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
  endif

  if (cs%biharmonic) then
    cs%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, time,    &
        'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T, &
        cmor_field_name='difmxybo',                                             &
        cmor_long_name='Ocean lateral biharmonic viscosity',                     &
        cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')
    cs%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, time, &
        'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)
    cs%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, time, &
        'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim')

    if (cs%id_grid_Re_Ah > 0) &
      ! Compute the smallest biharmonic viscosity capable of modifying the
      ! velocity at floating point precision.
      cs%min_grid_Ah = spacing(1.) * min_grid_sp_h4 * idt
  endif
  if (cs%Laplacian) then
    cs%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, time,   &
        'Laplacian Horizontal Viscosity at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
        cmor_field_name='difmxylo',                                             &
        cmor_long_name='Ocean lateral Laplacian viscosity',                     &
        cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')
    cs%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, time, &
        'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
    cs%id_grid_Re_Kh = register_diag_field('ocean_model', 'grid_Re_Kh', diag%axesTL, time, &
        'Grid Reynolds number for the Laplacian horizontal viscosity at h points', 'nondim')
    cs%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, time, &
      'Vertical vorticity at q Points', 's-1', conversion=us%s_to_T)
    cs%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, time, &
      'Horizontal divergence at h Points', 's-1', conversion=us%s_to_T)
    cs%id_sh_xy_q = register_diag_field('ocean_model', 'sh_xy_q', diag%axesBL, time, &
      'Shearing strain at q Points', 's-1', conversion=us%s_to_T)
    cs%id_sh_xx_h = register_diag_field('ocean_model', 'sh_xx_h', diag%axesTL, time, &
      'Horizontal tension at h Points', 's-1', conversion=us%s_to_T)

    if (cs%id_grid_Re_Kh > 0) &
      ! Compute a smallest Laplacian viscosity capable of modifying the
      ! velocity at floating point precision.
      cs%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * idt
  endif
  if (cs%use_GME) then
      cs%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, time, &
        'GME coefficient at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
      cs%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, time, &
        'GME coefficient at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
      cs%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,time,&
      'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', &
      'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
  endif
  cs%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,time,&
      'Integral work done by lateral friction terms', &
      'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
  cs%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,time,      &
      'Depth integrated work done by lateral friction', &
      'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2, &
      cmor_field_name='dispkexyfo',                                                              &
      cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',&
      cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
  if (cs%Laplacian .or. get_all) then
  endif
end subroutine hor_visc_init
!> Calculates factors in the anisotropic orientation tensor to be align with the grid.
!! With n1=1 and n2=0, this recovers the approach of Large et al, 2001.
subroutine align_aniso_tensor_to_grid(CS, n1, n2)
  type(hor_visc_cs), pointer :: CS !< Control structure for horizontal viscosity
  real,              intent(in) :: n1 !< i-component of direction vector [nondim]
  real,              intent(in) :: n2 !< j-component of direction vector [nondim]
  ! Local variables
  real :: recip_n2_norm
  ! For normalizing n=(n1,n2) in case arguments are not a unit vector
  recip_n2_norm = n1**2 + n2**2
  if (recip_n2_norm > 0.) recip_n2_norm = 1./recip_n2_norm
  cs%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
  cs%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
  cs%n1n1_m_n2n2_h(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
  cs%n1n1_m_n2n2_q(:,:) = ( n1 * n1 - n2 * n2 ) * recip_n2_norm
end subroutine align_aniso_tensor_to_grid
!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any
!! horizontal two-grid-point noise
subroutine smooth_gme(CS,G,GME_flux_h,GME_flux_q)
  ! Arguments
  type(hor_visc_cs),                            pointer       :: CS        !< Control structure
  type(ocean_grid_type),                        intent(in)    :: G         !< Ocean grid
  real, dimension(SZI_(G),SZJ_(G)),   optional, intent(inout) :: GME_flux_h!< GME diffusive flux
                                                              !! at h points
  real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< GME diffusive flux
                                                              !! at q points
  ! local variables
  real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original
  real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original
  real :: wc, ww, we, wn, ws ! averaging weights for smoothing
  integer :: i, j, k, s
  do s=1,1
    ! Update halos
    if (present(gme_flux_h)) then
      call pass_var(gme_flux_h, g%Domain)
      gme_flux_h_original = gme_flux_h
      ! apply smoothing on GME
      do j = g%jsc, g%jec
        do i = g%isc, g%iec
          ! skip land points
          if (g%mask2dT(i,j)==0.) cycle
          ! compute weights
          ww = 0.125 * g%mask2dT(i-1,j)
          we = 0.125 * g%mask2dT(i+1,j)
          ws = 0.125 * g%mask2dT(i,j-1)
          wn = 0.125 * g%mask2dT(i,j+1)
          wc = 1.0 - (ww+we+wn+ws)
          gme_flux_h(i,j) =  wc * gme_flux_h_original(i,j)   &
                           + ww * gme_flux_h_original(i-1,j) &
                           + we * gme_flux_h_original(i+1,j) &
                           + ws * gme_flux_h_original(i,j-1) &
                           + wn * gme_flux_h_original(i,j+1)
      enddo; enddo
    endif
    ! Update halos
    if (present(gme_flux_q)) then
      call pass_var(gme_flux_q, g%Domain, position=corner, complete=.true.)
      gme_flux_q_original = gme_flux_q
      ! apply smoothing on GME
      do j = g%JscB, g%JecB
        do i = g%IscB, g%IecB
          ! skip land points
          if (g%mask2dBu(i,j)==0.) cycle
          ! compute weights
          ww = 0.125 * g%mask2dBu(i-1,j)
          we = 0.125 * g%mask2dBu(i+1,j)
          ws = 0.125 * g%mask2dBu(i,j-1)
          wn = 0.125 * g%mask2dBu(i,j+1)
          wc = 1.0 - (ww+we+wn+ws)
          gme_flux_q(i,j) =  wc * gme_flux_q_original(i,j)   &
                           + ww * gme_flux_q_original(i-1,j) &
                           + we * gme_flux_q_original(i+1,j) &
                           + ws * gme_flux_q_original(i,j-1) &
                           + wn * gme_flux_q_original(i,j+1)
      enddo; enddo
    endif
  enddo ! s-loop
end subroutine smooth_gme
!> Deallocates any variables allocated in hor_visc_init.
subroutine hor_visc_end(CS)
  type(hor_visc_cs), pointer :: CS !< The control structure returned by a
                                   !! previous call to hor_visc_init.
  if (cs%Laplacian .or. cs%biharmonic) then
    dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)
    dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)
    dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)
  endif
  if (cs%Laplacian) then
    dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)
    dealloc_(cs%grid_sp_h2)
    if (cs%bound_Kh) then
      dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
    endif
    if (cs%Smagorinsky_Kh) then
      dealloc_(cs%Laplac2_const_xx) ; dealloc_(cs%Laplac2_const_xy)
    endif
    if (cs%Leith_Kh) then
      dealloc_(cs%Laplac3_const_xx) ; dealloc_(cs%Laplac3_const_xy)
    endif
  endif
  if (cs%biharmonic) then
    dealloc_(cs%grid_sp_h3)
    dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)
    dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)
    dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)
    if (cs%bound_Ah) then
      dealloc_(cs%Ah_Max_xx) ; dealloc_(cs%Ah_Max_xy)
    endif
    if (cs%Smagorinsky_Ah) then
      dealloc_(cs%Biharm6_const_xx) ; dealloc_(cs%Biharm6_const_xy)
    endif
    if (cs%Leith_Ah) then
      dealloc_(cs%Biharm_const_xx) ; dealloc_(cs%Biharm_const_xy)
    endif
    if (cs%Re_Ah > 0.0) then
      dealloc_(cs%Re_Ah_const_xx) ; dealloc_(cs%Re_Ah_const_xy)
    endif
  endif
  if (cs%anisotropic) then
    dealloc_(cs%n1n2_h)
    dealloc_(cs%n1n2_q)
    dealloc_(cs%n1n1_m_n2n2_h)
    dealloc_(cs%n1n1_m_n2n2_q)
  endif
  deallocate(cs)
end subroutine hor_visc_end
!> \namespace mom_hor_visc
!!
!! This module contains the subroutine horizontal_viscosity() that calculates the
!! effects of horizontal viscosity, including parameterizations of the value of
!! the viscosity itself. horizontal_viscosity() calculates the acceleration due to
!! some combination of a biharmonic viscosity and a Laplacian viscosity. Either or
!! both may use a coefficient that depends on the shear and strain of the flow.
!! All metric terms are retained. The Laplacian is calculated as the divergence of
!! a stress tensor, using the form suggested by Smagorinsky (1993). The biharmonic
!! is calculated by twice applying the divergence of the stress tensor that is
!! used to calculate the Laplacian, but without the dependence on thickness in the
!! first pass. This form permits a variable viscosity, and indicates no
!! acceleration for either resting fluid or solid body rotation.
!!
!! The form of the viscous accelerations is discussed extensively in Griffies and
!! Hallberg (2000), and the implementation here follows that discussion closely.
!! We use the notation of Smith and McWilliams (2003) with the exception that the
!! isotropic viscosity is \f$\kappa_h\f$.
!!
!! \section section_horizontal_viscosity Horizontal viscosity in MOM
!!
!! In general, the horizontal stress tensor can be written as
!! \f[
!! {\bf \sigma} =
!! \begin{pmatrix}
!! \frac{1}{2} \left( \sigma_D + \sigma_T \right) & \frac{1}{2} \sigma_S \\\\
!! \frac{1}{2} \sigma_S & \frac{1}{2} \left( \sigma_D - \sigma_T \right)
!! \end{pmatrix}
!! \f]
!! where \f$\sigma_D\f$, \f$\sigma_T\f$ and \f$\sigma_S\f$ are stresses associated with
!! invariant factors in the strain-rate tensor. For a Newtonian fluid, the stress
!! tensor is usually linearly related to the strain-rate tensor. The horizontal
!! strain-rate tensor is
!! \f[
!! \dot{\bf e} =
!! \begin{pmatrix}
!! \frac{1}{2} \left( \dot{e}_D + \dot{e}_T \right) & \frac{1}{2} \dot{e}_S \\\\
!! \frac{1}{2} \dot{e}_S & \frac{1}{2} \left( \dot{e}_D - \dot{e}_T \right)
!! \end{pmatrix}
!! \f]
!! where \f$\dot{e}_D = \partial_x u + \partial_y v\f$ is the horizontal divergence,
!! \f$\dot{e}_T = \partial_x u - \partial_y v\f$ is the horizontal tension, and
!! \f$\dot{e}_S = \partial_y u + \partial_x v\f$ is the horizontal shear strain.
!!
!! The trace of the stress tensor, \f$tr(\bf \sigma) = \sigma_D\f$, is usually
!! absorbed into the pressure and only the deviatoric stress tensor considered.
!! From here on, we drop \f$\sigma_D\f$. The trace of the strain tensor, \f$tr(\bf e) =
!! \dot{e}_D\f$ is non-zero for horizontally divergent flow but only enters the
!! stress tensor through \f$\sigma_D\f$ and so we will drop \f$\sigma_D\f$ from
!! calculations of the strain tensor in the code. Therefore the horizontal stress
!! tensor can be considered to be
!! \f[
!! {\bf \sigma} =
!! \begin{pmatrix}
!! \frac{1}{2} \sigma_T & \frac{1}{2} \sigma_S \\\\
!! \frac{1}{2} \sigma_S & - \frac{1}{2} \sigma_T
!! \end{pmatrix}
!! .\f]
!!
!! The stresses above are linearly related to the strain through a viscosity
!! coefficient, \f$\kappa_h\f$:
!! \f{eqnarray*}{
!! \sigma_T & = & 2 \kappa_h \dot{e}_T \\\\
!! \sigma_S & = & 2 \kappa_h \dot{e}_S
!! .
!! \f}
!!
!! The viscosity \f$\kappa_h\f$ may either be a constant or variable. For example,
!! \f$\kappa_h\f$ may vary with the shear, as proposed by Smagorinsky (1993).
!!
!! The accelerations resulting form the divergence of the stress tensor are
!! \f{eqnarray*}{
!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right)
!! & = &
!! \partial_x \left( \frac{1}{2} \sigma_T \right)
!! + \partial_y \left( \frac{1}{2} \sigma_S \right)
!! \\\\
!! & = &
!! \partial_x \left( \kappa_h \dot{e}_T \right)
!! + \partial_y \left( \kappa_h \dot{e}_S \right)
!! \\\\
!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right)
!! & = &
!! \partial_x \left( \frac{1}{2} \sigma_S \right)
!! + \partial_y \left( \frac{1}{2} \sigma_T \right)
!! \\\\
!! & = &
!! \partial_x \left( \kappa_h \dot{e}_S \right)
!! + \partial_y \left( - \kappa_h \dot{e}_T \right)
!! .
!! \f}
!!
!! The form of the Laplacian viscosity in general coordinates is:
!! \f{eqnarray*}{
!! \hat{\bf x} \cdot \left( \nabla \cdot \sigma \right)
!! & = &
!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_T \right)
!! + \partial_y \left( \kappa_h h \dot{e}_S \right) \right]
!! \\\\
!! \hat{\bf y} \cdot \left( \nabla \cdot \sigma \right)
!! & = &
!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_S \right)
!! - \partial_y \left( \kappa_h h \dot{e}_T \right) \right]
!! .
!! \f}
!!
!! \subsection section_laplacian_viscosity_coefficient Laplacian viscosity coefficient
!!
!! The horizontal viscosity coefficient, \f$\kappa_h\f$, can have multiple components.
!! The isotropic components are:
!!   - A uniform background component, \f$\kappa_{bg}\f$.
!!   - A constant but spatially variable 2D map, \f$\kappa_{2d}(x,y)\f$.
!!   - A ''MICOM'' viscosity, \f$U_\nu \Delta(x,y)\f$, which uses a constant
!! velocity scale, \f$U_\nu\f$ and a measure of the grid-spacing \f$\Delta(x,y)^2 =
!! \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}\f$.
!!   - A function of
!! latitude, \f$\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n\f$.
!!   - A dynamic Smagorinsky viscosity, \f$\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}\f$.
!!   - A dynamic Leith viscosity, \f$\kappa_{Lth}(x,y,t) =
!!                                    C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}\f$.
!!
!! A maximum stable viscosity, \f$\kappa_{max}(x,y)\f$ is calculated based on the
!! grid-spacing and time-step and used to clip calculated viscosities.
!!
!! The static components of \f$\kappa_h\f$ are first combined as follows:
!! \f[
!! \kappa_{static} = \min \left[ \max\left(
!! \kappa_{bg},
!! U_\nu \Delta(x,y),
!! \kappa_{2d}(x,y),
!! \kappa_\phi(x,y)
!! \right)
!! , \kappa_{max}(x,y) \right]
!! \f]
!! and stored in the module control structure as variables <code>Kh_bg_xx</code> and
!! <code>Kh_bg_xy</code> for the tension (h-points) and shear (q-points) components
!! respectively.
!!
!! The full viscosity includes the dynamic components as follows:
!! \f[
!! \kappa_h(x,y,t) = r(\Delta,L_d)
!! \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right)
!! \f]
!! where \f$r(\Delta,L_d)\f$ is a resolution function.
!!
!! The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each
!! other.
!!
!! \subsection section_viscous_boundary_conditions Viscous boundary conditions
!!
!! Free slip boundary conditions have been coded, although no slip boundary
!! conditions can be used with the Laplacian viscosity based on the 2D land-sea
!! mask. For a western boundary, for example, the boundary conditions with the
!! biharmonic operator would be written as:
!! \f[
!!   \partial_x v = 0 ; \partial_x^3 v = 0 ; u = 0 ; \partial_x^2 u = 0 ,
!! \f]
!! while for a Laplacian operator, they are simply
!! \f[
!!   \partial_x v = 0 ; u = 0 .
!! \f]
!! These boundary conditions are largely dictated by the use of an Arakawa
!! C-grid and by the varying layer thickness.
!!
!! \subsection section_anisotropic_viscosity Anisotropic viscosity
!!
!! Large et al., 2001, proposed enhancing viscosity in a particular direction and the
!! approach was generalized in Smith and McWilliams, 2003. We use the second form of their
!! two coefficient anisotropic viscosity (section 4.3). We also replace their
!! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and
!! \f$\kappa_a = D\f$ so that \f$\kappa_h\f$ can be considered the isotropic
!! viscosity and \f$\kappa_a=D\f$ can be consider the anisotropic viscosity. The
!! direction of anisotropy is defined by a unit vector \f$\hat{\bf
!! n}=(n_1,n_2)\f$.
!!
!! The contributions to the stress tensor are
!! \f[
!! \begin{pmatrix}
!! \sigma_T \\\\ \sigma_S
!! \end{pmatrix}
!! =
!! \left[
!! \begin{pmatrix}
!! 2 \kappa_h + \kappa_a & 0 \\\\
!! 0 & 2 \kappa_h
!! \end{pmatrix}
!! + 2 \kappa_a n_1 n_2
!! \begin{pmatrix}
!! - 2 n_1 n_2 & n_1^2 - n_2^2 \\\\
!! n_1^2 - n_2^2 & 2 n_1 n_2
!! \end{pmatrix}
!! \right]
!! \begin{pmatrix}
!! \dot{e}_T \\\\ \dot{e}_S
!! \end{pmatrix}
!! \f]
!! Dissipation of kinetic energy requires \f$\kappa_h \geq 0\f$ and \f$2 \kappa_h + \kappa_a \geq 0\f$.
!! Note that when anisotropy is aligned with the x-direction, \f$n_1 = \pm 1\f$, then
!! \f$n_2 = 0\f$ and the cross terms vanish. The accelerations in this aligned limit
!! with constant coefficients become
!! \f{eqnarray*}{
!! \hat{\bf x} \cdot \nabla \cdot {\bf \sigma}
!! & = &
!! \partial_x \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)
!! + \partial_y \left( \kappa_h \dot{e}_S \right)
!! \\\\
!! & = &
!! \left( \kappa_h + \kappa_a \right) \partial_{xx} u
!! + \kappa_h \partial_{yy} u
!! - \frac{1}{2} \kappa_a \partial_x \left( \partial_x u + \partial_y v \right)
!! \\\\
!! \hat{\bf y} \cdot \nabla \cdot {\bf \sigma}
!! & = &
!! \partial_x \left( \kappa_h \dot{e}_S \right)
!! - \partial_y \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)
!! \\\\
!! & = &
!! \kappa_h \partial_{xx} v
!! + \left( \kappa_h + \kappa_a \right) \partial_{yy} v
!! - \frac{1}{2} \kappa_a \partial_y \left( \partial_x u + \partial_y v \right)
!! \f}
!! which has contributions akin to a negative divergence damping (a divergence
!! enhancement?) but which is weaker than the enhanced tension terms by half.
!!
!! \subsection section_viscous_discretization Discretization
!!
!! The horizontal tension, \f$\dot{e}_T\f$, is stored in variable <code>sh_xx</code> and
!! discretized as
!! \f[
!! \dot{e}_T
!! = \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} u \right)
!! - \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} v \right)
!! .
!! \f]
!! The horizontal divergent strain, \f$\dot{e}_D\f$, is stored in variable
!! <code>div_xx</code> and discretized as
!! \f[
!! \dot{e}_D
!! = \frac{1}{h A} \left( \delta_i \left( \overline{h}^i \Delta y \, u \right)
!! + \delta_j \left( \overline{h}^j\Delta x \, v \right) \right)
!! .
!! \f]
!! Note that for expediency this is the exact discretization used in the
!! continuity equation.
!!
!! The horizontal shear strain, \f$\dot{e}_S\f$, is stored in variable <code>sh_xy</code>
!! and discretized as
!! \f[
!! \dot{e}_S = v_x + u_y
!! \f]
!! where
!! \f{align*}{
!! v_x &= \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} v \right) \\\\
!! u_y &= \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} u \right)
!! \f}
!! which are calculated separately so that no-slip or free-slip boundary
!! conditions can be applied to \f$v_x\f$ and \f$u_y\f$ where appropriate.
!!
!! The tendency for the x-component of the divergence of stress is stored in
!! variable <code>diffu</code> and discretized as
!! \f[
!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) =
!! \frac{1}{A \overline{h}^i} \left(
!! \frac{1}{\Delta y} \delta_i \left( h \Delta y^2 \kappa_h \dot{e}_T \right) +
!! \frac{1}{\Delta x} \delta_j \left( \tilde{h}^{ij} \Delta x^2 \kappa_h \dot{e}_S \right)
!! \right)
!! .
!! \f]
!!
!! The tendency for the y-component of the divergence of stress is stored in
!! variable <code>diffv</code> and discretized as
!! \f[
!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) =
!! \frac{1}{A \overline{h}^j} \left(
!! \frac{1}{\Delta y} \delta_i \left( \tilde{h}^{ij} \Delta y^2 A_M \dot{e}_S \right)
!! - \frac{1}{\Delta x} \delta_j \left( h \Delta x^2 A_M \dot{e}_T \right)
!! \right)
!! .
!! \f]
!!
!! \subsection section_viscous_refs References
!!
!! Griffies, S.M., and Hallberg, R.W., 2000: Biharmonic friction with a
!! Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models.
!! Monthly Weather Review, 128(8), 2935-2946.
!! https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2
!!
!! Large, W.G., Danabasoglu, G., McWilliams, J.C., Gent, P.R. and Bryan, F.O.,
!! 2001: Equatorial circulation of a global ocean climate model with
!! anisotropic horizontal viscosity.
!! Journal of Physical Oceanography, 31(2), pp.518-536.
!! https://doi.org/10.1175/1520-0485(2001)031%3C0518:ECOAGO%3E2.0.CO;2
!!
!! Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear
!! viscosities. Large eddy simulation of complex engineering and geophysical
!! flows, 1, 69-106.
!!
!! Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for
!! ocean models. Ocean Modelling, 5(2), 129-156.
!! https://doi.org/10.1016/S1463-5003(02)00016-1
end module mom_hor_visc