MOM_lateral_mixing_coeffs.F90

!> Variable mixing coefficients
module mom_lateral_mixing_coeffs

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

use mom_debugging,         only : hchksum, uvchksum
use mom_error_handler,     only : mom_error, fatal, warning, mom_mesg
use mom_diag_mediator,     only : register_diag_field, safe_alloc_ptr, post_data
use mom_diag_mediator,     only : diag_ctrl, time_type, query_averaging_enabled
use mom_domains,           only : create_group_pass, do_group_pass
use mom_domains,           only : group_pass_type, pass_var, pass_vector
use mom_file_parser,       only : get_param, log_version, param_file_type
use mom_interface_heights, only : find_eta
use mom_isopycnal_slopes,  only : calc_isoneutral_slopes
use mom_grid,              only : ocean_grid_type
use mom_unit_scaling,      only : unit_scale_type
use mom_variables,         only : thermo_var_ptrs
use mom_verticalgrid,      only : verticalgrid_type
use mom_wave_speed,        only : wave_speed, wave_speed_cs, wave_speed_init
use mom_open_boundary,     only : ocean_obc_type, obc_none

implicit none ; private

#include <MOM_memory.h>

!> Variable mixing coefficients
type, public :: varmix_cs
  logical :: use_variable_mixing  !< If true, use the variable mixing.
  logical :: resoln_scaled_kh     !< If true, scale away the Laplacian viscosity
                                  !! when the deformation radius is well resolved.
  logical :: resoln_scaled_khth   !< If true, scale away the thickness diffusivity
                                  !! when the deformation radius is well resolved.
  logical :: depth_scaled_khth    !< If true, KHTH is scaled away when the depth is
                                  !! shallower than a reference depth.
  logical :: resoln_scaled_khtr   !< If true, scale away the tracer diffusivity
                                  !! when the deformation radius is well resolved.
  logical :: interpolate_res_fn   !< If true, interpolate the resolution function
                                  !! to the velocity points from the thickness
                                  !! points; otherwise interpolate the wave
                                  !! speed and calculate the resolution function
                                  !! independently at each point.
  logical :: use_stored_slopes    !< If true, stores isopycnal slopes in this structure.
  logical :: resoln_use_ebt       !< If true, uses the equivalent barotropic wave speed instead
                                  !! of first baroclinic wave for calculating the resolution fn.
  logical :: khth_use_ebt_struct  !< If true, uses the equivalent barotropic structure
                                  !! as the vertical structure of thickness diffusivity.
  logical :: calculate_cg1        !< If true, calls wave_speed() to calculate the first
                                  !! baroclinic wave speed and populate CS%cg1.
                                  !! This parameter is set depending on other parameters.
  logical :: calculate_rd_dx      !< If true, calculates Rd/dx and populate CS%Rd_dx_h.
                                  !! This parameter is set depending on other parameters.
  logical :: calculate_res_fns    !< If true, calculate all the resolution factors.
                                  !! This parameter is set depending on other parameters.
  logical :: calculate_depth_fns !< If true, calculate all the depth factors.
                                  !! This parameter is set depending on other parameters.
  logical :: calculate_eady_growth_rate !< If true, calculate all the Eady growth rate.
                                  !! This parameter is set depending on other parameters.
  real, dimension(:,:), pointer :: &
    sn_u => null(), &     !< S*N at u-points [T-1 ~> s-1]
    sn_v => null(), &     !< S*N at v-points [T-1 ~> s-1]
    l2u => null(), &      !< Length scale^2 at u-points [L2 ~> m2]
    l2v => null(), &      !< Length scale^2 at v-points [L2 ~> m2]
    cg1 => null(), &      !< The first baroclinic gravity wave speed [L T-1 ~> m s-1].
    res_fn_h => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at h points [nondim].
    res_fn_q => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at q points [nondim].
    res_fn_u => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at u points [nondim].
    res_fn_v => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at v points [nondim].
    depth_fn_u => null(), & !< Non-dimensional function of the ratio of the depth to
                            !! a reference depth (maximum 1) at u points [nondim]
    depth_fn_v => null(), & !< Non-dimensional function of the ratio of the depth to
                            !! a reference depth (maximum 1) at v points [nondim]
    beta_dx2_h => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at h points [L T-1 ~> m s-1].
    beta_dx2_q => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at q points [L T-1 ~> m s-1].
    beta_dx2_u => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at u points [L T-1 ~> m s-1].
    beta_dx2_v => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at v points [L T-1 ~> m s-1].
    f2_dx2_h => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at h [L2 T-2 ~> m2 s-2].
    f2_dx2_q => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at q [L2 T-2 ~> m2 s-2].
    f2_dx2_u => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at u [L2 T-2 ~> m2 s-2].
    f2_dx2_v => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at v [L2 T-2 ~> m2 s-2].
    rd_dx_h => null()     !< Deformation radius over grid spacing [nondim]

  real, dimension(:,:,:), pointer :: &
    slope_x => null(), &  !< Zonal isopycnal slope [nondim]
    slope_y => null(), &  !< Meridional isopycnal slope [nondim]
    ebt_struct => null()  !< Vertical structure function to scale diffusivities with [nondim]
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
    laplac3_const_u       !< Laplacian metric-dependent constants [L3 ~> m3]

  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
    laplac3_const_v       !< Laplacian metric-dependent constants [L3 ~> m3]

  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
    kh_u_qg               !< QG Leith GM coefficient at u-points [L2 T-1 ~> m2 s-1]

  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
    kh_v_qg               !< QG Leith GM coefficient at v-points [L2 T-1 ~> m2 s-1]

  ! Parameters
  logical :: use_visbeck  !< Use Visbeck formulation for thickness diffusivity
  integer :: varmix_ktop  !< Top layer to start downward integrals
  real :: visbeck_l_scale !< Fixed length scale in Visbeck formula
  real :: res_coef_khth   !< A non-dimensional number that determines the function
                          !! of resolution, used for thickness and tracer mixing, as:
                          !!  F = 1 / (1 + (Res_coef_khth*Ld/dx)^Res_fn_power)
  real :: res_coef_visc   !< A non-dimensional number that determines the function
                          !! of resolution, used for lateral viscosity, as:
                          !!  F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power)
  real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m]
  real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim]
  real :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1]
  integer :: res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function.  Any
                               !! positive integer power may be used, but even powers
                               !! and especially 2 are coded to be more efficient.
  integer :: res_fn_power_visc !< The power of dx/Ld in the Kh resolution function.  Any
                               !! positive integer power may be used, but even powers
                               !! and especially 2 are coded to be more efficient.
  real :: visbeck_s_max   !< Upper bound on slope used in Eady growth rate [nondim].

  ! Leith parameters
  logical :: use_qg_leith_gm      !< If true, uses the QG Leith viscosity as the GM coefficient
  logical :: use_beta_in_qg_leith !< If true, includes the beta term in the QG Leith GM coefficient

  ! Diagnostics
  !>@{
  !! Diagnostic identifier
  integer :: id_sn_u=-1, id_sn_v=-1, id_l2u=-1, id_l2v=-1, id_res_fn = -1
  integer :: id_n2_u=-1, id_n2_v=-1, id_s2_u=-1, id_s2_v=-1
  integer :: id_rd_dx=-1, id_kh_u_qg = -1, id_kh_v_qg = -1
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                   !! timing of diagnostic output.
  !>@}

  type(wave_speed_cs), pointer :: wave_speed_csp => null() !< Wave speed control structure
  type(group_pass_type) :: pass_cg1 !< For group halo pass
  logical :: debug      !< If true, write out checksums of data for debugging
end type varmix_cs

public varmix_init, calc_slope_functions, calc_resoln_function
public calc_qg_leith_viscosity, calc_depth_function

contains

!> Calculates the non-dimensional depth functions.
subroutine calc_depth_function(G, CS)
  type(ocean_grid_type),                    intent(in) :: g  !< Ocean grid structure
  type(varmix_cs),                          pointer       :: cs !< Variable mixing coefficients

  ! Local variables
  integer :: is, ie, js, je, isq, ieq, jsq, jeq
  integer :: i, j
  real    :: h0 ! local variable for reference depth
  real    :: expo ! exponent used in the depth dependent scaling
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  if (.not. associated(cs)) call mom_error(fatal, "calc_depth_function:"// &
         "Module must be initialized before it is used.")
  if (.not. cs%calculate_depth_fns) return
  if (.not. associated(cs%Depth_fn_u)) call mom_error(fatal, &
    "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
  if (.not. associated(cs%Depth_fn_v)) call mom_error(fatal, &
    "calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.")

  h0 = cs%depth_scaled_khth_h0
  expo = cs%depth_scaled_khth_exp
!$OMP do
  do j=js,je ; do i=is-1,ieq
    cs%Depth_fn_u(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))/h0))**expo
  enddo ; enddo
!$OMP do
  do j=js-1,jeq ; do i=is,ie
    cs%Depth_fn_v(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))/h0))**expo
  enddo ; enddo

end subroutine calc_depth_function

!> Calculates and stores the non-dimensional resolution functions
subroutine calc_resoln_function(h, tv, G, GV, US, CS)
  type(ocean_grid_type),                    intent(inout) :: g  !< Ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h  !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                    intent(in)    :: tv !< Thermodynamic variables
  type(verticalgrid_type),                  intent(in)    :: gv !< Vertical grid structure
  type(unit_scale_type),                    intent(in)    :: us !< A dimensional unit scaling type
  type(varmix_cs),                          pointer       :: cs !< Variable mixing coefficients

  ! Local variables
  ! Depending on the power-function being used, dimensional rescaling may be limited, so some
  ! of the following variables have units that depend on that power.
  real :: cg1_q  ! The gravity wave speed interpolated to q points [L T-1 ~> m s-1] or [m s-1].
  real :: cg1_u  ! The gravity wave speed interpolated to u points [L T-1 ~> m s-1] or [m s-1].
  real :: cg1_v  ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
  real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
  integer :: power_2
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
  integer :: i, j, k
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  if (.not. associated(cs)) call mom_error(fatal, "calc_resoln_function:"// &
         "Module must be initialized before it is used.")
  if (cs%calculate_cg1) then
    if (.not. associated(cs%cg1)) call mom_error(fatal, &
      "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
    if (cs%khth_use_ebt_struct) then
      if (.not. associated(cs%ebt_struct)) call mom_error(fatal, &
        "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
      if (cs%Resoln_use_ebt) then
        ! Both resolution fn and vertical structure are using EBT
        call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct)
      else
        ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode
        call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct, &
                        use_ebt_mode=.true.)
        call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
      endif
      call pass_var(cs%ebt_struct, g%Domain)
    else
      call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
    endif

    call create_group_pass(cs%pass_cg1, cs%cg1, g%Domain)
    call do_group_pass(cs%pass_cg1, g%Domain)
  endif

  ! Calculate and store the ratio between deformation radius and grid-spacing
  ! at h-points [nondim].
  if (cs%calculate_rd_dx) then
    if (.not. associated(cs%Rd_dx_h)) call mom_error(fatal, &
      "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
    !$OMP parallel do default(shared)
    do j=js-1,je+1 ; do i=is-1,ie+1
      cs%Rd_dx_h(i,j) = cs%cg1(i,j) / &
            (sqrt(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))
    enddo ; enddo
    if (query_averaging_enabled(cs%diag)) then
      if (cs%id_Rd_dx > 0) call post_data(cs%id_Rd_dx, cs%Rd_dx_h, cs%diag)
    endif
  endif

  if (.not. cs%calculate_res_fns) return

  if (.not. associated(cs%Res_fn_h)) call mom_error(fatal, &
    "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%Res_fn_q)) call mom_error(fatal, &
    "calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%Res_fn_u)) call mom_error(fatal, &
    "calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%Res_fn_v)) call mom_error(fatal, &
    "calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%f2_dx2_h)) call mom_error(fatal, &
    "calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%f2_dx2_q)) call mom_error(fatal, &
    "calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%f2_dx2_u)) call mom_error(fatal, &
    "calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%f2_dx2_v)) call mom_error(fatal, &
    "calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%beta_dx2_h)) call mom_error(fatal, &
    "calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%beta_dx2_q)) call mom_error(fatal, &
    "calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%beta_dx2_u)) call mom_error(fatal, &
    "calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
  if (.not. associated(cs%beta_dx2_v)) call mom_error(fatal, &
    "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")

  !   Do this calculation on the extent used in MOM_hor_visc.F90, and
  ! MOM_tracer.F90 so that no halo update is needed.

!$OMP parallel default(none) shared(is,ie,js,je,Ieq,Jeq,CS,US) &
!$OMP                       private(dx_term,cg1_q,power_2,cg1_u,cg1_v)
  if (cs%Res_fn_power_visc >= 100) then
!$OMP do
    do j=js-1,je+1 ; do i=is-1,ie+1
      dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
      if ((cs%Res_coef_visc * cs%cg1(i,j))**2 > dx_term) then
        cs%Res_fn_h(i,j) = 0.0
      else
        cs%Res_fn_h(i,j) = 1.0
      endif
    enddo ; enddo
!$OMP do
    do j=js-1,jeq ; do i=is-1,ieq
      cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
      dx_term = cs%f2_dx2_q(i,j) +  cg1_q * cs%beta_dx2_q(i,j)
      if ((cs%Res_coef_visc * cg1_q)**2 > dx_term) then
        cs%Res_fn_q(i,j) = 0.0
      else
        cs%Res_fn_q(i,j) = 1.0
      endif
    enddo ; enddo
  elseif (cs%Res_fn_power_visc == 2) then
!$OMP do
    do j=js-1,je+1 ; do i=is-1,ie+1
      dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
      cs%Res_fn_h(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**2)
    enddo ; enddo
!$OMP do
    do j=js-1,jeq ; do i=is-1,ieq
      cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
      dx_term = cs%f2_dx2_q(i,j) +  cg1_q * cs%beta_dx2_q(i,j)
      cs%Res_fn_q(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cg1_q)**2)
    enddo ; enddo
  elseif (mod(cs%Res_fn_power_visc, 2) == 0) then
    power_2 = cs%Res_fn_power_visc / 2
!$OMP do
    do j=js-1,je+1 ; do i=is-1,ie+1
      dx_term = (us%L_T_to_m_s**2*(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**power_2
      cs%Res_fn_h(i,j) = dx_term / &
          (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
    enddo ; enddo
!$OMP do
    do j=js-1,jeq ; do i=is-1,ieq
      cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
      dx_term = (us%L_T_to_m_s**2*(cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)))**power_2
      cs%Res_fn_q(i,j) = dx_term / &
          (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
    enddo ; enddo
  else
!$OMP do
    do j=js-1,je+1 ; do i=is-1,ie+1
      dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_h(i,j) + &
                                    cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**cs%Res_fn_power_visc
      cs%Res_fn_h(i,j) = dx_term / &
         (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
    enddo ; enddo
!$OMP do
    do j=js-1,jeq ; do i=is-1,ieq
      cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
      dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_q(i,j) + &
                                    cg1_q * cs%beta_dx2_q(i,j)))**cs%Res_fn_power_visc
      cs%Res_fn_q(i,j) = dx_term / &
          (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
    enddo ; enddo
  endif

  if (cs%interpolate_Res_fn) then
    do j=js,je ; do i=is-1,ieq
      cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
    enddo ; enddo
    do j=js-1,jeq ; do i=is,ie
      cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
    enddo ; enddo
  else ! .not.CS%interpolate_Res_fn
    if (cs%Res_fn_power_khth >= 100) then
!$OMP do
      do j=js,je ; do i=is-1,ieq
        cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
        dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
        if ((cs%Res_coef_khth * cg1_u)**2 > dx_term) then
          cs%Res_fn_u(i,j) = 0.0
        else
          cs%Res_fn_u(i,j) = 1.0
        endif
      enddo ; enddo
!$OMP do
      do j=js-1,jeq ; do i=is,ie
        cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
        dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
        if ((cs%Res_coef_khth * cg1_v)**2 > dx_term) then
          cs%Res_fn_v(i,j) = 0.0
        else
          cs%Res_fn_v(i,j) = 1.0
        endif
      enddo ; enddo
    elseif (cs%Res_fn_power_khth == 2) then
!$OMP do
      do j=js,je ; do i=is-1,ieq
        cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
        dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
        cs%Res_fn_u(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_u)**2)
      enddo ; enddo
!$OMP do
      do j=js-1,jeq ; do i=is,ie
        cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
        dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
        cs%Res_fn_v(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_v)**2)
      enddo ; enddo
    elseif (mod(cs%Res_fn_power_khth, 2) == 0) then
      power_2 = cs%Res_fn_power_khth / 2
!$OMP do
      do j=js,je ; do i=is-1,ieq
        cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
        dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)))**power_2
        cs%Res_fn_u(i,j) = dx_term / &
            (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
      enddo ; enddo
!$OMP do
      do j=js-1,jeq ; do i=is,ie
        cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
        dx_term = (us%L_T_to_m_s**2 * (cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)))**power_2
        cs%Res_fn_v(i,j) = dx_term / &
            (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
      enddo ; enddo
    else
!$OMP do
      do j=js,je ; do i=is-1,ieq
        cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
        dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_u(i,j) + &
                                      cg1_u * cs%beta_dx2_u(i,j)))**cs%Res_fn_power_khth
        cs%Res_fn_u(i,j) = dx_term / &
            (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
      enddo ; enddo
!$OMP do
      do j=js-1,jeq ; do i=is,ie
        cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
        dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_v(i,j) + &
                                      cg1_v * cs%beta_dx2_v(i,j)))**cs%Res_fn_power_khth
        cs%Res_fn_v(i,j) = dx_term / &
            (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
      enddo ; enddo
    endif
  endif
!$OMP end parallel

  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_Res_fn > 0) call post_data(cs%id_Res_fn, cs%Res_fn_h, cs%diag)
  endif

end subroutine calc_resoln_function

!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
!! style scaling of diffusivity
subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
  type(ocean_grid_type),                    intent(inout) :: 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]
  type(thermo_var_ptrs),                    intent(in)    :: tv !< Thermodynamic variables
  real,                                     intent(in)    :: dt !< Time increment [T ~> s]
  type(varmix_cs),                          pointer       :: cs !< Variable mixing coefficients
  type(ocean_obc_type),           optional, pointer       :: obc !< Open boundaries control structure.
  ! Local variables
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    e             ! The interface heights relative to mean sea level [Z ~> m].
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: n2_u ! Square of Brunt-Vaisala freq at u-points [T-2 ~> s-2]
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: n2_v ! Square of Brunt-Vaisala freq at v-points [T-2 ~> s-2]

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

  if (cs%calculate_Eady_growth_rate) then
    call find_eta(h, tv, g, gv, us, e, halo_size=2)
    if (cs%use_stored_slopes) then
      call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, &
                                  cs%slope_x, cs%slope_y, n2_u, n2_v, 1, obc=obc)
      call calc_visbeck_coeffs(h, cs%slope_x, cs%slope_y, n2_u, n2_v, g, gv, us, cs, obc=obc)
!     call calc_slope_functions_using_just_e(h, G, CS, e, .false.)
    else
      !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y)
      call calc_slope_functions_using_just_e(h, g, gv, us, cs, e, .true., obc=obc)
    endif
  endif

  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_SN_u > 0) call post_data(cs%id_SN_u, cs%SN_u, cs%diag)
    if (cs%id_SN_v > 0) call post_data(cs%id_SN_v, cs%SN_v, cs%diag)
    if (cs%id_L2u > 0)  call post_data(cs%id_L2u, cs%L2u, cs%diag)
    if (cs%id_L2v > 0)  call post_data(cs%id_L2v, cs%L2v, cs%diag)
    if (cs%calculate_Eady_growth_rate .and. cs%use_stored_slopes) then
      if (cs%id_N2_u > 0) call post_data(cs%id_N2_u, n2_u, cs%diag)
      if (cs%id_N2_v > 0) call post_data(cs%id_N2_v, n2_v, cs%diag)
    endif
  endif

end subroutine calc_slope_functions

!> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al.
subroutine calc_visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
  type(ocean_grid_type),                       intent(inout) :: G  !< Ocean grid structure
  type(verticalGrid_type),                     intent(in)    :: GV !< Vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)    :: h  !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: slope_x !< Zonal isoneutral slope
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: N2_u    !< Buoyancy (Brunt-Vaisala) frequency
                                                                        !! at u-points [T-2 ~> s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in)    :: slope_y !< Meridional isoneutral slope
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in)    :: N2_v    !< Buoyancy (Brunt-Vaisala) frequency
                                                                        !! at v-points [T-2 ~> s-2]
  type(unit_scale_type),                       intent(in)    :: US !< A dimensional unit scaling type
  type(VarMix_CS),                             pointer       :: CS !< Variable mixing coefficients
  type(ocean_OBC_type),              optional, pointer       :: OBC !< Open boundaries control structure.

  ! Local variables
  real :: S2            ! Interface slope squared [nondim]
  real :: N2            ! Positive buoyancy frequency or zero [T-2 ~> s-2]
  real :: Hup, Hdn      ! Thickness from above, below [H ~> m or kg m-2]
  real :: H_geom        ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
  integer :: is, ie, js, je, nz
  integer :: i, j, k, kb_max
  integer :: l_seg
  real :: S2max, wNE, wSE, wSW, wNW
  real :: H_u(SZIB_(G)), H_v(SZI_(G))
  real :: S2_u(SZIB_(G), SZJ_(G))
  real :: S2_v(SZI_(G), SZJB_(G))
  logical :: local_open_u_BC, local_open_v_BC

  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
         "Module must be initialized before it is used.")
  if (.not. cs%calculate_Eady_growth_rate) return
  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
         "%SN_u is not associated with use_variable_mixing.")
  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
         "%SN_v is not associated with use_variable_mixing.")

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

  local_open_u_bc = .false.
  local_open_v_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_u_bc = obc%open_u_BCs_exist_globally
    local_open_v_bc = obc%open_v_BCs_exist_globally
  endif ; endif

  s2max = cs%Visbeck_S_max**2

  !$OMP parallel do default(shared)
  do j=js-1,je+1 ; do i=is-1,ie+1
    cs%SN_u(i,j) = 0.0
    cs%SN_v(i,j) = 0.0
  enddo ; enddo

  ! To set the length scale based on the deformation radius, use wave_speed to
  ! calculate the first-mode gravity wave speed and then blend the equatorial
  ! and midlatitude deformation radii, using calc_resoln_function as a template.

  !$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
  do j = js,je
    do i=is-1,ie
      cs%SN_u(i,j) = 0. ; h_u(i) = 0. ; s2_u(i,j) = 0.
    enddo
    do k=2,nz ; do i=is-1,ie
      hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
      hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
      h_geom = sqrt( hdn * hup )
     !H_geom = H_geom * sqrt(N2) ! WKB-ish
     !H_geom = H_geom * N2       ! WKB-ish
      wse = g%mask2dCv(i+1,j-1) * ( (h(i+1,j,k)*h(i+1,j-1,k)) * (h(i+1,j,k-1)*h(i+1,j-1,k-1)) )
      wnw = g%mask2dCv(i  ,j  ) * ( (h(i  ,j,k)*h(i  ,j+1,k)) * (h(i  ,j,k-1)*h(i  ,j+1,k-1)) )
      wne = g%mask2dCv(i+1,j  ) * ( (h(i+1,j,k)*h(i+1,j+1,k)) * (h(i+1,j,k-1)*h(i+1,j+1,k-1)) )
      wsw = g%mask2dCv(i  ,j-1) * ( (h(i  ,j,k)*h(i  ,j-1,k)) * (h(i  ,j,k-1)*h(i  ,j-1,k-1)) )
      s2 =  slope_x(i,j,k)**2 + &
              ((wnw*slope_y(i,j,k)**2 + wse*slope_y(i+1,j-1,k)**2) + &
               (wne*slope_y(i+1,j,k)**2 + wsw*slope_y(i,j-1,k)**2) ) / &
              ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
      if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2

      n2 = max(0., n2_u(i,j,k))
      cs%SN_u(i,j) = cs%SN_u(i,j) + sqrt( s2*n2 )*h_geom
      s2_u(i,j) = s2_u(i,j) + s2*h_geom
      h_u(i) = h_u(i) + h_geom
    enddo ; enddo
    do i=is-1,ie
      if (h_u(i)>0.) then
        cs%SN_u(i,j) = g%mask2dCu(i,j) * cs%SN_u(i,j) / h_u(i)
        s2_u(i,j) =  g%mask2dCu(i,j) * s2_u(i,j) / h_u(i)
      else
        cs%SN_u(i,j) = 0.
      endif
      if (local_open_u_bc) then
        l_seg = obc%segnum_u(i,j)

        if (l_seg /= obc_none) then
          if (obc%segment(l_seg)%open) then
            cs%SN_u(i,j) = 0.
          endif
        endif
      endif
    enddo
  enddo

  !$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
  do j = js-1,je
    do i=is,ie
      cs%SN_v(i,j) = 0. ; h_v(i) = 0. ; s2_v(i,j) = 0.
    enddo
    do k=2,nz ; do i=is,ie
      hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
      hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
      h_geom = sqrt( hdn * hup )
     !H_geom = H_geom * sqrt(N2) ! WKB-ish
     !H_geom = H_geom * N2       ! WKB-ish
      wse = g%mask2dCu(i,j)     * ( (h(i,j  ,k)*h(i+1,j  ,k)) * (h(i,j  ,k-1)*h(i+1,j  ,k-1)) )
      wnw = g%mask2dCu(i-1,j+1) * ( (h(i,j+1,k)*h(i-1,j+1,k)) * (h(i,j+1,k-1)*h(i-1,j+1,k-1)) )
      wne = g%mask2dCu(i,j+1)   * ( (h(i,j+1,k)*h(i+1,j+1,k)) * (h(i,j+1,k-1)*h(i+1,j+1,k-1)) )
      wsw = g%mask2dCu(i-1,j)   * ( (h(i,j  ,k)*h(i-1,j  ,k)) * (h(i,j  ,k-1)*h(i-1,j  ,k-1)) )
      s2 = slope_y(i,j,k)**2 + &
             ((wse*slope_x(i,j,k)**2 + wnw*slope_x(i-1,j+1,k)**2) + &
              (wne*slope_x(i,j+1,k)**2 + wsw*slope_x(i-1,j,k)**2) ) / &
             ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
      if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2

      n2 = max(0., n2_v(i,j,k))
      cs%SN_v(i,j) = cs%SN_v(i,j) + sqrt( s2*n2 )*h_geom
      s2_v(i,j) = s2_v(i,j) + s2*h_geom
      h_v(i) = h_v(i) + h_geom
    enddo ; enddo
    do i=is,ie
      if (h_v(i)>0.) then
        cs%SN_v(i,j) = g%mask2dCv(i,j) * cs%SN_v(i,j) / h_v(i)
        s2_v(i,j) = g%mask2dCv(i,j) * s2_v(i,j) / h_v(i)
      else
        cs%SN_v(i,j) = 0.
      endif
      if (local_open_v_bc) then
        l_seg = obc%segnum_v(i,j)

        if (l_seg /= obc_none) then
          if (obc%segment(obc%segnum_v(i,j))%open) then
            cs%SN_v(i,j) = 0.
          endif
        endif
      endif
    enddo
  enddo

! Offer diagnostic fields for averaging.
  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_S2_u > 0) call post_data(cs%id_S2_u, s2_u, cs%diag)
    if (cs%id_S2_v > 0) call post_data(cs%id_S2_v, s2_v, cs%diag)
  endif

  if (cs%debug) then
    call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, g%HI, haloshift=1)
    call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", n2_u, n2_v, g%HI, &
                  scale=us%s_to_T**2, scalar_pair=.true.)
    call uvchksum("calc_Visbeck_coeffs SN_[uv]", cs%SN_u, cs%SN_v, g%HI, &
                  scale=us%s_to_T, scalar_pair=.true.)
  endif

end subroutine calc_visbeck_coeffs

!> The original calc_slope_function() that calculated slopes using
!! interface positions only, not accounting for density variations.
subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slopes, OBC)
  type(ocean_grid_type),                      intent(inout) :: G  !< Ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2]
  type(verticalGrid_type),                    intent(in)    :: GV !< Vertical grid structure
  type(unit_scale_type),                      intent(in)    :: US !< A dimensional unit scaling type
  type(VarMix_CS),                            pointer       :: CS !< Variable mixing coefficients
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: e  !< Interface position [Z ~> m]
  logical,                                    intent(in)    :: calculate_slopes !< If true, calculate slopes internally
                                                                  !! otherwise use slopes stored in CS
  type(ocean_OBC_type),             optional, pointer       :: OBC !< Open boundaries control structure.
  ! Local variables
  real :: E_x(SZIB_(G), SZJ_(G))  ! X-slope of interface at u points [nondim] (for diagnostics)
  real :: E_y(SZI_(G), SZJB_(G))  ! Y-slope of interface at v points [nondim] (for diagnostics)
  real :: H_cutoff      ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
  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 :: S2            ! Interface slope squared [nondim]
  real :: N2            ! Brunt-Vaisala frequency squared [T-2 ~> s-2]
  real :: Hup, Hdn      ! Thickness from above, below [H ~> m or kg m-2]
  real :: H_geom        ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
  real :: one_meter     ! One meter in thickness units [H ~> m or kg m-2].
  integer :: is, ie, js, je, nz
  integer :: i, j, k, kb_max
  integer :: l_seg
  real    :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(G))
  real    :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(G))
  logical :: local_open_u_BC, local_open_v_BC

  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
         "Module must be initialized before it is used.")
  if (.not. cs%calculate_Eady_growth_rate) return
  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
         "%SN_u is not associated with use_variable_mixing.")
  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
         "%SN_v is not associated with use_variable_mixing.")

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

  local_open_u_bc = .false.
  local_open_v_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_u_bc = obc%open_u_BCs_exist_globally
    local_open_v_bc = obc%open_v_BCs_exist_globally
  endif ; endif

  one_meter = 1.0 * gv%m_to_H
  h_neglect = gv%H_subroundoff
  h_cutoff = real(2*nz) * (GV%Angstrom_H + h_neglect)

  ! To set the length scale based on the deformation radius, use wave_speed to
  ! calculate the first-mode gravity wave speed and then blend the equatorial
  ! and midlatitude deformation radii, using calc_resoln_function as a template.

  !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
  do k=nz,cs%VarMix_Ktop,-1

    if (calculate_slopes) then
      ! Calculate the interface slopes E_x and E_y and u- and v- points respectively
      do j=js-1,je+1 ; do i=is-1,ie
        e_x(i,j) = us%Z_to_L*(e(i+1,j,k)-e(i,j,k))*g%IdxCu(i,j)
        ! Mask slopes where interface intersects topography
        if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        e_y(i,j) = us%Z_to_L*(e(i,j+1,k)-e(i,j,k))*g%IdyCv(i,j)
        ! Mask slopes where interface intersects topography
        if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
      enddo ; enddo
    else
      do j=js-1,je+1 ; do i=is-1,ie
        e_x(i,j) = cs%slope_x(i,j,k)
        if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
      enddo ; enddo
      do j=js-1,je ; do i=is-1,ie+1
        e_y(i,j) = cs%slope_y(i,j,k)
        if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
      enddo ; enddo
    endif

    ! Calculate N*S*h from this layer and add to the sum
    do j=js,je ; do i=is-1,ie
      s2 = ( e_x(i,j)**2  + 0.25*( &
            (e_y(i,j)**2+e_y(i+1,j-1)**2)+(e_y(i+1,j)**2+e_y(i,j-1)**2) ) )
      hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
      hup = 2.*h(i+1,j,k)*h(i+1,j,k-1) / (h(i+1,j,k) + h(i+1,j,k-1) + h_neglect)
      h_geom = sqrt(hdn*hup)
      n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
      if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < h_cutoff) &
        s2 = 0.0
      s2n2_u_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
    enddo ; enddo
    do j=js-1,je ; do i=is,ie
      s2 = ( e_y(i,j)**2  + 0.25*( &
            (e_x(i,j)**2+e_x(i-1,j+1)**2)+(e_x(i,j+1)**2+e_x(i-1,j)**2) ) )
      hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
      hup = 2.*h(i,j+1,k)*h(i,j+1,k-1) / (h(i,j+1,k) + h(i,j+1,k-1) + h_neglect)
      h_geom = sqrt(hdn*hup)
      n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
      if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < h_cutoff) &
        s2 = 0.0
      s2n2_v_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
    enddo ; enddo

  enddo ! k
  !$OMP parallel do default(shared)
  do j=js,je
    do i=is-1,ie ; cs%SN_u(i,j) = 0.0 ; enddo
    do k=nz,cs%VarMix_Ktop,-1 ; do i=is-1,ie
      cs%SN_u(i,j) = cs%SN_u(i,j) + s2n2_u_local(i,j,k)
    enddo ; enddo
    ! SN above contains S^2*N^2*H, convert to vertical average of S*N
    do i=is-1,ie
      !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom_Z ) ))
      !The code below behaves better than the line above. Not sure why? AJA
      if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
        cs%SN_u(i,j) = g%mask2dCu(i,j) * sqrt( cs%SN_u(i,j) / &
                                               (max(g%bathyT(i,j), g%bathyT(i+1,j))) )
      else
        cs%SN_u(i,j) = 0.0
      endif
      if (local_open_u_bc) then
        l_seg = obc%segnum_u(i,j)

        if (l_seg /= obc_none) then
          if (obc%segment(l_seg)%open) then
            cs%SN_u(i,j) = 0.
          endif
        endif
      endif
    enddo
  enddo
  !$OMP parallel do default(shared)
  do j=js-1,je
    do i=is,ie ; cs%SN_v(i,j) = 0.0 ; enddo
    do k=nz,cs%VarMix_Ktop,-1 ; do i=is,ie
      cs%SN_v(i,j) = cs%SN_v(i,j) + s2n2_v_local(i,j,k)
    enddo ; enddo
    do i=is,ie
      !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom_Z ) ))
      !The code below behaves better than the line above. Not sure why? AJA
      if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
        cs%SN_v(i,j) = g%mask2dCv(i,j) * sqrt( cs%SN_v(i,j) / &
                                               (max(g%bathyT(i,j), g%bathyT(i,j+1))) )
      else
        cs%SN_v(i,j) = 0.0
      endif
      if (local_open_v_bc) then
        l_seg = obc%segnum_v(i,j)

        if (l_seg /= obc_none) then
          if (obc%segment(obc%segnum_v(i,j))%open) then
            cs%SN_v(i,j) = 0.
          endif
        endif
      endif
    enddo
  enddo

end subroutine calc_slope_functions_using_just_e

!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
subroutine calc_qg_leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy)
  type(varmix_cs),                           pointer     :: cs !< Variable mixing coefficients
  type(ocean_grid_type),                     intent(in)  :: g  !< Ocean grid structure
  type(verticalgrid_type),                   intent(in)  :: gv !< The ocean's vertical grid structure.
  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]
  integer,                                   intent(in)  :: k  !< Layer for which to calculate vorticity magnitude
  real, dimension(SZIB_(G),SZJ_(G)),         intent(in)  :: div_xx_dx  !< x-derivative of horizontal divergence
                                                                 !! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
  real, dimension(SZI_(G),SZJB_(G)),         intent(in)  :: div_xx_dy  !< y-derivative of horizontal divergence
                                                                 !! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
  real, dimension(SZI_(G),SZJB_(G)),         intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
                                                                 !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
  real, dimension(SZIB_(G),SZJ_(G)),         intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
                                                                 !! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
  ! Local variables
  real, dimension(SZI_(G),SZJB_(G)) :: &
    dslopey_dz, & ! z-derivative of y-slope at v-points [Z-1 ~> m-1]
    h_at_v,     & ! Thickness at v-points [H ~> m or kg m-2]
    beta_v,     & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
    grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
    grad_div_mag_v     ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]

  real, dimension(SZIB_(G),SZJ_(G)) :: &
    dslopex_dz, & ! z-derivative of x-slope at u-points [Z-1 ~> m-1]
    h_at_u,     & ! Thickness at u-points [H ~> m or kg m-2]
    beta_u,     & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
    grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
    grad_div_mag_u     ! Magnitude of divergence gradient at u-points [T-1 L-1 ~> s-1 m-1]
  real :: h_at_slope_above ! The thickness above [H ~> m or kg m-2]
  real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2]
  real :: ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
  real :: f  ! A copy of the Coriolis parameter [T-1 ~> s-1]
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq,nz
  real :: inv_pi3

  is  = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  nz = g%ke

  inv_pi3 = 1.0/((4.0*atan(1.0))**3)

  if ((k > 1) .and. (k < nz)) then

    do j=js-1,je+1 ; do i=is-2,ieq+1
      h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / &
                         ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
                         + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + gv%H_subroundoff )
      h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / &
                         ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
                         + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + gv%H_subroundoff )
      ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
      dslopex_dz(i,j) = 2. * ( cs%slope_x(i,j,k) - cs%slope_x(i,j,k+1) ) * ih
      h_at_u(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
    enddo ; enddo

    do j=js-2,jeq+1 ; do i=is-1,ie+1
      h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / &
                         ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
                         + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + gv%H_subroundoff )
      h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / &
                         ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
                         + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + gv%H_subroundoff )
      ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
      dslopey_dz(i,j) = 2. * ( cs%slope_y(i,j,k) - cs%slope_y(i,j,k+1) ) * ih
      h_at_v(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
    enddo ; enddo

    do j=js-1,je ; do i=is-1,ieq+1
      f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
      vort_xy_dx(i,j) = vort_xy_dx(i,j) - f * us%L_to_Z * &
            ( ( h_at_u(i,j) * dslopex_dz(i,j) + h_at_u(i-1,j+1) * dslopex_dz(i-1,j+1) ) &
            + ( h_at_u(i-1,j) * dslopex_dz(i-1,j) + h_at_u(i,j+1) * dslopex_dz(i,j+1) ) ) / &
              ( ( h_at_u(i,j) + h_at_u(i-1,j+1) ) + ( h_at_u(i-1,j) + h_at_u(i,j+1) ) + gv%H_subroundoff)
    enddo ; enddo

    do j=js-1,jeq+1 ; do i=is-1,ie
      f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
      vort_xy_dy(i,j) = vort_xy_dy(i,j) - f * us%L_to_Z * &
            ( ( h_at_v(i,j) * dslopey_dz(i,j) + h_at_v(i+1,j-1) * dslopey_dz(i+1,j-1) ) &
            + ( h_at_v(i,j-1) * dslopey_dz(i,j-1) + h_at_v(i+1,j) * dslopey_dz(i+1,j) ) ) / &
              ( ( h_at_v(i,j) + h_at_v(i+1,j-1) ) + ( h_at_v(i,j-1) + h_at_v(i+1,j) ) + gv%H_subroundoff)
    enddo ; enddo
  endif ! k > 1

  if (cs%use_QG_Leith_GM) then

    do j=js,je ; do i=is-1,ieq
      grad_vort_mag_u(i,j) = sqrt(vort_xy_dy(i,j)**2  + (0.25*((vort_xy_dx(i,j) + vort_xy_dx(i+1,j-1)) &
                                                             + (vort_xy_dx(i+1,j) + vort_xy_dx(i,j-1))))**2)
      grad_div_mag_u(i,j) = sqrt(div_xx_dx(i,j)**2  + (0.25*((div_xx_dy(i,j) + div_xx_dy(i+1,j-1)) &
                                                           + (div_xx_dy(i+1,j) + div_xx_dy(i,j-1))))**2)
      if (cs%use_beta_in_QG_Leith) then
        beta_u(i,j) = sqrt((0.5*(g%dF_dx(i,j)+g%dF_dx(i+1,j))**2) + &
                          (0.5*(g%dF_dy(i,j)+g%dF_dy(i+1,j))**2))
        cs%KH_u_QG(i,j,k) = min(grad_vort_mag_u(i,j) + grad_div_mag_u(i,j), 3.0*beta_u(i,j)) * &
                            cs%Laplac3_const_u(i,j) * inv_pi3
      else
        cs%KH_u_QG(i,j,k) = (grad_vort_mag_u(i,j) + grad_div_mag_u(i,j)) * &
                            cs%Laplac3_const_u(i,j) * inv_pi3
      endif
    enddo ; enddo

    do j=js-1,jeq ; do i=is,ie
      grad_vort_mag_v(i,j) = sqrt(vort_xy_dx(i,j)**2  + (0.25*((vort_xy_dy(i,j) + vort_xy_dy(i-1,j+1)) &
                                                             + (vort_xy_dy(i,j+1) + vort_xy_dy(i-1,j))))**2)
      grad_div_mag_v(i,j) = sqrt(div_xx_dy(i,j)**2  + (0.25*((div_xx_dx(i,j) + div_xx_dx(i-1,j+1)) &
                                                           + (div_xx_dx(i,j+1) + div_xx_dx(i-1,j))))**2)
      if (cs%use_beta_in_QG_Leith) then
        beta_v(i,j) = sqrt((0.5*(g%dF_dx(i,j)+g%dF_dx(i,j+1))**2) + &
                          (0.5*(g%dF_dy(i,j)+g%dF_dy(i,j+1))**2))
        cs%KH_v_QG(i,j,k) = min(grad_vort_mag_v(i,j) + grad_div_mag_v(i,j), 3.0*beta_v(i,j)) * &
                            cs%Laplac3_const_v(i,j) * inv_pi3
      else
        cs%KH_v_QG(i,j,k) = (grad_vort_mag_v(i,j) + grad_div_mag_v(i,j)) * &
                            cs%Laplac3_const_v(i,j) * inv_pi3
      endif
    enddo ; enddo
    ! post diagnostics

    if (k==nz) then
      if (cs%id_KH_v_QG > 0)  call post_data(cs%id_KH_v_QG, cs%KH_v_QG, cs%diag)
      if (cs%id_KH_u_QG > 0)  call post_data(cs%id_KH_u_QG, cs%KH_u_QG, cs%diag)
    endif
  endif

end subroutine calc_qg_leith_viscosity

!> Initializes the variables mixing coefficients container
subroutine varmix_init(Time, G, GV, US, param_file, diag, 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   !< The ocean's 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(varmix_cs),               pointer :: cs   !< Variable mixing coefficients
  ! Local variables
  real :: khtr_slope_cff, khth_slope_cff, oneortwo
  real :: n2_filter_depth  ! A depth below which stratification is treated as monotonic when
                           ! calculating the first-mode wave speed [Z ~> m]
  real :: khtr_passivity_coeff
  real :: absurdly_small_freq  ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1].  The
             ! default value is roughly (pi / (the age of the universe)).
  logical :: gill_equatorial_ld, use_fgnv_streamfn, use_meke, in_use
  logical :: default_2018_answers, remap_answers_2018
  real :: mle_front_length
  real :: leith_lap_const      ! The non-dimensional coefficient in the Leith viscosity
  real :: grid_sp_u2, grid_sp_v2 ! Intermediate quantities for Leith metrics [L2 ~> m2]
  real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
  real :: wave_speed_min      ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
  real :: wave_speed_tol      ! The fractional tolerance for finding the wave speeds [nondim]
  logical :: better_speed_est ! If true, use a more robust estimate of the first
                              ! mode wave speed as the starting point for iterations.
! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
  integer :: is, ie, js, je, isq, ieq, jsq, jeq, i, j
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (associated(cs)) then
    call mom_error(warning, "VarMix_init called with an associated "// &
                             "control structure.")
    return
  endif

  allocate(cs)
  in_use = .false. ! Set to true to avoid deallocating
  cs%diag => diag ! Diagnostics pointer
  cs%calculate_cg1 = .false.
  cs%calculate_Rd_dx = .false.
  cs%calculate_res_fns = .false.
  cs%calculate_Eady_growth_rate = .false.
  cs%calculate_depth_fns = .false.
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "USE_VARIABLE_MIXING", cs%use_variable_mixing,&
                 "If true, the variable mixing code will be called.  This "//&
                 "allows diagnostics to be created even if the scheme is "//&
                 "not used.  If KHTR_SLOPE_CFF>0 or  KhTh_Slope_Cff>0, "//&
                 "this is set to true regardless of what is in the "//&
                 "parameter file.", default=.false.)
  call get_param(param_file, mdl, "USE_VISBECK", cs%use_Visbeck,&
                 "If true, use the Visbeck et al. (1997) formulation for \n"//&
                 "thickness diffusivity.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
                 "If true, the Laplacian lateral viscosity is scaled away "//&
                 "when the first baroclinic deformation radius is well "//&
                 "resolved.", default=.false.)
  call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", cs%Depth_scaled_KhTh, &
                 "If true, KHTH is scaled away when the depth is shallower"//&
                 "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//&
                 "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
                 "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
                 default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
                 "If true, the interface depth diffusivity is scaled away "//&
                 "when the first baroclinic deformation radius is well "//&
                 "resolved.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
                 "If true, the epipycnal tracer diffusivity is scaled "//&
                 "away when the first baroclinic deformation radius is "//&
                 "well resolved.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_USE_EBT", cs%Resoln_use_ebt, &
                 "If true, uses the equivalent barotropic wave speed instead "//&
                 "of first baroclinic wave for calculating the resolution fn.",&
                 default=.false.)
  call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
                 "If true, uses the equivalent barotropic structure "//&
                 "as the vertical structure of thickness diffusivity.",&
                 default=.false.)
  call get_param(param_file, mdl, "KHTH_SLOPE_CFF", 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, "KHTR_SLOPE_CFF", khtr_slope_cff, &
                 "The nondimensional coefficient in the Visbeck formula "//&
                 "for the epipycnal tracer diffusivity", units="nondim", &
                 default=0.0)
  call get_param(param_file, mdl, "USE_STORED_SLOPES", cs%use_stored_slopes,&
                 "If true, the isopycnal slopes are calculated once and "//&
                 "stored for re-use. This uses more memory but avoids calling "//&
                 "the equation of state more times than should be necessary.", &
                 default=.false.)
  call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
                 "A miniscule frequency that is used to avoid division by 0.  The default "//&
                 "value is roughly (pi / (the age of the universe)).", &
                 default=1.0e-17, units="s-1", scale=us%T_to_s)
  call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
                 default=.false., do_not_log=.true.)
  cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn
  call get_param(param_file, mdl, "USE_MEKE", use_meke, &
                 default=.false., do_not_log=.true.)
  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
  cs%calculate_Eady_growth_rate = cs%calculate_Eady_growth_rate .or. use_meke
  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
                 default=0., do_not_log=.true.)
  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
  call get_param(param_file, mdl, "MLE_FRONT_LENGTH", mle_front_length, &
                 default=0., do_not_log=.true.)
  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)

  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)


  if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct) then
    in_use = .true.
    call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
                 "The depth below which N2 is monotonized to avoid stratification "//&
                 "artifacts from altering the equivalent barotropic mode structure.",&
                 units="m", default=2000., scale=us%m_to_Z)
    allocate(cs%ebt_struct(isd:ied,jsd:jed,g%ke)) ; cs%ebt_struct(:,:,:) = 0.0
  endif

  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
    cs%calculate_Eady_growth_rate = .true.
    call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
          "If non-zero, is an upper bound on slopes used in the "//&
          "Visbeck formula for diffusivity. This does not affect the "//&
          "isopycnal slope calculation used within thickness diffusion.",  &
          units="nondim", default=0.0)
  endif

  if (cs%use_stored_slopes) then
    in_use = .true.
    allocate(cs%slope_x(isdb:iedb,jsd:jed,g%ke+1)) ; cs%slope_x(:,:,:) = 0.0
    allocate(cs%slope_y(isd:ied,jsdb:jedb,g%ke+1)) ; cs%slope_y(:,:,:) = 0.0
    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)
  endif

  if (cs%calculate_Eady_growth_rate) then
    in_use = .true.
    allocate(cs%SN_u(isdb:iedb,jsd:jed)) ; cs%SN_u(:,:) = 0.0
    allocate(cs%SN_v(isd:ied,jsdb:jedb)) ; cs%SN_v(:,:) = 0.0
    cs%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, time, &
       'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=us%s_to_T)
    cs%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, time, &
       'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=us%s_to_T)
    call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
                 "The layer number at which to start vertical integration "//&
                 "of S*N for purposes of finding the Eady growth rate.", &
                 units="nondim", default=2)
  endif

  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
    in_use = .true.
    call get_param(param_file, mdl, "VISBECK_L_SCALE", cs%Visbeck_L_scale, &
                 "The fixed length scale in the Visbeck formula.", units="m", &
                 default=0.0)
    allocate(cs%L2u(isdb:iedb,jsd:jed)) ; cs%L2u(:,:) = 0.0
    allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = 0.0
    if (cs%Visbeck_L_scale<0) then
      do j=js,je ; do i=is-1,ieq
        cs%L2u(i,j) = cs%Visbeck_L_scale**2 * g%areaCu(i,j)
      enddo; enddo
      do j=js-1,jeq ; do i=is,ie
        cs%L2v(i,j) = cs%Visbeck_L_scale**2 * g%areaCv(i,j)
      enddo; enddo
    else
      cs%L2u(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
      cs%L2v(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
    endif

    cs%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, time, &
       'Length scale squared for mixing coefficient, at u-points', &
       'm2', conversion=us%L_to_m**2)
    cs%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, time, &
       'Length scale squared for mixing coefficient, at v-points', &
       'm2', conversion=us%L_to_m**2)
  endif

  if (cs%calculate_Eady_growth_rate .and. cs%use_stored_slopes) then
    cs%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, time, &
         'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
         's-2', conversion=us%s_to_T**2)
    cs%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, time, &
         'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', &
         's-2', conversion=us%s_to_T**2)
  endif
  if (cs%use_stored_slopes) then
    cs%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, time, &
         'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 'nondim')
    cs%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, time, &
         'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 'nondim')
  endif

  oneortwo = 1.0
  if (cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. cs%Resoln_scaled_KhTr) then
    cs%calculate_Rd_dx = .true.
    cs%calculate_res_fns = .true.
    allocate(cs%Res_fn_h(isd:ied,jsd:jed))       ; cs%Res_fn_h(:,:) = 0.0
    allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb))   ; cs%Res_fn_q(:,:) = 0.0
    allocate(cs%Res_fn_u(isdb:iedb,jsd:jed))     ; cs%Res_fn_u(:,:) = 0.0
    allocate(cs%Res_fn_v(isd:ied,jsdb:jedb))     ; cs%Res_fn_v(:,:) = 0.0
    allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%beta_dx2_q(:,:) = 0.0
    allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed))   ; cs%beta_dx2_u(:,:) = 0.0
    allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb))   ; cs%beta_dx2_v(:,:) = 0.0
    allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb))   ; cs%f2_dx2_q(:,:) = 0.0
    allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed))     ; cs%f2_dx2_u(:,:) = 0.0
    allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb))     ; cs%f2_dx2_v(:,:) = 0.0

    cs%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, time, &
       'Resolution function for scaling diffusivities', 'nondim')

    call get_param(param_file, mdl, "KH_RES_SCALE_COEF", cs%Res_coef_khth, &
                 "A coefficient that determines how KhTh is scaled away if "//&
                 "RESOLN_SCALED_... is true, as "//&
                 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
                 units="nondim", default=1.0)
    call get_param(param_file, mdl, "KH_RES_FN_POWER", cs%Res_fn_power_khth, &
                 "The power of dx/Ld in the Kh resolution function.  Any "//&
                 "positive integer may be used, although even integers "//&
                 "are more efficient to calculate.  Setting this greater "//&
                 "than 100 results in a step-function being used.", &
                 units="nondim", default=2)
    call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
                 "A coefficient that determines how Kh is scaled away if "//&
                 "RESOLN_SCALED_... is true, as "//&
                 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
                 "This function affects lateral viscosity, Kh, and not KhTh.", &
                 units="nondim", default=cs%Res_coef_khth)
    call get_param(param_file, mdl, "VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
                 "The power of dx/Ld in the Kh resolution function.  Any "//&
                 "positive integer may be used, although even integers "//&
                 "are more efficient to calculate.  Setting this greater "//&
                 "than 100 results in a step-function being used. "//&
                 "This function affects lateral viscosity, Kh, and not KhTh.", &
                 units="nondim", default=cs%Res_fn_power_khth)
    call get_param(param_file, mdl, "INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
                 "If true, interpolate the resolution function to the "//&
                 "velocity points from the thickness points; otherwise "//&
                 "interpolate the wave speed and calculate the resolution "//&
                 "function independently at each point.", default=.false.)
    if (cs%interpolate_Res_fn) then
      if (cs%Res_coef_visc /= cs%Res_coef_khth) call mom_error(fatal, &
           "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
           "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
      if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth) call mom_error(fatal, &
           "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
           "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
    endif
    call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", gill_equatorial_ld, &
                 "If true, uses Gill's definition of the baroclinic "//&
                 "equatorial deformation radius, otherwise, if false, use "//&
                 "Pedlosky's definition. These definitions differ by a factor "//&
                 "of 2 in front of the beta term in the denominator. Gill's "//&
                 "is the more appropriate definition.", default=.true.)
    if (gill_equatorial_ld) then
      oneortwo = 2.0
    endif

    do j=js-1,jeq ; do i=is-1,ieq
      cs%f2_dx2_q(i,j) = (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * &
                         max(g%CoriolisBu(i,j)**2, absurdly_small_freq**2)
      cs%beta_dx2_q(i,j) = oneortwo * ((g%dxBu(i,j))**2 + (g%dyBu(i,j))**2) * (sqrt(0.5 * &
          ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) ) ))
    enddo ; enddo

    do j=js,je ; do i=is-1,ieq
      cs%f2_dx2_u(i,j) = (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * &
          max(0.5* (g%CoriolisBu(i,j)**2+g%CoriolisBu(i,j-1)**2), absurdly_small_freq**2)
      cs%beta_dx2_u(i,j) = oneortwo * ((g%dxCu(i,j))**2 + (g%dyCu(i,j))**2) * (sqrt( &
          0.25*( (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2 + &
                  ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
                 (((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2 + &
                  ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) ) + &
                  ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 ))
    enddo ; enddo

    do j=js-1,jeq ; do i=is,ie
      cs%f2_dx2_v(i,j) = ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * &
          max(0.5*(g%CoriolisBu(i,j)**2+g%CoriolisBu(i-1,j)**2), absurdly_small_freq**2)
      cs%beta_dx2_v(i,j) = oneortwo * ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * (sqrt( &
          ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
          0.25*( (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
                  ((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2) + &
                 (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2 + &
                  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
    enddo ; enddo

  endif

  if (cs%Depth_scaled_KhTh) then
    cs%calculate_depth_fns = .true.
    allocate(cs%Depth_fn_u(isdb:iedb,jsd:jed))     ; cs%Depth_fn_u(:,:) = 0.0
    allocate(cs%Depth_fn_v(isd:ied,jsdb:jedb))     ; cs%Depth_fn_v(:,:) = 0.0
    call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", cs%depth_scaled_khth_h0, &
    "The depth above which KHTH is scaled away.",&
    units="m", default=1000.)
    call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", cs%depth_scaled_khth_exp, &
    "The exponent used in the depth dependent scaling function for KHTH.",&
    units="nondim", default=3.0)
  endif

  ! Resolution %Rd_dx_h
  cs%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, time, &
       'Ratio between deformation radius and grid spacing', 'm m-1')
  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)

  if (cs%calculate_Rd_dx) then
    cs%calculate_cg1 = .true. ! We will need %cg1
    allocate(cs%Rd_dx_h(isd:ied,jsd:jed))   ; cs%Rd_dx_h(:,:) = 0.0
    allocate(cs%beta_dx2_h(isd:ied,jsd:jed)); cs%beta_dx2_h(:,:) = 0.0
    allocate(cs%f2_dx2_h(isd:ied,jsd:jed))  ; cs%f2_dx2_h(:,:) = 0.0
    do j=js-1,je+1 ; do i=is-1,ie+1
      cs%f2_dx2_h(i,j) = (g%dxT(i,j)**2 + g%dyT(i,j)**2) * &
          max(0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
                      (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)), &
              absurdly_small_freq**2)
      cs%beta_dx2_h(i,j) = oneortwo * ((g%dxT(i,j))**2 + (g%dyT(i,j))**2) * (sqrt(0.5 * &
          ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
    enddo ; enddo
  endif

  if (cs%calculate_cg1) then
    in_use = .true.
    allocate(cs%cg1(isd:ied,jsd:jed)) ; cs%cg1(:,:) = 0.0
    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, "REMAPPING_2018_ANSWERS", remap_answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
    call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
                 "The fractional tolerance for finding the wave speeds.", &
                 units="nondim", default=0.001)
    !### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
    call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
                 "A floor in the first mode speed below which 0 used instead.", &
                 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
    call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
                 "If true, use a more robust estimate of the first mode wave speed as the "//&
                 "starting point for iterations.", default=.false.) !### Change the default.
    call wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, &
                         mono_n2_depth=n2_filter_depth, remap_answers_2018=remap_answers_2018, &
                         better_speed_est=better_speed_est, min_speed=wave_speed_min, &
                         wave_speed_tol=wave_speed_tol)
  endif

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

  if (cs%Use_QG_Leith_GM) then
    call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
               "The nondimensional Laplacian Leith constant, \n"//&
               "often set to 1.0", units="nondim", default=0.0)

    call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
               "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
               default=.true.)

    alloc_(cs%Laplac3_const_u(isdb:iedb,jsd:jed)) ; cs%Laplac3_const_u(:,:) = 0.0
    alloc_(cs%Laplac3_const_v(isd:ied,jsdb:jedb)) ; cs%Laplac3_const_v(:,:) = 0.0
    alloc_(cs%KH_u_QG(isdb:iedb,jsd:jed,g%ke)) ; cs%KH_u_QG(:,:,:) = 0.0
    alloc_(cs%KH_v_QG(isd:ied,jsdb:jedb,g%ke)) ; cs%KH_v_QG(:,:,:) = 0.0
    ! register diagnostics

    cs%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, time, &
       'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
    cs%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, time, &
       'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

    do j=jsq,jeq+1 ; do i=is-1,ieq
      ! Static factors in the Leith schemes
      grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
      grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
      cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
    enddo ; enddo
    do j=js-1,jeq ; do i=isq,ieq+1
      ! Static factors in the Leith schemes
      grid_sp_v2 = g%dyCv(i,j)*g%dxCv(i,j)
      grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
      cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
    enddo ; enddo

    if (.not. cs%use_stored_slopes) call mom_error(fatal, &
           "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
           "USE_STORED_SLOPES must be True when using QG Leith.")
  endif

  ! If nothing is being stored in this class then deallocate
  if (in_use) then
    cs%use_variable_mixing = .true.
  else
    deallocate(cs)
    return
  endif

end subroutine varmix_init

!> \namespace mom_lateral_mixing_coeffs
!!
!! This module provides a container for various factors used in prescribing diffusivities, that are
!! a function of the state (in particular the stratification and isoneutral slopes).
!!
!! \section section_Resolution_Function The resolution function
!!
!! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
!! The square of the resolution parameter is
!!
!! \f[
!! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
!! \f]
!!
!! where the grid spacing is calculated as
!!
!! \f[
!! \Delta^2 = \Delta x^2 + \Delta y^2 .
!! \f]
!!
!! \todo Check this reference to Bob on/off paper.
!! The resolution function used in scaling diffusivities (Hallberg, 2010) is
!!
!! \f[
!! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
!! \f]
!!
!! The resolution function can be applied independently to thickness diffusion (module mom_thickness_diffuse),
!! tracer diffusion (mom_tracer_hordiff) lateral viscosity (mom_hor_visc).
!!
!! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
!! Ocean Modelling, 71, pp 92-103.  http://dx.doi.org/10.1016/j.ocemod.2013.08.007
!!
!! | Symbol                | Module parameter |
!! | ------                | --------------- |
!! | -                     | <code>USE_VARIABLE_MIXING</code> |
!! | -                     | <code>RESOLN_SCALED_KH</code> |
!! | -                     | <code>RESOLN_SCALED_KHTH</code> |
!! | -                     | <code>RESOLN_SCALED_KHTR</code> |
!! | \f$ \alpha \f$        | <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) |
!! | \f$ p \f$             | <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) |
!! | \f$ \alpha \f$        | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) |
!! | \f$ p \f$             | <code>VISC_RES_FN_POWER</code> (for lateral viscosity) |
!! | -                     | <code>GILL_EQUATORIAL_LD</code> |
!!
!!
!!
!! \section section_Vicbeck Visbeck diffusivity
!!
!! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
!! scheme.  The factors are combined in mom_thickness_diffuse::thickness_diffuse() but calculated in this module.
!!
!! \f[
!! \kappa_h = \alpha_s L_s^2 S N
!! \f]
!!
!! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
!!
!! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
!! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
!!
!! | Symbol                | Module parameter |
!! | ------                | --------------- |
!! | -                     | <code>USE_VARIABLE_MIXING</code> |
!! | \f$ \alpha_s \f$      | <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)|
!! | \f$ \alpha_s \f$      | <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)|
!! | \f$ L_{s} \f$         | <code>VISBECK_L_SCALE</code> |
!! | \f$ S_{max} \f$       | <code>VISBECK_MAX_SLOPE</code> |
!!
!!
!! \section section_vertical_structure_khth Vertical structure function for KhTh
!!
!! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
!! velocity mode.  The structure function is stored in the control structure for thie module (varmix_cs) but is
!! calculated using subroutines in mom_wave_speed.
!!
!! | Symbol                | Module parameter |
!! | ------                | --------------- |
!! | -                     | <code>KHTH_USE_EBT_STRUCT</code> |

end module mom_lateral_mixing_coeffs