MOM_hor_visc.F90

Go to the documentation of this file.

module mom_hor_visc
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - June 2002.                        *
!*                                                                     *
!*    This program contains the subroutine that calculates the         *
!*  effects of horizontal viscosity, including parameterizations of    *
!*  the value of the viscosity itself. horizontal_viscosity calc-      *
!*  ulates 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.                                                          *
!*                                                                     *
!*    set_up_hor_visc calculates and stores the values of a number of  *
!*  metric functions that are used in horizontal_viscosity.  It is     *
!*  called by horizontal_viscosity the first time that the latter is   *
!*  called.                                                            *
!*                                                                     *
!*    The form of the Laplacian viscosity is:                          *
!*                                                                     *
!*    diffu = 1/h * {d/dx[KH*h*sh_xx] + d/dy[KH*h*sh_xy]}              *
!*    diffv = 1/h * {d/dx[KH*h*sh_xy] - d/dy[KH*h*sh_xx]}              *
!*                                                                     *
!*    sh_xx = du/dx - dv/dy         sh_xy = du/dy + dv/dx              *
!*                                                                     *
!*  with appropriate metric terms thrown in.  KH may either be a       *
!*  constant or may vary with the shear, as proposed by Smagorinsky.   *
!*  The form of this term is discussed extensively in Griffies and     *
!*  Hallberg (MWR, 2000), and the implementation here follows that     *
!*  discussion closely.                                                *
!*                                                                     *
!*    Only free slip boundary conditions have been coded, although     *
!*  no slip boundary conditions could be used with the Laplacian       *
!*  viscosity.  For a western boundary, for example, the boundary      *
!*  conditions with the biharmonic operator would be written as:       *
!*    dv/dx = 0, d^3v/dx^3 = 0, u = 0, d^2u/dx^2 = 0 ,                 *
!*  while for a Laplacian operator, they are simply:                   *
!*    dv/dx = 0, u = 0 .                                               *
!*  These boundary conditions are largely dictated by the use of       *
!*  a an Arakawa C-grid and by the varying layer thickness.            *
!*                                                                     *
!*                                                                     *
!*                                                                     *
!*                                                                     *
!*                                                                     *
!* Macros written all in capital letters are defined in MOM_memory.h.  *
!*                                                                     *
!*     A small fragment of the C-grid is shown below:                  *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu, hq, str_xy, sh_xy         *
!*    j+1  > o > o >   At ^:  v, diffv, v0                             *
!*    j    x ^ x ^ x   At >:  u, diffu, u0                             *
!*    j    > o > o >   At o:  h, str_xx, sh_xx                         *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

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
use mom_error_handler,         only : mom_error, fatal, warning
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
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_verticalgrid,          only : verticalgrid_type
use mom_io,                    only : read_data, slasher

implicit none ; private

#include <MOM_memory.h>

public horizontal_viscosity, hor_visc_init, hor_visc_end

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 :: 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    :: bound_coef      ! The nondimensional coefficient of the ratio of
                             ! the viscosity bounds to the theoretical maximum
                             ! for stability without considering other terms.
                             ! 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 :: leith_ah        ! If true, use a biharmonic form of 2D Leith
                             ! nonlinear eddy viscosity. AH is the background.
  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. The default is 0.0

  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    kh_bg_xx,        &! The background Laplacian viscosity at h points, in units
                      ! of m2 s-1. The actual viscosity may be the larger of this
                      ! viscosity and the Smagorinsky and Leith viscosities.
    kh_bg_2d,        &! The background Laplacian viscosity at h points, in units
                      ! of m2 s-1. The actual viscosity may be the larger of this
                      ! viscosity and the Smagorinsky and Leith viscosities.
    ah_bg_xx,        &! The background biharmonic viscosity at h points, in units
                      ! of m4 s-1. The actual viscosity may be the larger of this
                      ! viscosity and the Smagorinsky and Leith viscosities.
    kh_max_xx,       &! The maximum permitted Laplacian viscosity, m2 s-1.
    ah_max_xx,       &! The maximum permitted biharmonic viscosity, m4 s-1.
    biharm_const2_xx,&! A constant relating the biharmonic viscosity to the
                      ! square of the velocity shear, in m4 s.  This value is
                      ! set to be the magnitude of the Coriolis terms once the
                      ! velocity differences reach a value of order 1/2 MAXVEL.

    reduction_xx      ! The amount by which stresses through h points are reduced
                      ! due to partial barriers. Nondimensional.

  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    kh_bg_xy,        &! The background Laplacian viscosity at q points, in units
                      ! of m2 s-1. The actual viscosity may be the larger of this
                      ! viscosity and the Smagorinsky and Leith viscosities.
    ah_bg_xy,        &! The background biharmonic viscosity at q points, in units
                      ! of m4 s-1. The actual viscosity may be the larger of this
                      ! viscosity and the Smagorinsky and Leith viscosities.
    kh_max_xy,       &! The maximum permitted Laplacian viscosity, m2 s-1.
    ah_max_xy,       &! The maximum permitted biharmonic viscosity, m4 s-1.
    biharm_const2_xy,&! A constant relating the biharmonic viscosity to the
                      ! square of the velocity shear, in m4 s.  This value is
                      ! set to be the magnitude of the Coriolis terms once the
                      ! velocity differences reach a value of order 1/2 MAXVEL.
    reduction_xy      ! The amount by which stresses through q points are reduced
                      ! due to partial barriers. Nondimensional.

! The following variables are precalculated combinations of metric terms.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    dx2h, dy2h, &       ! dx^2  and dy^2  at h points, in m2
    dx_dyt, dy_dxt      ! dx/dy and dy/dx at h points, nondim
  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    dx2q, dy2q, &       ! dx^2  and dy^2  at q points, in m2
    dx_dybu, dy_dxbu    ! dx/dy and dy/dx at q points, nondim
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    idx2dycu, idxdy2u   ! 1/(dx^2 dy) and 1/(dx dy^2) at u points, in m-3
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    idx2dycv, idxdy2v   ! 1/(dx^2 dy) and 1/(dx dy^2) at v points, in m-3

! The following variables are precalculated time-invariant combinations of
! parameters and metric terms.
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    laplac_const_xx, &  ! Laplacian  metric-dependent constants (nondim)
    biharm_const_xx, &  ! Biharmonic metric-dependent constants (nondim)
    laplac3_const_xx, & ! Laplacian  metric-dependent constants (nondim)
    biharm5_const_xx    ! Biharmonic metric-dependent constants (nondim)

  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
    laplac_const_xy, &  ! Laplacian  metric-dependent constants (nondim)
    biharm_const_xy, &  ! Biharmonic metric-dependent constants (nondim)
    laplac3_const_xy, & ! Laplacian  metric-dependent constants (nondim)
    biharm5_const_xy    ! Biharmonic metric-dependent constants (nondim)

  type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic timing

  ! diagnostic ids
  integer :: id_diffu     = -1, id_diffv         = -1
  integer :: id_ah_h      = -1, id_ah_q          = -1
  integer :: id_kh_h      = -1, id_kh_q          = -1
  integer :: id_frictwork = -1, id_frictworkintz = -1

end type hor_visc_cs

contains

!>    This subroutine determines 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 to ie range before this subroutine is called:
!!  v[is-2,is-1,ie+1,ie+2], u[is-2,is-1,ie+1,ie+2], and h[is-1,ie+1],
!! with a similarly sized halo in the y-direction.
subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, OBC)
  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, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                 intent(in)  :: v      !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                 intent(in)  :: h      !< Layer thicknesses, in H
                                                       !! (usually 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 (m/s2)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                 intent(out) :: diffv  !< Meridional acceleration due to convergence
                                                       !! of along-coordinate stress tensor (m/s2).
  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(hor_visc_cs),             pointer     :: CS     !< Pontrol structure returned by a previous
                                                       !! call to hor_visc_init.
  type(ocean_obc_type), pointer, optional    :: OBC    !< Pointer to an open boundary condition type

! Arguments:
!  (in)      u      - zonal velocity (m/s)
!  (in)      v      - meridional velocity (m/s)
!  (in)      h      - layer thickness (m or kg m-2); h units are referred to as H.
!  (out)     diffu  - zonal acceleration due to convergence of
!                     along-coordinate stress tensor (m/s2)
!  (out)     diffv  - meridional acceleration due to convergence of
!                     along-coordinate stress tensor (m/s2)
!  (inout)   MEKE   - pointer to a structure containing fields related to
!                     Mesoscale Eddy Kinetic Energy
!  (in)      VarMix - pointer to a structure with fields that specify the
!                     spatially variable viscosities
!  (in)      G      - ocean grid structure
!  (in)      GV     - The ocean's vertical grid structure.
!  (in)      CS     - control structure returned by a previous call to
!                     hor_visc_init
!  (in)      OBC    - pointer to an open boundary condition type

!  By R. Hallberg, August 1998 - November 1998.
!    This subroutine determines 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 to ie range before this subroutine is called:
!   v[is-2,is-1,ie+1,ie+2], u[is-2,is-1,ie+1,ie+2], and h[is-1,ie+1],
!  with a similarly sized halo in the y-direction.

  real, dimension(SZIB_(G),SZJ_(G)) :: &
    u0, &   ! Laplacian of u (m-1 s-1)
    h_u     ! Thickness interpolated to u points, in H.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    v0, &   ! Laplacian of v (m-1 s-1)
    h_v     ! Thickness interpolated to v points, in H.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    sh_xx, &      ! horizontal tension (du/dx - dv/dy) (1/sec) including metric terms
    str_xx,&      ! str_xx is the diagonal term in the stress tensor (H m2 s-2)
    bhstr_xx,&    ! A copy of str_xx that only contains the biharmonic contribution (H m2 s-2)
    div_xx, &     ! horizontal divergence (du/dx + dv/dy) (1/sec) including metric terms
    FrictWorkIntz ! depth integrated energy dissipated by lateral friction (W/m2)

  real, dimension(SZIB_(G),SZJB_(G)) :: &
    dvdx, dudy, & ! components in the shearing strain (s-1)
    sh_xy,  &     ! horizontal shearing strain (du/dy + dv/dx) (1/sec) including metric terms
    str_xy, &     ! str_xy is the cross term in the stress tensor (H m2 s-2)
    bhstr_xy, &   ! A copy of str_xy that only contains the biharmonic contribution (H m2 s-2)
    vort_xy       ! vertical vorticity (dv/dx - du/dy) (1/sec) including metric terms

  real, dimension(SZI_(G),SZJB_(G)) :: &
    vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) (m-1 sec-1) including metric terms
    div_xx_dy     ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) (m-1 sec-1) including metric terms

  real, dimension(SZIB_(G),SZJ_(G)) :: &
    vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) (m-1 sec-1) including metric terms
    div_xx_dx     ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) (m-1 sec-1) including metric terms

  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
    Ah_q, &   ! biharmonic viscosity at corner points (m4/s)
    Kh_q      ! Laplacian viscosity at corner points (m2/s)

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    Ah_h, &          ! biharmonic viscosity at thickness points (m4/s)
    Kh_h, &          ! Laplacian viscosity at thickness points (m2/s)
    FrictWork        ! energy dissipated by lateral friction (W/m2)

  real :: Ah         ! biharmonic viscosity (m4/s)
  real :: Kh         ! Laplacian  viscosity (m2/s)
  real :: AhSm       ! Smagorinsky biharmonic viscosity (m4/s)
  real :: KhSm       ! Smagorinsky Laplacian viscosity  (m2/s)
  real :: AhLth      ! 2D Leith biharmonic viscosity (m4/s)
  real :: KhLth      ! 2D Leith Laplacian viscosity  (m2/s)
  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 (1/s)
  real :: Vort_mag   ! magnitude of the vorticity (1/s)
  real :: h2uq, h2vq ! temporary variables in units of H^2 (i.e. 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, in units of H (i.e. m or kg m-2).
  real :: hq         ! harmonic mean of the harmonic means of the u- & v-
                     ! point thicknesses, in H; This form guarantees that hq/hu < 4.
  real :: h_neglect  ! thickness so small it can be lost in roundoff and so neglected (H)
  real :: h_neglect3 ! h_neglect^3, in H3
  real :: hrat_min   ! minimum thicknesses at the 4 neighboring
                     ! velocity points divided by the thickness at the stress
                     ! point (h or q point) (nondimensional)
  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
  real :: RoScl     ! The scaling function for MEKE source term
  real :: FatH      ! abs(f) at h-point for MEKE source term (s-1)

  logical :: rescale_Kh
  logical :: find_FrictWork
  logical :: apply_OBC = .false.
  logical :: use_MEKE_Ku
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: i, j, k, n

  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

  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.
  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

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

!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h,  &
!$OMP                                  rescale_Kh,VarMix,h_neglect,h_neglect3,        &
!$OMP                                  Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
!$OMP                                  find_FrictWork,FrictWork,use_MEKE_Ku,MEKE)     &
!$OMP                          private(u0, v0, sh_xx, str_xx, visc_bound_rem,         &
!$OMP                                  sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
!$OMP                                  bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
!$OMP                                  vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
!$OMP                                  div_xx, div_xx_dx, div_xx_dy, mod_Leith,       &
!$OMP                                  Shear_mag, h2uq, h2vq, hq, Kh_scale, hrat_min)
  do k=1,nz

!    This code uses boundary conditions that are consistent with
!  free slip and no normal flow boundary conditions.  The boundary
!  conditions for the western boundary, for example, are:
!    dv/dx = 0,  d^3v/dx^3 = 0,    u = 0,     d^2u/dx^2 = 0 .
!  The overall scheme is second order accurate.
!    All of the metric terms are retained, and the repeated use of
!  the symmetric stress tensor insures that no stress is applied with
!  no flow or solid-body rotation, even with non-constant values of
!  of the biharmonic viscosity.

!  The following are the forms of the horizontal tension and hori-
!  shearing strain advocated by Smagorinsky (1993) and discussed in
!  Griffies and Hallberg (MWR, 2000).
    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
      sh_xx(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)) - &
                    cs%DX_dyT(i,j)*(g%IdxCv(i,j) * v(i,j,k) - &
                                    g%IdxCv(i,j-1)*v(i,j-1,k)))
      div_xx(i,j) = 0.5*((g%dyCu(i,j) * u(i,j,k) * (h(i+1,j,k)+h(i,j,k)) - &
                          g%dyCu(i-1,j) * u(i-1,j,k) * (h(i-1,j,k)+h(i,j,k)) ) + &
                         (g%dxCv(i,j) * v(i,j,k) * (h(i,j,k)+h(i,j+1,k)) - &
                          g%dxCv(i,j-1)*v(i,j-1,k)*(h(i,j,k)+h(i,j-1,k))))*g%IareaT(i,j)/ &
                                    (h(i,j,k) + h_neglect)
    enddo ; enddo
    ! Components for the shearing strain
    do j=js-2,jeq+1 ; do i=is-2,ieq+1
      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 accomodate 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
    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

    ! 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) 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.
            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.
            endif
          enddo
        endif
      endif
      if (obc%segment(n)%direction == obc_direction_n) then
        ! There are extra wide halos here to accomodate 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

    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

    if (cs%no_slip) then
      do j=js-2,jeq+1 ; do i=is-2,ieq+1
        vort_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
        vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )
      enddo ; enddo
    endif

! Vorticity gradient
    do j=js-2,jeq+1 ; do i=is-1,ieq+1
      vort_xy_dx(i,j) = cs%DY_dxBu(i,j)*(vort_xy(i,j)*g%IdyCu(i,j) - vort_xy(i-1,j)*g%IdyCu(i-1,j))
    enddo ; enddo

    do j=js-1,jeq+1 ; do i=is-2,ieq+1
      vort_xy_dy(i,j) = cs%DX_dyBu(i,j)*(vort_xy(i,j)*g%IdxCv(i,j) - vort_xy(i,j-1)*g%IdxCv(i,j-1))
    enddo ; enddo

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

    do j=js-2,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

! Coefficient for modified Leith
    if (cs%Modified_Leith) then
      mod_leith = 1.0
    else
      mod_leith = 0.0
    endif

!  Evaluate u0 = x.Div(Grad u) and v0 = y.Div( Grad u)
    if (cs%biharmonic) then
      do j=js-1,jeq+1 ; do i=isq-1,ieq+1
        u0(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
        v0(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 .and. 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
              v0(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
              u0(i,j) = 0.
            enddo
          endif
        enddo
      endif
    endif

    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
        vort_mag = sqrt( &
          0.5*((vort_xy_dx(i,j-1)*vort_xy_dx(i,j-1) + vort_xy_dx(i,j)*vort_xy_dx(i,j)) + &
                (vort_xy_dy(i-1,j)*vort_xy_dy(i-1,j) + vort_xy_dy(i,j)*vort_xy_dy(i,j))) + &
          mod_leith*0.5*((div_xx_dx(i,j)*div_xx_dx(i,j) + div_xx_dx(i-1,j)*div_xx_dx(i-1,j)) + &
                (div_xx_dy(i,j)*div_xx_dy(i,j) + div_xx_dy(i,j-1)*div_xx_dy(i,j-1))))
      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_scale = 1.0 ; if (rescale_kh) kh_scale = varmix%Res_fn_h(i,j)
        khsm = 0.0; khlth = 0.0
        if ((cs%Smagorinsky_Kh) .or. (cs%Leith_Kh)) then
          if (cs%Smagorinsky_Kh) &
            khsm = cs%LAPLAC_CONST_xx(i,j) * shear_mag
          if (cs%Leith_Kh) &
            khlth = cs%LAPLAC3_CONST_xx(i,j) * vort_mag
          kh = kh_scale * max(khlth, max(cs%Kh_bg_xx(i,j), khsm))
          if (cs%bound_Kh .and. .not.cs%better_bound_Kh) &
            kh = min(kh, cs%Kh_Max_xx(i,j))
        else
          kh = kh_scale * cs%Kh_bg_xx(i,j)
        endif

        if (use_meke_ku) then
          kh = kh + meke%Ku(i,j)
        endif
        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 - abs(Kh) / (hrat_min*CS%Kh_Max_xx(i,j))
            visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xx(i,j))
          endif
        endif

        if (cs%id_Kh_h>0) kh_h(i,j,k) = kh

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

      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 = vort_mag * (cs%BIHARM_CONST_xx(i,j))
          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 (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) ah_h(i,j,k) = ah

        str_xx(i,j) = str_xx(i,j) + ah * &
          (cs%DY_dxT(i,j)*(g%IdyCu(i,j)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
           cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(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)*u0(i,j) - g%IdyCu(i-1,j)*u0(i-1,j)) - &
           cs%DX_dyT(i,j) *(g%IdxCv(i,j)*v0(i,j) - g%IdxCv(i,j-1)*v0(i,j-1)))
        bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))

      endif  ! biharmonic

      str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))
    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
        dvdx(i,j) = cs%DY_dxBu(i,j)*(v0(i+1,j)*g%IdyCv(i+1,j) - v0(i,j)*g%IdyCv(i,j))
        dudy(i,j) = cs%DX_dyBu(i,j)*(u0(i,j+1)*g%IdxCu(i,j+1) - u0(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
                dvdx(i,j) = 0. ; dudy(i,j) = 0.
              elseif (obc%freeslip_strain) then
                dudy(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
                dvdx(i,j) = 0. ; dudy(i,j) = 0.
              elseif (obc%freeslip_strain) then
                dvdx(i,j) = 0.
              endif
            enddo
          endif
        enddo
      endif ; endif
    endif

    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)) &
        vort_mag = sqrt( &
          0.5*((vort_xy_dx(i,j)*vort_xy_dx(i,j) + vort_xy_dx(i+1,j)*vort_xy_dx(i+1,j)) + &
                (vort_xy_dy(i,j)*vort_xy_dy(i,j) + vort_xy_dy(i,j+1)*vort_xy_dy(i,j+1))) + &
          mod_leith*0.5*((div_xx_dx(i,j)*div_xx_dx(i,j) + div_xx_dx(i,j+1)*div_xx_dx(i,j+1)) + &
                (div_xx_dy(i,j)*div_xx_dy(i,j) + div_xx_dy(i+1,j)*div_xx_dy(i+1,j))))

      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 = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * &
!          ((h(i,j,k) + h(i+1,j+1,k)) + (h(i,j+1,k) + h(i+1,j,k))))
      hq = 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 + 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 = hu + hv
            hrat_min = 1.0
          else
            ! Both hu and hv are nonzero, so take the harmonic mean.
            hq = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
            hrat_min = min(1.0, min(hu, hv) / (hq + h_neglect) )
          endif
        endif
      endif

      if (cs%Laplacian) then
        ! Determine the Laplacian viscosity at q points, using the
        ! largest value from several parameterizations.
        kh_scale = 1.0 ; if (rescale_kh) kh_scale = varmix%Res_fn_q(i,j)
        khsm = 0.0; khlth = 0.0
        if ((cs%Smagorinsky_Kh) .or. (cs%Leith_Kh)) then
          if (cs%Smagorinsky_Kh) &
            khsm = cs%LAPLAC_CONST_xy(i,j) * shear_mag
          if (cs%Leith_Kh) &
            khlth = cs%LAPLAC3_CONST_xy(i,j) * vort_mag
          kh = kh_scale * max(max(cs%Kh_bg_xy(i,j), khsm), khlth)
          if (cs%bound_Kh .and. .not.cs%better_bound_Kh) &
            kh = min(kh, cs%Kh_Max_xy(i,j))
        else
          kh = kh_scale * cs%Kh_bg_xy(i,j)
        endif
        if (use_meke_ku) then
          kh = kh + 0.25*( (meke%Ku(i,j)+meke%Ku(i+1,j+1))    &
                          +(meke%Ku(i+1,j)+meke%Ku(i,j+1)) )
        endif
        ! Place a floor on the viscosity, if desired.
        kh = max(kh,cs%Kh_bg_min)

        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 - abs(Kh) / (hrat_min*CS%Kh_Max_xy(I,J))
            visc_bound_rem = 1.0 - kh / (hrat_min*cs%Kh_Max_xy(i,j))
          endif
        endif

        if (cs%id_Kh_q>0) kh_q(i,j,k) = kh

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

      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 =  vort_mag * (cs%BIHARM5_CONST_xy(i,j))
          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 (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) ah_q(i,j,k) = ah

        str_xy(i,j) = str_xy(i,j) + ah * ( dvdx(i,j) + dudy(i,j) )

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

      endif  ! biharmonic

      if (cs%no_slip) then
        str_xy(i,j) = str_xy(i,j) * (hq * cs%reduction_xy(i,j))
      else
        str_xy(i,j) = str_xy(i,j) * (hq * g%mask2dBu(i,j) * cs%reduction_xy(i,j))
      endif
    enddo ; enddo

    do j=js,je ; do i=isq,ieq
!  Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
      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)
      frictwork(i,j,k) = gv%H_to_kg_m2 * ( &
              (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
      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))))
          fath = fath ** meke%backscatter_Ro_pow ! f^n
          shear_mag = ( ( 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)
          meke%mom_src(i,j) = meke%mom_src(i,j) + gv%H_to_kg_m2 * (                   &
                ((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
      else
        do j=js,je ; do i=is,ie
         meke%mom_src(i,j) = meke%mom_src(i,j) + frictwork(i,j,k)
        enddo ; enddo
      endif
    endif ; endif

  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_Ah_h>0)      call post_data(cs%id_Ah_h, ah_h, 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_Kh_q>0)      call post_data(cs%id_Kh_q, kh_q, cs%diag)

  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


end subroutine horizontal_viscosity

!> This subroutine allocates space for and calculates static variables
!! used by this module. The metrics may be 0, 1, or 2-D arrays,
!! while fields like the background viscosities are 2-D arrays.
!! ALLOC is a macro defined in MOM_memory.h to either allocate
!! for dynamic memory, or do nothing when using static memory.
subroutine hor_visc_init(Time, G, param_file, diag, CS)
  type(time_type),         intent(in)    :: Time !< current model time.
  type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure.
  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

! This subroutine allocates space for and calculates static variables
! used by this module. The metrics may be 0, 1, or 2-D arrays,
! while fields like the background viscosities are 2-D arrays.
! ALLOC is a macro defined in MOM_memory.h to either allocate
! for dynamic memory, or do nothing when using static memory.
!
! Arguments:
!  (in)      Time       - current model time
!  (in)      G          - ocean grid structure
!  (in)      param_file - structure to parse for model parameter values
!  (in)      diag       - structure to regulate diagnostic output
!  (in/out)  CS         - pointer to the control structure for this module

  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, in m-2, with u0u, v0u, and v0v defined similarly.
  real :: grid_sp_h2       ! Harmonic mean of the squares of the grid
  real :: grid_sp_h3       ! Harmonic mean of the squares of the grid^(3/2)
  real :: grid_sp_q2       ! spacings at h and q points (m2)
  real :: grid_sp_q3       ! spacings at h and q points^(3/2) (m3)
  real :: Kh_Limit         ! A coefficient (1/s) 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 (1/s)
  real :: BoundCorConst    ! constant (s2/m2)
  real :: Ah_Limit         ! coefficient (1/s) used, along with the
                           ! grid spacing, to limit biharmonic viscosity
  real :: Kh               ! Lapacian horizontal viscosity (m2/s)
  real :: Ah               ! biharmonic horizontal viscosity (m4/s)
  real :: Kh_vel_scale     ! this speed (m/s) times grid spacing gives Lap visc
  real :: Ah_vel_scale     ! this speed (m/s) times grid spacing cubed gives bih visc
  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               ! dynamics time step (sec)
  real :: Idt              ! inverse of dt (1/s)
  real :: denom            ! work variable; the denominator of a fraction
  real :: maxvel           ! largest permitted velocity components (m/s)
  real :: bound_Cor_vel    ! grid-scale velocity variations at which value
                           ! the quadratically varying biharmonic viscosity
                           ! balances Coriolis acceleration (m/s)
  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.
  character(len=64) :: inputdir, filename
  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%bound_Coriolis = .false.
  cs%Modified_Leith = .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, "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)
    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)
    call get_param(param_file, mdl, "KH_VEL_SCALE", kh_vel_scale, &
                 "The velocity scale which is multiplied by the grid \n"//&
                 "spacing to calculate the Laplacian viscosity. \n"//&
                 "The final viscosity is the largest of this scaled \n"//&
                 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
                 units="m s-1", default=0.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, \n"//&
                 "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 \n"//&
                 "proportional to the gradient of divergence.", &
                 default=.false.)

    if (cs%Leith_Kh .or. get_all) &
      call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
                 "The nondimensional Laplacian Leith constant, \n"//&
                 "often ??", units="nondim", default=0.0, &
                  fail_if_missing = cs%Leith_Kh)

    call get_param(param_file, mdl, "BOUND_KH", cs%bound_Kh, &
                 "If true, the Laplacian coefficient is locally limited \n"//&
                 "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 \n"//&
                 "to be stable with a better bounding than just BOUND_KH.", &
                 default=cs%bound_Kh)
  endif

  call get_param(param_file, mdl, "BIHARMONIC", cs%biharmonic, &
                 "If true, use a biharmonic horizontal viscosity. \n"//&
                 "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)
    call get_param(param_file, mdl, "AH_VEL_SCALE", ah_vel_scale, &
                 "The velocity scale which is multiplied by the cube of \n"//&
                 "the grid spacing to calculate the biharmonic viscosity. \n"//&
                 "The final viscosity is the largest of this scaled \n"//&
                 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
                 units="m s-1", default=0.0)
    call get_param(param_file, mdl, "SMAGORINSKY_AH", cs%Smagorinsky_Ah, &
                 "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//&
                 "viscosity.", default=.false.)
    call get_param(param_file, mdl, "LEITH_AH", cs%Leith_Ah, &
                 "If true, use a biharmonic Leith nonlinear eddy \n"//&
                 "viscosity.", default=.false.)

    call get_param(param_file, mdl, "BOUND_AH", cs%bound_Ah, &
                 "If true, the biharmonic coefficient is locally limited \n"//&
                 "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 \n"//&
                 "to be stable with a better bounding than just BOUND_AH.", &
                 default=cs%bound_Ah)

    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, \n"//&
                 "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 \n"//&
                 "of the velocity shears, so that the resulting viscous \n"//&
                 "drag is of comparable magnitude to the Coriolis terms \n"//&
                 "when the velocity differences between adjacent grid \n"//&
                 "points is 0.5*BOUND_CORIOLIS_VEL.  The default is the \n"//&
                 "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 \n"//&
                 "the biharmonic drag to have comparable magnitude to the \n"//&
                 "Coriolis acceleration.  The default is set by MAXVEL.", &
                 units="m s-1", default=maxvel)
      endif
    endif

    if (cs%Leith_Ah .or. get_all) then
      call get_param(param_file, mdl, "LEITH_BI_CONST",leith_bi_const, &
                 "The nondimensional biharmonic Leith constant, \n"//&
                 "typical values are thus far undetermined", units="nondim", default=0.0, &
                 fail_if_missing = cs%Leith_Ah)
    endif

  endif

  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 \n"//&
                 "viscosity bounds to the theoretical maximum for \n"//&
                 "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 \n"//&
                 "free slip boundary conditions are assumed. The \n"//&
                 "implementation of the free slip BCs on a C-grid is much \n"//&
                 "cleaner than the no slip BCs. The use of free slip BCs \n"//&
                 "is strongly encouraged, and no slip BCs are not used with \n"//&
                 "the biharmonic viscosity.", default=.false.)

  call get_param(param_file, mdl, "USE_KH_BG_2D", cs%use_Kh_bg_2d, &
                 "If true, read a file containing 2-d background harmonic  \n"//&
                 "viscosities. The final viscosity is the maximum of the other "//&
                 "terms and this background value.", default=.false.)


  if (cs%bound_Kh .or. cs%bound_Ah .or. cs%better_bound_Kh .or. cs%better_bound_Ah) &
    call get_param(param_file, mdl, "DT", dt, &
                 "The (baroclinic) dynamics time step.", units = "s", &
                 fail_if_missing=.true.)

  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

  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%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(isdb:iedb,jsdb:jedb)) ; 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%Laplac_Const_xx(isd:ied,jsd:jed))     ; cs%Laplac_Const_xx(:,:) = 0.0
      alloc_(cs%Laplac_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Laplac_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%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 read_data(trim(inputdir)//trim(filename), 'Kh', cs%Kh_bg_2d, &
                   domain=g%domain%mpp_domain, timelevel=1)
    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
    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%Biharm5_Const_xx(isd:ied,jsd:jed))     ; cs%Biharm5_Const_xx(:,:) = 0.0
      alloc_(cs%Biharm5_Const_xy(isdb:iedb,jsdb:jedb)) ; cs%Biharm5_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)
    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)
      if (cs%Smagorinsky_Kh) cs%LAPLAC_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

      cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))

      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))

      if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
        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
    enddo ; enddo
    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_Kh) cs%LAPLAC_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

      cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))

      if (cs%use_Kh_bg_2d) cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xy(i,j))

      if (cs%bound_Kh .and. .not.cs%better_bound_Kh) then
        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.
    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))
    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)

      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%BIHARM5_CONST_xx(i,j) = leith_bi_const * (grid_sp_h2 * grid_sp_h3)
      endif

      cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
      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
    enddo ; enddo
    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%BIHARM5_CONST_xy(i,j) = leith_bi_const * (grid_sp_q2 * grid_sp_q3)
      endif

      cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
      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
    idt = 1.0 / dt
    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
  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
    idt = 1.0 / dt
    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
  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', 'meter second-2')

  cs%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, time, &
      'Meridional Acceleration from Horizontal Viscosity', 'meter second-2')

  if (cs%biharmonic) then
    cs%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, time,    &
        'Biharmonic Horizontal Viscosity at h Points', 'meter4 second-1',        &
        cmor_field_name='difmxybo', cmor_units='m4 s-1',                        &
        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', 'meter4 second-1')
  endif

  if (cs%Laplacian) then
    cs%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, time,   &
        'Laplacian Horizontal Viscosity at h Points', 'meter2 second-1',        &
        cmor_field_name='difmxylo', cmor_units='m2 s-1',                        &
        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', 'meter2 second-1')
  endif

  cs%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,time,&
      'Integral work done by lateral friction terms', 'Watt meter-2')

  cs%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,time,      &
      'Depth integrated work done by lateral friction', 'Watt meter-2',                          &
      cmor_field_name='dispkexyfo', cmor_units='W m-2',                                          &
      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

subroutine hor_visc_end(CS)
! This subroutine deallocates any variables allocated in hor_visc_init.
! Argument:  CS - The control structure returned by a previous call to
!                 hor_visc_init.
  type(hor_visc_cs), pointer :: CS

  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)
    if (cs%bound_Kh) then
      dealloc_(cs%Kh_Max_xx) ; dealloc_(cs%Kh_Max_xy)
    endif
    if (cs%Smagorinsky_Kh) then
      dealloc_(cs%Laplac_Const_xx) ; dealloc_(cs%Laplac_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%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%Biharm_Const_xx) ; dealloc_(cs%Biharm_Const_xy)
      if (cs%bound_Coriolis) then
        dealloc_(cs%Biharm_Const2_xx) ; dealloc_(cs%Biharm_Const2_xy)
      endif
    endif
    if (cs%Leith_Ah) then
      dealloc_(cs%Biharm5_Const_xx) ; dealloc_(cs%Biharm5_Const_xy)
    endif
  endif
  deallocate(cs)

end subroutine hor_visc_end

end module mom_hor_visc