MOM_lateral_mixing_coeffs.F90

Go to the documentation of this file.

!> 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
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_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_wave_speed, only : wave_speed, wave_speed_cs, wave_speed_init

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 :: 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 :: use_visbeck_slope_bug !< If true, then retain a legacy bug in the calculation of weights
                                   !! applied to isoneutral slopes. There was an erroneous k-indexing
                                   !! for layer thicknesses. In addition, masking at coastlines was not
                                   !! used which introduced potential restart issues.  This flag will be
                                   !! deprecated in a future release.
  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_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 (s^-1)
    sn_v => null(), &  !< S*N at v-points (s^-1)
    l2u => null(), &   !< Length scale^2 at u-points (m^2)
    l2v => null(), &   !< Length scale^2 at v-points (m^2)
    cg1 => null(), &   !< The first baroclinic gravity wave speed in 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.
    res_fn_q => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at q points.
    res_fn_u => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at u points.
    res_fn_v => null(), & !< Non-dimensional function of the ratio the first baroclinic
                          !! deformation radius to the grid spacing at v points.
    beta_dx2_h => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at h points.
    beta_dx2_q => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at q points.
    beta_dx2_u => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at u points.
    beta_dx2_v => null(), & !< The magnitude of the gradient of the Coriolis parameter
                            !! times the grid spacing squared at v points.
    f2_dx2_h => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at h, in m2 s-2.
    f2_dx2_q => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at q, in m2 s-2.
    f2_dx2_u => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at u, in m2 s-2.
    f2_dx2_v => null(), & !< The Coriolis parameter squared times the grid
                          !! spacing squared at v, in m2 s-2.
    rd_dx_h => null()     !< Deformation radius over grid spacing (non-dim.)

  real, dimension(:,:,:), pointer :: &
    slope_x => null(), &  !< Zonal isopycnal slope (non-dimensional)
    slope_y => null(), &  !< Meridional isopycnal slope (non-dimensional)
    ebt_struct => null()  !< Vertical structure function to scale diffusivities with (non-dim)

  ! Parameters
  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 :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers (m2/s)
  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).

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

contains

!> Calculates and stores the non-dimensional resolution functions
subroutine calc_resoln_function(h, tv, G, GV, CS)
  type(ocean_grid_type),                    intent(inout) :: G  !< Ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h  !< Layer thickness (m or kg/m2)
  type(thermo_var_ptrs),                    intent(in)    :: tv !< Thermodynamic variables
  type(verticalgrid_type),                  intent(in)    :: GV !< Vertical grid structure
  type(varmix_cs),                          pointer       :: CS !< Variable mixing coefficients
  ! Local variables
  real :: cg1_q  ! The gravity wave speed interpolated to q points, in m s-1.
  real :: cg1_u  ! The gravity wave speed interpolated to u points, in m s-1.
  real :: cg1_v  ! The gravity wave speed interpolated to v points, in m s-1.
  real :: dx_term
  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, 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, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct, use_ebt_mode=.true.)
        call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp)
      endif
      call pass_var(cs%ebt_struct, g%Domain)
    else
      call wave_speed(h, tv, g, gv, 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 (non-dimensional).
  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 default(none) shared(is,ie,js,je,CS)
!$OMP do
    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
!$OMP end parallel
    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) &
!$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 = (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 * 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 = (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 * 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 = (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 * 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 = (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 * 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 = (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 * 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 = (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 * 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 = (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 * 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 = (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 * 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, CS)
  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(inout) :: h  !< Layer thickness (m or kg/m2)
  type(thermo_var_ptrs),                    intent(in)    :: tv !< Thermodynamic variables
  real,                                     intent(in)    :: dt !< Time increment (s)
  type(varmix_cs),                          pointer       :: CS !< Variable mixing coefficients
  ! Local variables
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    e             ! The interface heights relative to mean sea level, in m.
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: N2_v ! Square of Brunt-Vaisala freq at u-points

  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, gv%g_Earth, g, gv, e, halo_size=2)
    if (cs%use_stored_slopes) then
      call calc_isoneutral_slopes(g, gv, h, e, tv, dt*cs%kappa_smooth, &
                                  cs%slope_x, cs%slope_y, n2_u, n2_v, 1)
      call calc_visbeck_coeffs(h, e, cs%slope_x, cs%slope_y, n2_u, n2_v, g, gv, cs)
!     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, cs, e, .true.)
    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%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

end subroutine calc_slope_functions

!> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al.
subroutine calc_visbeck_coeffs(h, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS)
  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 (m or kg/m2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1),  intent(in)    :: e  !< Interface position (m)
  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    !< Brunt-Vaisala frequency at u-points (1/s2)
  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    !< Brunt-Vaisala frequency at v-points (1/s2)
  type(varmix_cs),                             pointer       :: CS !< Variable mixing coefficients
  ! Local variables
  real :: E_x(szib_(g), szj_(g))  ! X-slope of interface at u points (for diagnostics)
  real :: E_y(szi_(g), szjb_(g))  ! Y-slope of interface at u points (for diagnostics)
  real :: Khth_Loc      ! Locally calculated thickness mixing coefficient (m2/s)
  real :: S2            ! Interface slope squared (non-dim)
  real :: N2            ! Brunt-Vaisala frequency (1/s)
  real :: Hup, Hdn      ! Thickness from above, below (m or kg m-2)
  real :: H_geom        ! The geometric mean of Hup*Hdn, in m or kg m-2.
  integer :: is, ie, js, je, nz
  integer :: i, j, k, kb_max
  real :: S2max, wNE, wSE, wSW, wNW
  real :: SN_u_local(szib_(g), szj_(g),szk_(g))
  real :: SN_v_local(szi_(g), szjb_(g),szk_(g))
  real :: H_u(szib_(g)), H_v(szi_(g))
  real :: S2_u(szib_(g), szj_(g))
  real :: S2_v(szi_(g), szjb_(g))

  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

  s2max = cs%Visbeck_S_max**2

!$OMP parallel default(none) shared(is,ie,js,je,CS,nz,e,G,GV,h, &
!$OMP                               S2_u,S2_v,slope_x,slope_y,   &
!$OMP                               SN_u_local,SN_v_local,N2_u,N2_v, S2max)   &
!$OMP                       private(E_x,E_y,S2,H_u,H_v,Hdn,Hup,H_geom,N2, &
!$OMP                       wNE, wSE, wSW, wNW)
!$OMP do
  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 do
  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
      if (cs%use_Visbeck_slope_bug) then
        wse = h(i+1,j,k)*h(i+1,j-1,k) * h(i+1,j,k)*h(i+1,j-1,k-1)
        wnw = h(i  ,j,k)*h(i  ,j+1,k) * h(i  ,j,k)*h(i  ,j+1,k-1)
        wne = h(i+1,j,k)*h(i+1,j+1,k) * h(i+1,j,k)*h(i+1,j+1,k-1)
        wsw = h(i  ,j,k)*h(i  ,j-1,k) * h(i  ,j,k)*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**2 ) !### This should be **4 for consistent units.
      else
        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 )
      endif
      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
    enddo
  enddo

!$OMP do
  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
      if (cs%use_Visbeck_slope_bug) then
        wse = h(i,j  ,k)*h(i+1,j  ,k) * h(i,j  ,k)*h(i+1,j  ,k-1)
        wnw = h(i,j+1,k)*h(i-1,j+1,k) * h(i,j+1,k)*h(i-1,j+1,k-1)
        wne = h(i,j+1,k)*h(i+1,j+1,k) * h(i,j+1,k)*h(i+1,j+1,k-1)
        wsw = h(i,j  ,k)*h(i-1,j  ,k) * h(i,j  ,k)*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**2 ) !### This should be **4 for consistent units.
      else
        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 ) !### This should be **4 for consistent units.
      endif
      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
    enddo
  enddo

!$OMP end parallel

! 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)
    call uvchksum("calc_Visbeck_coeffs SN_[uv]", cs%SN_u, cs%SN_v, g%HI)
  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, CS, e, calculate_slopes)
  type(ocean_grid_type),                      intent(inout) :: G  !< Ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: h  !< Layer thickness (m or kg/m2)
  type(verticalgrid_type),                    intent(in)    :: GV !< Vertical grid structure
  type(varmix_cs),                            pointer       :: CS !< Variable mixing coefficients
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: e  !< Interface position (m)
  logical,                                    intent(in)    :: calculate_slopes !< If true, calculate slopes internally
                                                                  !! otherwise use slopes stored in CS
  ! Local variables
  real :: E_x(szib_(g), szj_(g))  ! X-slope of interface at u points (for diagnostics)
  real :: E_y(szi_(g), szjb_(g))  ! Y-slope of interface at u points (for diagnostics)
  real :: Khth_Loc      ! Locally calculated thickness mixing coefficient (m2/s)
  real :: H_cutoff      ! Local estimate of a minimum thickness for masking (m)
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected, in H.
  real :: S2            ! Interface slope squared (non-dim)
  real :: N2            ! Brunt-Vaisala frequency (1/s)
  real :: Hup, Hdn      ! Thickness from above, below (m or kg m-2)
  real :: H_geom        ! The geometric mean of Hup*Hdn, in m or kg m-2.
  real :: one_meter     ! One meter in thickness units of m or kg m-2.
  integer :: is, ie, js, je, nz
  integer :: i, j, k, kb_max
  real    :: SN_u_local(szib_(g), szj_(g),szk_(g))
  real    :: SN_v_local(szi_(g), szjb_(g),szk_(g))

  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

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

!$OMP parallel default(none) shared(is,ie,js,je,CS,nz,e,G,GV,h,H_cutoff,h_neglect, &
!$OMP                               one_meter,SN_u_local,SN_v_local,calculate_slopes)   &
!$OMP                       private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
!$OMP do
  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 do
  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) = (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) = (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) / (gv%H_to_m * 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
      sn_u_local(i,j,k) = (h_geom * gv%H_to_m) * 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) / (gv%H_to_m * 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
      sn_v_local(i,j,k) = (h_geom * gv%H_to_m) * s2 * n2
    enddo ; enddo

  enddo ! k
!$OMP do
  do j = js,je;
    do k=nz,cs%VarMix_Ktop,-1 ; do i=is-1,ie
      cs%SN_u(i,j) = cs%SN_u(i,j) + sn_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 ) )
      !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 ) 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
    enddo
  enddo
!$OMP do
  do j=js-1,je
    do k=nz,cs%VarMix_Ktop,-1 ; do i=is,ie
      cs%SN_v(i,j) = cs%SN_v(i,j) + sn_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 ) )
      !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 ) 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
    enddo
  enddo
!$OMP end parallel

end subroutine calc_slope_functions_using_just_e

!> Initializes the variables mixing coefficients container
subroutine varmix_init(Time, G, param_file, diag, CS)
  type(time_type),            intent(in) :: Time !< Current model time
  type(ocean_grid_type),      intent(in) :: G    !< Ocean grid structure
  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, N2_filter_depth
  real :: KhTr_passivity_coeff
  real, parameter :: absurdly_small_freq2 = 1e-34  ! A miniscule frequency
             ! squared that is used to avoid division by 0, in s-2.  This
             ! value is roughly (pi / (the age of the universe) )^2.
  logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
  real :: MLE_front_length
! 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.

  ! 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 \n"//&
                 "allows diagnostics to be created even if the scheme is \n"//&
                 "not used.  If KHTR_SLOPE_CFF>0 or  KhTh_Slope_Cff>0, \n"//&
                 "this is set to true regardless of what is in the \n"//&
                 "parameter file.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
                 "If true, the Laplacian lateral viscosity is scaled away \n"//&
                 "when the first baroclinic deformation radius is well \n"//&
                 "resolved.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
                 "If true, the interface depth diffusivity is scaled away \n"//&
                 "when the first baroclinic deformation radius is well \n"//&
                 "resolved.", default=.false.)
  call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
                 "If true, the epipycnal tracer diffusivity is scaled \n"//&
                 "away when the first baroclinic deformation radius is \n"//&
                 "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\n"//&
                 "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\n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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\n"//&
                 "stored for re-use. This uses more memory but avoids calling\n"//&
                 "the equation of state more times than should be necessary.", &
                 default=.false.)
  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_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\n"//&
                 "artifacts from altering the equivalent barotropic mode structure.",&
                 units='m', default=2000.)
    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\n"//       &
          "Visbeck formula for diffusivity. This does not affect the\n"//  &
          "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 \n"//&
                 "more sensible values of T & S into thin layers.", &
                 default=1.0e-6)
  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')
    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')
    call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
                 "The layer number at which to start vertical integration \n"//&
                 "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(:,:) = cs%Visbeck_L_scale**2
    allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = cs%Visbeck_L_scale**2

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

  if (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')
    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')
    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.', 's^-2')
    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.', 's^-2')
  endif

  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 \n"//&
                 "RESOLN_SCALED_... is true, as \n"//&
                 "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 \n"//&
                 "positive integer may be used, although even integers \n"//&
                 "are more efficient to calculate.  Setting this greater \n"//&
                 "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 \n"//&
                 "RESOLN_SCALED_... is true, as \n"//&
                 "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).\n"//&
                 "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 \n"//&
                 "positive integer may be used, although even integers \n"//&
                 "are more efficient to calculate.  Setting this greater \n"//&
                 "than 100 results in a step-function being used.\n"//&
                 "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 \n"//&
                 "velocity points from the thickness points; otherwise \n"//&
                 "interpolate the wave speed and calculate the resolution \n"//&
                 "function independently at each point.", default=.true.)
    call get_param(param_file, mdl, "USE_VISBECK_SLOPE_BUG", cs%use_Visbeck_slope_bug, &
                 "If true, then retain a legacy bug in the calculation of weights \n"//&
                 "applied to isoneutral slopes. There was an erroneous k-indexing \n"//&
                 "for layer thicknesses. In addition, masking at coastlines was not \n"//&
                 "used which introduced potential restart issues.  This flag will be \n"//&
                 "deprecated in a future release.", default=.false.)
    if (cs%interpolate_Res_fn) then
      if (cs%Res_coef_visc .ne. 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 .ne. 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\n"//&
                 "equatorial deformation radius, otherwise, if false, use\n"//&
                 "Pedlosky's definition. These definitions differ by a factor\n"//&
                 "of 2 infront of the beta term in the denominator. Gill's"//&
                 "is the more appropriate definition.\n", default=.false.)
    if (gill_equatorial_ld) then
      oneortwo = 2.0
    else
      oneortwo = 1.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_freq2)
      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_freq2)
      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_freq2)
      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

  ! 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', 'Nondim')
  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_freq2)
      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 wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, mono_n2_depth=n2_filter_depth)
  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 use subroutines in
!! mom_wave_speed.
!!
!! | Symbol                | Module parameter |
!! | ------                | --------------- |
!! | -                     | <code>KHTH_USE_EBT_STRUCT</code> |

end module mom_lateral_mixing_coeffs