mom_lateral_mixing_coeffs Module Reference

Detailed Description

Variable mixing coefficients.

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

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

\[ R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 } \]

where the grid spacing is calculated as

\[ \Delta^2 = \Delta x^2 + \Delta y^2 . \]

Todo:

Check this reference to Bob on/off paper. The resolution function used in scaling diffusivities (Hallberg, 2010) is

\[ r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p} \]

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
-	USE_VARIABLE_MIXING
-	RESOLN_SCALED_KH
-	RESOLN_SCALED_KHTH
-	RESOLN_SCALED_KHTR
\( \alpha \)	KH_RES_SCALE_COEF (for thickness and tracer diffusivity)
\( p \)	KH_RES_FN_POWER (for thickness and tracer diffusivity)
\( \alpha \)	VISC_RES_SCALE_COEF (for lateral viscosity)
\( p \)	VISC_RES_FN_POWER (for lateral viscosity)
-	GILL_EQUATORIAL_LD

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.

\[ \kappa_h = \alpha_s L_s^2 S N \]

where \(S\) is the magnitude of the isoneutral slope and \(N\) 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
-	USE_VARIABLE_MIXING
\( \alpha_s \)	KHTH_SLOPE_CFF (for mom_thickness_diffuse module)
\( \alpha_s \)	KHTR_SLOPE_CFF (for mom_tracer_hordiff module)
\( L_{s} \)	VISBECK_L_SCALE
\( S_{max} \)	VISBECK_MAX_SLOPE

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

Data Types
type	varmix_cs
	Variable mixing coefficients. More...

Functions/Subroutines
subroutine, public	calc_depth_function (G, CS)
	Calculates the non-dimensional depth functions. More...
subroutine, public	calc_resoln_function (h, tv, G, GV, US, CS)
	Calculates and stores the non-dimensional resolution functions. More...
subroutine, public	calc_slope_functions (h, tv, dt, G, GV, US, CS, OBC)
	Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al. style scaling of diffusivity. More...
subroutine	calc_visbeck_coeffs (h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
	Calculates factors used when setting diffusivity coefficients similar to Visbeck et al. More...
subroutine	calc_slope_functions_using_just_e (h, G, GV, US, CS, e, calculate_slopes, OBC)
	The original calc_slope_function() that calculated slopes using interface positions only, not accounting for density variations. More...
subroutine, public	calc_qg_leith_viscosity (CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy)
	Calculates the Leith Laplacian and bi-harmonic viscosity coefficients. More...
subroutine, public	varmix_init (Time, G, GV, US, param_file, diag, CS)
	Initializes the variables mixing coefficients container. More...

Function/Subroutine Documentation

calc_depth_function()

subroutine, public mom_lateral_mixing_coeffs::calc_depth_function	(	type(ocean_grid_type), intent(in)	G,
		type(varmix_cs), pointer	CS
	)

Calculates the non-dimensional depth functions.

Parameters

[in]	g	Ocean grid structure
	cs	Variable mixing coefficients

Definition at line 156 of file MOM_lateral_mixing_coeffs.F90.

  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
 

calc_qg_leith_viscosity()

subroutine, public mom_lateral_mixing_coeffs::calc_qg_leith_viscosity	(	type(varmix_cs), pointer	CS,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h,
		integer, intent(in)	k,
		real, dimension(szib_(g),szj_(g)), intent(in)	div_xx_dx,
		real, dimension(szi_(g),szjb_(g)), intent(in)	div_xx_dy,
		real, dimension(szi_(g),szjb_(g)), intent(inout)	vort_xy_dx,
		real, dimension(szib_(g),szj_(g)), intent(inout)	vort_xy_dy
	)

Calculates the Leith Laplacian and bi-harmonic viscosity coefficients.

Parameters

	cs	Variable mixing coefficients
[in]	g	Ocean grid structure
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in,out]	h	Layer thickness [H ~> m or kg m-2]
[in]	k	Layer for which to calculate vorticity magnitude
[in]	div_xx_dx	x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
[in]	div_xx_dy	y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
[in,out]	vort_xy_dx	x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
[in,out]	vort_xy_dy	y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]

Definition at line 807 of file MOM_lateral_mixing_coeffs.F90.

  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
 

calc_resoln_function()

subroutine, public mom_lateral_mixing_coeffs::calc_resoln_function	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(varmix_cs), pointer	CS
	)

Calculates and stores the non-dimensional resolution functions.

Parameters

[in,out]	g	Ocean grid structure
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	tv	Thermodynamic variables
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
	cs	Variable mixing coefficients

Definition at line 190 of file MOM_lateral_mixing_coeffs.F90.

  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
 

calc_slope_functions()

subroutine, public mom_lateral_mixing_coeffs::calc_slope_functions	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(varmix_cs), pointer	CS,
		type(ocean_obc_type), optional, pointer	OBC
	)

Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al. style scaling of diffusivity.

Parameters

[in,out]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	h	Layer thickness [H ~> m or kg m-2]
[in]	tv	Thermodynamic variables
[in]	dt	Time increment [T ~> s]
	cs	Variable mixing coefficients
	obc	Open boundaries control structure.

Definition at line 436 of file MOM_lateral_mixing_coeffs.F90.

  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
 

calc_slope_functions_using_just_e()

subroutine mom_lateral_mixing_coeffs::calc_slope_functions_using_just_e ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  h, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(varmix_cs), pointer  CS, real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)  e, logical, intent(in)  calculate_slopes, type(ocean_obc_type), optional, pointer  OBC  )

private

The original calc_slope_function() that calculated slopes using interface positions only, not accounting for density variations.

Parameters

[in,out]	g	Ocean grid structure
[in,out]	h	Layer thickness [H ~> m or kg m-2]
[in]	gv	Vertical grid structure
[in]	us	A dimensional unit scaling type
	cs	Variable mixing coefficients
[in]	e	Interface position [Z ~> m]
[in]	calculate_slopes	If true, calculate slopes internally otherwise use slopes stored in CS
	obc	Open boundaries control structure.

Definition at line 646 of file MOM_lateral_mixing_coeffs.F90.

  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
 

calc_visbeck_coeffs()

subroutine mom_lateral_mixing_coeffs::calc_visbeck_coeffs ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke+1), intent(in)  slope_x, real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke+1), intent(in)  slope_y, real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke+1), intent(in)  N2_u, real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke+1), intent(in)  N2_v, type(ocean_grid_type), intent(inout)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(varmix_cs), pointer  CS, type(ocean_obc_type), optional, pointer  OBC  )

private

Calculates factors used when setting diffusivity coefficients similar to Visbeck et al.

Parameters

[in,out]	g	Ocean grid structure
[in]	gv	Vertical grid structure
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	slope_x	Zonal isoneutral slope
[in]	n2_u	Buoyancy (Brunt-Vaisala) frequency at u-points [T-2 ~> s-2]
[in]	slope_y	Meridional isoneutral slope
[in]	n2_v	Buoyancy (Brunt-Vaisala) frequency at v-points [T-2 ~> s-2]
[in]	us	A dimensional unit scaling type
	cs	Variable mixing coefficients
	obc	Open boundaries control structure.

Definition at line 481 of file MOM_lateral_mixing_coeffs.F90.

  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
 

varmix_init()

subroutine, public mom_lateral_mixing_coeffs::varmix_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(varmix_cs), pointer	CS
	)

Initializes the variables mixing coefficients container.

Parameters

[in]	time	Current model time
[in]	g	Ocean grid structure
[in]	gv	The ocean's vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	param_file	Parameter file handles
[in,out]	diag	Diagnostics control structure
	cs	Variable mixing coefficients

Definition at line 936 of file MOM_lateral_mixing_coeffs.F90.

  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