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 use 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_resoln_function (h, tv, G, GV, CS)
	Calculates and stores the non-dimensional resolution functions. More...
subroutine, public	calc_slope_functions (h, tv, dt, G, GV, CS)
	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, e, slope_x, slope_y, N2_u, N2_v, G, GV, CS)
	Calculates factors used when setting diffusivity coefficients similar to Visbeck et al. More...
subroutine	calc_slope_functions_using_just_e (h, G, GV, CS, e, calculate_slopes)
	The original calc_slope_function() that calculated slopes using interface positions only, not accounting for density variations. More...
subroutine, public	varmix_init (Time, G, param_file, diag, CS)
	Initializes the variables mixing coefficients container. More...

Function/Subroutine Documentation

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(varmix_cs), pointer	CS
	)

Calculates and stores the non-dimensional resolution functions.

Parameters

[in,out]	g	Ocean grid structure
[in]	h	Layer thickness (m or kg/m2)
[in]	tv	Thermodynamic variables
[in]	gv	Vertical grid structure
	cs	Variable mixing coefficients

Definition at line 133 of file MOM_lateral_mixing_coeffs.F90.

References mom_error_handler::mom_error(), mom_diag_mediator::query_averaging_enabled(), and mom_wave_speed::wave_speed().

Referenced by mom::step_mom(), and mom::step_offline().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

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(varmix_cs), pointer	CS
	)

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,out]	h	Layer thickness (m or kg/m2)
[in]	tv	Thermodynamic variables
[in]	dt	Time increment (s)
	cs	Variable mixing coefficients

Definition at line 380 of file MOM_lateral_mixing_coeffs.F90.

References mom_isopycnal_slopes::calc_isoneutral_slopes(), calc_slope_functions_using_just_e(), calc_visbeck_coeffs(), mom_error_handler::mom_error(), and mom_diag_mediator::query_averaging_enabled().

  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

Here is the call graph for this function:

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

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 (m or kg/m2)
[in]	gv	Vertical grid structure
	cs	Variable mixing coefficients
[in]	e	Interface position (m)
[in]	calculate_slopes	If true, calculate slopes internally otherwise use slopes stored in CS

Definition at line 585 of file MOM_lateral_mixing_coeffs.F90.

References mom_error_handler::mom_error().

Referenced by calc_slope_functions().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

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 %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)  e, 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(varmix_cs), pointer  CS  )

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 (m or kg/m2)
[in]	e	Interface position (m)
[in]	slope_x	Zonal isoneutral slope
[in]	n2_u	Brunt-Vaisala frequency at u-points (1/s2)
[in]	slope_y	Meridional isoneutral slope
[in]	n2_v	Brunt-Vaisala frequency at v-points (1/s2)
	cs	Variable mixing coefficients

Definition at line 421 of file MOM_lateral_mixing_coeffs.F90.

References mom_error_handler::mom_error(), and mom_diag_mediator::query_averaging_enabled().

Referenced by calc_slope_functions().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

varmix_init()

subroutine, public mom_lateral_mixing_coeffs::varmix_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		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]	param_file	Parameter file handles
[in,out]	diag	Diagnostics control structure
	cs	Variable mixing coefficients

Definition at line 724 of file MOM_lateral_mixing_coeffs.F90.

References mom_error_handler::mom_error(), and mom_wave_speed::wave_speed_init().

  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

Here is the call graph for this function: