MOM_thickness_diffuse.F90

!> Thickness diffusion (or Gent McWilliams)
module mom_thickness_diffuse

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

use mom_debugging,             only : hchksum, uvchksum
use mom_diag_mediator,         only : post_data, query_averaging_enabled, diag_ctrl
use mom_diag_mediator,         only : register_diag_field, safe_alloc_ptr, time_type
use mom_diag_mediator,         only : diag_update_remap_grids
use mom_domains,               only : pass_var, corner, pass_vector
use mom_error_handler,         only : mom_error, fatal, warning, is_root_pe
use mom_eos,                   only : calculate_density, calculate_density_derivs, eos_domain
use mom_eos,                   only : calculate_density_second_derivs
use mom_file_parser,           only : get_param, log_version, param_file_type
use mom_grid,                  only : ocean_grid_type
use mom_interface_heights,     only : find_eta
use mom_isopycnal_slopes,      only : vert_fill_ts
use mom_lateral_mixing_coeffs, only : varmix_cs
use mom_meke_types,            only : meke_type
use mom_unit_scaling,          only : unit_scale_type
use mom_variables,             only : thermo_var_ptrs, cont_diag_ptrs
use mom_verticalgrid,          only : verticalgrid_type
implicit none ; private

#include <MOM_memory.h>

public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
! public vert_fill_TS
public thickness_diffuse_get_kh

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Control structure for thickness diffusion
type, public :: thickness_diffuse_cs ; private
  real    :: khth                !< Background interface depth diffusivity [L2 T-1 ~> m2 s-1]
  real    :: khth_slope_cff      !< Slope dependence coefficient of Khth [nondim]
  real    :: max_khth_cfl        !< Maximum value of the diffusive CFL for thickness diffusion
  real    :: khth_min            !< Minimum value of Khth [L2 T-1 ~> m2 s-1]
  real    :: khth_max            !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max
  real    :: slope_max           !< Slopes steeper than slope_max are limited in some way [nondim].
  real    :: kappa_smooth        !< Vertical diffusivity used to interpolate more
                                 !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1].
  logical :: thickness_diffuse   !< If true, interfaces heights are diffused.
  logical :: use_fgnv_streamfn   !< If true, use the streamfunction formulation of
                                 !! Ferrari et al., 2010, which effectively emphasizes
                                 !! graver vertical modes by smoothing in the vertical.
  real    :: fgnv_scale          !< A coefficient scaling the vertical smoothing term in the
                                 !! Ferrari et al., 2010, streamfunction formulation [nondim].
  real    :: fgnv_c_min          !< A minimum wave speed used in the Ferrari et al., 2010,
                                 !! streamfunction formulation [L T-1 ~> m s-1].
  real    :: n2_floor            !< A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010,
                                 !! streamfunction formulation [T-2 ~> s-2].
  logical :: detangle_interfaces !< If true, add 3-d structured interface height
                                 !! diffusivities to horizontally smooth jagged layers.
  real    :: detangle_time       !< If detangle_interfaces is true, this is the
                                 !! timescale over which maximally jagged grid-scale
                                 !! thickness variations are suppressed [T ~> s].  This must be
                                 !! longer than DT, or 0 (the default) to use DT.
  integer :: nkml                !< number of layers within mixed layer
  logical :: debug               !< write verbose checksums for debugging purposes
  logical :: use_gme_thickness_diffuse !< If true, passes GM coefficients to MOM_hor_visc for use
                                 !! with GME closure.
  logical :: meke_geometric      !< If true, uses the GM coefficient formulation from the GEOMETRIC
                                 !! framework (Marshall et al., 2012)
  real    :: meke_geometric_alpha!< The nondimensional coefficient governing the efficiency of
                                 !! the GEOMETRIC thickness difussion [nondim]
  real    :: meke_geometric_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness
                                 !! diffusivity [T-1 ~> s-1].
  logical :: meke_geom_answers_2018  !< If true, use expressions in the MEKE_GEOMETRIC calculation
                                 !! that recover the answers from the original implementation.
                                 !! Otherwise, use expressions that satisfy rotational symmetry.
  logical :: use_kh_in_meke      !< If true, uses the thickness diffusivity calculated here to diffuse MEKE.
  logical :: gm_src_alt          !< If true, use the GM energy conversion form S^2*N^2*kappa rather
                                 !! than the streamfunction for the GM source term.
  logical :: use_gm_work_bug     !< If true, use the incorrect sign for the
                                 !! top-level work tendency on the top layer.
  real :: stanley_det_coeff      !< The coefficient correlating SGS temperature variance with the mean
                                 !! temperature gradient in the deterministic part of the Stanley parameterization.
                                 !! Negative values disable the scheme." [nondim]

  type(diag_ctrl), pointer :: diag => null() !< structure used to regulate timing of diagnostics
  real, pointer :: gmwork(:,:)       => null()  !< Work by thickness diffusivity [R Z L2 T-3 ~> W m-2]
  real, pointer :: diagslopex(:,:,:) => null()  !< Diagnostic: zonal neutral slope [nondim]
  real, pointer :: diagslopey(:,:,:) => null()  !< Diagnostic: zonal neutral slope [nondim]

  real, dimension(:,:,:), pointer :: &
    kh_u_gme => null(), &        !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
    kh_v_gme => null()           !< interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]

  !>@{
  !! Diagnostic identifier
  integer :: id_uhgm    = -1, id_vhgm    = -1, id_gmwork = -1
  integer :: id_kh_u    = -1, id_kh_v    = -1, id_kh_t   = -1
  integer :: id_kh_u1   = -1, id_kh_v1   = -1, id_kh_t1  = -1
  integer :: id_slope_x = -1, id_slope_y = -1
  integer :: id_sfn_unlim_x = -1, id_sfn_unlim_y = -1, id_sfn_x = -1, id_sfn_y = -1
  !>@}
end type thickness_diffuse_cs

contains

!> Calculates thickness diffusion coefficients and applies thickness diffusion to layer
!! thicknesses, h. Diffusivities are limited to ensure stability.
!! Also returns along-layer mass fluxes used in the continuity equation.
subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS)
  type(ocean_grid_type),                     intent(in)    :: G      !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: GV     !< Vertical grid structure
  type(unit_scale_type),                     intent(in)    :: US     !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h      !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr   !< Accumulated zonal mass flux
                                                                     !! [L2 H ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr   !< Accumulated meridional mass flux
                                                                     !! [L2 H ~> m3 or kg]
  type(thermo_var_ptrs),                     intent(in)    :: tv     !< Thermodynamics structure
  real,                                      intent(in)    :: dt     !< Time increment [T ~> s]
  type(meke_type),                           pointer       :: MEKE   !< MEKE control structure
  type(varmix_cs),                           pointer       :: VarMix !< Variable mixing coefficients
  type(cont_diag_ptrs),                      intent(inout) :: CDp    !< Diagnostics for the continuity equation
  type(thickness_diffuse_cs),                pointer       :: CS     !< Control structure for thickness diffusion
  ! Local variables
  real :: e(szi_(g), szj_(g), szk_(g)+1) ! heights of interfaces, relative to mean
                                         ! sea level [Z ~> m], positive up.
  real :: uhD(szib_(g), szj_(g), szk_(g)) ! Diffusive u*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]
  real :: vhD(szi_(g), szjb_(g), szk_(g)) ! Diffusive v*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]

  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
    KH_u, &       ! interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
    int_slope_u   ! A nondimensional ratio from 0 to 1 that gives the relative
                  ! weighting of the interface slopes to that calculated also
                  ! using density gradients at u points.  The physically correct
                  ! slopes occur at 0, while 1 is used for numerical closures [nondim].
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
    KH_v, &       ! interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
    int_slope_v   ! A nondimensional ratio from 0 to 1 that gives the relative
                  ! weighting of the interface slopes to that calculated also
                  ! using density gradients at v points.  The physically correct
                  ! slopes occur at 0, while 1 is used for numerical closures [nondim].
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
    KH_t          ! diagnosed diffusivity at tracer points [L2 T-1 ~> m2 s-1]

  real, dimension(SZIB_(G), SZJ_(G)) :: &
    KH_u_CFL      ! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G), SZJB_(G)) :: &
    KH_v_CFL      ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1]
  real :: Khth_Loc_u(szib_(g), szj_(g))
  real :: Khth_Loc_v(szi_(g), szjb_(g))
  real :: Khth_Loc(szib_(g), szjb_(g))  ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1]
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected [H ~> m or kg m-2].
  real, dimension(:,:), pointer :: cg1 => null() !< Wave speed [L T-1 ~> m s-1]
  logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_ebt_struct, use_Visbeck
  logical :: use_QG_Leith
  integer :: i, j, k, is, ie, js, je, nz
  real :: hu(szi_(g), szj_(g))       ! u-thickness [H ~> m or kg m-2]
  real :: hv(szi_(g), szj_(g))       ! v-thickness [H ~> m or kg m-2]
  real :: KH_u_lay(szi_(g), szj_(g)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]
  real :: KH_v_lay(szi_(g), szj_(g)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]

  if (.not. associated(cs)) call mom_error(fatal, "MOM_thickness_diffuse: "//&
         "Module must be initialized before it is used.")

  if ((.not.cs%thickness_diffuse) .or. &
       .not.( cs%Khth > 0.0 .or. associated(varmix) .or. associated(meke) ) ) return

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  h_neglect = gv%H_subroundoff

  if (associated(meke)) then
    if (associated(meke%GM_src)) then
      do j=js,je ; do i=is,ie ; meke%GM_src(i,j) = 0. ; enddo ; enddo
    endif
  endif

  use_varmix = .false. ; resoln_scaled = .false. ; use_stored_slopes = .false.
  khth_use_ebt_struct = .false. ; use_visbeck = .false. ; use_qg_leith = .false.
  depth_scaled = .false.

  if (associated(varmix)) then
    use_varmix = varmix%use_variable_mixing .and. (cs%KHTH_Slope_Cff > 0.)
    resoln_scaled = varmix%Resoln_scaled_KhTh
    depth_scaled = varmix%Depth_scaled_KhTh
    use_stored_slopes = varmix%use_stored_slopes
    khth_use_ebt_struct = varmix%khth_use_ebt_struct
    use_visbeck = varmix%use_Visbeck
    use_qg_leith = varmix%use_QG_Leith_GM
    if (associated(varmix%cg1)) cg1 => varmix%cg1
  else
    cg1 => null()
  endif


!$OMP parallel do default(none) shared(is,ie,js,je,KH_u_CFL,dt,G,CS)
  do j=js,je ; do i=is-1,ie
    kh_u_cfl(i,j) = (0.25*cs%max_Khth_CFL) /  &
      (dt * (g%IdxCu(i,j)*g%IdxCu(i,j) + g%IdyCu(i,j)*g%IdyCu(i,j)))
  enddo ; enddo
!$OMP parallel do default(none) shared(is,ie,js,je,KH_v_CFL,dt,G,CS)
  do j=js-1,je ; do i=is,ie
    kh_v_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
      (dt * (g%IdxCv(i,j)*g%IdxCv(i,j) + g%IdyCv(i,j)*g%IdyCv(i,j)))
  enddo ; enddo

  ! Calculates interface heights, e, in [Z ~> m].
  call find_eta(h, tv, g, gv, us, e, halo_size=1)

  ! Set the diffusivities.
!$OMP parallel default(none) shared(is,ie,js,je,Khth_Loc_u,CS,use_VarMix,VarMix,    &
!$OMP                               MEKE,Resoln_scaled,KH_u,G,use_QG_Leith,use_Visbeck,&
!$OMP                               KH_u_CFL,nz,Khth_Loc,KH_v,KH_v_CFL,int_slope_u, &
!$OMP                               int_slope_v,khth_use_ebt_struct, Depth_scaled, &
!$OMP                               Khth_loc_v)
!$OMP do
  do j=js,je; do i=is-1,ie
    khth_loc_u(i,j) = cs%Khth
  enddo ; enddo

  if (use_varmix) then
    if (use_visbeck) then
!$OMP do
      do j=js,je ; do i=is-1,ie
        khth_loc_u(i,j) = khth_loc_u(i,j) + &
          cs%KHTH_Slope_Cff*varmix%L2u(i,j) * varmix%SN_u(i,j)
      enddo ; enddo
    endif
  endif

  if (associated(meke)) then ; if (associated(meke%Kh)) then
    if (cs%MEKE_GEOMETRIC) then
!$OMP do
      do j=js,je ; do i=is-1,ie
        khth_loc_u(i,j) = khth_loc_u(i,j) + g%mask2dCu(i,j) * cs%MEKE_GEOMETRIC_alpha * &
                          0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)) / &
                          (varmix%SN_u(i,j) + cs%MEKE_GEOMETRIC_epsilon)
      enddo ; enddo
    else
      do j=js,je ; do i=is-1,ie
        khth_loc_u(i,j) = khth_loc_u(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
      enddo ; enddo
    endif
  endif ; endif

  if (resoln_scaled) then
!$OMP do
    do j=js,je; do i=is-1,ie
      khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Res_fn_u(i,j)
    enddo ; enddo
  endif

  if (depth_scaled) then
!$OMP do
    do j=js,je; do i=is-1,ie
      khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Depth_fn_u(i,j)
    enddo ; enddo
  endif

  if (cs%Khth_Max > 0) then
!$OMP do
    do j=js,je; do i=is-1,ie
      khth_loc_u(i,j) = max(cs%Khth_Min, min(khth_loc_u(i,j), cs%Khth_Max))
    enddo ; enddo
  else
!$OMP do
    do j=js,je; do i=is-1,ie
      khth_loc_u(i,j) = max(cs%Khth_Min, khth_loc_u(i,j))
    enddo ; enddo
  endif
!$OMP do
  do j=js,je; do i=is-1,ie
    kh_u(i,j,1) = min(kh_u_cfl(i,j), khth_loc_u(i,j))
  enddo ; enddo

  if (khth_use_ebt_struct) then
!$OMP do
    do k=2,nz+1 ; do j=js,je ; do i=is-1,ie
      kh_u(i,j,k) = kh_u(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i+1,j,k-1) )
    enddo ; enddo ; enddo
  else
!$OMP do
    do k=2,nz+1 ; do j=js,je ; do i=is-1,ie
      kh_u(i,j,k) = kh_u(i,j,1)
    enddo ; enddo ; enddo
  endif

  if (use_varmix) then
    if (use_qg_leith) then
!$OMP do
      do k=1,nz ; do j=js,je ; do i=is-1,ie
        kh_u(i,j,k) = varmix%KH_u_QG(i,j,k)
      enddo ; enddo ; enddo
    endif
  endif

  if (cs%use_GME_thickness_diffuse) then
!$OMP do
    do k=1,nz+1 ; do j=js,je ; do i=is-1,ie
      cs%KH_u_GME(i,j,k) = kh_u(i,j,k)
    enddo ; enddo ; enddo
  endif

!$OMP do
  do j=js-1,je ; do i=is,ie
    khth_loc_v(i,j) = cs%Khth
  enddo ; enddo

  if (use_varmix) then
    if (use_visbeck) then
!$OMP do
      do j=js-1,je ; do i=is,ie
        khth_loc_v(i,j) = khth_loc_v(i,j) + cs%KHTH_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
      enddo ; enddo
    endif
  endif
  if (associated(meke)) then ; if (associated(meke%Kh)) then
    if (cs%MEKE_GEOMETRIC) then
!$OMP do
      do j=js-1,je ; do i=is,ie
        khth_loc_v(i,j) = khth_loc_v(i,j) + g%mask2dCv(i,j) * cs%MEKE_GEOMETRIC_alpha * &
                        0.5*(meke%MEKE(i,j)+meke%MEKE(i,j+1)) / &
                        (varmix%SN_v(i,j) + cs%MEKE_GEOMETRIC_epsilon)
      enddo ; enddo
    else
      do j=js-1,je ; do i=is,ie
        khth_loc_v(i,j) = khth_loc_v(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
      enddo ; enddo
    endif
  endif ; endif

  if (resoln_scaled) then
!$OMP do
    do j=js-1,je; do i=is,ie
      khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Res_fn_v(i,j)
    enddo ; enddo
  endif

  if (depth_scaled) then
!$OMP do
    do j=js-1,je ; do i=is,ie
      khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Depth_fn_v(i,j)
    enddo ; enddo
  endif

  if (cs%Khth_Max > 0) then
!$OMP do
    do j=js-1,je ; do i=is,ie
      khth_loc_v(i,j) = max(cs%Khth_Min, min(khth_loc_v(i,j), cs%Khth_Max))
    enddo ; enddo
  else
!$OMP do
    do j=js-1,je ; do i=is,ie
      khth_loc_v(i,j) = max(cs%Khth_Min, khth_loc_v(i,j))
    enddo ; enddo
  endif

  if (cs%max_Khth_CFL > 0.0) then
!$OMP do
    do j=js-1,je ; do i=is,ie
      kh_v(i,j,1) = min(kh_v_cfl(i,j), khth_loc_v(i,j))
    enddo ; enddo
  endif

  if (khth_use_ebt_struct) then
!$OMP do
    do k=2,nz+1 ; do j=js-1,je ; do i=is,ie
      kh_v(i,j,k) = kh_v(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i,j+1,k-1) )
    enddo ; enddo ; enddo
  else
!$OMP do
    do k=2,nz+1 ; do j=js-1,je ; do i=is,ie
      kh_v(i,j,k) = kh_v(i,j,1)
    enddo ; enddo ; enddo
  endif

  if (use_varmix) then
    if (use_qg_leith) then
!$OMP do
      do k=1,nz ; do j=js-1,je ; do i=is,ie
        kh_v(i,j,k) = varmix%KH_v_QG(i,j,k)
      enddo ; enddo ; enddo
    endif
  endif

  if (cs%use_GME_thickness_diffuse) then
!$OMP do
    do k=1,nz+1 ; do j=js-1,je ; do i=is,ie
      cs%KH_v_GME(i,j,k) = kh_v(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (associated(meke)) then ; if (associated(meke%Kh)) then
    if (cs%MEKE_GEOMETRIC) then
      if (cs%MEKE_GEOM_answers_2018) then
        !$OMP do
        do j=js,je ; do i=is,ie
          ! This does not give bitwise rotational symmetry.
          meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
                         (0.25*(varmix%SN_u(i,j)+varmix%SN_u(i-1,j) + &
                                varmix%SN_v(i,j)+varmix%SN_v(i,j-1)) + &
                          cs%MEKE_GEOMETRIC_epsilon)
        enddo ; enddo
      else
        !$OMP do
        do j=js,je ; do i=is,ie
          ! With the additional parentheses this gives bitwise rotational symmetry.
          meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
                         (0.25*((varmix%SN_u(i,j)+varmix%SN_u(i-1,j)) + &
                                (varmix%SN_v(i,j)+varmix%SN_v(i,j-1))) + &
                          cs%MEKE_GEOMETRIC_epsilon)
        enddo ; enddo
      endif
    endif
  endif ; endif


!$OMP do
  do k=1,nz+1 ; do j=js,je ; do i=is-1,ie ; int_slope_u(i,j,k) = 0.0 ; enddo ; enddo ; enddo
!$OMP do
  do k=1,nz+1 ; do j=js-1,je ; do i=is,ie ; int_slope_v(i,j,k) = 0.0 ; enddo ; enddo ; enddo
!$OMP end parallel

  if (cs%detangle_interfaces) then
    call add_detangling_kh(h, e, kh_u, kh_v, kh_u_cfl, kh_v_cfl, tv, dt, g, gv, us, &
                           cs, int_slope_u, int_slope_v)
  endif

  if (cs%debug) then
    call uvchksum("Kh_[uv]", kh_u, kh_v, g%HI, haloshift=0, &
                  scale=(us%L_to_m**2)*us%s_to_T, scalar_pair=.true.)
    call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, g%HI, haloshift=0)
    call hchksum(h, "thickness_diffuse_1 h", g%HI, haloshift=1, scale=gv%H_to_m)
    call hchksum(e, "thickness_diffuse_1 e", g%HI, haloshift=1, scale=us%Z_to_m)
    if (use_stored_slopes) then
      call uvchksum("VarMix%slope_[xy]", varmix%slope_x, varmix%slope_y, &
                    g%HI, haloshift=0)
    endif
    if (associated(tv%eqn_of_state)) then
      call hchksum(tv%T, "thickness_diffuse T", g%HI, haloshift=1)
      call hchksum(tv%S, "thickness_diffuse S", g%HI, haloshift=1)
    endif
  endif

  ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S
  if (use_stored_slopes) then
    call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
                                int_slope_u, int_slope_v, varmix%slope_x, varmix%slope_y)
  else
    call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
                                int_slope_u, int_slope_v)
  endif

  if (associated(meke) .AND. associated(varmix)) then
    if (associated(meke%Rd_dx_h) .and. associated(varmix%Rd_dx_h)) then
!$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix)
      do j=js,je ; do i=is,ie
        meke%Rd_dx_h(i,j) = varmix%Rd_dx_h(i,j)
      enddo ; enddo
    endif
  endif

  ! offer diagnostic fields for averaging
  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_uhGM > 0)   call post_data(cs%id_uhGM, uhd, cs%diag)
    if (cs%id_vhGM > 0)   call post_data(cs%id_vhGM, vhd, cs%diag)
    if (cs%id_GMwork > 0) call post_data(cs%id_GMwork, cs%GMwork, cs%diag)
    if (cs%id_KH_u > 0)   call post_data(cs%id_KH_u, kh_u, cs%diag)
    if (cs%id_KH_v > 0)   call post_data(cs%id_KH_v, kh_v, cs%diag)
    if (cs%id_KH_u1 > 0)  call post_data(cs%id_KH_u1, kh_u(:,:,1), cs%diag)
    if (cs%id_KH_v1 > 0)  call post_data(cs%id_KH_v1, kh_v(:,:,1), cs%diag)

    ! Diagnose diffusivity at T-cell point.  Do simple average, rather than
    ! thickness-weighted average, in order that KH_t is depth-independent
    ! in the case where KH_u and KH_v are depth independent.  Otherwise,
    ! if use thickness weighted average, the variations of thickness with
    ! depth will place a spurious depth dependence to the diagnosed KH_t.
    if (cs%id_KH_t > 0 .or. cs%id_KH_t1 > 0 .or. cs%Use_KH_in_MEKE) then
      do k=1,nz
        ! thicknesses across u and v faces, converted to 0/1 mask
        ! layer average of the interface diffusivities KH_u and KH_v
        do j=js,je ; do i=is-1,ie
          hu(i,j)       = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect)
          if (hu(i,j) /= 0.0) hu(i,j) = 1.0
          kh_u_lay(i,j) = 0.5*(kh_u(i,j,k)+kh_u(i,j,k+1))
        enddo ; enddo
        do j=js-1,je ; do i=is,ie
          hv(i,j)       = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
          if (hv(i,j) /= 0.0) hv(i,j) = 1.0
          kh_v_lay(i,j) = 0.5*(kh_v(i,j,k)+kh_v(i,j,k+1))
        enddo ; enddo
        ! diagnose diffusivity at T-point
        do j=js,je ; do i=is,ie
          kh_t(i,j,k) = ((hu(i-1,j)*kh_u_lay(i-1,j)+hu(i,j)*kh_u_lay(i,j))  &
                        +(hv(i,j-1)*kh_v_lay(i,j-1)+hv(i,j)*kh_v_lay(i,j))) &
                       / (hu(i-1,j)+hu(i,j)+hv(i,j-1)+hv(i,j)+h_neglect)
        enddo ; enddo
      enddo

      if (cs%Use_KH_in_MEKE) then
        meke%Kh_diff(:,:) = 0.0
        do k=1,nz
          do j=js,je ; do i=is,ie
            meke%Kh_diff(i,j) = meke%Kh_diff(i,j) + kh_t(i,j,k) * h(i,j,k)
          enddo; enddo
        enddo

        do j=js,je ; do i=is,ie
          meke%Kh_diff(i,j) = meke%Kh_diff(i,j) / max(1.0,g%bathyT(i,j))
        enddo ; enddo
      endif

      if (cs%id_KH_t  > 0) call post_data(cs%id_KH_t,  kh_t,        cs%diag)
      if (cs%id_KH_t1 > 0) call post_data(cs%id_KH_t1, kh_t(:,:,1), cs%diag)
    endif

  endif

  !$OMP parallel do default(none) shared(is,ie,js,je,nz,uhtr,uhD,dt,vhtr,CDp,vhD,h,G,GV)
  do k=1,nz
    do j=js,je ; do i=is-1,ie
      uhtr(i,j,k) = uhtr(i,j,k) + uhd(i,j,k) * dt
      if (associated(cdp%uhGM)) cdp%uhGM(i,j,k) = uhd(i,j,k)
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      vhtr(i,j,k) = vhtr(i,j,k) + vhd(i,j,k) * dt
      if (associated(cdp%vhGM)) cdp%vhGM(i,j,k) = vhd(i,j,k)
    enddo ; enddo
    do j=js,je ; do i=is,ie
      h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * &
          ((uhd(i,j,k) - uhd(i-1,j,k)) + (vhd(i,j,k) - vhd(i,j-1,k)))
      if (h(i,j,k) < gv%Angstrom_H) h(i,j,k) = gv%Angstrom_H
    enddo ; enddo
  enddo

  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  ! This needs to happen after the H update and before the next post_data.
  call diag_update_remap_grids(cs%diag)

  if (cs%debug) then
    call uvchksum("thickness_diffuse [uv]hD", uhd, vhd, &
                  g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
    call uvchksum("thickness_diffuse [uv]htr", uhtr, vhtr, &
                  g%HI, haloshift=0, scale=us%L_to_m**2*gv%H_to_m)
    call hchksum(h, "thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
  endif

end subroutine thickness_diffuse

!> Calculates parameterized layer transports for use in the continuity equation.
!! Fluxes are limited to give positive definite thicknesses.
!! Called by thickness_diffuse().
subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, &
                                  CS, int_slope_u, int_slope_v, slope_x, slope_y)
  type(ocean_grid_type),                       intent(in)  :: G      !< Ocean grid structure
  type(verticalgrid_type),                     intent(in)  :: GV     !< Vertical grid structure
  type(unit_scale_type),                       intent(in)  :: US     !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)  :: h      !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1),  intent(in)  :: e      !< Interface positions [Z ~> m]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in)  :: Kh_u   !< Thickness diffusivity on interfaces
                                                                     !! at u points [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in)  :: Kh_v   !< Thickness diffusivity on interfaces
                                                                     !! at v points [L2 T-1 ~> m2 s-1]
  type(thermo_var_ptrs),                       intent(in)  :: tv     !< Thermodynamics structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)),   intent(out) :: uhD    !< Zonal mass fluxes
                                                                     !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)),   intent(out) :: vhD    !< Meridional mass fluxes
                                                                     !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(:,:),                        pointer     :: cg1    !< Wave speed [L T-1 ~> m s-1]
  real,                                        intent(in)  :: dt     !< Time increment [T ~> s]
  type(meke_type),                             pointer     :: MEKE   !< MEKE control structure
  type(thickness_diffuse_cs),                  pointer     :: CS     !< Control structure for thickness diffusion
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in)  :: int_slope_u !< Ratio that determine how much of
                                                                     !! the isopycnal slopes are taken directly from
                                                                     !! the interface slopes without consideration of
                                                                     !! density gradients [nondim].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in)  :: int_slope_v !< Ratio that determine how much of
                                                                     !! the isopycnal slopes are taken directly from
                                                                     !! the interface slopes without consideration of
                                                                     !! density gradients [nondim].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in)  :: slope_x !< Isopycnal slope at u-points
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in)  :: slope_y !< Isopycnal slope at v-points
  ! Local variables
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
    T, &          ! The temperature (or density) [degC], with the values in
                  ! in massless layers filled vertically by diffusion.
    s, &          ! The filled salinity [ppt], with the values in
                  ! in massless layers filled vertically by diffusion.
    h_avail, &    ! The mass available for diffusion out of each face, divided
                  ! by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
    h_frac        ! The fraction of the mass in the column above the bottom
                  ! interface of a layer that is within a layer [nondim]. 0<h_frac<=1
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
    Slope_y_PE, &  ! 3D array of neutral slopes at v-points, set equal to Slope (below, nondim)
    hN2_y_PE       ! thickness in m times Brunt-Vaisala freqeuncy at v-points [L2 Z-1 T-2 ~> m s-2],
                   ! used for calculating PE release
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
    Slope_x_PE, &  ! 3D array of neutral slopes at u-points, set equal to Slope (below, nondim)
    hN2_x_PE       ! thickness in m times Brunt-Vaisala freqeuncy at u-points [L2 Z-1 T-2 ~> m s-2],
                   ! used for calculating PE release
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    pres, &       ! The pressure at an interface [R L2 T-2 ~> Pa].
    h_avail_rsum  ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)) :: &
    drho_dT_u, &  ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1]
    drho_dS_u, &  ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1].
    drho_dT_dT_u  ! The second derivative of density with temperature at u points [R degC-2 ~> kg m-3 degC-2]
  real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored.
  real, dimension(SZI_(G)) :: &
    drho_dT_v, &  ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1]
    drho_dS_v, &  ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1].
    drho_dT_dT_v  ! The second derivative of density with temperature at v points [R degC-2 ~> kg m-3 degC-2]
  real :: uhtot(szib_(g), szj_(g))  ! The vertical sum of uhD [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: vhtot(szi_(g), szjb_(g))  ! The vertical sum of vhD [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G)) :: &
    T_u, &        ! Temperature on the interface at the u-point [degC].
    S_u, &        ! Salinity on the interface at the u-point [ppt].
    pres_u        ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
  real, dimension(SZI_(G)) :: &
    T_v, &        ! Temperature on the interface at the v-point [degC].
    S_v, &        ! Salinity on the interface at the v-point [ppt].
    pres_v        ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
  real :: Work_u(szib_(g), szj_(g)) ! The work being done by the thickness
  real :: Work_v(szi_(g), szjb_(g)) ! diffusion integrated over a cell [R Z L4 T-3  ~> W ]
  real :: Work_h        ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2].
  real :: PE_release_h  ! The amount of potential energy released by GM averaged over an h-cell [L4 Z-1 T-3 ~> m3 s-3]
                        ! The calculation is equal to h * S^2 * N^2 * kappa_GM.
  real :: I4dt          ! 1 / 4 dt [T-1 ~> s-1].
  real :: drdiA, drdiB  ! Along layer zonal- and meridional- potential density
  real :: drdjA, drdjB  ! gradients in the layers above (A) and below(B) the
                        ! interface times the grid spacing [R ~> kg m-3].
  real :: drdkL, drdkR  ! Vertical density differences across an interface [R ~> kg m-3].
  real :: drdi_u(szib_(g), szk_(g)) ! Copy of drdi at u-points [R ~> kg m-3].
  real :: drdj_v(szi_(g), szk_(g))  ! Copy of drdj at v-points [R ~> kg m-3].
  real :: drdkDe_u(szib_(g),szk_(g)+1) ! Lateral difference of product of drdk and e at u-points
                                       ! [Z R ~> kg m-2].
  real :: drdkDe_v(szi_(g),szk_(g)+1)  ! Lateral difference of product of drdk and e at v-points
                                       ! [Z R ~> kg m-2].
  real :: hg2A, hg2B, hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
  real :: haA, haB, haL, haR     ! Arithmetic mean thicknesses [H ~> m or kg m-2].
  real :: dzaL, dzaR    ! Temporary thicknesses [Z ~> m].
  real :: wtA, wtB, wtL, wtR  ! Unscaled weights, with various units.
  real :: drdx, drdy    ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
  real :: drdz          ! Vertical density gradient [R Z-1 ~> kg m-4].
  real :: h_harm        ! Harmonic mean layer thickness [H ~> m or kg m-2].
  real :: c2_h_u(szib_(g), szk_(g)+1) ! Wave speed squared divided by h at u-points [L2 Z-1 T-2 ~> m s-2].
  real :: c2_h_v(szi_(g), szk_(g)+1)  ! Wave speed squared divided by h at v-points [L2 Z-1 T-2 ~> m s-2].
  real :: hN2_u(szib_(g), szk_(g)+1)  ! Thickness in m times N2 at interfaces above u-points [L2 Z-1 T-2 ~> m s-2].
  real :: hN2_v(szi_(g), szk_(g)+1)   ! Thickness in m times N2 at interfaces above v-points [L2 Z-1 T-2 ~> m s-2].
  real :: Sfn_est       ! A preliminary estimate (before limiting) of the overturning
                        ! streamfunction [Z L2 T-1 ~> m3 s-1].
  real :: Sfn_unlim_u(szib_(g), szk_(g)+1) ! Streamfunction for u-points [Z L2 T-1 ~> m3 s-1].
  real :: Sfn_unlim_v(szi_(g), szk_(g)+1)  ! Streamfunction for v-points [Z L2 T-1 ~> m3 s-1].
  real :: slope2_Ratio_u(szib_(g), szk_(g)+1) ! The ratio of the slope squared to slope_max squared.
  real :: slope2_Ratio_v(szi_(g), szk_(g)+1)  ! The ratio of the slope squared to slope_max squared.
  real :: Sfn_in_h      ! The overturning streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] (note that
                        ! the units are different from other Sfn vars).
  real :: Sfn_safe      ! The streamfunction that goes linearly back to 0 at the surface.  This is a
                        ! good thing to use when the slope is so large as to be meaningless [Z L2 T-1 ~> m3 s-1].
  real :: Slope         ! The slope of density surfaces, calculated in a way
                        ! that is always between -1 and 1, nondimensional.
  real :: mag_grad2     ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
  real :: I_slope_max2  ! The inverse of slope_max squared, nondimensional.
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected [H ~> m or kg m-2].
  real :: h_neglect2    ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
  real :: dz_neglect    ! A thickness [Z ~> m], that is so small it is usually lost
                        ! in roundoff and can be neglected [Z ~> m].
  real :: G_scale       ! The gravitational acceleration times a unit conversion
                        ! factor [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  logical :: find_work  ! If true, find the change in energy due to the fluxes.
  integer :: nk_linear  ! The number of layers over which the streamfunction goes to 0.
  real :: G_rho0        ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2].
  real :: N2_floor      ! A floor for N2 to avoid degeneracy in the elliptic solver
                        ! times unit conversion factors [T-2 L2 Z-2 ~> s-2]
  real :: Tl(5)         ! copy and T in local stencil [degC]
  real :: mn_T          ! mean of T in local stencil [degC]
  real :: mn_T2         ! mean of T**2 in local stencil [degC]
  real :: hl(5)         ! Copy of local stencil of H [H ~> m]
  real :: r_sm_H        ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tsgs2 ! Sub-grid temperature variance [degC2]

  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: diag_sfn_x, diag_sfn_unlim_x ! Diagnostics
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics
  logical :: present_int_slope_u, present_int_slope_v
  logical :: present_slope_x, present_slope_y, calc_derivatives
  integer, dimension(2) ::  EOSdom_u ! The shifted i-computational domain to use for equation of
                                     ! state calculations at u-points.
  integer, dimension(2) ::  EOSdom_v ! The shifted I-computational domain to use for equation of
                                     ! state calculations at v-points.
  logical :: use_Stanley
  integer :: is, ie, js, je, nz, IsdB, halo
  integer :: i, j, k
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke ; isdb = g%IsdB

  i4dt = 0.25 / dt
  i_slope_max2 = 1.0 / (cs%slope_max**2)
  g_scale = gv%g_Earth * gv%H_to_Z

  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
  dz_neglect = gv%H_subroundoff*gv%H_to_Z
  g_rho0 = gv%g_Earth / gv%Rho0
  n2_floor = cs%N2_floor*us%Z_to_L**2

  use_eos = associated(tv%eqn_of_state)
  present_int_slope_u = PRESENT(int_slope_u)
  present_int_slope_v = PRESENT(int_slope_v)
  present_slope_x = PRESENT(slope_x)
  present_slope_y = PRESENT(slope_y)
  use_stanley = cs%Stanley_det_coeff >= 0.

  nk_linear = max(gv%nkml, 1)

  slope_x_pe(:,:,:) = 0.0
  slope_y_pe(:,:,:) = 0.0
  hn2_x_pe(:,:,:) = 0.0
  hn2_y_pe(:,:,:) = 0.0

  find_work = .false.
  if (associated(meke)) find_work = associated(meke%GM_src)
  find_work = (associated(cs%GMwork) .or. find_work)

  if (use_eos) then
    halo = 1 ! Default halo to fill is 1
    if (use_stanley) halo = 2 ! Need wider valid halo for gradients of T
    call vert_fill_ts(h, tv%T, tv%S, cs%kappa_smooth*dt, t, s, g, gv, halo, larger_h_denom=.true.)
  endif

  if (cs%use_FGNV_streamfn .and. .not. associated(cg1)) call mom_error(fatal, &
       "cg1 must be associated when using FGNV streamfunction.")

!$OMP parallel default(none) shared(is,ie,js,je,h_avail_rsum,pres,h_avail,I4dt, use_Stanley, &
!$OMP                               CS,G,GV,tv,h,h_frac,nz,uhtot,Work_u,vhtot,Work_v,Tsgs2,T, &
!$OMP                               diag_sfn_x, diag_sfn_y, diag_sfn_unlim_x, diag_sfn_unlim_y ) &
!$OMP          private(hl,r_sm_H,Tl,mn_T,mn_T2)
  ! Find the maximum and minimum permitted streamfunction.
!$OMP do
  do j=js-1,je+1 ; do i=is-1,ie+1
    h_avail_rsum(i,j,1) = 0.0
    pres(i,j,1) = 0.0
    if (associated(tv%p_surf)) then ; pres(i,j,1) = tv%p_surf(i,j) ; endif

    h_avail(i,j,1) = max(i4dt*g%areaT(i,j)*(h(i,j,1)-gv%Angstrom_H),0.0)
    h_avail_rsum(i,j,2) = h_avail(i,j,1)
    h_frac(i,j,1) = 1.0
    pres(i,j,2) = pres(i,j,1) + (gv%g_Earth*gv%H_to_RZ) * h(i,j,1)
  enddo ; enddo
  if (use_stanley) then
!$OMP do
    do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1
      !! SGS variance in i-direction [degC2]
      !dTdi2 = ( ( G%mask2dCu(I  ,j) * G%IdxCu(I  ,j) * ( T(i+1,j,k) - T(i,j,k) ) &
      !          + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) &
      !          ) * G%dxT(i,j) * 0.5 )**2
      !! SGS variance in j-direction [degC2]
      !dTdj2 = ( ( G%mask2dCv(i,J  ) * G%IdyCv(i,J  ) * ( T(i,j+1,k) - T(i,j,k) ) &
      !          + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) &
      !          ) * G%dyT(i,j) * 0.5 )**2
      !Tsgs2(i,j,k) = CS%Stanley_det_coeff * 0.5 * ( dTdi2 + dTdj2 )
      ! This block does a thickness weighted variance calculation and helps control for
      ! extreme gradients along layers which are vanished against topography. It is
      ! still a poor approximation in the interior when coordinates are strongly tilted.
      hl(1) = h(i,j,k) * g%mask2dT(i,j)
      hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
      hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
      hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
      hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
      r_sm_h = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + gv%H_subroundoff )
      ! Mean of T
      tl(1) = t(i,j,k) ; tl(2) = t(i-1,j,k) ; tl(3) = t(i+1,j,k)
      tl(4) = t(i,j-1,k) ; tl(5) = t(i,j+1,k)
      mn_t = ( hl(1)*tl(1) + ( ( hl(2)*tl(2) + hl(3)*tl(3) ) + ( hl(4)*tl(4) + hl(5)*tl(5) ) ) ) * r_sm_h
      ! Adjust T vectors to have zero mean
      tl(:) = tl(:) - mn_t ; mn_t = 0.
      ! Variance of T
      mn_t2 = ( hl(1)*tl(1)*tl(1) + ( ( hl(2)*tl(2)*tl(2) + hl(3)*tl(3)*tl(3) ) &
                                    + ( hl(4)*tl(4)*tl(4) + hl(5)*tl(5)*tl(5) ) ) ) * r_sm_h
      ! Variance should be positive but round-off can violate this. Calculating
      ! variance directly would fix this but requires more operations.
      tsgs2(i,j,k) = cs%Stanley_det_coeff * max(0., mn_t2)
     enddo ; enddo ; enddo
  endif
!$OMP do
  do j=js-1,je+1
    do k=2,nz ; do i=is-1,ie+1
      h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
      h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
      h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) &
        h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
      pres(i,j,k+1) = pres(i,j,k) + (gv%g_Earth*gv%H_to_RZ) * h(i,j,k)
    enddo ; enddo
  enddo
!$OMP do
  do j=js,je ; do i=is-1,ie
    uhtot(i,j) = 0.0 ; work_u(i,j) = 0.0
    diag_sfn_x(i,j,1) = 0.0 ; diag_sfn_unlim_x(i,j,1) = 0.0
    diag_sfn_x(i,j,nz+1) = 0.0 ; diag_sfn_unlim_x(i,j,nz+1) = 0.0
  enddo ; enddo
!$OMP do
  do j=js-1,je ; do i=is,ie
    vhtot(i,j) = 0.0 ; work_v(i,j) = 0.0
    diag_sfn_y(i,j,1) = 0.0 ; diag_sfn_unlim_y(i,j,1) = 0.0
    diag_sfn_y(i,j,nz+1) = 0.0 ; diag_sfn_unlim_y(i,j,nz+1) = 0.0
  enddo ; enddo
!$OMP end parallel

    eosdom_u(1) = (is-1) - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
!$OMP                                  nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,  &
!$OMP                                  I_slope_max2,h_neglect2,present_int_slope_u, &
!$OMP                                  int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum,  &
!$OMP                                  uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1,   &
!$OMP                                  diag_sfn_x, diag_sfn_unlim_x,N2_floor,EOSdom_u, &
!$OMP                                  use_stanley, Tsgs2,                          &
!$OMP                                  present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE)  &
!$OMP                          private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u,      &
!$OMP                                  drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP                                  drho_dT_dT_u,scrap,                          &
!$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &
!$OMP                                  drdx,mag_grad2,Slope,slope2_Ratio_u,hN2_u,   &
!$OMP                                  Sfn_unlim_u,drdi_u,drdkDe_u,h_harm,c2_h_u,   &
!$OMP                                  Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
  do j=js,je
    do i=is-1,ie ; hn2_u(i,1) = 0. ; hn2_u(i,nz+1) = 0. ; enddo
    do k=nz,2,-1
      if (find_work .and. .not.(use_eos)) then
        drdia = 0.0 ; drdib = 0.0
        drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
      endif

      calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
         (find_work .or. .not. present_slope_x .or. cs%use_FGNV_streamfn .or. use_stanley)

      ! Calculate the zonal fluxes and gradients.
      if (calc_derivatives) then
        do i=is-1,ie
          pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
          t_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
          s_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
        enddo
        call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
                                      tv%eqn_of_state, eosdom_u)
      endif
      if (use_stanley) then
        ! The second line below would correspond to arguments
        !            drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
        call calculate_density_second_derivs(t_u, s_u, pres_u, &
                     scrap, scrap, drho_dt_dt_u, scrap, scrap, &
                     (is-isdb+1)-1, ie-is+2, tv%eqn_of_state)
      endif

      do i=is-1,ie
        if (calc_derivatives) then
          ! Estimate the horizontal density gradients along layers.
          drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
                  drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
          drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
                  drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))

          ! Estimate the vertical density gradients times the grid spacing.
          drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
                   drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
          drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
                   drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
          drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
        elseif (find_work) then ! This is used in pure stacked SW mode
          drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
        endif
        if (use_stanley) then
          ! Correction to the horizontal density gradient due to nonlinearity in
          ! the EOS rectifying SGS temperature anomalies
          drdia = drdia + drho_dt_dt_u(i) * 0.5 * ( tsgs2(i+1,j,k-1)-tsgs2(i,j,k-1) )
          drdib = drdib + drho_dt_dt_u(i) * 0.5 * ( tsgs2(i+1,j,k)-tsgs2(i,j,k) )
        endif
        if (find_work) drdi_u(i,k) = drdib

        if (k > nk_linear) then
          if (use_eos) then
            if (cs%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x) then
              hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
              hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
              hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
              har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
              if (gv%Boussinesq) then
                dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
              else
                dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
                dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
              endif
              ! Use the harmonic mean thicknesses to weight the horizontal gradients.
              ! These unnormalized weights have been rearranged to minimize divisions.
              wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)

              drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
              ! The expression for drdz above is mathematically equivalent to:
              !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
              !          ((hg2L/haL) + (hg2R/haR))
              hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
              hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
              haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
              hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect

              ! hN2_u is used with the FGNV streamfunction formulation
              hn2_u(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
                           max(drdz*g_rho0, n2_floor)
            endif
            if (present_slope_x) then
              slope = slope_x(i,j,k)
              slope2_ratio_u(i,k) = slope**2 * i_slope_max2
            else
              ! Use the harmonic mean thicknesses to weight the horizontal gradients.
              ! These unnormalized weights have been rearranged to minimize divisions.
              wta = hg2a*hab ; wtb = hg2b*haa
              ! This is the gradient of density along geopotentials.
              drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
                      drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)

              ! This estimate of slope is accurate for small slopes, but bounded
              ! to be between -1 and 1.
              mag_grad2 = drdx**2 + (us%L_to_Z*drdz)**2
              if (mag_grad2 > 0.0) then
                slope = drdx / sqrt(mag_grad2)
                slope2_ratio_u(i,k) = slope**2 * i_slope_max2
              else ! Just in case mag_grad2 = 0 ever.
                slope = 0.0
                slope2_ratio_u(i,k) = 1.0e20  ! Force the use of the safe streamfunction.
              endif
            endif

            ! Adjust real slope by weights that bias towards slope of interfaces
            ! that ignore density gradients along layers.
            if (present_int_slope_u) then
              slope = (1.0 - int_slope_u(i,j,k)) * slope + &
                      int_slope_u(i,j,k) * us%Z_to_L*((e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j))
              slope2_ratio_u(i,k) = (1.0 - int_slope_u(i,j,k)) * slope2_ratio_u(i,k)
            endif

            slope_x_pe(i,j,k) = min(slope,cs%slope_max)
            hn2_x_pe(i,j,k) = hn2_u(i,k)
            if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope

            ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].
            sfn_unlim_u(i,k) = -((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)

            ! Avoid moving dense water upslope from below the level of
            ! the bottom on the receiving side.
            if (sfn_unlim_u(i,k) > 0.0) then ! The flow below this interface is positive.
              if (e(i,j,k) < e(i+1,j,nz+1)) then
                sfn_unlim_u(i,k) = 0.0 ! This is not uhtot, because it may compensate for
                                ! deeper flow in very unusual cases.
              elseif (e(i+1,j,nz+1) > e(i,j,k+1)) then
                ! Scale the transport with the fraction of the donor layer above
                ! the bottom on the receiving side.
                sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i,j,k) - e(i+1,j,nz+1)) / &
                                         ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
              endif
            else
              if (e(i+1,j,k) < e(i,j,nz+1)) then ; sfn_unlim_u(i,k) = 0.0
              elseif (e(i,j,nz+1) > e(i+1,j,k+1)) then
                sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i+1,j,k) - e(i,j,nz+1)) / &
                                       ((e(i+1,j,k) - e(i+1,j,k+1)) + dz_neglect))
              endif
            endif

          else ! .not. use_EOS
            if (present_slope_x) then
              slope = slope_x(i,j,k)
            else
              slope = us%Z_to_L*((e(i,j,k)-e(i+1,j,k))*g%IdxCu(i,j)) * g%mask2dCu(i,j)
            endif
            if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
            sfn_unlim_u(i,k) = ((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)
            hn2_u(i,k) = gv%g_prime(k)
          endif ! if (use_EOS)
        else ! if (k > nk_linear)
          hn2_u(i,k) = n2_floor * dz_neglect
          sfn_unlim_u(i,k) = 0.
        endif ! if (k > nk_linear)
        if (cs%id_sfn_unlim_x>0) diag_sfn_unlim_x(i,j,k) = sfn_unlim_u(i,k)
      enddo ! i-loop
    enddo ! k-loop

    if (cs%use_FGNV_streamfn) then
      do k=1,nz ; do i=is-1,ie ; if (g%mask2dCu(i,j)>0.) then
        h_harm = max( h_neglect, &
              2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) )
        c2_h_u(i,k) = cs%FGNV_scale * &
            ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (gv%H_to_Z*h_harm)
      endif ; enddo ; enddo

      ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
      do i=is-1,ie
        if (g%mask2dCu(i,j)>0.) then
          sfn_unlim_u(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_u(i,:)
          call streamfn_solver(nz, c2_h_u(i,:), hn2_u(i,:), sfn_unlim_u(i,:))
        else
          sfn_unlim_u(i,:) = 0.
        endif
      enddo
    endif

    do k=nz,2,-1
      do i=is-1,ie
        if (k > nk_linear) then
          if (use_eos) then

            if (uhtot(i,j) <= 0.0) then
              ! The transport that must balance the transport below is positive.
              sfn_safe = uhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
            else !  (uhtot(I,j) > 0.0)
              sfn_safe = uhtot(i,j) * (1.0 - h_frac(i+1,j,k)) * gv%H_to_Z
            endif

            ! The actual streamfunction at each interface.
            sfn_est = (sfn_unlim_u(i,k) + slope2_ratio_u(i,k)*sfn_safe) / (1.0 + slope2_ratio_u(i,k))
          else  ! With .not.use_EOS, the layers are constant density.
            sfn_est = sfn_unlim_u(i,k)
          endif

          ! Make sure that there is enough mass above to allow the streamfunction
          ! to satisfy the boundary condition of 0 at the surface.
          sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i+1,j,k))

          ! The actual transport is limited by the mass available in the two
          ! neighboring grid cells.
          uhd(i,j,k) = max(min((sfn_in_h - uhtot(i,j)), h_avail(i,j,k)), &
                           -h_avail(i+1,j,k))

          if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
!         sfn_x(I,j,K) = max(min(Sfn_in_h, uhtot(I,j)+h_avail(i,j,k)), &
!                            uhtot(I,j)-h_avail(i+1,j,K))
!         sfn_slope_x(I,j,K) = max(uhtot(I,j)-h_avail(i+1,j,k), &
!                                  min(uhtot(I,j)+h_avail(i,j,k), &
!               min(h_avail_rsum(i+1,j,K), max(-h_avail_rsum(i,j,K), &
!               (KH_u(I,j,K)*G%dy_Cu(I,j)) * ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) )) ))
        else ! k <= nk_linear
          ! Balance the deeper flow with a return flow uniformly distributed
          ! though the remaining near-surface layers.  This is the same as
          ! using Sfn_safe above.  There is no need to apply the limiters in
          ! this case.
          if (uhtot(i,j) <= 0.0) then
            uhd(i,j,k) = -uhtot(i,j) * h_frac(i,j,k)
          else !  (uhtot(I,j) > 0.0)
            uhd(i,j,k) = -uhtot(i,j) * h_frac(i+1,j,k)
          endif

          if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
!         if (sfn_slope_x(I,j,K+1) <= 0.0) then
!           sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i,j,k))
!         else
!           sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i+1,j,k))
!         endif
        endif

        uhtot(i,j) = uhtot(i,j) + uhd(i,j,k)

        if (find_work) then
          !   This is the energy tendency based on the original profiles, and does
          ! not include any nonlinear terms due to a finite time step (which would
          ! involve interactions between the fluxes through the different faces.
          !   A second order centered estimate is used for the density transfered
          ! between water columns.

          work_u(i,j) = work_u(i,j) + g_scale * &
            ( uhtot(i,j) * drdkde_u(i,k) - &
              (uhd(i,j,k) * drdi_u(i,k)) * 0.25 * &
              ((e(i,j,k) + e(i,j,k+1)) + (e(i+1,j,k) + e(i+1,j,k+1))) )
        endif

      enddo
    enddo ! end of k-loop
  enddo ! end of j-loop

    ! Calculate the meridional fluxes and gradients.
    eosdom_v(:) = eos_domain(g%HI)
!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
!$OMP                                  nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,  &
!$OMP                                  I_slope_max2,h_neglect2,present_int_slope_v, &
!$OMP                                  int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum,  &
!$OMP                                  vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1,   &
!$OMP                                  diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,&
!$OMP                                  use_stanley, Tsgs2,                          &
!$OMP                                  present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE)  &
!$OMP                          private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v,      &
!$OMP                                  drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP                                  drho_dT_dT_v,scrap,                          &
!$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &
!$OMP                                  drdy,mag_grad2,Slope,slope2_Ratio_v,hN2_v,   &
!$OMP                                  Sfn_unlim_v,drdj_v,drdkDe_v,h_harm,c2_h_v,   &
!$OMP                                  Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
  do j=js-1,je
    do k=nz,2,-1
      if (find_work .and. .not.(use_eos)) then
        drdja = 0.0 ; drdjb = 0.0
        drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
      endif

      calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
         (find_work .or. .not. present_slope_y .or. cs%use_FGNV_streamfn .or. use_stanley)

      if (calc_derivatives) then
        do i=is,ie
          pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
          t_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
          s_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
        enddo
        call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, &
                                      tv%eqn_of_state, eosdom_v)
      endif
      if (use_stanley) then
        ! The second line below would correspond to arguments
        !            drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
        call calculate_density_second_derivs(t_v, s_v, pres_v, &
                     scrap, scrap, drho_dt_dt_v, scrap, scrap, &
                     is, ie-is+1, tv%eqn_of_state)
      endif
      do i=is,ie
        if (calc_derivatives) then
          ! Estimate the horizontal density gradients along layers.
          drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
                  drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
          drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
                  drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))

          ! Estimate the vertical density gradients times the grid spacing.
          drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
                   drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
          drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
                   drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
          drdkde_v(i,k) =  drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
        elseif (find_work) then ! This is used in pure stacked SW mode
          drdkde_v(i,k) =  drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
        endif
        if (use_stanley) then
          ! Correction to the horizontal density gradient due to nonlinearity in
          ! the EOS rectifying SGS temperature anomalies
          drdja = drdja + drho_dt_dt_v(i) * 0.5 * ( tsgs2(i,j+1,k-1)-tsgs2(i,j,k-1) )
          drdjb = drdjb + drho_dt_dt_v(i) * 0.5 * ( tsgs2(i,j+1,k)-tsgs2(i,j,k) )
        endif

        if (find_work) drdj_v(i,k) = drdjb

        if (k > nk_linear) then
          if (use_eos) then
            if (cs%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y) then
              hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
              hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
              hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
              har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
              if (gv%Boussinesq) then
                dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
              else
                dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
                dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
              endif
              ! Use the harmonic mean thicknesses to weight the horizontal gradients.
              ! These unnormalized weights have been rearranged to minimize divisions.
              wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)

              drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
              ! The expression for drdz above is mathematically equivalent to:
              !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
              !          ((hg2L/haL) + (hg2R/haR))
              hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
              hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
              haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
              hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect

              ! hN2_v is used with the FGNV streamfunction formulation
              hn2_v(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
                           max(drdz*g_rho0, n2_floor)
            endif
            if (present_slope_y) then
              slope = slope_y(i,j,k)
              slope2_ratio_v(i,k) = slope**2 * i_slope_max2
            else
              ! Use the harmonic mean thicknesses to weight the horizontal gradients.
              ! These unnormalized weights have been rearranged to minimize divisions.
              wta = hg2a*hab ; wtb = hg2b*haa
              ! This is the gradient of density along geopotentials.
              drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
                      drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)

              ! This estimate of slope is accurate for small slopes, but bounded
              ! to be between -1 and 1.
              mag_grad2 = drdy**2 + (us%L_to_Z*drdz)**2
              if (mag_grad2 > 0.0) then
                slope = drdy / sqrt(mag_grad2)
                slope2_ratio_v(i,k) = slope**2 * i_slope_max2
              else ! Just in case mag_grad2 = 0 ever.
                slope = 0.0
                slope2_ratio_v(i,k) = 1.0e20  ! Force the use of the safe streamfunction.
              endif
            endif

            ! Adjust real slope by weights that bias towards slope of interfaces
            ! that ignore density gradients along layers.
            if (present_int_slope_v) then
              slope = (1.0 - int_slope_v(i,j,k)) * slope + &
                      int_slope_v(i,j,k) * us%Z_to_L*((e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j))
              slope2_ratio_v(i,k) = (1.0 - int_slope_v(i,j,k)) * slope2_ratio_v(i,k)
            endif

            slope_y_pe(i,j,k) = min(slope,cs%slope_max)
            hn2_y_pe(i,j,k) = hn2_v(i,k)
            if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope

            ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].
            sfn_unlim_v(i,k) = -((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)

            ! Avoid moving dense water upslope from below the level of
            ! the bottom on the receiving side.
            if (sfn_unlim_v(i,k) > 0.0) then ! The flow below this interface is positive.
              if (e(i,j,k) < e(i,j+1,nz+1)) then
                sfn_unlim_v(i,k) = 0.0 ! This is not vhtot, because it may compensate for
                                ! deeper flow in very unusual cases.
              elseif (e(i,j+1,nz+1) > e(i,j,k+1)) then
                ! Scale the transport with the fraction of the donor layer above
                ! the bottom on the receiving side.
                sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j,k) - e(i,j+1,nz+1)) / &
                                         ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
              endif
            else
              if (e(i,j+1,k) < e(i,j,nz+1)) then ; sfn_unlim_v(i,k) = 0.0
              elseif (e(i,j,nz+1) > e(i,j+1,k+1)) then
                sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j+1,k) - e(i,j,nz+1)) / &
                                       ((e(i,j+1,k) - e(i,j+1,k+1)) + dz_neglect))
              endif
            endif

          else ! .not. use_EOS
            if (present_slope_y) then
              slope = slope_y(i,j,k)
            else
              slope = us%Z_to_L*((e(i,j,k)-e(i,j+1,k))*g%IdyCv(i,j)) * g%mask2dCv(i,j)
            endif
            if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
            sfn_unlim_v(i,k) = ((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)
            hn2_v(i,k) = gv%g_prime(k)
          endif ! if (use_EOS)
        else ! if (k > nk_linear)
          hn2_v(i,k) = n2_floor * dz_neglect
          sfn_unlim_v(i,k) = 0.
        endif ! if (k > nk_linear)
        if (cs%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,j,k) = sfn_unlim_v(i,k)
      enddo ! i-loop
    enddo ! k-loop

    if (cs%use_FGNV_streamfn) then
      do k=1,nz ; do i=is,ie ; if (g%mask2dCv(i,j)>0.) then
        h_harm = max( h_neglect, &
              2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) )
        c2_h_v(i,k) = cs%FGNV_scale * &
            ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (gv%H_to_Z*h_harm)
      endif ; enddo ; enddo

      ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
      do i=is,ie
        if (g%mask2dCv(i,j)>0.) then
          sfn_unlim_v(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_v(i,:)
          call streamfn_solver(nz, c2_h_v(i,:), hn2_v(i,:), sfn_unlim_v(i,:))
        else
          sfn_unlim_v(i,:) = 0.
        endif
      enddo
    endif

    do k=nz,2,-1
      do i=is,ie
        if (k > nk_linear) then
          if (use_eos) then

            if (vhtot(i,j) <= 0.0) then
              ! The transport that must balance the transport below is positive.
              sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
            else !  (vhtot(I,j) > 0.0)
              sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j+1,k)) * gv%H_to_Z
            endif

            ! The actual streamfunction at each interface.
            sfn_est = (sfn_unlim_v(i,k) + slope2_ratio_v(i,k)*sfn_safe) / (1.0 + slope2_ratio_v(i,k))
          else      ! With .not.use_EOS, the layers are constant density.
            sfn_est = sfn_unlim_v(i,k)
          endif

          ! Make sure that there is enough mass above to allow the streamfunction
          ! to satisfy the boundary condition of 0 at the surface.
          sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i,j+1,k))

          ! The actual transport is limited by the mass available in the two
          ! neighboring grid cells.
          vhd(i,j,k) = max(min((sfn_in_h - vhtot(i,j)), h_avail(i,j,k)), -h_avail(i,j+1,k))

          if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
!         sfn_y(i,J,K) = max(min(Sfn_in_h, vhtot(i,J)+h_avail(i,j,k)), &
!                            vhtot(i,J)-h_avail(i,j+1,k))
!         sfn_slope_y(i,J,K) = max(vhtot(i,J)-h_avail(i,j+1,k), &
!                                  min(vhtot(i,J)+h_avail(i,j,k), &
!               min(h_avail_rsum(i,j+1,K), max(-h_avail_rsum(i,j,K), &
!               (KH_v(i,J,K)*G%dx_Cv(i,J)) * ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) )) ))
        else ! k <= nk_linear
          ! Balance the deeper flow with a return flow uniformly distributed
          ! though the remaining near-surface layers.  This is the same as
          ! using Sfn_safe above.  There is no need to apply the limiters in
          ! this case.
          if (vhtot(i,j) <= 0.0) then
            vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j,k)
          else !  (vhtot(i,J) > 0.0)
            vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j+1,k)
          endif

          if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
!         if (sfn_slope_y(i,J,K+1) <= 0.0) then
!           sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j,k))
!         else
!           sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j+1,k))
!         endif
        endif

        vhtot(i,j) = vhtot(i,j)  + vhd(i,j,k)

        if (find_work) then
          !   This is the energy tendency based on the original profiles, and does
          ! not include any nonlinear terms due to a finite time step (which would
          ! involve interactions between the fluxes through the different faces.
          !   A second order centered estimate is used for the density transfered
          ! between water columns.

          work_v(i,j) = work_v(i,j) + g_scale * &
            ( vhtot(i,j) * drdkde_v(i,k) - &
             (vhd(i,j,k) * drdj_v(i,k)) * 0.25 * &
             ((e(i,j,k) + e(i,j,k+1)) + (e(i,j+1,k) + e(i,j+1,k+1))) )
        endif

      enddo
    enddo ! end of k-loop
  enddo ! end of j-loop

  ! In layer 1, enforce the boundary conditions that Sfn(z=0) = 0.0
  if (.not.find_work .or. .not.(use_eos)) then
    do j=js,je ; do i=is-1,ie ; uhd(i,j,1) = -uhtot(i,j) ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; vhd(i,j,1) = -vhtot(i,j) ; enddo ; enddo
  else
    eosdom_u(1) = (is-1) - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
    !$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB)
    do j=js,je
      if (use_eos) then
        do i=is-1,ie
          pres_u(i) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
          t_u(i) = 0.5*(t(i,j,1) + t(i+1,j,1))
          s_u(i) = 0.5*(s(i,j,1) + s(i+1,j,1))
        enddo
        call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
                                      tv%eqn_of_state, eosdom_u )
      endif
      do i=is-1,ie
        uhd(i,j,1) = -uhtot(i,j)

        if (use_eos) then
          drdib = drho_dt_u(i) * (t(i+1,j,1)-t(i,j,1)) + &
                  drho_ds_u(i) * (s(i+1,j,1)-s(i,j,1))
        endif
        if (cs%use_GM_work_bug) then
          work_u(i,j) = work_u(i,j) + g_scale * &
              ( (uhd(i,j,1) * drdib) * 0.25 * &
                ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
        else
          work_u(i,j) = work_u(i,j) - g_scale * &
              ( (uhd(i,j,1) * drdib) * 0.25 * &
                ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
        endif
      enddo
    enddo

    eosdom_v(:) = eos_domain(g%HI)
    !$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB)
    do j=js-1,je
      if (use_eos) then
        do i=is,ie
          pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
          t_v(i) = 0.5*(t(i,j,1) + t(i,j+1,1))
          s_v(i) = 0.5*(s(i,j,1) + s(i,j+1,1))
        enddo
        call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, &
                                      tv%eqn_of_state, eosdom_v)
      endif
      do i=is,ie
        vhd(i,j,1) = -vhtot(i,j)

        if (use_eos) then
          drdjb = drho_dt_v(i) * (t(i,j+1,1)-t(i,j,1)) + &
                  drho_ds_v(i) * (s(i,j+1,1)-s(i,j,1))
        endif
        work_v(i,j) = work_v(i,j) - g_scale * &
            ( (vhd(i,j,1) * drdjb) * 0.25 * &
              ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
      enddo
    enddo
  endif

  if (find_work) then ; do j=js,je ; do i=is,ie
    ! Note that the units of Work_v and Work_u are W, while Work_h is W m-2.
    work_h = 0.5 * g%IareaT(i,j) * &
      ((work_u(i-1,j) + work_u(i,j)) + (work_v(i,j-1) + work_v(i,j)))
    if (associated(cs%GMwork)) cs%GMwork(i,j) = work_h
    if (associated(meke) .and. .not.cs%GM_src_alt) then ; if (associated(meke%GM_src)) then
      meke%GM_src(i,j) = meke%GM_src(i,j) + work_h
    endif ; endif
  enddo ; enddo ; endif

  if (find_work .and. cs%GM_src_alt .and. associated(meke)) then ; if (associated(meke%GM_src)) then
    do j=js,je ; do i=is,ie ; do k=nz,1,-1
      pe_release_h = -0.25*(kh_u(i,j,k)*(slope_x_pe(i,j,k)**2) * hn2_x_pe(i,j,k) + &
                            kh_u(i-1,j,k)*(slope_x_pe(i-1,j,k)**2) * hn2_x_pe(i-1,j,k) + &
                            kh_v(i,j,k)*(slope_y_pe(i,j,k)**2) * hn2_y_pe(i,j,k) + &
                            kh_v(i,j-1,k)*(slope_y_pe(i,j-1,k)**2) * hn2_y_pe(i,j-1,k))
      meke%GM_src(i,j) = meke%GM_src(i,j) + us%L_to_Z**2 * gv%Rho0 * pe_release_h
    enddo ; enddo ; enddo
  endif ; endif

  if (cs%id_slope_x > 0) call post_data(cs%id_slope_x, cs%diagSlopeX, cs%diag)
  if (cs%id_slope_y > 0) call post_data(cs%id_slope_y, cs%diagSlopeY, cs%diag)
  if (cs%id_sfn_x > 0) call post_data(cs%id_sfn_x, diag_sfn_x, cs%diag)
  if (cs%id_sfn_y > 0) call post_data(cs%id_sfn_y, diag_sfn_y, cs%diag)
  if (cs%id_sfn_unlim_x > 0) call post_data(cs%id_sfn_unlim_x, diag_sfn_unlim_x, cs%diag)
  if (cs%id_sfn_unlim_y > 0) call post_data(cs%id_sfn_unlim_y, diag_sfn_unlim_y, cs%diag)

end subroutine thickness_diffuse_full

!> Tridiagonal solver for streamfunction at interfaces
subroutine streamfn_solver(nk, c2_h, hN2, sfn)
  integer,               intent(in)    :: nk   !< Number of layers
  real, dimension(nk),   intent(in)    :: c2_h !< Wave speed squared over thickness in layers [L2 Z-1 T-2 ~> m s-2]
  real, dimension(nk+1), intent(in)    :: hN2  !< Thickness times N2 at interfaces [L2 Z-1 T-2 ~> m s-2]
  real, dimension(nk+1), intent(inout) :: sfn  !< Streamfunction [Z L2 T-1 ~> m3 s-1] or arbitrary units
                                               !! On entry, equals diffusivity times slope.
                                               !! On exit, equals the streamfunction.
  ! Local variables
  integer :: k

  real :: b_denom, beta, d1, c1(nk)

  sfn(1) = 0.
  b_denom = hn2(2) + c2_h(1)
  beta = 1.0 / ( b_denom + c2_h(2) )
  d1 = beta * b_denom
  sfn(2) = ( beta * hn2(2) )*sfn(2)
  do k=3,nk
    c1(k-1) = beta * c2_h(k-1)
    b_denom = hn2(k) + d1*c2_h(k-1)
    beta = 1.0 / (b_denom + c2_h(k))
    d1 = beta * b_denom
    sfn(k) = beta * (hn2(k)*sfn(k) + c2_h(k-1)*sfn(k-1))
  enddo
  c1(nk) = beta * c2_h(nk)
  sfn(nk+1) = 0.
  do k=nk,2,-1
    sfn(k) = sfn(k) + c1(k)*sfn(k+1)
  enddo

end subroutine streamfn_solver

!> Modifies thickness diffusivities to untangle layer structures
subroutine add_detangling_kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, CS, &
                             int_slope_u, int_slope_v)
  type(ocean_grid_type),                       intent(in)    :: G    !< Ocean grid structure
  type(verticalgrid_type),                     intent(in)    :: GV   !< Vertical grid structure
  type(unit_scale_type),                       intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)    :: h    !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1),  intent(in)    :: e    !< Interface positions [Z ~> m]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kh_u !< Thickness diffusivity on interfaces
                                                                     !! at u points [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: Kh_v !< Thickness diffusivity on interfaces
                                                                     !! at v points [L2 T-1 ~> m2 s-1]
  real, dimension(SZIB_(G),SZJ_(G)),           intent(in)    :: Kh_u_CFL !< Maximum stable thickness diffusivity
                                                                     !! at u points [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJB_(G)),           intent(in)    :: Kh_v_CFL !< Maximum stable thickness diffusivity
                                                                     !! at v points [L2 T-1 ~> m2 s-1]
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< Thermodynamics structure
  real,                                        intent(in)    :: dt   !< Time increment [T ~> s]
  type(thickness_diffuse_cs),                  pointer       :: CS   !< Control structure for thickness diffusion
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: int_slope_u !< Ratio that determine how much of
                                                                     !! the isopycnal slopes are taken directly from
                                                                     !! the interface slopes without consideration
                                                                     !! of density gradients.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: int_slope_v !< Ratio that determine how much of
                                                                     !! the isopycnal slopes are taken directly from
                                                                     !! the interface slopes without consideration
                                                                     !! of density gradients.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    de_top     ! The distances between the top of a layer and the top of the
               ! region where the detangling is applied [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
    Kh_lay_u   ! The tentative interface height diffusivity for each layer at
               ! u points [L2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
    Kh_lay_v   ! The tentative interface height diffusivity for each layer at
               ! v points [L2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    de_bot     ! The distances from the bottom of the region where the
               ! detangling is applied [H ~> m or kg m-2].
  real :: h1, h2    ! The thinner and thicker surrounding thicknesses [H ~> m or kg m-2],
                    ! with the thinner modified near the boundaries to mask out
                    ! thickness variations due to topography, etc.
  real :: jag_Rat   ! The nondimensional jaggedness ratio for a layer, going
                    ! from 0 (smooth) to 1 (jagged).  This is the difference
                    ! between the arithmetic and harmonic mean thicknesses
                    ! normalized by the arithmetic mean thickness.
  real :: Kh_scale  ! A ratio by which Kh_u_CFL is scaled for maximally jagged
                    ! layers [nondim].
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected [H ~> m or kg m-2].

  real :: I_sl      ! The absolute value of the larger in magnitude of the slopes
                    ! above and below [L Z-1 ~> nondim].
  real :: Rsl       ! The ratio of the smaller magnitude slope to the larger
                    ! magnitude one [nondim]. 0 <= Rsl <1.
  real :: IRsl      ! The (limited) inverse of Rsl [nondim]. 1 < IRsl <= 1e9.
  real :: dH        ! The thickness gradient divided by the damping timescale
                    ! and the ratio of the face length to the adjacent cell
                    ! areas for comparability with the diffusivities [L Z T-1 ~> m2 s-1].
  real :: adH       ! The absolute value of dH [L Z T-1 ~> m2 s-1].
  real :: sign      ! 1 or -1, with the same sign as the layer thickness gradient.
  real :: sl_K      ! The sign-corrected slope of the interface above [Z L-1 ~> nondim].
  real :: sl_Kp1    ! The sign-corrected slope of the interface below [Z L-1 ~> nondim].
  real :: I_sl_K    ! The (limited) inverse of sl_K [L Z-1 ~> nondim].
  real :: I_sl_Kp1  ! The (limited) inverse of sl_Kp1 [L Z-1 ~> nondim].
  real :: I_4t      ! A quarter of a flux scaling factor divided by
                    ! the damping timescale [T-1 ~> s-1].
  real :: Fn_R      ! A function of Rsl, such that Rsl < Fn_R < 1.
  real :: denom, I_denom ! A denominator and its inverse, various units.
  real :: Kh_max    ! A local ceiling on the diffusivity [L2 T-1 ~> m2 s-1].
  real :: wt1, wt2  ! Nondimensional weights.
  !   Variables used only in testing code.
  ! real, dimension(SZK_(G)) :: uh_here
  ! real, dimension(SZK_(G)+1) :: Sfn
  real :: dKh       ! An increment in the diffusivity [L2 T-1 ~> m2 s-1].

  real, dimension(SZIB_(G),SZK_(G)+1) :: &
    Kh_bg, &        ! The background (floor) value of Kh [L2 T-1 ~> m2 s-1].
    Kh, &           ! The tentative value of Kh [L2 T-1 ~> m2 s-1].
    Kh_detangle, &  ! The detangling diffusivity that could be used [L2 T-1 ~> m2 s-1].
    Kh_min_max_p, & ! The smallest ceiling that can be placed on Kh(I,K)
                    ! based on the value of Kh(I,K+1) [L2 T-1 ~> m2 s-1].
    kh_min_max_m, & ! The smallest ceiling that can be placed on Kh(I,K)
                    ! based on the value of Kh(I,K-1) [L2 T-1 ~> m2 s-1].
    ! The following are variables that define the relationships between
    ! successive values of Kh.
    ! Search for Kh that satisfy...
    !    Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
    !    Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
    !    Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
    !    Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
    kh_min_m , &   ! See above [nondim].
    kh0_min_m , &  ! See above [L2 T-1 ~> m2 s-1].
    kh_max_m , &   ! See above [nondim].
    kh0_max_m, &   ! See above [L2 T-1 ~> m2 s-1].
    kh_min_p , &   ! See above [nondim].
    kh0_min_p , &  ! See above [L2 T-1 ~> m2 s-1].
    kh_max_p , &   ! See above [nondim].
    kh0_max_p      ! See above [L2 T-1 ~> m2 s-1].
  real, dimension(SZIB_(G)) :: &
    Kh_max_max  ! The maximum diffusivity permitted in a column [L2 T-1 ~> m2 s-1]..
  logical, dimension(SZIB_(G)) :: &
    do_i        ! If true, work on a column.
  integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  k_top = gv%nk_rho_varies + 1
  h_neglect = gv%H_subroundoff
  !   The 0.5 is because we are not using uniform weightings, but are
  ! distributing the diffusivities more effectively (with wt1 & wt2), but this
  ! means that the additions to a single interface can be up to twice as large.
  kh_scale = 0.5
  if (cs%detangle_time > dt) kh_scale = 0.5 * dt / cs%detangle_time

  do j=js-1,je+1 ; do i=is-1,ie+1
    de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
  enddo ; enddo
  do k=k_top+1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
    de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
  enddo ; enddo ; enddo

  do j=js,je ; do i=is-1,ie
    kh_lay_u(i,j,nz) = 0.0 ; kh_lay_u(i,j,k_top) = 0.0
  enddo ; enddo
  do j=js-1,je ; do i=is,ie
    kh_lay_v(i,j,nz) = 0.0 ; kh_lay_v(i,j,k_top) = 0.0
  enddo ; enddo

  do k=nz-1,k_top+1,-1
    ! Find the diffusivities associated with each layer.
    do j=js-1,je+1 ; do i=is-1,ie+1
      de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
    enddo ; enddo

    do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.0) then
      if (h(i,j,k) > h(i+1,j,k)) then
        h2 = h(i,j,k)
        h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
      else
        h2 = h(i+1,j,k)
        h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
      endif
      jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
      kh_lay_u(i,j,k) = (kh_scale * kh_u_cfl(i,j)) * jag_rat**2
    endif ; enddo ; enddo

    do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.0) then
      if (h(i,j,k) > h(i,j+1,k)) then
        h2 = h(i,j,k)
        h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
      else
        h2 = h(i,j+1,k)
        h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
      endif
      jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
      kh_lay_v(i,j,k) = (kh_scale * kh_v_cfl(i,j)) * jag_rat**2
    endif ; enddo ; enddo
  enddo

  ! Limit the diffusivities

  i_4t = kh_scale / (4.0 * dt)

  do n=1,2
    if (n==1) then ; jsh = js ; ish = is-1
    else ; jsh = js-1 ; ish = is ; endif

    do j=jsh,je

      ! First, populate the diffusivities
      if (n==1) then ! This is a u-column.
        do i=ish,ie
          do_i(i) = (g%mask2dCu(i,j) > 0.0)
          kh_max_max(i) = kh_u_cfl(i,j)
        enddo
        do k=1,nz+1 ; do i=ish,ie
          kh_bg(i,k) = kh_u(i,j,k) ; kh(i,k) = kh_bg(i,k)
          kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
          kh_detangle(i,k) = 0.0
        enddo ; enddo
      else ! This is a v-column.
        do i=ish,ie
          do_i(i) = (g%mask2dCv(i,j) > 0.0) ; kh_max_max(i) = kh_v_cfl(i,j)
        enddo
        do k=1,nz+1 ; do i=ish,ie
          kh_bg(i,k) = kh_v(i,j,k) ; kh(i,k) = kh_bg(i,k)
          kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
          kh_detangle(i,k) = 0.0
        enddo ; enddo
      endif

      ! Determine the limits on the diffusivities.
      do k=k_top,nz ; do i=ish,ie ; if (do_i(i)) then
        if (n==1) then ! This is a u-column.
          dh = 0.0
          denom = ((g%IareaT(i+1,j) + g%IareaT(i,j)) * g%dy_Cu(i,j))
          !   This expression uses differences in e in place of h for better
          ! consistency with the slopes.
          if (denom > 0.0) &
            dh = i_4t * ((e(i+1,j,k) - e(i+1,j,k+1)) - &
                         (e(i,j,k) - e(i,j,k+1))) / denom
           ! dH = I_4t * (h(i+1,j,k) - h(i,j,k)) / denom

          adh = abs(dh)
          sign = 1.0 ; if (dh < 0) sign = -1.0
          sl_k = sign * (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
          sl_kp1 = sign * (e(i+1,j,k+1)-e(i,j,k+1)) * g%IdxCu(i,j)

          ! Add the incremental diffusivites to the surrounding interfaces.
          ! Adding more to the more steeply sloping layers (as below) makes
          ! the diffusivities more than twice as effective.
          denom = (sl_k**2 + sl_kp1**2)
          wt1 = 0.5 ; wt2 = 0.5
          if (denom > 0.0) then
            wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
          endif
          kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_u(i,j,k)
          kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_u(i,j,k)
        else ! This is a v-column.
          dh = 0.0
          denom = ((g%IareaT(i,j+1) + g%IareaT(i,j)) * g%dx_Cv(i,j))
          if (denom > 0.0) &
            dh = i_4t * ((e(i,j+1,k) - e(i,j+1,k+1)) - &
                         (e(i,j,k) - e(i,j,k+1))) / denom
           ! dH = I_4t * (h(i,j+1,k) - h(i,j,k)) / denom

          adh = abs(dh)
          sign = 1.0 ; if (dh < 0) sign = -1.0
          sl_k = sign * (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
          sl_kp1 = sign * (e(i,j+1,k+1)-e(i,j,k+1)) * g%IdyCv(i,j)

          ! Add the incremental diffusviites to the surrounding interfaces.
          ! Adding more to the more steeply sloping layers (as below) makes
          ! the diffusivities more than twice as effective.
          denom = (sl_k**2 + sl_kp1**2)
          wt1 = 0.5 ; wt2 = 0.5
          if (denom > 0.0) then
            wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
          endif
          kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_v(i,j,k)
          kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_v(i,j,k)
        endif

        if (adh == 0.0) then
          kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
          kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
          kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
          kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
        elseif (adh > 0.0) then
          if (sl_k <= sl_kp1) then
            ! This case should only arise from nonlinearities in the equation of state.
            ! Treat it as though dedx(K) = dedx(K+1) & dH = 0.
            kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
            kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
            kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
            kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
          elseif (sl_k <= 0.0) then   ! Both slopes are opposite to dH
            i_sl = -1.0 / sl_kp1
            rsl = -sl_k * i_sl                            ! 0 <= Rsl < 1
            irsl = 1e9 ; if (rsl > 1e-9) irsl = 1.0/rsl   ! 1 < IRsl <= 1e9

            fn_r = rsl
            if (kh_max_max(i) > 0) &
              fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / (kh_max_max(i)))

            kh_min_m(i,k+1) = fn_r ; kh0_min_m(i,k+1) = 0.0
            kh_max_m(i,k+1) = rsl ; kh0_max_m(i,k+1) = adh * i_sl
            kh_min_p(i,k) = irsl ; kh0_min_p(i,k) = -adh * (i_sl*irsl)
            kh_max_p(i,k) = 1.0/(fn_r + 1.0e-30) ; kh0_max_p(i,k) = 0.0
          elseif (sl_kp1 < 0.0) then  ! Opposite (nonzero) signs of slopes.
            i_sl_k = 1e18*us%Z_to_L ; if (sl_k > 1e-18*us%L_to_Z) i_sl_k = 1.0 / sl_k
            i_sl_kp1 = 1e18*us%Z_to_L ; if (-sl_kp1 > 1e-18*us%L_to_Z) i_sl_kp1 = -1.0 / sl_kp1

            kh_min_m(i,k+1) = 0.0 ; kh0_min_m(i,k+1) = 0.0
            kh_max_m(i,k+1) = - sl_k*i_sl_kp1 ; kh0_max_m(i,k+1) = adh*i_sl_kp1
            kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
            kh_max_p(i,k) = sl_kp1*i_sl_k ; kh0_max_p(i,k) = adh*i_sl_k

            ! This limit does not use the slope weighting so that potentially
            ! sharp gradients in diffusivities are not forced to occur.
            kh_max = adh / (sl_k - sl_kp1)
            kh_min_max_p(i,k) = max(kh_min_max_p(i,k), kh_max)
            kh_min_max_m(i,k+1) = max(kh_min_max_m(i,k+1), kh_max)
          else ! Both slopes are of the same sign as dH.
            i_sl = 1.0 / sl_k
            rsl = sl_kp1 * i_sl                           ! 0 <= Rsl < 1
            irsl = 1e9 ; if (rsl > 1e-9) irsl = 1.0/rsl   ! 1 < IRsl <= 1e9

            ! Rsl <= Fn_R <= 1
            fn_r = rsl
            if (kh_max_max(i) > 0) &
              fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / kh_max_max(i))

            kh_min_m(i,k+1) = irsl ; kh0_min_m(i,k+1) = -adh * (i_sl*irsl)
            kh_max_m(i,k+1) = 1.0/(fn_r + 1.0e-30) ; kh0_max_m(i,k+1) = 0.0
            kh_min_p(i,k) = fn_r ; kh0_min_p(i,k) = 0.0
            kh_max_p(i,k) = rsl ; kh0_max_p(i,k) = adh * i_sl
          endif
        endif
      endif ; enddo ; enddo ! I-loop & k-loop

      do k=k_top,nz+1,nz+1-k_top ; do i=ish,ie ; if (do_i(i)) then
        ! The diffusivities at k_top and nz+1 are both fixed.
        kh_min_m(i,k) = 0.0 ; kh0_min_m(i,k) = 0.0
        kh_max_m(i,k) = 0.0 ; kh0_max_m(i,k) = 0.0
        kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
        kh_max_p(i,k) = 0.0 ; kh0_max_p(i,k) = 0.0
        kh_min_max_p(i,k) = kh_bg(i,k)
        kh_min_max_m(i,k) = kh_bg(i,k)
      endif ; enddo ; enddo ! I-loop and k_top/nz+1 loop

      ! Search for Kh that satisfy...
      !    Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
      !    Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
      !    Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
      !    Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)

      ! Increase the diffusivies to satisfy the min constraints.
      ! All non-zero min constraints on one diffusivity are max constraints on another.
      do k=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
        kh(i,k) = max(kh_bg(i,k), kh_detangle(i,k), &
                      min(kh_min_m(i,k)*kh(i,k-1) + kh0_min_m(i,k), kh(i,k-1)))

        if (kh0_max_m(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_m(i,k))
        if (kh0_max_p(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_p(i,k))
      endif ; enddo ; enddo ! I-loop & k-loop
      ! This is still true... do i=ish,ie ; Kh(I,nz+1) = Kh_bg(I,nz+1) ; enddo
      do k=nz,k_top+1,-1 ; do i=ish,ie ; if (do_i(i)) then
        kh(i,k) = max(kh(i,k), min(kh_min_p(i,k)*kh(i,k+1) + kh0_min_p(i,k), kh(i,k+1)))

        kh_max = max(kh_min_max_p(i,k), kh_max_p(i,k)*kh(i,k+1) + kh0_max_p(i,k))
        kh(i,k) = min(kh(i,k), kh_max)
      endif ; enddo ; enddo ! I-loop & k-loop
      !  All non-zero min constraints on one diffusivity are max constraints on
      ! another layer, so the min constraints can now be discounted.

      ! Decrease the diffusivities to satisfy the max constraints.
        do k=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
          kh_max = max(kh_min_max_m(i,k), kh_max_m(i,k)*kh(i,k-1) + kh0_max_m(i,k))
          if (kh(i,k) > kh_max) kh(i,k) = kh_max
        endif ; enddo ; enddo  ! i- and K-loops

      ! This code tests the solutions...
!     do i=ish,ie
!       Sfn(:) = 0.0 ; uh_here(:) = 0.0
!       do K=k_top,nz
!         if ((Kh(i,K) > Kh_bg(i,K)) .or. (Kh(i,K+1) > Kh_bg(i,K+1))) then
!           if (n==1) then ! u-point.
!             if ((h(i+1,j,k) - h(i,j,k)) * &
!                 ((e(i+1,j,K)-e(i+1,j,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
!               Sfn(K) = -Kh(i,K) * (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
!               Sfn(K+1) = -Kh(i,K+1) * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IdxCu(I,j)
!               uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dy_Cu(I,j)
!               if (abs(uh_here(k)) * min(G%IareaT(i,j), G%IareaT(i+1,j)) > &
!                   (1e-10*GV%m_to_H)) then
!                 if (uh_here(k) * (h(i+1,j,k) - h(i,j,k)) > 0.0) then
!                   call MOM_error(WARNING, "Corrective u-transport is up the thickness gradient.", .true.)
!                 endif
!                 if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh_here(k)) - &
!                      (h(i+1,j,k) + 4.0*dt*G%IareaT(i+1,j)*uh_here(k))) * &
!                     (h(i,j,k) - h(i+1,j,k)) < 0.0) then
!                   call MOM_error(WARNING, "Corrective u-transport is too large.", .true.)
!                 endif
!               endif
!             endif
!           else ! v-point
!             if ((h(i,j+1,k) - h(i,j,k)) * &
!                 ((e(i,j+1,K)-e(i,j+1,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
!               Sfn(K) = -Kh(i,K) * (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)
!               Sfn(K+1) = -Kh(i,K+1) * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IdyCv(i,J)
!               uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dx_Cv(i,J)
!               if (abs(uh_here(K)) * min(G%IareaT(i,j), G%IareaT(i,j+1)) > &
!                   (1e-10*GV%m_to_H)) then
!                 if (uh_here(K) * (h(i,j+1,k) - h(i,j,k)) > 0.0) then
!                   call MOM_error(WARNING, &
!                          "Corrective v-transport is up the thickness gradient.", .true.)
!                 endif
!                 if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh_here(K)) - &
!                      (h(i,j+1,k) + 4.0*dt*G%IareaT(i,j+1)*uh_here(K))) * &
!                     (h(i,j,k) - h(i,j+1,k)) < 0.0) then
!                   call MOM_error(WARNING, &
!                          "Corrective v-transport is too large.", .true.)
!                 endif
!               endif
!             endif
!           endif ! u- or v- selection.
!          !  de_dx(I,K) = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
!         endif
!       enddo
!     enddo

      if (n==1) then ! This is a u-column.
        do k=k_top+1,nz ; do i=ish,ie
          if (kh(i,k) > kh_u(i,j,k)) then
            dkh = (kh(i,k) - kh_u(i,j,k))
            int_slope_u(i,j,k) = dkh / kh(i,k)
            kh_u(i,j,k) = kh(i,k)
          endif
        enddo ; enddo
      else ! This is a v-column.
        do k=k_top+1,nz ; do i=ish,ie
          if (kh(i,k) > kh_v(i,j,k)) then
            dkh = kh(i,k) - kh_v(i,j,k)
            int_slope_v(i,j,k) = dkh / kh(i,k)
            kh_v(i,j,k) = kh(i,k)
          endif
        enddo ; enddo
      endif

    enddo ! j-loop
  enddo  ! n-loop over u- and v- directions.

end subroutine add_detangling_kh

!> Initialize the thickness diffusion module/structure
subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
  type(time_type),         intent(in) :: Time    !< Current model time
  type(ocean_grid_type),   intent(in) :: G       !< Ocean grid structure
  type(verticalgrid_type), intent(in) :: GV      !< Vertical grid structure
  type(unit_scale_type),   intent(in) :: US      !< A dimensional unit scaling type
  type(param_file_type),   intent(in) :: param_file !< Parameter file handles
  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
  type(cont_diag_ptrs),    intent(inout) :: CDp  !< Continuity equation diagnostics
  type(thickness_diffuse_cs), pointer    :: CS   !< Control structure for thickness diffusion

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_thickness_diffuse" ! This module's name.
  real :: omega        ! The Earth's rotation rate [T-1 ~> s-1]
  real :: strat_floor  ! A floor for Brunt-Vasaila frequency in the Ferrari et al. 2010,
                       ! streamfunction formulation, expressed as a fraction of planetary
                       ! rotation [nondim].
  logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags.

  if (associated(cs)) then
    call mom_error(warning, &
      "Thickness_diffuse_init called with an associated control structure.")
    return
  else ; allocate(cs) ; endif

  cs%diag => diag

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "THICKNESSDIFFUSE", cs%thickness_diffuse, &
                 "If true, interface heights are diffused with a "//&
                 "coefficient of KHTH.", default=.false.)
  call get_param(param_file, mdl, "KHTH", cs%Khth, &
                 "The background horizontal thickness diffusivity.", &
                 default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
  call get_param(param_file, mdl, "KHTH_SLOPE_CFF", cs%KHTH_Slope_Cff, &
                 "The nondimensional coefficient in the Visbeck formula "//&
                 "for the interface depth diffusivity", units="nondim", &
                 default=0.0)
  call get_param(param_file, mdl, "KHTH_MIN", cs%KHTH_Min, &
                 "The minimum horizontal thickness diffusivity.", &
                 default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
  call get_param(param_file, mdl, "KHTH_MAX", cs%KHTH_Max, &
                 "The maximum horizontal thickness diffusivity.", &
                 default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
  call get_param(param_file, mdl, "KHTH_MAX_CFL", cs%max_Khth_CFL, &
                 "The maximum value of the local diffusive CFL ratio that "//&
                 "is permitted for the thickness diffusivity. 1.0 is the "//&
                 "marginally unstable value in a pure layered model, but "//&
                 "much smaller numbers (e.g. 0.1) seem to work better for "//&
                 "ALE-based models.", units = "nondimensional", default=0.8)

!  call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%QG_Leith_GM, &
!               "If true, use the QG Leith viscosity as the GM coefficient.", &
!               default=.false.)

  if (cs%max_Khth_CFL < 0.0) cs%max_Khth_CFL = 0.0
  call get_param(param_file, mdl, "DETANGLE_INTERFACES", cs%detangle_interfaces, &
                 "If defined add 3-d structured enhanced interface height "//&
                 "diffusivities to horizontally smooth jagged layers.", &
                 default=.false.)
  cs%detangle_time = 0.0
  if (cs%detangle_interfaces) &
    call get_param(param_file, mdl, "DETANGLE_TIMESCALE", cs%detangle_time, &
                 "A timescale over which maximally jagged grid-scale "//&
                 "thickness variations are suppressed.  This must be "//&
                 "longer than DT, or 0 to use DT.", units="s", default=0.0, scale=us%s_to_T)
  call get_param(param_file, mdl, "KHTH_SLOPE_MAX", cs%slope_max, &
                 "A slope beyond which the calculated isopycnal slope is "//&
                 "not reliable and is scaled away.", units="nondim", default=0.01)
  call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
                 "A diapycnal diffusivity that is used to interpolate "//&
                 "more sensible values of T & S into thin layers.", &
                 units="m2 s-1", default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
  call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", cs%use_FGNV_streamfn, &
                 "If true, use the streamfunction formulation of "//&
                 "Ferrari et al., 2010, which effectively emphasizes "//&
                 "graver vertical modes by smoothing in the vertical.",  &
                 default=.false.)
  call get_param(param_file, mdl, "FGNV_FILTER_SCALE", cs%FGNV_scale, &
                 "A coefficient scaling the vertical smoothing term in the "//&
                 "Ferrari et al., 2010, streamfunction formulation.", &
                 units="nondim", default=1., do_not_log=.not.cs%use_FGNV_streamfn)
  call get_param(param_file, mdl, "FGNV_C_MIN", cs%FGNV_c_min, &
                 "A minium wave speed used in the Ferrari et al., 2010, "//&
                 "streamfunction formulation.", &
                 default=0., units="m s-1", scale=us%m_s_to_L_T, do_not_log=.not.cs%use_FGNV_streamfn)
  call get_param(param_file, mdl, "FGNV_STRAT_FLOOR", strat_floor, &
                 "A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "//&
                 "streamfunction formulation, expressed as a fraction of planetary "//&
                 "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
                 default=1.e-15, units="nondim", do_not_log=.not.cs%use_FGNV_streamfn)
  call get_param(param_file, mdl, "STANLEY_PRM_DET_COEFF", cs%Stanley_det_coeff, &
                 "The coefficient correlating SGS temperature variance with the mean "//&
                 "temperature gradient in the deterministic part of the Stanley parameterization. "//&
                 "Negative values disable the scheme.", units="nondim", default=-1.0)
  call get_param(param_file, mdl, "OMEGA", omega, &
                 "The rotation rate of the earth.", &
                 default=7.2921e-5, units="s-1", scale=us%T_to_s, do_not_log=.not.cs%use_FGNV_streamfn)
  if (cs%use_FGNV_streamfn) cs%N2_floor = (strat_floor*omega)**2
  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)

  call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", cs%GM_src_alt, &
                 "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&
                 "than the streamfunction for the GM source term.", default=.false.)
  call get_param(param_file, mdl, "MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &
                 "If true, uses the GM coefficient formulation from the GEOMETRIC "//&
                 "framework (Marshall et al., 2012).", default=.false.)
  if (cs%MEKE_GEOMETRIC) then
    call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", cs%MEKE_GEOMETRIC_epsilon, &
                 "Minimum Eady growth rate used in the calculation of GEOMETRIC "//&
                 "thickness diffusivity.", units="s-1", default=1.0e-7, scale=us%T_to_s)
    call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &
                 "The nondimensional coefficient governing the efficiency of the GEOMETRIC "//&
                 "thickness diffusion.", units="nondim", default=0.05)

    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, "MEKE_GEOMETRIC_2018_ANSWERS", cs%MEKE_GEOM_answers_2018, &
                 "If true, use expressions in the MEKE_GEOMETRIC calculation that recover the "//&
                 "answers from the original implementation.  Otherwise, use expressions that "//&
                 "satisfy rotational symmetry.", default=default_2018_answers)
  endif

  call get_param(param_file, mdl, "USE_KH_IN_MEKE", cs%Use_KH_in_MEKE, &
                 "If true, uses the thickness diffusivity calculated here to diffuse MEKE.", &
                 default=.false.)

  call get_param(param_file, mdl, "USE_GME", cs%use_GME_thickness_diffuse, &
                 "If true, use the GM+E backscatter scheme in association "//&
                 "with the Gent and McWilliams parameterization.", default=.false.)

  call get_param(param_file, mdl, "USE_GM_WORK_BUG", cs%use_GM_work_bug, &
                 "If true, compute the top-layer work tendency on the u-grid "//&
                 "with the incorrect sign, for legacy reproducibility.", &
                 default=.false.)

  if (cs%use_GME_thickness_diffuse) then
    call safe_alloc_ptr(cs%KH_u_GME,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
    call safe_alloc_ptr(cs%KH_v_GME,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
  endif

  cs%id_uhGM = register_diag_field('ocean_model', 'uhGM', diag%axesCuL, time, &
           'Time Mean Diffusive Zonal Thickness Flux', &
           'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
           y_cell_method='sum', v_extensive=.true.)
  if (cs%id_uhGM > 0) call safe_alloc_ptr(cdp%uhGM,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
  cs%id_vhGM = register_diag_field('ocean_model', 'vhGM', diag%axesCvL, time, &
           'Time Mean Diffusive Meridional Thickness Flux', &
           'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
           x_cell_method='sum', v_extensive=.true.)
  if (cs%id_vhGM > 0) call safe_alloc_ptr(cdp%vhGM,g%isd,g%ied,g%JsdB,g%JedB,g%ke)

  cs%id_GMwork = register_diag_field('ocean_model', 'GMwork', diag%axesT1, time, &
          'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
          'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2, cmor_field_name='tnkebto', &
          cmor_long_name='Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
          cmor_standard_name='tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection')
  if (cs%id_GMwork > 0) call safe_alloc_ptr(cs%GMwork,g%isd,g%ied,g%jsd,g%jed)

  cs%id_KH_u = register_diag_field('ocean_model', 'KHTH_u', diag%axesCui, time, &
           'Parameterized mesoscale eddy advection diffusivity at U-point', &
           'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  cs%id_KH_v = register_diag_field('ocean_model', 'KHTH_v', diag%axesCvi, time, &
           'Parameterized mesoscale eddy advection diffusivity at V-point', &
           'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  cs%id_KH_t = register_diag_field('ocean_model', 'KHTH_t', diag%axesTL, time, &
          'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
          'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
          cmor_field_name='diftrblo', &
          cmor_long_name='Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
          cmor_units='m2 s-1', &
          cmor_standard_name='ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection')

  cs%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', diag%axesCu1, time,         &
           'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)', &
           'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  cs%id_KH_v1 = register_diag_field('ocean_model', 'KHTH_v1', diag%axesCv1, time,         &
           'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)', &
           'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
  cs%id_KH_t1 = register_diag_field('ocean_model', 'KHTH_t1', diag%axesT1, time, &
           'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)', &
           'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

  cs%id_slope_x =  register_diag_field('ocean_model', 'neutral_slope_x', diag%axesCui, time, &
           'Zonal slope of neutral surface', 'nondim')
  if (cs%id_slope_x > 0) call safe_alloc_ptr(cs%diagSlopeX,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
  cs%id_slope_y =  register_diag_field('ocean_model', 'neutral_slope_y', diag%axesCvi, time, &
           'Meridional slope of neutral surface', 'nondim')
  if (cs%id_slope_y > 0) call safe_alloc_ptr(cs%diagSlopeY,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
  cs%id_sfn_x =  register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, time, &
           'Parameterized Zonal Overturning Streamfunction', &
           'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
  cs%id_sfn_y =  register_diag_field('ocean_model', 'GM_sfn_y', diag%axesCvi, time, &
           'Parameterized Meridional Overturning Streamfunction', &
           'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
  cs%id_sfn_unlim_x =  register_diag_field('ocean_model', 'GM_sfn_unlim_x', diag%axesCui, time, &
           'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
           'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
  cs%id_sfn_unlim_y =  register_diag_field('ocean_model', 'GM_sfn_unlim_y', diag%axesCvi, time, &
           'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
           'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)

end subroutine thickness_diffuse_init

!> Copies ubtav and vbtav from private type into arrays
subroutine thickness_diffuse_get_kh(CS, KH_u_GME, KH_v_GME, G)
  type(thickness_diffuse_cs),          pointer     :: CS   !< Control structure for
                                                   !! this module
  type(ocean_grid_type),               intent(in)  :: G    !< Grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: KH_u_GME !< interface height
                                                   !! diffusivities at u-faces [L2 T-1 ~> m2 s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: KH_v_GME !< interface height
                                                   !! diffusivities at v-faces [L2 T-1 ~> m2 s-1]
  ! Local variables
  integer :: i,j,k

  do k=1,g%ke+1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
    kh_u_gme(i,j,k) = cs%KH_u_GME(i,j,k)
  enddo ; enddo ; enddo

  do k=1,g%ke+1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
    kh_v_gme(i,j,k) = cs%KH_v_GME(i,j,k)
  enddo ; enddo ; enddo

end subroutine thickness_diffuse_get_kh

!> Deallocate the thickness diffusion control structure
subroutine thickness_diffuse_end(CS)
  type(thickness_diffuse_cs), pointer :: CS !< Control structure for thickness diffusion

  if (associated(cs)) deallocate(cs)
end subroutine thickness_diffuse_end

!> \namespace mom_thickness_diffuse
!!
!! \section section_gm Thickness diffusion (aka Gent-McWilliams)
!!
!! Thickness diffusion is implemented via along-layer mass fluxes
!! \f[
!! h^\dagger \leftarrow h^n - \Delta t \nabla \cdot ( \vec{uh}^* )
!! \f]
!! where the mass fluxes are cast as the difference in vector streamfunction
!!
!! \f[
!! \vec{uh}^* = \delta_k \vec{\psi} .
!! \f]
!!
!! The GM implementation of thickness diffusion made the streamfunction proportional to the potential density slope
!! \f[
!! \vec{\psi} = - \kappa_h \frac{\nabla_z \rho}{\partial_z \rho}
!! = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = \kappa_h \frac{M^2}{N^2}
!! \f]
!! but for robustness the scheme is implemented as
!! \f[
!! \vec{\psi} = \kappa_h \frac{M^2}{\sqrt{N^4 + M^4}}
!! \f]
!! since the quantity \f$\frac{M^2}{\sqrt{N^2 + M^2}}\f$ is bounded between $-1$ and $1$ and does not change sign
!! if \f$N^2<0\f$.
!!
!! Optionally, the method of Ferrari et al, 2010, can be used to obtain the streamfunction which solves the
!! vertically elliptic equation:
!! \f[
!! \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi  = ( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}}
!! \f]
!! which recovers the previous streamfunction relation in the limit that \f$ c \rightarrow 0 \f$.
!! Here, \f$c=\max(c_{min},c_g)\f$ is the maximum of either \f$c_{min}\f$ and either the first baroclinic mode
!! wave-speed or the equivalent barotropic mode wave-speed.
!! \f$N_*^2 = \max(N^2,0)\f$ is a non-negative form of the square of the Brunt-Vaisala frequency.
!! The parameter \f$\gamma_F\f$ is used to reduce the vertical smoothing length scale.
!! \f[
!! \kappa_h = \left( \kappa_o + \alpha_{s} L_{s}^2 < S N > + \alpha_{M} \kappa_{M} \right) r(\Delta x,L_d)
!! \f]
!! where \f$ S \f$ is the isoneutral slope magnitude, \f$ N \f$ is the square root of Brunt-Vaisala frequency,
!! \f$\kappa_{M}\f$ is the diffusivity calculated by the MEKE parameterization (mom_meke module) and
!! \f$ r(\Delta x,L_d) \f$ is a function of the local resolution (ratio of grid-spacing, \f$\Delta x\f$,
!! to deformation radius, \f$L_d\f$). The length \f$L_s\f$ is provided by the mom_lateral_mixing_coeffs module
!! (enabled with <code>USE_VARIABLE_MIXING=True</code> and the term \f$<SN>\f$ is the vertical average slope
!! times the Brunt-Vaisala frequency prescribed by Visbeck et al., 1996.
!!
!! The result of the above expression is subsequently bounded by minimum and maximum values, including an upper
!! diffusivity consistent with numerical stability (\f$ \kappa_{cfl} \f$ is calculated internally).
!! \f[
!! \kappa_h \leftarrow \min{\left( \kappa_{max}, \kappa_{cfl}, \max{\left( \kappa_{min}, \kappa_h \right)} \right)}
!!                      f(c_g,z)
!! \f]
!!
!! where \f$f(c_g,z)\f$ is a vertical structure function.
!! \f$f(c_g,z)\f$ is calculated in module mom_lateral_mixing_coeffs.
!! If <code>KHTH_USE_EBT_STRUCT=True</code> then \f$f(c_g,z)\f$ is set to look like the equivalent barotropic
!! modal velocity structure. Otherwise \f$f(c_g,z)=1\f$ and the diffusivity is independent of depth.
!!
!! In order to calculate meaningful slopes in vanished layers, temporary copies of the thermodynamic variables
!! are passed through a vertical smoother, function vert_fill_ts():
!! \f{eqnarray*}{
!! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] \theta & \leftarrow & \theta \\
!! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] s & \leftarrow & s
!! \f}
!!
!! \subsection section_khth_module_parameters Module mom_thickness_diffuse parameters
!!
!! | Symbol                | Module parameter |
!! | ------                | --------------- |
!! | -                     | <code>THICKNESSDIFFUSE</code> |
!! | \f$ \kappa_o \f$      | <code>KHTH</code> |
!! | \f$ \alpha_{s} \f$    | <code>KHTH_SLOPE_CFF</code> |
!! | \f$ \kappa_{min} \f$  | <code>KHTH_MIN</code> |
!! | \f$ \kappa_{max} \f$  | <code>KHTH_MAX</code> |
!! | -                     | <code>KHTH_MAX_CFL</code> |
!! | \f$ \kappa_{smth} \f$ | <code>KD_SMOOTH</code> |
!! | \f$ \alpha_{M} \f$    | <code>MEKE_KHTH_FAC</code> (from mom_meke module) |
!! | -                     | <code>KHTH_USE_EBT_STRUCT</code> (from mom_lateral_mixing_coeffs module) |
!! | -                     | <code>KHTH_USE_FGNV_STREAMFUNCTION</code> |
!! | \f$ \gamma_F \f$      | <code>FGNV_FILTER_SCALE</code> |
!! | \f$ c_{min} \f$       | <code>FGNV_C_MIN</code> |
!!
!! \subsection section_khth_module_reference References
!!
!! Ferrari, R., S.M. Griffies, A.J.G. Nurser and G.K. Vallis, 2010:
!! A boundary-value problem for the parameterized mesoscale eddy transport.
!! Ocean Modelling, 32, 143-156. http://doi.org/10.1016/j.ocemod.2010.01.004
!!
!! Viscbeck, M., J.C. Marshall, H. Jones, 1996:
!! Dynamics of isolated convective regions in the ocean. J. Phys. Oceangr., 26, 1721-1734.
!! http://dx.doi.org/10.1175/1520-0485(1996)026%3C1721:DOICRI%3E2.0.CO;2

end module mom_thickness_diffuse