MOM6
Data Types | Functions/Subroutines
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)
 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)
 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)
 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]gOcean grid structure
csVariable mixing coefficients

Definition at line 159 of file MOM_lateral_mixing_coeffs.F90.

159  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
160  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
161 
162  ! Local variables
163  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq
164  integer :: i, j
165  real :: H0 ! local variable for reference depth
166  real :: expo ! exponent used in the depth dependent scaling
167  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
168  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
169 
170  if (.not. associated(cs)) call mom_error(fatal, "calc_depth_function:"// &
171  "Module must be initialized before it is used.")
172  if (.not. cs%calculate_depth_fns) return
173  if (.not. associated(cs%Depth_fn_u)) call mom_error(fatal, &
174  "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
175  if (.not. associated(cs%Depth_fn_v)) call mom_error(fatal, &
176  "calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.")
177 
178  h0 = cs%depth_scaled_khth_h0
179  expo = cs%depth_scaled_khth_exp
180 !$OMP do
181  do j=js,je ; do i=is-1,ieq
182  cs%Depth_fn_u(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))/h0))**expo
183  enddo ; enddo
184 !$OMP do
185  do j=js-1,jeq ; do i=is,ie
186  cs%Depth_fn_v(i,j) = (min(1.0, 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))/h0))**expo
187  enddo ; enddo
188 

References mom_error_handler::mom_error().

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

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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
csVariable mixing coefficients
[in]gOcean grid structure
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in,out]hLayer thickness [H ~> m or kg m-2]
[in]kLayer for which to calculate vorticity magnitude
[in]div_xx_dxx-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
[in]div_xx_dyy-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
[in,out]vort_xy_dxx-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
[in,out]vort_xy_dyy-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]

Definition at line 751 of file MOM_lateral_mixing_coeffs.F90.

751  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
752  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
753  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
754  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
755 ! real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal flow [L T-1 ~> m s-1]
756 ! real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional flow [L T-1 ~> m s-1]
757  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
758  integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
759  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence
760  !! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
761  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence
762  !! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
763  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
764  !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
765  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
766  !! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
767 ! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Kh_h !< Leith Laplacian viscosity
768  !! at h-points [L2 T-1 ~> m2 s-1]
769 ! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Kh_q !< Leith Laplacian viscosity
770  !! at q-points [L2 T-1 ~> m2 s-1]
771 ! real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Leith_Ah_h !< Leith bi-harmonic viscosity
772  !! at h-points [L4 T-1 ~> m4 s-1]
773 ! real, dimension(SZIB_(G),SZJB_(G)), intent(out) :: Leith_Ah_q !< Leith bi-harmonic viscosity
774  !! at q-points [L4 T-1 ~> m4 s-1]
775 
776  ! Local variables
777  real, dimension(SZI_(G),SZJB_(G)) :: &
778  dslopey_dz, & ! z-derivative of y-slope at v-points [Z-1 ~> m-1]
779  h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
780  beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
781  grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
782  grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]
783 
784  real, dimension(SZIB_(G),SZJ_(G)) :: &
785  dslopex_dz, & ! z-derivative of x-slope at u-points [Z-1 ~> m-1]
786  h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
787  beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
788  grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
789  grad_div_mag_u ! Magnitude of divergence gradient at u-points [T-1 L-1 ~> s-1 m-1]
790  real :: h_at_slope_above ! The thickness above [H ~> m or kg m-2]
791  real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2]
792  real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
793  real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1]
794  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq,nz
795  real :: inv_PI3
796 
797  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
798  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
799  nz = g%ke
800 
801  inv_pi3 = 1.0/((4.0*atan(1.0))**3)
802 
803  !### I believe this halo update to be unnecessary. -RWH
804  call pass_var(h, g%Domain)
805 
806  if ((k > 1) .and. (k < nz)) then
807 
808  ! Add in stretching term for the QG Leith vsicosity
809 ! if (CS%use_QG_Leith) then
810 
811  !### do j=js-1,je+1 ; do I=is-2,Ieq+1
812  do j=js-2,jeq+2 ; do i=is-2,ieq+1
813  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) ) / &
814  ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
815  + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + gv%H_subroundoff )
816  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) ) / &
817  ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
818  + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + gv%H_subroundoff )
819  ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
820  dslopex_dz(i,j) = 2. * ( cs%slope_x(i,j,k) - cs%slope_x(i,j,k+1) ) * ih
821  h_at_u(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
822  enddo ; enddo
823 
824  !### do J=js-2,Jeq+1 ; do i=is-1,ie+1
825  do j=js-2,jeq+1 ; do i=is-2,ieq+2
826  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) ) / &
827  ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
828  + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + gv%H_subroundoff )
829  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) ) / &
830  ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
831  + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + gv%H_subroundoff )
832  ih = 1./ ( ( h_at_slope_above + h_at_slope_below + gv%H_subroundoff ) * gv%H_to_Z )
833  dslopey_dz(i,j) = 2. * ( cs%slope_y(i,j,k) - cs%slope_y(i,j,k+1) ) * ih
834  h_at_v(i,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * ih
835  enddo ; enddo
836 
837  !### do J=js-1,je ; do i=is-1,Ieq+1
838  do j=js-2,jeq+1 ; do i=is-1,ieq+1
839  f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
840  vort_xy_dx(i,j) = vort_xy_dx(i,j) - f * us%L_to_Z * &
841  ( ( h_at_u(i,j) * dslopex_dz(i,j) + h_at_u(i-1,j+1) * dslopex_dz(i-1,j+1) ) &
842  + ( h_at_u(i-1,j) * dslopex_dz(i-1,j) + h_at_u(i,j+1) * dslopex_dz(i,j+1) ) ) / &
843  ( ( 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)
844  enddo ; enddo
845 
846  !### do j=js-1,Jeq+1 ; do I=is-1,ie
847  do j=js-1,jeq+1 ; do i=is-2,ieq+1
848  f = 0.5 * ( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
849  !### I think that this should be vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * &
850  vort_xy_dy(i,j) = vort_xy_dx(i,j) - f * us%L_to_Z * &
851  ( ( h_at_v(i,j) * dslopey_dz(i,j) + h_at_v(i+1,j-1) * dslopey_dz(i+1,j-1) ) &
852  + ( h_at_v(i,j-1) * dslopey_dz(i,j-1) + h_at_v(i+1,j) * dslopey_dz(i+1,j) ) ) / &
853  ( ( 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)
854  enddo ; enddo
855  endif ! k > 1
856 
857  !### I believe this halo update to be unnecessary. -RWH
858  call pass_vector(vort_xy_dy,vort_xy_dx,g%Domain)
859 
860  if (cs%use_QG_Leith_GM) then
861 
862  do j=js,je ; do i=is-1,ieq
863  !### These expressions are not rotationally symmetric. Add parentheses and regroup, as in:
864  ! 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)) +
865  ! (vort_xy_dx(i+1,J) + vort_xy_dx(i,J-1))))**2 )
866  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) &
867  + vort_xy_dx(i,j-1) + vort_xy_dx(i+1,j-1)))**2)
868  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) &
869  + div_xx_dy(i,j-1) + div_xx_dy(i+1,j-1)))**2)
870  if (cs%use_beta_in_QG_Leith) then
871  beta_u(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i+1,j))**2) + &
872  (0.5*(g%dF_dy(i,j)+g%dF_dy(i+1,j))**2) )
873  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)) * &
874  cs%Laplac3_const_u(i,j) * inv_pi3
875  else
876  cs%KH_u_QG(i,j,k) = (grad_vort_mag_u(i,j) + grad_div_mag_u(i,j)) * &
877  cs%Laplac3_const_u(i,j) * inv_pi3
878  endif
879  enddo ; enddo
880 
881  do j=js-1,jeq ; do i=is,ie
882  !### These expressions are not rotationally symmetric. Add parentheses and regroup.
883  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) &
884  + vort_xy_dy(i,j+1) + vort_xy_dy(i-1,j+1)))**2)
885  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) &
886  + div_xx_dx(i,j+1) + div_xx_dx(i-1,j+1)))**2)
887  if (cs%use_beta_in_QG_Leith) then
888  beta_v(i,j) = sqrt( (0.5*(g%dF_dx(i,j)+g%dF_dx(i,j+1))**2) + &
889  (0.5*(g%dF_dy(i,j)+g%dF_dy(i,j+1))**2) )
890  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)) * &
891  cs%Laplac3_const_v(i,j) * inv_pi3
892  else
893  cs%KH_v_QG(i,j,k) = (grad_vort_mag_v(i,j) + grad_div_mag_v(i,j)) * &
894  cs%Laplac3_const_v(i,j) * inv_pi3
895  endif
896  enddo ; enddo
897  ! post diagnostics
898 
899  if (k==nz) then
900  if (cs%id_KH_v_QG > 0) call post_data(cs%id_KH_v_QG, cs%KH_v_QG, cs%diag)
901  if (cs%id_KH_u_QG > 0) call post_data(cs%id_KH_u_QG, cs%KH_u_QG, cs%diag)
902  endif
903  endif
904 

Referenced by mom_hor_visc::horizontal_viscosity().

Here is the caller graph for this function:

◆ 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]gOcean grid structure
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvThermodynamic variables
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
csVariable mixing coefficients

Definition at line 193 of file MOM_lateral_mixing_coeffs.F90.

193  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
194  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
195  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
196  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
197  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
198  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
199 
200  ! Local variables
201  ! Depending on the power-function being used, dimensional rescaling may be limited, so some
202  ! of the following variables have units that depend on that power.
203  real :: cg1_q ! The gravity wave speed interpolated to q points [L T-1 ~> m s-1] or [m s-1].
204  real :: cg1_u ! The gravity wave speed interpolated to u points [L T-1 ~> m s-1] or [m s-1].
205  real :: cg1_v ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
206  real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
207  integer :: power_2
208  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
209  integer :: i, j, k
210  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
211  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
212 
213  if (.not. associated(cs)) call mom_error(fatal, "calc_resoln_function:"// &
214  "Module must be initialized before it is used.")
215  if (cs%calculate_cg1) then
216  if (.not. associated(cs%cg1)) call mom_error(fatal, &
217  "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
218  if (cs%khth_use_ebt_struct) then
219  if (.not. associated(cs%ebt_struct)) call mom_error(fatal, &
220  "calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
221  if (cs%Resoln_use_ebt) then
222  ! Both resolution fn and vertical structure are using EBT
223  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct)
224  else
225  ! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode
226  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp, modal_structure=cs%ebt_struct, &
227  use_ebt_mode=.true.)
228  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
229  endif
230  call pass_var(cs%ebt_struct, g%Domain)
231  else
232  call wave_speed(h, tv, g, gv, us, cs%cg1, cs%wave_speed_CSp)
233  endif
234 
235  call create_group_pass(cs%pass_cg1, cs%cg1, g%Domain)
236  call do_group_pass(cs%pass_cg1, g%Domain)
237  endif
238 
239  ! Calculate and store the ratio between deformation radius and grid-spacing
240  ! at h-points [nondim].
241  if (cs%calculate_rd_dx) then
242  if (.not. associated(cs%Rd_dx_h)) call mom_error(fatal, &
243  "calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
244  !$OMP parallel do default(shared)
245  do j=js-1,je+1 ; do i=is-1,ie+1
246  cs%Rd_dx_h(i,j) = cs%cg1(i,j) / &
247  (sqrt(cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)))
248  enddo ; enddo
249  if (query_averaging_enabled(cs%diag)) then
250  if (cs%id_Rd_dx > 0) call post_data(cs%id_Rd_dx, cs%Rd_dx_h, cs%diag)
251  endif
252  endif
253 
254  if (.not. cs%calculate_res_fns) return
255 
256  if (.not. associated(cs%Res_fn_h)) call mom_error(fatal, &
257  "calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
258  if (.not. associated(cs%Res_fn_q)) call mom_error(fatal, &
259  "calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
260  if (.not. associated(cs%Res_fn_u)) call mom_error(fatal, &
261  "calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
262  if (.not. associated(cs%Res_fn_v)) call mom_error(fatal, &
263  "calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
264  if (.not. associated(cs%f2_dx2_h)) call mom_error(fatal, &
265  "calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
266  if (.not. associated(cs%f2_dx2_q)) call mom_error(fatal, &
267  "calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
268  if (.not. associated(cs%f2_dx2_u)) call mom_error(fatal, &
269  "calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
270  if (.not. associated(cs%f2_dx2_v)) call mom_error(fatal, &
271  "calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
272  if (.not. associated(cs%beta_dx2_h)) call mom_error(fatal, &
273  "calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
274  if (.not. associated(cs%beta_dx2_q)) call mom_error(fatal, &
275  "calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
276  if (.not. associated(cs%beta_dx2_u)) call mom_error(fatal, &
277  "calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
278  if (.not. associated(cs%beta_dx2_v)) call mom_error(fatal, &
279  "calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")
280 
281  ! Do this calculation on the extent used in MOM_hor_visc.F90, and
282  ! MOM_tracer.F90 so that no halo update is needed.
283 
284 !$OMP parallel default(none) shared(is,ie,js,je,Ieq,Jeq,CS,US) &
285 !$OMP private(dx_term,cg1_q,power_2,cg1_u,cg1_v)
286  if (cs%Res_fn_power_visc >= 100) then
287 !$OMP do
288  do j=js-1,je+1 ; do i=is-1,ie+1
289  dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
290  if ((cs%Res_coef_visc * cs%cg1(i,j))**2 > dx_term) then
291  cs%Res_fn_h(i,j) = 0.0
292  else
293  cs%Res_fn_h(i,j) = 1.0
294  endif
295  enddo ; enddo
296 !$OMP do
297  do j=js-1,jeq ; do i=is-1,ieq
298  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
299  dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
300  if ((cs%Res_coef_visc * cg1_q)**2 > dx_term) then
301  cs%Res_fn_q(i,j) = 0.0
302  else
303  cs%Res_fn_q(i,j) = 1.0
304  endif
305  enddo ; enddo
306  elseif (cs%Res_fn_power_visc == 2) then
307 !$OMP do
308  do j=js-1,je+1 ; do i=is-1,ie+1
309  dx_term = cs%f2_dx2_h(i,j) + cs%cg1(i,j)*cs%beta_dx2_h(i,j)
310  cs%Res_fn_h(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cs%cg1(i,j))**2)
311  enddo ; enddo
312 !$OMP do
313  do j=js-1,jeq ; do i=is-1,ieq
314  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
315  dx_term = cs%f2_dx2_q(i,j) + cg1_q * cs%beta_dx2_q(i,j)
316  cs%Res_fn_q(i,j) = dx_term / (dx_term + (cs%Res_coef_visc * cg1_q)**2)
317  enddo ; enddo
318  elseif (mod(cs%Res_fn_power_visc, 2) == 0) then
319  power_2 = cs%Res_fn_power_visc / 2
320 !$OMP do
321  do j=js-1,je+1 ; do i=is-1,ie+1
322  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
323  cs%Res_fn_h(i,j) = dx_term / &
324  (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
325  enddo ; enddo
326 !$OMP do
327  do j=js-1,jeq ; do i=is-1,ieq
328  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
329  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
330  cs%Res_fn_q(i,j) = dx_term / &
331  (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
332  enddo ; enddo
333  else
334 !$OMP do
335  do j=js-1,je+1 ; do i=is-1,ie+1
336  dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_h(i,j) + &
337  cs%cg1(i,j)*cs%beta_dx2_h(i,j)))**cs%Res_fn_power_visc
338  cs%Res_fn_h(i,j) = dx_term / &
339  (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cs%cg1(i,j))**cs%Res_fn_power_visc)
340  enddo ; enddo
341 !$OMP do
342  do j=js-1,jeq ; do i=is-1,ieq
343  cg1_q = 0.25 * ((cs%cg1(i,j) + cs%cg1(i+1,j+1)) + (cs%cg1(i+1,j) + cs%cg1(i,j+1)))
344  dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_q(i,j) + &
345  cg1_q * cs%beta_dx2_q(i,j)))**cs%Res_fn_power_visc
346  cs%Res_fn_q(i,j) = dx_term / &
347  (dx_term + (cs%Res_coef_visc * us%L_T_to_m_s*cg1_q)**cs%Res_fn_power_visc)
348  enddo ; enddo
349  endif
350 
351  if (cs%interpolate_Res_fn) then
352  do j=js,je ; do i=is-1,ieq
353  cs%Res_fn_u(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i+1,j))
354  enddo ; enddo
355  do j=js-1,jeq ; do i=is,ie
356  cs%Res_fn_v(i,j) = 0.5*(cs%Res_fn_h(i,j) + cs%Res_fn_h(i,j+1))
357  enddo ; enddo
358  else ! .not.CS%interpolate_Res_fn
359  if (cs%Res_fn_power_khth >= 100) then
360 !$OMP do
361  do j=js,je ; do i=is-1,ieq
362  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
363  dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
364  if ((cs%Res_coef_khth * cg1_u)**2 > dx_term) then
365  cs%Res_fn_u(i,j) = 0.0
366  else
367  cs%Res_fn_u(i,j) = 1.0
368  endif
369  enddo ; enddo
370 !$OMP do
371  do j=js-1,jeq ; do i=is,ie
372  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
373  dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
374  if ((cs%Res_coef_khth * cg1_v)**2 > dx_term) then
375  cs%Res_fn_v(i,j) = 0.0
376  else
377  cs%Res_fn_v(i,j) = 1.0
378  endif
379  enddo ; enddo
380  elseif (cs%Res_fn_power_khth == 2) then
381 !$OMP do
382  do j=js,je ; do i=is-1,ieq
383  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
384  dx_term = cs%f2_dx2_u(i,j) + cg1_u * cs%beta_dx2_u(i,j)
385  cs%Res_fn_u(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_u)**2)
386  enddo ; enddo
387 !$OMP do
388  do j=js-1,jeq ; do i=is,ie
389  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
390  dx_term = cs%f2_dx2_v(i,j) + cg1_v * cs%beta_dx2_v(i,j)
391  cs%Res_fn_v(i,j) = dx_term / (dx_term + (cs%Res_coef_khth * cg1_v)**2)
392  enddo ; enddo
393  elseif (mod(cs%Res_fn_power_khth, 2) == 0) then
394  power_2 = cs%Res_fn_power_khth / 2
395 !$OMP do
396  do j=js,je ; do i=is-1,ieq
397  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
398  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
399  cs%Res_fn_u(i,j) = dx_term / &
400  (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
401  enddo ; enddo
402 !$OMP do
403  do j=js-1,jeq ; do i=is,ie
404  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
405  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
406  cs%Res_fn_v(i,j) = dx_term / &
407  (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
408  enddo ; enddo
409  else
410 !$OMP do
411  do j=js,je ; do i=is-1,ieq
412  cg1_u = 0.5 * (cs%cg1(i,j) + cs%cg1(i+1,j))
413  dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_u(i,j) + &
414  cg1_u * cs%beta_dx2_u(i,j)))**cs%Res_fn_power_khth
415  cs%Res_fn_u(i,j) = dx_term / &
416  (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_u)**cs%Res_fn_power_khth)
417  enddo ; enddo
418 !$OMP do
419  do j=js-1,jeq ; do i=is,ie
420  cg1_v = 0.5 * (cs%cg1(i,j) + cs%cg1(i,j+1))
421  dx_term = (us%L_T_to_m_s*sqrt(cs%f2_dx2_v(i,j) + &
422  cg1_v * cs%beta_dx2_v(i,j)))**cs%Res_fn_power_khth
423  cs%Res_fn_v(i,j) = dx_term / &
424  (dx_term + (cs%Res_coef_khth * us%L_T_to_m_s*cg1_v)**cs%Res_fn_power_khth)
425  enddo ; enddo
426  endif
427  endif
428 !$OMP end parallel
429 
430  if (query_averaging_enabled(cs%diag)) then
431  if (cs%id_Res_fn > 0) call post_data(cs%id_Res_fn, cs%Res_fn_h, cs%diag)
432  endif
433 

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

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(unit_scale_type), intent(in)  US,
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]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in,out]hLayer thickness [H ~> m or kg m-2]
[in]tvThermodynamic variables
[in]dtTime increment [T ~> s]
csVariable mixing coefficients

Definition at line 439 of file MOM_lateral_mixing_coeffs.F90.

439  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
440  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
441  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
442  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
443  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
444  real, intent(in) :: dt !< Time increment [T ~> s]
445  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
446  ! Local variables
447  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
448  e ! The interface heights relative to mean sea level [Z ~> m].
449  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [T-2 ~> s-2]
450  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [T-2 ~> s-2]
451 
452  if (.not. associated(cs)) call mom_error(fatal, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//&
453  "Module must be initialized before it is used.")
454 
455  if (cs%calculate_Eady_growth_rate) then
456  call find_eta(h, tv, g, gv, us, e, halo_size=2)
457  if (cs%use_stored_slopes) then
458  call calc_isoneutral_slopes(g, gv, us, h, e, tv, dt*cs%kappa_smooth, &
459  cs%slope_x, cs%slope_y, n2_u, n2_v, 1)
460  call calc_visbeck_coeffs(h, cs%slope_x, cs%slope_y, n2_u, n2_v, g, gv, us, cs)
461 ! call calc_slope_functions_using_just_e(h, G, CS, e, .false.)
462  else
463  !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y)
464  call calc_slope_functions_using_just_e(h, g, gv, us, cs, e, .true.)
465  endif
466  endif
467 
468  if (query_averaging_enabled(cs%diag)) then
469  if (cs%id_SN_u > 0) call post_data(cs%id_SN_u, cs%SN_u, cs%diag)
470  if (cs%id_SN_v > 0) call post_data(cs%id_SN_v, cs%SN_v, cs%diag)
471  if (cs%id_L2u > 0) call post_data(cs%id_L2u, cs%L2u, cs%diag)
472  if (cs%id_L2v > 0) call post_data(cs%id_L2v, cs%L2v, cs%diag)
473  !### I do not believe that CS%N2_u and CS%N2_v are ever set, but because the contents
474  ! of CS are public, they might be set somewhere outside of this module.
475  if (cs%id_N2_u > 0) call post_data(cs%id_N2_u, cs%N2_u, cs%diag)
476  if (cs%id_N2_v > 0) call post_data(cs%id_N2_v, cs%N2_v, cs%diag)
477  endif
478 

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

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

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

Parameters
[in,out]gOcean grid structure
[in,out]hLayer thickness [H ~> m or kg m-2]
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
csVariable mixing coefficients
[in]eInterface position [Z ~> m]
[in]calculate_slopesIf true, calculate slopes internally otherwise use slopes stored in CS

Definition at line 618 of file MOM_lateral_mixing_coeffs.F90.

618  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
619  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
620  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
621  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
622  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
623  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface position [Z ~> m]
624  logical, intent(in) :: calculate_slopes !< If true, calculate slopes internally
625  !! otherwise use slopes stored in CS
626  ! Local variables
627  real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics)
628  real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics)
629  real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
630  real :: h_neglect ! A thickness that is so small it is usually lost
631  ! in roundoff and can be neglected [H ~> m or kg m-2].
632  real :: S2 ! Interface slope squared [nondim]
633  real :: N2 ! Brunt-Vaisala frequency squared [T-2 ~> s-2]
634  real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
635  real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
636  real :: one_meter ! One meter in thickness units [H ~> m or kg m-2].
637  integer :: is, ie, js, je, nz
638  integer :: i, j, k, kb_max
639  real :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(G))
640  real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(G))
641 
642  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
643  "Module must be initialized before it is used.")
644  if (.not. cs%calculate_Eady_growth_rate) return
645  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
646  "%SN_u is not associated with use_variable_mixing.")
647  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
648  "%SN_v is not associated with use_variable_mixing.")
649 
650  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
651 
652  one_meter = 1.0 * gv%m_to_H
653  h_neglect = gv%H_subroundoff
654  h_cutoff = real(2*nz) * (gv%Angstrom_H + h_neglect)
655 
656  ! To set the length scale based on the deformation radius, use wave_speed to
657  ! calculate the first-mode gravity wave speed and then blend the equatorial
658  ! and midlatitude deformation radii, using calc_resoln_function as a template.
659 
660  !$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
661  do k=nz,cs%VarMix_Ktop,-1
662 
663  if (calculate_slopes) then
664  ! Calculate the interface slopes E_x and E_y and u- and v- points respectively
665  do j=js-1,je+1 ; do i=is-1,ie
666  e_x(i,j) = us%Z_to_L*(e(i+1,j,k)-e(i,j,k))*g%IdxCu(i,j)
667  ! Mask slopes where interface intersects topography
668  if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
669  enddo ; enddo
670  do j=js-1,je ; do i=is-1,ie+1
671  e_y(i,j) = us%Z_to_L*(e(i,j+1,k)-e(i,j,k))*g%IdyCv(i,j)
672  ! Mask slopes where interface intersects topography
673  if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
674  enddo ; enddo
675  else
676  do j=js-1,je+1 ; do i=is-1,ie
677  e_x(i,j) = cs%slope_x(i,j,k)
678  if (min(h(i,j,k),h(i+1,j,k)) < h_cutoff) e_x(i,j) = 0.
679  enddo ; enddo
680  do j=js-1,je ; do i=is-1,ie+1
681  e_y(i,j) = cs%slope_y(i,j,k)
682  if (min(h(i,j,k),h(i,j+1,k)) < h_cutoff) e_y(i,j) = 0.
683  enddo ; enddo
684  endif
685 
686  ! Calculate N*S*h from this layer and add to the sum
687  do j=js,je ; do i=is-1,ie
688  s2 = ( e_x(i,j)**2 + 0.25*( &
689  (e_y(i,j)**2+e_y(i+1,j-1)**2)+(e_y(i+1,j)**2+e_y(i,j-1)**2) ) )
690  hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
691  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)
692  h_geom = sqrt(hdn*hup)
693  n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
694  if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < h_cutoff) &
695  s2 = 0.0
696  s2n2_u_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
697  enddo ; enddo
698  do j=js-1,je ; do i=is,ie
699  s2 = ( e_y(i,j)**2 + 0.25*( &
700  (e_x(i,j)**2+e_x(i-1,j+1)**2)+(e_x(i,j+1)**2+e_x(i-1,j)**2) ) )
701  hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
702  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)
703  h_geom = sqrt(hdn*hup)
704  n2 = gv%g_prime(k)*us%L_to_Z**2 / (gv%H_to_Z * max(hdn,hup,one_meter))
705  if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < h_cutoff) &
706  s2 = 0.0
707  s2n2_v_local(i,j,k) = (h_geom * gv%H_to_Z) * s2 * n2
708  enddo ; enddo
709 
710  enddo ! k
711  !$OMP parallel do default(shared)
712  do j=js,je
713  do i=is-1,ie ; cs%SN_u(i,j) = 0.0 ; enddo
714  do k=nz,cs%VarMix_Ktop,-1 ; do i=is-1,ie
715  cs%SN_u(i,j) = cs%SN_u(i,j) + s2n2_u_local(i,j,k)
716  enddo ; enddo
717  ! SN above contains S^2*N^2*H, convert to vertical average of S*N
718  do i=is-1,ie
719  !SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(I,j), G%bathyT(I+1,j)) + GV%Angstrom_Z ) ))
720  !The code below behaves better than the line above. Not sure why? AJA
721  if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
722  cs%SN_u(i,j) = g%mask2dCu(i,j) * sqrt( cs%SN_u(i,j) / &
723  (max(g%bathyT(i,j), g%bathyT(i+1,j))) )
724  else
725  cs%SN_u(i,j) = 0.0
726  endif
727  enddo
728  enddo
729  !$OMP parallel do default(shared)
730  do j=js-1,je
731  do i=is,ie ; cs%SN_v(i,j) = 0.0 ; enddo
732  do k=nz,cs%VarMix_Ktop,-1 ; do i=is,ie
733  cs%SN_v(i,j) = cs%SN_v(i,j) + s2n2_v_local(i,j,k)
734  enddo ; enddo
735  do i=is,ie
736  !SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + GV%Angstrom_Z ) ))
737  !The code below behaves better than the line above. Not sure why? AJA
738  if ( min(g%bathyT(i,j), g%bathyT(i+1,j)) > h_cutoff*gv%H_to_Z ) then
739  cs%SN_v(i,j) = g%mask2dCv(i,j) * sqrt( cs%SN_v(i,j) / &
740  (max(g%bathyT(i,j), g%bathyT(i,j+1))) )
741  else
742  cs%SN_v(i,j) = 0.0
743  endif
744  enddo
745  enddo
746 

References mom_error_handler::mom_error().

Referenced by calc_slope_functions().

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 %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 
)
private

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

Parameters
[in,out]gOcean grid structure
[in]gvVertical grid structure
[in]hLayer thickness [H ~> m or kg m-2]
[in]slope_xZonal isoneutral slope
[in]n2_uBuoyancy (Brunt-Vaisala) frequency at u-points [T-2 ~> s-2]
[in]slope_yMeridional isoneutral slope
[in]n2_vBuoyancy (Brunt-Vaisala) frequency at v-points [T-2 ~> s-2]
[in]usA dimensional unit scaling type
csVariable mixing coefficients

Definition at line 483 of file MOM_lateral_mixing_coeffs.F90.

483  type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
484  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
485  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
486  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: slope_x !< Zonal isoneutral slope
487  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency
488  !! at u-points [T-2 ~> s-2]
489  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: slope_y !< Meridional isoneutral slope
490  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency
491  !! at v-points [T-2 ~> s-2]
492  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
493  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
494 
495  ! Local variables
496  real :: S2 ! Interface slope squared [nondim]
497  real :: N2 ! Positive buoyancy frequency or zero [T-2 ~> s-2]
498  real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
499  real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
500  integer :: is, ie, js, je, nz
501  integer :: i, j, k, kb_max
502  real :: S2max, wNE, wSE, wSW, wNW
503  real :: H_u(SZIB_(G)), H_v(SZI_(G))
504  real :: S2_u(SZIB_(G), SZJ_(G))
505  real :: S2_v(SZI_(G), SZJB_(G))
506 
507  if (.not. associated(cs)) call mom_error(fatal, "calc_slope_function:"// &
508  "Module must be initialized before it is used.")
509  if (.not. cs%calculate_Eady_growth_rate) return
510  if (.not. associated(cs%SN_u)) call mom_error(fatal, "calc_slope_function:"// &
511  "%SN_u is not associated with use_variable_mixing.")
512  if (.not. associated(cs%SN_v)) call mom_error(fatal, "calc_slope_function:"// &
513  "%SN_v is not associated with use_variable_mixing.")
514 
515  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
516 
517  s2max = cs%Visbeck_S_max**2
518 
519  !$OMP parallel do default(shared)
520  do j=js-1,je+1 ; do i=is-1,ie+1
521  cs%SN_u(i,j) = 0.0
522  cs%SN_v(i,j) = 0.0
523  enddo ; enddo
524 
525  ! To set the length scale based on the deformation radius, use wave_speed to
526  ! calculate the first-mode gravity wave speed and then blend the equatorial
527  ! and midlatitude deformation radii, using calc_resoln_function as a template.
528 
529  !$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
530  do j = js,je
531  do i=is-1,ie
532  cs%SN_u(i,j) = 0. ; h_u(i) = 0. ; s2_u(i,j) = 0.
533  enddo
534  do k=2,nz ; do i=is-1,ie
535  hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
536  hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
537  h_geom = sqrt( hdn * hup )
538  !H_geom = H_geom * sqrt(N2) ! WKB-ish
539  !H_geom = H_geom * N2 ! WKB-ish
540  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)) )
541  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)) )
542  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)) )
543  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)) )
544  s2 = slope_x(i,j,k)**2 + &
545  ((wnw*slope_y(i,j,k)**2 + wse*slope_y(i+1,j-1,k)**2) + &
546  (wne*slope_y(i+1,j,k)**2 + wsw*slope_y(i,j-1,k)**2) ) / &
547  ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
548  if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
549 
550  n2 = max(0., n2_u(i,j,k))
551  cs%SN_u(i,j) = cs%SN_u(i,j) + sqrt( s2*n2 )*h_geom
552  s2_u(i,j) = s2_u(i,j) + s2*h_geom
553  h_u(i) = h_u(i) + h_geom
554  enddo ; enddo
555  do i=is-1,ie
556  if (h_u(i)>0.) then
557  cs%SN_u(i,j) = g%mask2dCu(i,j) * cs%SN_u(i,j) / h_u(i)
558  s2_u(i,j) = g%mask2dCu(i,j) * s2_u(i,j) / h_u(i)
559  else
560  cs%SN_u(i,j) = 0.
561  endif
562  enddo
563  enddo
564 
565  !$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
566  do j = js-1,je
567  do i=is,ie
568  cs%SN_v(i,j) = 0. ; h_v(i) = 0. ; s2_v(i,j) = 0.
569  enddo
570  do k=2,nz ; do i=is,ie
571  hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
572  hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
573  h_geom = sqrt( hdn * hup )
574  !H_geom = H_geom * sqrt(N2) ! WKB-ish
575  !H_geom = H_geom * N2 ! WKB-ish
576  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)) )
577  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)) )
578  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)) )
579  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)) )
580  s2 = slope_y(i,j,k)**2 + &
581  ((wse*slope_x(i,j,k)**2 + wnw*slope_x(i-1,j+1,k)**2) + &
582  (wne*slope_x(i,j+1,k)**2 + wsw*slope_x(i-1,j,k)**2) ) / &
583  ( ((wse+wnw) + (wne+wsw)) + gv%H_subroundoff**4 )
584  if (s2max>0.) s2 = s2 * s2max / (s2 + s2max) ! Limit S2
585 
586  n2 = max(0., n2_v(i,j,k))
587  cs%SN_v(i,j) = cs%SN_v(i,j) + sqrt( s2*n2 )*h_geom
588  s2_v(i,j) = s2_v(i,j) + s2*h_geom
589  h_v(i) = h_v(i) + h_geom
590  enddo ; enddo
591  do i=is,ie
592  if (h_v(i)>0.) then
593  cs%SN_v(i,j) = g%mask2dCv(i,j) * cs%SN_v(i,j) / h_v(i)
594  s2_v(i,j) = g%mask2dCv(i,j) * s2_v(i,j) / h_v(i)
595  else
596  cs%SN_v(i,j) = 0.
597  endif
598  enddo
599  enddo
600 
601 ! Offer diagnostic fields for averaging.
602  if (query_averaging_enabled(cs%diag)) then
603  if (cs%id_S2_u > 0) call post_data(cs%id_S2_u, s2_u, cs%diag)
604  if (cs%id_S2_v > 0) call post_data(cs%id_S2_v, s2_v, cs%diag)
605  endif
606 
607  if (cs%debug) then
608  call uvchksum("calc_Visbeck_coeffs slope_[xy]", slope_x, slope_y, g%HI, haloshift=1)
609  call uvchksum("calc_Visbeck_coeffs N2_u, N2_v", n2_u, n2_v, g%HI, scale=us%s_to_T**2)
610  call uvchksum("calc_Visbeck_coeffs SN_[uv]", cs%SN_u, cs%SN_v, g%HI, scale=us%s_to_T)
611  endif
612 

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

Referenced by calc_slope_functions().

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(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]timeCurrent model time
[in]gOcean grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileParameter file handles
[in,out]diagDiagnostics control structure
csVariable mixing coefficients

Definition at line 909 of file MOM_lateral_mixing_coeffs.F90.

909  type(time_type), intent(in) :: Time !< Current model time
910  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
911  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
912  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
913  type(param_file_type), intent(in) :: param_file !< Parameter file handles
914  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
915  type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
916  ! Local variables
917  real :: KhTr_Slope_Cff, KhTh_Slope_Cff, oneOrTwo
918  real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when
919  ! calculating the first-mode wave speed [Z ~> m]
920  real :: KhTr_passivity_coeff
921  real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
922  ! default value is roughly (pi / (the age of the universe)).
923  logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
924  real :: MLE_front_length
925  real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity
926  real :: grid_sp_u2, grid_sp_v2 ! Intermediate quantities for Leith metrics [L2 ~> m2]
927  real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
928 ! This include declares and sets the variable "version".
929 #include "version_variable.h"
930  character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
931  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j
932  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
933  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
934  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
935  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
936  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
937 
938  if (associated(cs)) then
939  call mom_error(warning, "VarMix_init called with an associated "// &
940  "control structure.")
941  return
942  endif
943 
944  allocate(cs)
945  in_use = .false. ! Set to true to avoid deallocating
946  cs%diag => diag ! Diagnostics pointer
947  cs%calculate_cg1 = .false.
948  cs%calculate_Rd_dx = .false.
949  cs%calculate_res_fns = .false.
950  cs%calculate_Eady_growth_rate = .false.
951  cs%calculate_depth_fns = .false.
952  ! Read all relevant parameters and write them to the model log.
953  call log_version(param_file, mdl, version, "")
954  call get_param(param_file, mdl, "USE_VARIABLE_MIXING", cs%use_variable_mixing,&
955  "If true, the variable mixing code will be called. This "//&
956  "allows diagnostics to be created even if the scheme is "//&
957  "not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
958  "this is set to true regardless of what is in the "//&
959  "parameter file.", default=.false.)
960  call get_param(param_file, mdl, "USE_VISBECK", cs%use_Visbeck,&
961  "If true, use the Visbeck et al. (1997) formulation for \n"//&
962  "thickness diffusivity.", default=.false.)
963  call get_param(param_file, mdl, "RESOLN_SCALED_KH", cs%Resoln_scaled_Kh, &
964  "If true, the Laplacian lateral viscosity is scaled away "//&
965  "when the first baroclinic deformation radius is well "//&
966  "resolved.", default=.false.)
967  call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", cs%Depth_scaled_KhTh, &
968  "If true, KHTH is scaled away when the depth is shallower"//&
969  "than a reference depth: KHTH = MIN(1,H/H0)**N * KHTH, "//&
970  "where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
971  "and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
972  default=.false.)
973  call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", cs%Resoln_scaled_KhTh, &
974  "If true, the interface depth diffusivity is scaled away "//&
975  "when the first baroclinic deformation radius is well "//&
976  "resolved.", default=.false.)
977  call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", cs%Resoln_scaled_KhTr, &
978  "If true, the epipycnal tracer diffusivity is scaled "//&
979  "away when the first baroclinic deformation radius is "//&
980  "well resolved.", default=.false.)
981  call get_param(param_file, mdl, "RESOLN_USE_EBT", cs%Resoln_use_ebt, &
982  "If true, uses the equivalent barotropic wave speed instead "//&
983  "of first baroclinic wave for calculating the resolution fn.",&
984  default=.false.)
985  call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", cs%khth_use_ebt_struct, &
986  "If true, uses the equivalent barotropic structure "//&
987  "as the vertical structure of thickness diffusivity.",&
988  default=.false.)
989  call get_param(param_file, mdl, "KHTH_SLOPE_CFF", khth_slope_cff, &
990  "The nondimensional coefficient in the Visbeck formula "//&
991  "for the interface depth diffusivity", units="nondim", &
992  default=0.0)
993  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", khtr_slope_cff, &
994  "The nondimensional coefficient in the Visbeck formula "//&
995  "for the epipycnal tracer diffusivity", units="nondim", &
996  default=0.0)
997  call get_param(param_file, mdl, "USE_STORED_SLOPES", cs%use_stored_slopes,&
998  "If true, the isopycnal slopes are calculated once and "//&
999  "stored for re-use. This uses more memory but avoids calling "//&
1000  "the equation of state more times than should be necessary.", &
1001  default=.false.)
1002  call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1003  "A miniscule frequency that is used to avoid division by 0. The default "//&
1004  "value is roughly (pi / (the age of the universe)).", &
1005  default=1.0e-17, units="s-1", scale=us%T_to_s)
1006  call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_fgnv_streamfn, &
1007  default=.false., do_not_log=.true.)
1008  cs%calculate_cg1 = cs%calculate_cg1 .or. use_fgnv_streamfn
1009  call get_param(param_file, mdl, "USE_MEKE", use_meke, &
1010  default=.false., do_not_log=.true.)
1011  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. use_meke
1012  cs%calculate_Eady_growth_rate = cs%calculate_Eady_growth_rate .or. use_meke
1013  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", khtr_passivity_coeff, &
1014  default=0., do_not_log=.true.)
1015  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (khtr_passivity_coeff>0.)
1016  call get_param(param_file, mdl, "MLE_FRONT_LENGTH", mle_front_length, &
1017  default=0., do_not_log=.true.)
1018  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (mle_front_length>0.)
1019 
1020  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
1021 
1022 
1023  if (cs%Resoln_use_ebt .or. cs%khth_use_ebt_struct) then
1024  in_use = .true.
1025  call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", n2_filter_depth, &
1026  "The depth below which N2 is monotonized to avoid stratification "//&
1027  "artifacts from altering the equivalent barotropic mode structure.",&
1028  units="m", default=2000., scale=us%m_to_Z)
1029  allocate(cs%ebt_struct(isd:ied,jsd:jed,g%ke)) ; cs%ebt_struct(:,:,:) = 0.0
1030  endif
1031 
1032  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1033  cs%calculate_Eady_growth_rate = .true.
1034  call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", cs%Visbeck_S_max, &
1035  "If non-zero, is an upper bound on slopes used in the "//&
1036  "Visbeck formula for diffusivity. This does not affect the "//&
1037  "isopycnal slope calculation used within thickness diffusion.", &
1038  units="nondim", default=0.0)
1039  endif
1040 
1041  if (cs%use_stored_slopes) then
1042  in_use = .true.
1043  allocate(cs%slope_x(isdb:iedb,jsd:jed,g%ke+1)) ; cs%slope_x(:,:,:) = 0.0
1044  allocate(cs%slope_y(isd:ied,jsdb:jedb,g%ke+1)) ; cs%slope_y(:,:,:) = 0.0
1045  allocate(cs%N2_u(isdb:iedb,jsd:jed,g%ke+1)) ; cs%N2_u(:,:,:) = 0.0
1046  allocate(cs%N2_v(isd:ied,jsdb:jedb,g%ke+1)) ; cs%N2_v(:,:,:) = 0.0
1047  call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
1048  "A diapycnal diffusivity that is used to interpolate "//&
1049  "more sensible values of T & S into thin layers.", &
1050  units="m2 s-1", default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
1051  endif
1052 
1053  if (cs%calculate_Eady_growth_rate) then
1054  in_use = .true.
1055  allocate(cs%SN_u(isdb:iedb,jsd:jed)) ; cs%SN_u(:,:) = 0.0
1056  allocate(cs%SN_v(isd:ied,jsdb:jedb)) ; cs%SN_v(:,:) = 0.0
1057  cs%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, time, &
1058  'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=us%s_to_T)
1059  cs%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, time, &
1060  'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=us%s_to_T)
1061  call get_param(param_file, mdl, "VARMIX_KTOP", cs%VarMix_Ktop, &
1062  "The layer number at which to start vertical integration "//&
1063  "of S*N for purposes of finding the Eady growth rate.", &
1064  units="nondim", default=2)
1065  endif
1066 
1067  if (khtr_slope_cff>0. .or. khth_slope_cff>0.) then
1068  in_use = .true.
1069  call get_param(param_file, mdl, "VISBECK_L_SCALE", cs%Visbeck_L_scale, &
1070  "The fixed length scale in the Visbeck formula.", units="m", &
1071  default=0.0)
1072  allocate(cs%L2u(isdb:iedb,jsd:jed)) ; cs%L2u(:,:) = 0.0
1073  allocate(cs%L2v(isd:ied,jsdb:jedb)) ; cs%L2v(:,:) = 0.0
1074  if (cs%Visbeck_L_scale<0) then
1075  do j=js,je ; do i=is-1,ieq
1076  cs%L2u(i,j) = cs%Visbeck_L_scale**2 * g%areaCu(i,j)
1077  enddo; enddo
1078  do j=js-1,jeq ; do i=is,ie
1079  cs%L2v(i,j) = cs%Visbeck_L_scale**2 * g%areaCv(i,j)
1080  enddo; enddo
1081  else
1082  cs%L2u(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
1083  cs%L2v(:,:) = us%m_to_L**2*cs%Visbeck_L_scale**2
1084  endif
1085 
1086  cs%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, time, &
1087  'Length scale squared for mixing coefficient, at u-points', &
1088  'm2', conversion=us%L_to_m**2)
1089  cs%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, time, &
1090  'Length scale squared for mixing coefficient, at v-points', &
1091  'm2', conversion=us%L_to_m**2)
1092  endif
1093 
1094  if (cs%use_stored_slopes) then
1095  cs%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, time, &
1096  'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', 's-2')
1097  cs%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, time, &
1098  'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', 's-2')
1099  !### The units of the next two diagnostics should be 'nondim'.
1100  cs%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, time, &
1101  'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 's-2')
1102  cs%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, time, &
1103  'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 's-2')
1104  endif
1105 
1106  oneortwo = 1.0
1107  if (cs%Resoln_scaled_Kh .or. cs%Resoln_scaled_KhTh .or. cs%Resoln_scaled_KhTr) then
1108  cs%calculate_Rd_dx = .true.
1109  cs%calculate_res_fns = .true.
1110  allocate(cs%Res_fn_h(isd:ied,jsd:jed)) ; cs%Res_fn_h(:,:) = 0.0
1111  allocate(cs%Res_fn_q(isdb:iedb,jsdb:jedb)) ; cs%Res_fn_q(:,:) = 0.0
1112  allocate(cs%Res_fn_u(isdb:iedb,jsd:jed)) ; cs%Res_fn_u(:,:) = 0.0
1113  allocate(cs%Res_fn_v(isd:ied,jsdb:jedb)) ; cs%Res_fn_v(:,:) = 0.0
1114  allocate(cs%beta_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%beta_dx2_q(:,:) = 0.0
1115  allocate(cs%beta_dx2_u(isdb:iedb,jsd:jed)) ; cs%beta_dx2_u(:,:) = 0.0
1116  allocate(cs%beta_dx2_v(isd:ied,jsdb:jedb)) ; cs%beta_dx2_v(:,:) = 0.0
1117  allocate(cs%f2_dx2_q(isdb:iedb,jsdb:jedb)) ; cs%f2_dx2_q(:,:) = 0.0
1118  allocate(cs%f2_dx2_u(isdb:iedb,jsd:jed)) ; cs%f2_dx2_u(:,:) = 0.0
1119  allocate(cs%f2_dx2_v(isd:ied,jsdb:jedb)) ; cs%f2_dx2_v(:,:) = 0.0
1120 
1121  cs%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, time, &
1122  'Resolution function for scaling diffusivities', 'nondim')
1123 
1124  call get_param(param_file, mdl, "KH_RES_SCALE_COEF", cs%Res_coef_khth, &
1125  "A coefficient that determines how KhTh is scaled away if "//&
1126  "RESOLN_SCALED_... is true, as "//&
1127  "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
1128  units="nondim", default=1.0)
1129  call get_param(param_file, mdl, "KH_RES_FN_POWER", cs%Res_fn_power_khth, &
1130  "The power of dx/Ld in the Kh resolution function. Any "//&
1131  "positive integer may be used, although even integers "//&
1132  "are more efficient to calculate. Setting this greater "//&
1133  "than 100 results in a step-function being used.", &
1134  units="nondim", default=2)
1135  call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", cs%Res_coef_visc, &
1136  "A coefficient that determines how Kh is scaled away if "//&
1137  "RESOLN_SCALED_... is true, as "//&
1138  "F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
1139  "This function affects lateral viscosity, Kh, and not KhTh.", &
1140  units="nondim", default=cs%Res_coef_khth)
1141  call get_param(param_file, mdl, "VISC_RES_FN_POWER", cs%Res_fn_power_visc, &
1142  "The power of dx/Ld in the Kh resolution function. Any "//&
1143  "positive integer may be used, although even integers "//&
1144  "are more efficient to calculate. Setting this greater "//&
1145  "than 100 results in a step-function being used. "//&
1146  "This function affects lateral viscosity, Kh, and not KhTh.", &
1147  units="nondim", default=cs%Res_fn_power_khth)
1148  call get_param(param_file, mdl, "INTERPOLATE_RES_FN", cs%interpolate_Res_fn, &
1149  "If true, interpolate the resolution function to the "//&
1150  "velocity points from the thickness points; otherwise "//&
1151  "interpolate the wave speed and calculate the resolution "//&
1152  "function independently at each point.", default=.true.)
1153  if (cs%interpolate_Res_fn) then
1154  if (cs%Res_coef_visc /= cs%Res_coef_khth) call mom_error(fatal, &
1155  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1156  "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
1157  if (cs%Res_fn_power_visc /= cs%Res_fn_power_khth) call mom_error(fatal, &
1158  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1159  "When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
1160  endif
1161  !### Change the default of GILL_EQUATORIAL_LD to True.
1162  call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", gill_equatorial_ld, &
1163  "If true, uses Gill's definition of the baroclinic "//&
1164  "equatorial deformation radius, otherwise, if false, use "//&
1165  "Pedlosky's definition. These definitions differ by a factor "//&
1166  "of 2 in front of the beta term in the denominator. Gill's "//&
1167  "is the more appropriate definition.", default=.false.)
1168  if (gill_equatorial_ld) then
1169  oneortwo = 2.0
1170  endif
1171 
1172  do j=js-1,jeq ; do i=is-1,ieq
1173  cs%f2_dx2_q(i,j) = (g%dxBu(i,j)**2 + g%dyBu(i,j)**2) * &
1174  max(g%CoriolisBu(i,j)**2, absurdly_small_freq**2)
1175  cs%beta_dx2_q(i,j) = oneortwo * ((g%dxBu(i,j))**2 + (g%dyBu(i,j))**2) * (sqrt(0.5 * &
1176  ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1177  ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1178  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1179  ((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2) ) ))
1180  enddo ; enddo
1181 
1182  do j=js,je ; do i=is-1,ieq
1183  cs%f2_dx2_u(i,j) = (g%dxCu(i,j)**2 + g%dyCu(i,j)**2) * &
1184  max(0.5* (g%CoriolisBu(i,j)**2+g%CoriolisBu(i,j-1)**2), absurdly_small_freq**2)
1185  cs%beta_dx2_u(i,j) = oneortwo * ((g%dxCu(i,j))**2 + (g%dyCu(i,j))**2) * (sqrt( &
1186  0.25*( (((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2 + &
1187  ((g%CoriolisBu(i+1,j)-g%CoriolisBu(i,j)) * g%IdxCv(i+1,j))**2) + &
1188  (((g%CoriolisBu(i+1,j-1)-g%CoriolisBu(i,j-1)) * g%IdxCv(i+1,j-1))**2 + &
1189  ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2) ) + &
1190  ((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 ))
1191  enddo ; enddo
1192 
1193  do j=js-1,jeq ; do i=is,ie
1194  cs%f2_dx2_v(i,j) = ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * &
1195  max(0.5*(g%CoriolisBu(i,j)**2+g%CoriolisBu(i-1,j)**2), absurdly_small_freq**2)
1196  cs%beta_dx2_v(i,j) = oneortwo * ((g%dxCv(i,j))**2 + (g%dyCv(i,j))**2) * (sqrt( &
1197  ((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1198  0.25*( (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1199  ((g%CoriolisBu(i-1,j+1)-g%CoriolisBu(i-1,j)) * g%IdyCu(i-1,j+1))**2) + &
1200  (((g%CoriolisBu(i,j+1)-g%CoriolisBu(i,j)) * g%IdyCu(i,j+1))**2 + &
1201  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1202  enddo ; enddo
1203 
1204  endif
1205 
1206  if (cs%Depth_scaled_KhTh) then
1207  cs%calculate_depth_fns = .true.
1208  allocate(cs%Depth_fn_u(isdb:iedb,jsd:jed)) ; cs%Depth_fn_u(:,:) = 0.0
1209  allocate(cs%Depth_fn_v(isd:ied,jsdb:jedb)) ; cs%Depth_fn_v(:,:) = 0.0
1210  call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", cs%depth_scaled_khth_h0, &
1211  "The depth above which KHTH is scaled away.",&
1212  units="m", default=1000.)
1213  call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", cs%depth_scaled_khth_exp, &
1214  "The exponent used in the depth dependent scaling function for KHTH.",&
1215  units="nondim", default=3.0)
1216  endif
1217 
1218  ! Resolution %Rd_dx_h
1219  cs%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, time, &
1220  'Ratio between deformation radius and grid spacing', 'm m-1')
1221  cs%calculate_Rd_dx = cs%calculate_Rd_dx .or. (cs%id_Rd_dx>0)
1222 
1223  if (cs%calculate_Rd_dx) then
1224  cs%calculate_cg1 = .true. ! We will need %cg1
1225  allocate(cs%Rd_dx_h(isd:ied,jsd:jed)) ; cs%Rd_dx_h(:,:) = 0.0
1226  allocate(cs%beta_dx2_h(isd:ied,jsd:jed)); cs%beta_dx2_h(:,:) = 0.0
1227  allocate(cs%f2_dx2_h(isd:ied,jsd:jed)) ; cs%f2_dx2_h(:,:) = 0.0
1228  do j=js-1,je+1 ; do i=is-1,ie+1
1229  cs%f2_dx2_h(i,j) = (g%dxT(i,j)**2 + g%dyT(i,j)**2) * &
1230  max(0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1231  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)), &
1232  absurdly_small_freq**2)
1233  cs%beta_dx2_h(i,j) = oneortwo * ((g%dxT(i,j))**2 + (g%dyT(i,j))**2) * (sqrt(0.5 * &
1234  ( (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
1235  ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
1236  (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
1237  ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ) ))
1238  enddo ; enddo
1239  endif
1240 
1241  if (cs%calculate_cg1) then
1242  in_use = .true.
1243  allocate(cs%cg1(isd:ied,jsd:jed)) ; cs%cg1(:,:) = 0.0
1244  call wave_speed_init(cs%wave_speed_CSp, use_ebt_mode=cs%Resoln_use_ebt, mono_n2_depth=n2_filter_depth)
1245  endif
1246 
1247  ! Leith parameters
1248  call get_param(param_file, mdl, "USE_QG_LEITH_GM", cs%use_QG_Leith_GM, &
1249  "If true, use the QG Leith viscosity as the GM coefficient.", &
1250  default=.false.)
1251 
1252  if (cs%Use_QG_Leith_GM) then
1253  call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &
1254  "The nondimensional Laplacian Leith constant, \n"//&
1255  "often set to 1.0", units="nondim", default=0.0)
1256 
1257  call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_QG_Leith, &
1258  "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
1259  default=.true.)
1260 
1261  alloc_(cs%Laplac3_const_u(isdb:iedb,jsd:jed)) ; cs%Laplac3_const_u(:,:) = 0.0
1262  alloc_(cs%Laplac3_const_v(isd:ied,jsdb:jedb)) ; cs%Laplac3_const_v(:,:) = 0.0
1263  alloc_(cs%KH_u_QG(isdb:iedb,jsd:jed,g%ke)) ; cs%KH_u_QG(:,:,:) = 0.0
1264  alloc_(cs%KH_v_QG(isd:ied,jsdb:jedb,g%ke)) ; cs%KH_v_QG(:,:,:) = 0.0
1265  ! register diagnostics
1266 
1267  cs%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, time, &
1268  'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1269  cs%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, time, &
1270  'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1271 
1272  do j=jsq,jeq+1 ; do i=is-1,ieq
1273  ! Static factors in the Leith schemes
1274  grid_sp_u2 = g%dyCu(i,j)*g%dxCu(i,j)
1275  grid_sp_u3 = sqrt(grid_sp_u2)
1276  cs%Laplac3_const_u(i,j) = leith_lap_const * grid_sp_u3
1277  enddo ; enddo
1278  do j=js-1,jeq ; do i=isq,ieq+1
1279  ! Static factors in the Leith schemes
1280  !### The second factor here is wrong. It should be G%dxCv(i,J).
1281  grid_sp_v2 = g%dyCv(i,j)*g%dxCu(i,j)
1282  grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
1283  cs%Laplac3_const_v(i,j) = leith_lap_const * grid_sp_v3
1284  enddo ; enddo
1285 
1286  if (.not. cs%use_stored_slopes) call mom_error(fatal, &
1287  "MOM_lateral_mixing_coeffs.F90, VarMix_init:"//&
1288  "USE_STORED_SLOPES must be True when using QG Leith.")
1289  endif
1290 
1291  ! If nothing is being stored in this class then deallocate
1292  if (in_use) then
1293  cs%use_variable_mixing = .true.
1294  else
1295  deallocate(cs)
1296  return
1297  endif
1298 

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

Here is the call graph for this function: