src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Variable mixing coefficients
6		module MOM_lateral_mixing_coeffs
7
8		use MOM_debugging, only : hchksum, uvchksum
9		use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg
10		use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, post_data
11		use MOM_diag_mediator, only : diag_ctrl, time_type, query_averaging_enabled
12		use MOM_domains, only : create_group_pass, do_group_pass
13		use MOM_domains, only : group_pass_type, pass_var, pass_vector
14		use MOM_EOS, only : calculate_density_derivs, EOS_domain
15		use MOM_file_parser, only : get_param, log_version, param_file_type
16		use MOM_interface_heights, only : find_eta, thickness_to_dz
17		use MOM_isopycnal_slopes, only : calc_isoneutral_slopes
18		use MOM_grid, only : ocean_grid_type
19		use MOM_unit_scaling, only : unit_scale_type
20		use MOM_variables, only : thermo_var_ptrs
21		use MOM_verticalGrid, only : verticalGrid_type
22		use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init
23		use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE
24		use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
25		use MOM_MEKE_types, only : MEKE_type
26
27		!$ use omp_lib, only: omp_get_num_devices
28
29		implicit none ; private
30
31		#include <MOM_memory.h>
32
33		!> Variable mixing coefficients
34		type, public :: VarMix_CS
35		logical :: initialized = .false. !< True if this control structure has been initialized.
36		logical :: use_variable_mixing !< If true, use the variable mixing.
37		logical :: Resoln_scaling_used !< If true, a resolution function is used somewhere to scale
38		!! away one of the viscosities or diffusivities when the
39		!! deformation radius is well resolved.
40		logical :: Resoln_scaled_Kh !< If true, scale away the Laplacian viscosity
41		!! when the deformation radius is well resolved.
42		logical :: Resoln_scaled_KhTh !< If true, scale away the thickness diffusivity
43		!! when the deformation radius is well resolved.
44		logical :: Depth_scaled_KhTh !< If true, KHTH is scaled away when the depth is
45		!! shallower than a reference depth.
46		logical :: Resoln_scaled_KhTr !< If true, scale away the tracer diffusivity
47		!! when the deformation radius is well resolved.
48		logical :: interpolate_Res_fn !< If true, interpolate the resolution function
49		!! to the velocity points from the thickness
50		!! points; otherwise interpolate the wave
51		!! speed and calculate the resolution function
52		!! independently at each point.
53		logical :: use_stored_slopes !< If true, stores isopycnal slopes in this structure.
54		logical :: Resoln_use_ebt !< If true, use the equivalent barotropic wave speed instead of the
55		!! first baroclinic wave speed for calculating the resolution function.
56		logical :: khth_use_ebt_struct !< If true, uses the equivalent barotropic structure
57		!! as the vertical structure of thickness diffusivity.
58		logical :: kdgl90_use_ebt_struct !< If true, uses the equivalent barotropic structure
59		!! as the vertical structure of diffusivity in the GL90 scheme.
60		logical :: kdgl90_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
61		!! as the vertical structure of diffusivity in the GL90 scheme.
62		logical :: khth_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
63		!! as the vertical structure of thickness diffusivity.
64		logical :: khtr_use_ebt_struct !< If true, uses the equivalent barotropic structure
65		!! as the vertical structure of tracer diffusivity.
66		logical :: khtr_use_sqg_struct !< If true, uses the surface quasigeostrophic structure
67		!! as the vertical structure of tracer diffusivity.
68		logical :: calculate_cg1 !< If true, calls wave_speed() to calculate the first
69		!! baroclinic wave speed and populate CS%cg1.
70		!! This parameter is set depending on other parameters.
71		logical :: calculate_Rd_dx !< If true, calculates Rd/dx and populate CS%Rd_dx_h.
72		!! This parameter is set depending on other parameters.
73		logical :: calculate_res_fns !< If true, calculate all the resolution factors.
74		!! This parameter is set depending on other parameters.
75		logical :: calculate_depth_fns !< If true, calculate all the depth factors.
76		!! This parameter is set depending on other parameters.
77		logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rates.
78		!! This parameter is set depending on other parameters.
79		logical :: use_stanley_iso !< If true, use Stanley parameterization in MOM_isopycnal_slopes
80		logical :: use_simpler_Eady_growth_rate !< If true, use a simpler method to calculate the
81		!! Eady growth rate that avoids division by layer thickness.
82		!! This parameter is set depending on other parameters.
83		logical :: full_depth_Eady_growth_rate !< If true, calculate the Eady growth rate based on an
84		!! average that includes contributions from sea-level changes
85		!! in its denominator, rather than just the nominal depth of
86		!! the bathymetry. This only applies when using the model
87		!! interface heights as a proxy for isopycnal slopes.
88		logical :: OBC_friendly !< If true, use only interior data for thickness weighting and
89		!! to calculate stratification and other fields at open boundary
90		!! condition faces.
91		logical :: res_fn_OBC_bug !< If false, use only interior data for calculating the resolution
92		!! functions at open boundary condition faces and vertices.
93		real :: cropping_distance !< Distance from surface or bottom to filter out outcropped or
94		!! incropped interfaces for the Eady growth rate calc [Z ~> m]
95		real :: h_min_N2 !< The minimum vertical distance to use in the denominator of the
96		!! buoyancy frequency used in the slope calculation [H ~> m or kg m-2]
97
98		real, allocatable :: SN_u(:,:) !< S*N at u-points [T-1 ~> s-1]
99		real, allocatable :: SN_v(:,:) !< S*N at v-points [T-1 ~> s-1]
100		real, allocatable :: L2u(:,:) !< Length scale^2 at u-points [L2 ~> m2]
101		real, allocatable :: L2v(:,:) !< Length scale^2 at v-points [L2 ~> m2]
102		real, allocatable :: cg1(:,:) !< The first baroclinic gravity wave speed [L T-1 ~> m s-1].
103		real, allocatable :: Res_fn_h(:,:) !< Non-dimensional function of the ratio the first baroclinic
104		!! deformation radius to the grid spacing at h points [nondim].
105		real, allocatable :: Res_fn_q(:,:) !< Non-dimensional function of the ratio the first baroclinic
106		!! deformation radius to the grid spacing at q points [nondim].
107		real, allocatable :: Res_fn_u(:,:) !< Non-dimensional function of the ratio the first baroclinic
108		!! deformation radius to the grid spacing at u points [nondim].
109		real, allocatable :: Res_fn_v(:,:) !< Non-dimensional function of the ratio the first baroclinic
110		!! deformation radius to the grid spacing at v points [nondim].
111		real, allocatable :: Depth_fn_u(:,:) !< Non-dimensional function of the ratio of the depth to
112		!! a reference depth (maximum 1) at u points [nondim]
113		real, allocatable :: Depth_fn_v(:,:) !< Non-dimensional function of the ratio of the depth to
114		!! a reference depth (maximum 1) at v points [nondim]
115		real, allocatable :: beta_dx2_h(:,:) !< The magnitude of the gradient of the Coriolis parameter
116		!! times the grid spacing squared at h points [L T-1 ~> m s-1].
117		real, allocatable :: beta_dx2_q(:,:) !< The magnitude of the gradient of the Coriolis parameter
118		!! times the grid spacing squared at q points [L T-1 ~> m s-1].
119		real, allocatable :: beta_dx2_u(:,:) !< The magnitude of the gradient of the Coriolis parameter
120		!! times the grid spacing squared at u points [L T-1 ~> m s-1].
121		real, allocatable :: beta_dx2_v(:,:) !< The magnitude of the gradient of the Coriolis parameter
122		!! times the grid spacing squared at v points [L T-1 ~> m s-1].
123		real, allocatable :: f2_dx2_h(:,:) !< The Coriolis parameter squared times the grid
124		!! spacing squared at h [L2 T-2 ~> m2 s-2].
125		real, allocatable :: f2_dx2_q(:,:) !< The Coriolis parameter squared times the grid
126		!! spacing squared at q [L2 T-2 ~> m2 s-2].
127		real, allocatable :: f2_dx2_u(:,:) !< The Coriolis parameter squared times the grid
128		!! spacing squared at u [L2 T-2 ~> m2 s-2].
129		real, allocatable :: f2_dx2_v(:,:) !< The Coriolis parameter squared times the grid
130		!! spacing squared at v [L2 T-2 ~> m2 s-2].
131		real, allocatable :: Rd_dx_h(:,:) !< Deformation radius over grid spacing [nondim]
132
133		real, allocatable :: slope_x(:,:,:) !< Zonal isopycnal slope [Z L-1 ~> nondim]
134		real, allocatable :: slope_y(:,:,:) !< Meridional isopycnal slope [Z L-1 ~> nondim]
135		real, allocatable :: ebt_struct(:,:,:) !< EBT vertical structure to scale diffusivities with [nondim]
136		real, allocatable :: sqg_struct(:,:,:) !< SQG vertical structure to scale diffusivities with [nondim]
137		real, allocatable :: BS_struct(:,:,:) !< Vertical structure function used in backscatter [nondim]
138		real, allocatable :: khth_struct(:,:,:) !< Vertical structure function used in thickness diffusivity [nondim]
139		real, allocatable :: khtr_struct(:,:,:) !< Vertical structure function used in tracer diffusivity [nondim]
140		real, allocatable :: kdgl90_struct(:,:,:) !< Vertical structure function used in GL90 diffusivity [nondim]
141		real :: BS_EBT_power !< Power to raise EBT vertical structure to. Default 0.0.
142		real :: sqg_expo !< Exponent for SQG vertical structure [nondim]. Default 1.0
143		logical :: interpolated_sqg_struct !< If true, interpolate properties to velocity points and then
144		!! interpolate the buoyancy frequencies and layer thicknesses
145		!! back to tracer points when calculating the SQG vertical
146		!! structure.
147		logical :: BS_use_sqg_struct !< If true, use sqg_stuct for backscatter vertical structure.
148
149		real, allocatable :: Laplac3_const_u(:,:) !< Laplacian metric-dependent constants at u-points [L3 ~> m3]
150		real, allocatable :: Laplac3_const_v(:,:) !< Laplacian metric-dependent constants at u-points [L3 ~> m3]
151		real, allocatable :: KH_u_QG(:,:,:) !< QG Leith GM coefficient at u-points [L2 T-1 ~> m2 s-1]
152		real, allocatable :: KH_v_QG(:,:,:) !< QG Leith GM coefficient at v-points [L2 T-1 ~> m2 s-1]
153
154		! Parameters
155		logical :: use_Visbeck !< Use Visbeck formulation for thickness diffusivity
156		integer :: VarMix_Ktop !< Top layer to start downward integrals
157		real :: Visbeck_L_scale !< Fixed length scale in Visbeck formula [L ~> m], or if negative a scaling
158		!! factor [nondim] relating this length scale squared to the cell area
159		real :: Eady_GR_D_scale !< Depth over which to average SN [Z ~> m]
160		real :: Res_coef_khth !< A coefficient [nondim] that determines the function
161		!! of resolution, used for thickness and tracer mixing, as:
162		!! F = 1 / (1 + (Res_coef_khth*Ld/dx)^Res_fn_power)
163		real :: Res_coef_visc !< A coefficient [nondim] that determines the function
164		!! of resolution, used for lateral viscosity, as:
165		!! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power)
166		real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m]
167		real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim]
168		real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
169		integer :: Res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any
170		!! positive integer power may be used, but even powers
171		!! and especially 2 are coded to be more efficient.
172		integer :: Res_fn_power_visc !< The power of dx/Ld in the Kh resolution function. Any
173		!! positive integer power may be used, but even powers
174		!! and especially 2 are coded to be more efficient.
175		real :: Visbeck_S_max !< Upper bound on slope used in Eady growth rate [Z L-1 ~> nondim].
176
177		! Leith parameters
178		logical :: use_QG_Leith_GM !< If true, uses the QG Leith viscosity as the GM coefficient
179		logical :: use_beta_in_QG_Leith !< If true, includes the beta term in the QG Leith GM coefficient
180
181		! Isoneutral blocking parameters
182		integer :: niblock !< The i block size used in calc_isoneutral_slopes [nondim].
183		integer :: njblock !< The j block size used in calc_isoneutral_slopes [nondim].
184		integer :: nkblock !< The k block size used in calc_isoneutral_slopes [nondim].
185
186		! Diagnostics
187		!>@{
188		!! Diagnostic identifier
189		integer :: id_SN_u=-1, id_SN_v=-1, id_L2u=-1, id_L2v=-1, id_Res_fn = -1
190		integer :: id_N2_u=-1, id_N2_v=-1, id_S2_u=-1, id_S2_v=-1
191		integer :: id_dzu=-1, id_dzv=-1, id_dzSxN=-1, id_dzSyN=-1
192		integer :: id_Rd_dx=-1, id_KH_u_QG = -1, id_KH_v_QG = -1
193		integer :: id_sqg_struct=-1, id_BS_struct=-1, id_khth_struct=-1, id_khtr_struct=-1
194		integer :: id_kdgl90_struct=-1
195		type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
196		!! timing of diagnostic output.
197		!>@}
198
199		type(wave_speed_CS) :: wave_speed !< Wave speed control structure
200		type(group_pass_type) :: pass_cg1 !< For group halo pass
201		logical :: debug !< If true, write out checksums of data for debugging
202		end type VarMix_CS
203
204		public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function
205		public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function, calc_sqg_struct
206
207		contains
208
209		!> Calculates the non-dimensional depth functions.
210	12	subroutine calc_depth_function(G, CS)
211		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
212		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
213
214		! Local variables
215		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq
216		integer :: i, j
217		real :: H0 ! The depth above which KHTH is linearly scaled away [Z ~> m]
218		real :: h1, h2 ! Temporary total thicknesses [Z ~> m]
219		real :: expo ! exponent used in the depth dependent scaling [nondim]
220	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
221	12	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
222
223	12	if (.not. CS%initialized) call MOM_error(FATAL, "calc_depth_function: "// &
224	0	"Module must be initialized before it is used.")
225
226	12	if (.not. CS%calculate_depth_fns) return
227	0	if (.not. allocated(CS%Depth_fn_u)) call MOM_error(FATAL, &
228	0	"calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
229	0	if (.not. allocated(CS%Depth_fn_v)) call MOM_error(FATAL, &
230	0	"calc_depth_function: %Depth_fn_v is not associated with Depth_scaled_KhTh.")
231
232		! For efficiency, the reciprocal of H0 should be used instead.
233	0	H0 = CS%depth_scaled_khth_h0
234	0	expo = CS%depth_scaled_khth_exp
235		!$OMP do
236	0	do j=js,je ; do I=is-1,Ieq
237	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
238	0	h2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)
239	0	CS%Depth_fn_u(I,j) = (MIN(1.0, (0.5 * (h1 + h2)) / H0))**expo
240		enddo ; enddo
241		!$OMP do
242	0	do J=js-1,Jeq ; do i=is,ie
243	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
244	0	h2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)
245	0	CS%Depth_fn_v(i,J) = (MIN(1.0, (0.5 * (h1 + h2)) / H0))**expo
246		enddo ; enddo
247
248		end subroutine calc_depth_function
249
250		!> Calculates and stores the non-dimensional resolution functions
251	12	subroutine calc_resoln_function(h, tv, G, GV, US, CS, MEKE, OBC, dt)
252		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
253		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
254		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
255		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
256		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
257		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
258		type(MEKE_type), intent(in) :: MEKE !< MEKE struct
259		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
260		real, intent(in) :: dt !< Time increment [T ~> s]
261
262		! Local variables
263		! Depending on the power-function being used, dimensional rescaling may be limited, so some
264		! of the following variables have units that depend on that power.
265	12	real :: cg1_q(SZIB_(G),SZJB_(G)) ! The gravity wave speed interpolated to q points [L T-1 ~> m s-1] or [m s-1].
266	24	real :: cg1_u(SZIB_(G),SZJ_(G)) ! The gravity wave speed interpolated to u points [L T-1 ~> m s-1] or [m s-1].
267	12	real :: cg1_v(SZI_(G),SZJB_(G)) ! The gravity wave speed interpolated to v points [L T-1 ~> m s-1] or [m s-1].
268		real :: dx_term ! A term in the denominator [L2 T-2 ~> m2 s-2] or [m2 s-2]
269		logical :: apply_u_OBC, apply_v_OBC ! If true, OBCs will be used to set the wave speed at some points on this PE.
270		integer :: power_2
271		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
272		integer :: i, j, k
273	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
274	12	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
275
276	12	if (.not. CS%initialized) call MOM_error(FATAL, "calc_resoln_function: "// &
277	0	"Module must be initialized before it is used.")
278
279	12	if (CS%calculate_cg1) then
280		!$omp target update from(h)
281	12	if (.not. allocated(CS%cg1)) call MOM_error(FATAL, &
282	0	"calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
283		if (CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct &
284	12	.or. CS%khtr_use_ebt_struct .or. CS%BS_EBT_power>0.) then
285	0	if (.not. allocated(CS%ebt_struct)) call MOM_error(FATAL, &
286	0	"calc_resoln_function: %ebt_struct is not associated with RESOLN_USE_EBT.")
287	0	if (CS%Resoln_use_ebt) then
288		! Both resolution fn and vertical structure are using EBT
289	0	call wave_speed(h, tv, G, GV, US, CS%cg1, CS%wave_speed, modal_structure=CS%ebt_struct)
290		else
291		! Use EBT to get vertical structure first and then re-calculate cg1 using first baroclinic mode
292		call wave_speed(h, tv, G, GV, US, CS%cg1, CS%wave_speed, modal_structure=CS%ebt_struct, &
293	0	use_ebt_mode=.true.)
294	0	call wave_speed(h, tv, G, GV, US, CS%cg1, CS%wave_speed)
295		endif
296	0	call pass_var(CS%ebt_struct, G%Domain)
297		else
298	12	call wave_speed(h, tv, G, GV, US, CS%cg1, CS%wave_speed)
299		endif
300
301	12	call create_group_pass(CS%pass_cg1, CS%cg1, G%Domain)
302	12	call do_group_pass(CS%pass_cg1, G%Domain)
303		endif
304
305		if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
306	12	.or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then
307		!$omp target update from(h)
308	0	call calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE, OBC)
309	0	call pass_var(CS%sqg_struct, G%Domain)
310		endif
311
312	12	if (CS%BS_EBT_power>0.) then
313	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
314	0	CS%BS_struct(i,j,k) = CS%ebt_struct(i,j,k)**CS%BS_EBT_power
315		enddo ; enddo ; enddo
316	12	elseif (CS%BS_use_sqg_struct) then
317	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
318	0	CS%BS_struct(i,j,k) = CS%sqg_struct(i,j,k)
319		enddo ; enddo ; enddo
320		endif
321
322	12	if (CS%khth_use_ebt_struct) then
323	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
324	0	CS%khth_struct(i,j,k) = CS%ebt_struct(i,j,k)
325		enddo ; enddo ; enddo
326	12	elseif (CS%khth_use_sqg_struct) then
327	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
328	0	CS%khth_struct(i,j,k) = CS%sqg_struct(i,j,k)
329		enddo ; enddo ; enddo
330		endif
331
332	12	if (CS%khtr_use_ebt_struct) then
333	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
334	0	CS%khtr_struct(i,j,k) = CS%ebt_struct(i,j,k)
335		enddo ; enddo ; enddo
336	12	elseif (CS%khtr_use_sqg_struct) then
337	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
338	0	CS%khtr_struct(i,j,k) = CS%sqg_struct(i,j,k)
339		enddo ; enddo ; enddo
340		endif
341
342	12	if (CS%kdgl90_use_ebt_struct) then
343	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
344	0	CS%kdgl90_struct(i,j,k) = CS%ebt_struct(i,j,k)
345		enddo ; enddo ; enddo
346	12	elseif (CS%kdgl90_use_sqg_struct) then
347	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
348	0	CS%kdgl90_struct(i,j,k) = CS%sqg_struct(i,j,k)
349		enddo ; enddo ; enddo
350		endif
351
352		! Calculate and store the ratio between deformation radius and grid-spacing
353		! at h-points [nondim].
354	12	if (CS%calculate_rd_dx) then
355	12	if (.not. allocated(CS%Rd_dx_h)) call MOM_error(FATAL, &
356	0	"calc_resoln_function: %Rd_dx_h is not associated with calculate_rd_dx.")
357		!$OMP parallel do default(shared)
358	91524	do j=js-1,je+1 ; do i=is-1,ie+1
359		CS%Rd_dx_h(i,j) = CS%cg1(i,j) / &
360	91512	(sqrt(CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j)))
361		enddo ; enddo
362	12	if (query_averaging_enabled(CS%diag)) then
363	12	if (CS%id_Rd_dx > 0) call post_data(CS%id_Rd_dx, CS%Rd_dx_h, CS%diag)
364		endif
365		endif
366
367	12	if (.not. CS%calculate_res_fns) return
368
369	12	if (.not. allocated(CS%Res_fn_h)) call MOM_error(FATAL, &
370	0	"calc_resoln_function: %Res_fn_h is not associated with Resoln_scaled_Kh.")
371	12	if (.not. allocated(CS%Res_fn_q)) call MOM_error(FATAL, &
372	0	"calc_resoln_function: %Res_fn_q is not associated with Resoln_scaled_Kh.")
373	12	if (.not. allocated(CS%Res_fn_u)) call MOM_error(FATAL, &
374	0	"calc_resoln_function: %Res_fn_u is not associated with Resoln_scaled_Kh.")
375	12	if (.not. allocated(CS%Res_fn_v)) call MOM_error(FATAL, &
376	0	"calc_resoln_function: %Res_fn_v is not associated with Resoln_scaled_Kh.")
377	12	if (.not. allocated(CS%f2_dx2_h)) call MOM_error(FATAL, &
378	0	"calc_resoln_function: %f2_dx2_h is not associated with Resoln_scaled_Kh.")
379	12	if (.not. allocated(CS%f2_dx2_q)) call MOM_error(FATAL, &
380	0	"calc_resoln_function: %f2_dx2_q is not associated with Resoln_scaled_Kh.")
381	12	if (.not. allocated(CS%f2_dx2_u)) call MOM_error(FATAL, &
382	0	"calc_resoln_function: %f2_dx2_u is not associated with Resoln_scaled_Kh.")
383	12	if (.not. allocated(CS%f2_dx2_v)) call MOM_error(FATAL, &
384	0	"calc_resoln_function: %f2_dx2_v is not associated with Resoln_scaled_Kh.")
385	12	if (.not. allocated(CS%beta_dx2_h)) call MOM_error(FATAL, &
386	0	"calc_resoln_function: %beta_dx2_h is not associated with Resoln_scaled_Kh.")
387	12	if (.not. allocated(CS%beta_dx2_q)) call MOM_error(FATAL, &
388	0	"calc_resoln_function: %beta_dx2_q is not associated with Resoln_scaled_Kh.")
389	12	if (.not. allocated(CS%beta_dx2_u)) call MOM_error(FATAL, &
390	0	"calc_resoln_function: %beta_dx2_u is not associated with Resoln_scaled_Kh.")
391	12	if (.not. allocated(CS%beta_dx2_v)) call MOM_error(FATAL, &
392	0	"calc_resoln_function: %beta_dx2_v is not associated with Resoln_scaled_Kh.")
393
394	12	apply_u_OBC = .false. ; apply_v_OBC = .false.
395	12	if (associated(OBC) .and. (.not.CS%res_fn_OBC_bug)) then
396	0	apply_u_OBC = OBC%u_OBCs_on_PE
397	0	apply_v_OBC = OBC%v_OBCs_on_PE
398		endif
399
400		!$OMP parallel default(shared) private(dx_term,power_2)
401
402	12	if (apply_u_OBC .or. apply_v_OBC) then
403		!$OMP do
404	0	do J=js-1,Jeq ; do I=is-1,Ieq
405		if ((OBC%segnum_u(I,j) /= 0) .or. (OBC%segnum_u(I,j+1) /= 0) .or. &
406	0	(OBC%segnum_v(i,J) /= 0) .or. (OBC%segnum_u(i+1,J) /= 0)) then
407		! This is an OBC node, so use the fact that G%mask2dT is zero behind OBCs. The nondimensional
408		! constant 1e-20 in the denominator makes this a de facto implementation of Adcroft's reciprocal
409		! rule with a value that works for either 64-bit or 32-bit real numbers.
410		cg1_q(I,J) = ((G%mask2dT(i,j) * CS%cg1(i,j) + G%mask2dT(i+1,j+1) * CS%cg1(i+1,j+1)) + &
411		(G%mask2dT(i+1,j) * CS%cg1(i+1,j) + G%mask2dT(i,j+1) * CS%cg1(i,j+1))) / &
412	0	((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-20)
413		else
414	0	cg1_q(I,J) = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + (CS%cg1(i+1,j) + CS%cg1(i,j+1)))
415		endif
416		enddo ; enddo
417		else
418		!$OMP do
419	89316	do J=js-1,Jeq ; do I=is-1,Ieq
420	89304	cg1_q(I,J) = 0.25 * ((CS%cg1(i,j) + CS%cg1(i+1,j+1)) + (CS%cg1(i+1,j) + CS%cg1(i,j+1)))
421		enddo ; enddo
422		endif
423
424		! Do this calculation on the extent used in MOM_hor_visc.F90, and
425		! MOM_tracer.F90 so that no halo update is needed.
426	12	if (CS%Res_fn_power_visc >= 100) then
427		!$OMP do
428	0	do j=js-1,je+1 ; do i=is-1,ie+1
429	0	dx_term = CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j)
430	0	if ((CS%Res_coef_visc * CS%cg1(i,j))**2 > dx_term) then
431	0	CS%Res_fn_h(i,j) = 0.0
432		else
433	0	CS%Res_fn_h(i,j) = 1.0
434		endif
435		enddo ; enddo
436		!$OMP do
437	0	do J=js-1,Jeq ; do I=is-1,Ieq
438	0	dx_term = CS%f2_dx2_q(I,J) + cg1_q(I,J) * CS%beta_dx2_q(I,J)
439	0	if ((CS%Res_coef_visc * cg1_q(I,J))**2 > dx_term) then
440	0	CS%Res_fn_q(I,J) = 0.0
441		else
442	0	CS%Res_fn_q(I,J) = 1.0
443		endif
444		enddo ; enddo
445	12	elseif (CS%Res_fn_power_visc == 2) then
446		!$OMP do
447	91524	do j=js-1,je+1 ; do i=is-1,ie+1
448	90768	dx_term = CS%f2_dx2_h(i,j) + CS%cg1(i,j)*CS%beta_dx2_h(i,j)
449	91512	CS%Res_fn_h(i,j) = dx_term / (dx_term + (CS%Res_coef_visc * CS%cg1(i,j))**2)
450		enddo ; enddo
451		!$OMP do
452	89316	do J=js-1,Jeq ; do I=is-1,Ieq
453	88572	dx_term = CS%f2_dx2_q(I,J) + cg1_q(I,J) * CS%beta_dx2_q(I,J)
454	89304	CS%Res_fn_q(I,J) = dx_term / (dx_term + (CS%Res_coef_visc * cg1_q(I,J))**2)
455		enddo ; enddo
456	0	elseif (mod(CS%Res_fn_power_visc, 2) == 0) then
457	0	power_2 = CS%Res_fn_power_visc / 2
458		!$OMP do
459	0	do j=js-1,je+1 ; do i=is-1,ie+1
460	0	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
461		CS%Res_fn_h(i,j) = dx_term / &
462	0	(dx_term + (CS%Res_coef_visc * US%L_T_to_m_sCS%cg1(i,j))*CS%Res_fn_power_visc)
463		enddo ; enddo
464		!$OMP do
465	0	do J=js-1,Jeq ; do I=is-1,Ieq
466	0	dx_term = (US%L_T_to_m_s*2(CS%f2_dx2_q(I,J) + cg1_q(I,J) * CS%beta_dx2_q(I,J)))**power_2
467		CS%Res_fn_q(I,J) = dx_term / &
468	0	(dx_term + (CS%Res_coef_visc * US%L_T_to_m_scg1_q(I,J))*CS%Res_fn_power_visc)
469		enddo ; enddo
470		else
471		!$OMP do
472	0	do j=js-1,je+1 ; do i=is-1,ie+1
473		dx_term = (US%L_T_to_m_s*sqrt(CS%f2_dx2_h(i,j) + &
474	0	CS%cg1(i,j)CS%beta_dx2_h(i,j)))*CS%Res_fn_power_visc
475		CS%Res_fn_h(i,j) = dx_term / &
476	0	(dx_term + (CS%Res_coef_visc * US%L_T_to_m_sCS%cg1(i,j))*CS%Res_fn_power_visc)
477		enddo ; enddo
478		!$OMP do
479	0	do J=js-1,Jeq ; do I=is-1,Ieq
480		dx_term = (US%L_T_to_m_s*sqrt(CS%f2_dx2_q(I,J) + &
481	0	cg1_q(I,J) * CS%beta_dx2_q(I,J)))**CS%Res_fn_power_visc
482		CS%Res_fn_q(I,J) = dx_term / &
483	0	(dx_term + (CS%Res_coef_visc * US%L_T_to_m_scg1_q(I,J))*CS%Res_fn_power_visc)
484		enddo ; enddo
485		endif
486
487	12	if (CS%interpolate_Res_fn) then
488	0	if (apply_u_OBC) then
489	0	do j=js,je ; do I=is-1,Ieq
490	0	CS%Res_fn_u(I,j) = 0.5*(CS%Res_fn_h(i,j) + CS%Res_fn_h(i+1,j))
491	0	if (OBC%segnum_u(I,j) > 0) CS%Res_fn_u(I,j) = CS%Res_fn_h(i,j) ! Eastern OBC
492	0	if (OBC%segnum_u(I,j) < 0) CS%Res_fn_u(I,j) = CS%Res_fn_h(i+1,j) ! Western OBC
493		enddo ; enddo
494		else
495	0	do j=js,je ; do I=is-1,Ieq
496	0	CS%Res_fn_u(I,j) = 0.5*(CS%Res_fn_h(i,j) + CS%Res_fn_h(i+1,j))
497		enddo ; enddo
498		endif
499
500	0	if (apply_v_OBC) then
501	0	do J=js-1,Jeq ; do i=is,ie
502	0	CS%Res_fn_v(i,J) = 0.5*(CS%Res_fn_h(i,j) + CS%Res_fn_h(i,j+1))
503	0	if (OBC%segnum_v(i,J) > 0) CS%Res_fn_v(i,J) = CS%Res_fn_h(i,j) ! Northern OBC
504	0	if (OBC%segnum_v(i,J) < 0) CS%Res_fn_v(i,J) = CS%Res_fn_h(i,j+1) ! Southern OBC
505		enddo ; enddo
506		else
507	0	do J=js-1,Jeq ; do i=is,ie
508	0	CS%Res_fn_v(i,J) = 0.5*(CS%Res_fn_h(i,j) + CS%Res_fn_h(i,j+1))
509		enddo ; enddo
510		endif
511
512		else ! .not.CS%interpolate_Res_fn
513	12	if (apply_u_OBC) then
514		!$OMP do
515	0	do j=js,je ; do I=is-1,Ieq
516	0	cg1_u(I,j) = 0.5 * (CS%cg1(i,j) + CS%cg1(i+1,j))
517	0	if (OBC%segnum_u(I,j) > 0) cg1_u(I,j) = CS%cg1(i,j) ! Eastern OBC
518	0	if (OBC%segnum_u(I,j) < 0) cg1_u(I,j) = CS%cg1(i+1,j) ! Western OBC
519		enddo ; enddo
520		else
521		!$OMP do
522	87852	do j=js,je ; do I=is-1,Ieq
523	87840	cg1_u(I,j) = 0.5 * (CS%cg1(i,j) + CS%cg1(i+1,j))
524		enddo ; enddo
525		endif
526
527	12	if (apply_v_OBC) then
528		!$OMP do
529	0	do J=js-1,Jeq ; do i=is,ie
530	0	cg1_v(i,J) = 0.5 * (CS%cg1(i,j) + CS%cg1(i,j+1))
531	0	if (OBC%segnum_v(i,J) > 0) cg1_v(i,J) = CS%cg1(i,j) ! Northern OBC
532	0	if (OBC%segnum_v(i,J) < 0) cg1_v(i,J) = CS%cg1(i,j+1) ! Southern OBC
533		enddo ; enddo
534		else
535		!$OMP do
536	88584	do J=js-1,Jeq ; do i=is,ie
537	88572	cg1_v(i,J) = 0.5 * (CS%cg1(i,j) + CS%cg1(i,j+1))
538		enddo ; enddo
539		endif
540
541	12	if (CS%Res_fn_power_khth >= 100) then
542		!$OMP do
543	0	do j=js,je ; do I=is-1,Ieq
544	0	dx_term = CS%f2_dx2_u(I,j) + cg1_u(I,j) * CS%beta_dx2_u(I,j)
545	0	if ((CS%Res_coef_khth * cg1_u(I,j))**2 > dx_term) then
546	0	CS%Res_fn_u(I,j) = 0.0
547		else
548	0	CS%Res_fn_u(I,j) = 1.0
549		endif
550		enddo ; enddo
551		!$OMP do
552	0	do J=js-1,Jeq ; do i=is,ie
553	0	dx_term = CS%f2_dx2_v(i,J) + cg1_v(i,J) * CS%beta_dx2_v(i,J)
554	0	if ((CS%Res_coef_khth * cg1_v(i,J))**2 > dx_term) then
555	0	CS%Res_fn_v(i,J) = 0.0
556		else
557	0	CS%Res_fn_v(i,J) = 1.0
558		endif
559		enddo ; enddo
560	12	elseif (CS%Res_fn_power_khth == 2) then
561		!$OMP do
562	87852	do j=js,je ; do I=is-1,Ieq
563	87120	dx_term = CS%f2_dx2_u(I,j) + cg1_u(I,j) * CS%beta_dx2_u(I,j)
564	87840	CS%Res_fn_u(I,j) = dx_term / (dx_term + (CS%Res_coef_khth * cg1_u(I,j))**2)
565		enddo ; enddo
566		!$OMP do
567	88584	do J=js-1,Jeq ; do i=is,ie
568	87840	dx_term = CS%f2_dx2_v(i,J) + cg1_v(i,J) * CS%beta_dx2_v(i,J)
569	88572	CS%Res_fn_v(i,J) = dx_term / (dx_term + (CS%Res_coef_khth * cg1_v(i,J))**2)
570		enddo ; enddo
571	0	elseif (mod(CS%Res_fn_power_khth, 2) == 0) then
572	0	power_2 = CS%Res_fn_power_khth / 2
573		!$OMP do
574	0	do j=js,je ; do I=is-1,Ieq
575	0	dx_term = (US%L_T_to_m_s*2 (CS%f2_dx2_u(I,j) + cg1_u(I,j) * CS%beta_dx2_u(I,j)))**power_2
576		CS%Res_fn_u(I,j) = dx_term / &
577	0	(dx_term + (CS%Res_coef_khth * US%L_T_to_m_scg1_u(I,j))*CS%Res_fn_power_khth)
578		enddo ; enddo
579		!$OMP do
580	0	do J=js-1,Jeq ; do i=is,ie
581	0	dx_term = (US%L_T_to_m_s*2 (CS%f2_dx2_v(i,J) + cg1_v(i,J) * CS%beta_dx2_v(i,J)))**power_2
582		CS%Res_fn_v(i,J) = dx_term / &
583	0	(dx_term + (CS%Res_coef_khth * US%L_T_to_m_scg1_v(i,J))*CS%Res_fn_power_khth)
584		enddo ; enddo
585		else
586		!$OMP do
587	0	do j=js,je ; do I=is-1,Ieq
588		dx_term = (US%L_T_to_m_s*sqrt(CS%f2_dx2_u(I,j) + &
589	0	cg1_u(I,j) * CS%beta_dx2_u(I,j)))**CS%Res_fn_power_khth
590		CS%Res_fn_u(I,j) = dx_term / &
591	0	(dx_term + (CS%Res_coef_khth * US%L_T_to_m_scg1_u(I,j))*CS%Res_fn_power_khth)
592		enddo ; enddo
593		!$OMP do
594	0	do J=js-1,Jeq ; do i=is,ie
595		dx_term = (US%L_T_to_m_s*sqrt(CS%f2_dx2_v(i,J) + &
596	0	cg1_v(i,J) * CS%beta_dx2_v(i,J)))**CS%Res_fn_power_khth
597		CS%Res_fn_v(i,J) = dx_term / &
598	0	(dx_term + (CS%Res_coef_khth * US%L_T_to_m_scg1_v(i,J))*CS%Res_fn_power_khth)
599		enddo ; enddo
600		endif
601		endif
602		!$OMP end parallel
603
604	12	if (query_averaging_enabled(CS%diag)) then
605	12	if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag)
606	12	if (CS%id_BS_struct > 0) call post_data(CS%id_BS_struct, CS%BS_struct, CS%diag)
607	12	if (CS%id_khth_struct > 0) call post_data(CS%id_khth_struct, CS%khth_struct, CS%diag)
608	12	if (CS%id_khtr_struct > 0) call post_data(CS%id_khtr_struct, CS%khtr_struct, CS%diag)
609	12	if (CS%id_kdgl90_struct > 0) call post_data(CS%id_kdgl90_struct, CS%kdgl90_struct, CS%diag)
610		endif
611
612	12	if (CS%debug) then
613	0	call hchksum(CS%cg1, "calc_resoln_fn cg1", G%HI, haloshift=1, unscale=US%L_T_to_m_s)
614		call uvchksum("Res_fn_[uv]", CS%Res_fn_u, CS%Res_fn_v, G%HI, haloshift=0, &
615	0	unscale=1.0, scalar_pair=.true.)
616		endif
617		end subroutine calc_resoln_function
618
619		!> Calculates and stores functions of SQG mode
620	0	subroutine calc_sqg_struct(h, tv, G, GV, US, CS, dt, MEKE, OBC)
621		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
622		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
623		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
624		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
625		type(thermo_var_ptrs), intent(in) :: tv !<Thermodynamic variables
626		real, intent(in) :: dt !< Time increment [T ~> s]
627		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control struct
628		type(MEKE_type), intent(in) :: MEKE !< MEKE struct
629		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
630
631		! Local variables
632	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
633	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2]
634	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2]
635	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
636	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
637	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! \|Sx\| N times dz at u-points [Z T-1 ~> m s-1]
638	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzSyN ! \|Sy\| N times dz at v-points [Z T-1 ~> m s-1]
639	0	real, dimension(SZI_(G),SZJ_(G)) :: f ! Absolute value of the Coriolis parameter at h point [T-1 ~> s-1]
640		real :: N2 ! Positive buoyancy frequency square or zero [L2 Z-2 T-2 ~> s-2]
641		real :: dzc ! Spacing between two adjacent layers in stretched vertical coordinate [Z ~> m]
642		real :: f_subround ! The minimal resolved value of Coriolis parameter to prevent division by zero [T-1 ~> s-1]
643	0	real, dimension(SZI_(G),SZJ_(G)) :: Le ! Eddy length scale [L ~> m]
644
645	0	real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Geometric layer thicknesses in height units [Z ~> m]
646	0	real :: I_f_Le(SZI_(G),SZJ_(G)) ! The inverse of the absolute value of f times the Eddy
647		! length scale [T L-1 ~> s m-1]
648	0	real :: p_i(SZI_(G),SZJ_(G)) ! Pressure at the interface [R L2 T-2 ~> Pa]
649	0	real :: T_i(SZI_(G)) ! Temperature at the interface [C ~> degC]
650	0	real :: S_i(SZI_(G)) ! Salinity at the interface [S ~> ppt]
651	0	real :: dRho_dS(SZI_(G)) ! Local change in density with salinity using the model EOS and
652		! state interpolated to an interface [R C-1 ~> kg m-3 ppt-1]
653	0	real :: dRho_dT(SZI_(G)) ! Local change in density with salinity using the model EOS and
654		! state interpolated [R C-1 ~> kg m-3 degC-1]
655		real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
656		real :: GxSpV ! Gravitiational acceleration times the specific volume at an interface
657		! [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
658		real :: drdk ! Vertical density differences across an interface [R ~> kg m-3]
659		real :: dz_int ! Average of thicknesses around an interface in height units [Z ~> m]
660		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
661		integer :: i, j, k, is, ie, js, je, nz
662		integer :: niblock, njblock, nkblock
663
664	0	niblock = CS%niblock
665	0	njblock = CS%njblock
666	0	nkblock = CS%nkblock
667
668	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
669	0	f_subround = 1.0e-40 * US%s_to_T
670
671	0	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
672	0	"Module must be initialized before it is used.")
673
674	0	if (CS%sqg_expo <= 0.) then
675	0	CS%sqg_struct(:,:,:) = 1.
676		else
677	0	if (allocated(MEKE%Le)) then
678	0	do j=js,je ; do i=is,ie
679	0	Le(i,j) = MEKE%Le(i,j)
680		enddo ; enddo
681		else
682	0	do j=js,je ; do i=is,ie
683	0	Le(i,j) = sqrt(G%areaT(i,j))
684		enddo ; enddo
685		endif
686
687	0	do j=js,je ; do i=is,ie
688		! Setting the structure averaged over the top layer to 1 is consistent with it being well mixed.
689	0	CS%sqg_struct(i,j,1) = 1.0
690		enddo ; enddo
691
692	0	if (CS%interpolated_sqg_struct) then
693	0	do j=js,je ; do i=is,ie
694		f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
695	0	(G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround)
696		enddo ; enddo
697		!$omp target update to(h)
698		!$omp target enter data map(alloc: e)
699	0	call find_eta(h, tv, G, GV, US, e, halo_size=2) !### Could be halo_size=1?
700
701	0	if (niblock == 0) niblock = ie - is + 1
702	0	if (njblock == 0) njblock = je - js + 1
703	0	if (nkblock == 0) nkblock = nz
704
705		!$omp target enter data map(to: tv, tv%T, tv%S)
706		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
707		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
708		!$omp target enter data map(alloc: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
709		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
710		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
711		N2_v=N2_v, dzu=dzu, dzv=dzv, dzSxN=dzSxN, dzSyN=dzSyN, halo=1, &
712	0	OBC=OBC, OBC_N2=CS%OBC_friendly)
713		!$omp target exit data map(release: tv%T, tv%S, tv)
714		!$omp target exit data map(delete: e)
715		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
716		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
717		!$omp target exit data map(from: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
718	0	do k=2,nz ; do j=js,je ; do i=is,ie
719	0	N2 = max(0.25 * ((N2_u(I-1,j,K) + N2_u(I,j,K)) + (N2_v(i,J-1,K) + N2_v(i,J,K))), 0.0)
720	0	dzc = 0.25 * ((dzu(I-1,j,K) + dzu(I,j,K)) + (dzv(i,J-1,K) + dzv(i,J,K)))
721		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
722	0	exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j))))
723		enddo ; enddo ; enddo
724		else
725	0	do j=js,je ; do i=is,ie
726		I_f_Le(i,j) = 1.0 / &
727		(Le(i,j) * max(0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
728	0	(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1)))), f_subround))
729		enddo ; enddo
730
731	0	call thickness_to_dz(h, tv, dz, G, GV, US)
732
733	0	if (associated(tv%eqn_of_state)) then
734	0	EOSdom(:) = EOS_domain(G%HI)
735	0	H_to_pres = GV%H_to_RZ * GV%g_Earth
736		! Set the pressure at the topmost interior interface.
737	0	p_i(:,:) = 0.0
738	0	if (associated(tv%p_surf)) then
739	0	do j=js,je ; do i=is,ie ; p_i(i,j) = tv%p_surf(i,j) ; enddo ; enddo
740		endif
741	0	if (.not.allocated(tv%SpV_avg)) GxSpV = GV%g_Earth / GV%Rho0
742	0	do K=2,nz ; do j=js,je
743		! Find the derivatives of density with T and S at the interface.
744	0	do i=is,ie
745	0	p_i(i,j) = p_i(i,j) + H_to_pres * h(i,j,k-1)
746	0	T_i(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
747	0	S_i(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
748		enddo
749	0	call calculate_density_derivs(T_i, S_i, p_i(:,j), dRho_dT, dRho_dS, tv%eqn_of_state, EOSdom)
750
751	0	do i=is,ie
752	0	if (allocated(tv%SpV_avg)) & ! GxSpV is in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
753	0	GxSpV = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
754
755		drdk = max(dRho_dT(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
756	0	dRho_dS(i) * (tv%S(i,j,k)-tv%S(i,j,k-1)), 0.0) ! Density difference [R ~> kg m-3]
757	0	dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
758		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
759	0	exp(-CS%sqg_expo * (sqrt((GxSpV * drdk) * dz_int) * I_f_Le(i,j)) )
760		! To derive the expression above, note that
761		! N2 = GxSpV * drdk / dzh(i,j,K) ! Square of positive buoyancy freq. [L2 Z-2 T-2 ~> s-2]
762		! CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
763		! exp(-CS%sqg_expo * (dz_int(i,j,K) * sqrt(N2) * I_f_Le(i,j)) )
764		enddo
765		enddo ; enddo
766		else ! (GV%Boussinesq .and. .not.use_EOS) then
767	0	do K=2,nz ; do j=js,je ; do i=is,ie
768	0	dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
769		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
770	0	exp(-CS%sqg_expo * (sqrt(GV%g_prime(K) * dz_int) * I_f_Le(i,j)) )
771		enddo ; enddo ; enddo
772		endif
773		endif
774		endif
775
776	0	if (query_averaging_enabled(CS%diag)) then
777	0	if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag)
778	0	if (CS%interpolated_sqg_struct .and. (CS%sqg_expo > 0.)) then
779	0	if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
780	0	if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
781		endif
782		endif
783
784	0	end subroutine calc_sqg_struct
785
786		!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
787		!! style scaling of diffusivity
788	12	subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
789		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
790		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
791		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
792		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
793		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
794		real, intent(in) :: dt !< Time increment [T ~> s]
795		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
796		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
797
798		! Local variables
799	24	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
800	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of buoyancy frequency at u-points [L2 Z-2 T-2 ~> s-2]
801	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of buoyancy frequency at v-points [L2 Z-2 T-2 ~> s-2]
802	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
803	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
804	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! \|Sx\| N times dz at u-points [Z T-1 ~> m s-1]
805	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzSyN ! \|Sy\| N times dz at v-points [Z T-1 ~> m s-1]
806		integer :: niblock, njblock, nkblock
807
808	12	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
809	0	"Module must be initialized before it is used.")
810
811	12	niblock = CS%niblock
812	12	njblock = CS%njblock
813	12	nkblock = CS%nkblock
814
815	12	if (niblock == 0) niblock = G%iec - G%isc + 1
816	12	if (njblock == 0) njblock = G%jec - G%jsc + 1
817	12	if (nkblock == 0) nkblock = GV%ke
818
819	12	if (CS%calculate_Eady_growth_rate) then
820		!$omp target update to(h)
821		!$omp target enter data map(alloc: e)
822	12	call find_eta(h, tv, G, GV, US, e, halo_size=2)
823
824	12	if (CS%use_simpler_Eady_growth_rate) then
825		!$omp target enter data map(to: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
826		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
827		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
828		!$omp target enter data map(alloc: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
829		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
830		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
831		N2_v=N2_v, dzu=dzu, dzv=dzv, dzSxN=dzSxN, dzSyN=dzSyN, halo=1, &
832	0	OBC=OBC, OBC_N2=CS%OBC_friendly)
833		!$omp target update from(e, dzu, dzv, dzSxN, dzSyN)
834	0	call calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v)
835		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
836		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
837		!$omp target exit data map(release: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
838		!$omp target exit data map(delete: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
839	12	elseif (CS%use_stored_slopes) then
840		!$omp target enter data map(to: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
841		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
842		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
843		!$omp target enter data map(alloc: N2_u, N2_v)
844		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
845		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
846	12	N2_v=N2_v, halo=1, OBC=OBC, OBC_N2=CS%OBC_friendly)
847		!$omp target exit data map(from: CS%slope_x, CS%slope_y, N2_u, N2_v)
848	12	call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
849		!$omp target exit data map(release: tv%T, tv%S, tv, CS%slope_x, CS%slope_y)
850		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
851		!$omp target exit data map(delete: N2_u, N2_v )
852		else
853		!$omp target update from(e)
854	0	call calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
855		endif
856		!$omp target exit data map(delete: e)
857		endif
858
859	12	if (query_averaging_enabled(CS%diag)) then
860	12	if (CS%id_dzu > 0) call post_data(CS%id_dzu, dzu, CS%diag)
861	12	if (CS%id_dzv > 0) call post_data(CS%id_dzv, dzv, CS%diag)
862	12	if (CS%id_dzSxN > 0) call post_data(CS%id_dzSxN, dzSxN, CS%diag)
863	12	if (CS%id_dzSyN > 0) call post_data(CS%id_dzSyN, dzSyN, CS%diag)
864	12	if (CS%id_SN_u > 0) call post_data(CS%id_SN_u, CS%SN_u, CS%diag)
865	12	if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag)
866	12	if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag)
867	12	if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag)
868	12	if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
869	12	if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
870		endif
871
872	12	end subroutine calc_slope_functions
873
874		!> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al., 1997.
875		!! This is on older implementation that is susceptible to large values of Eady growth rate
876		!! for incropping layers.
877	12	subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
878		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
879		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
880		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
881		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [Z L-1 ~> nondim]
882		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency
883		!! at u-points [L2 Z-2 T-2 ~> s-2]
884		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope
885		!! [Z L-1 ~> nondim]
886		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency
887		!! at v-points [L2 Z-2 T-2 ~> s-2]
888		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
889		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
890		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
891
892		! Local variables
893		real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
894		real :: N2 ! Positive buoyancy frequency or zero [L2 Z-2 T-2 ~> s-2]
895		real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
896		real :: H_geom ! The geometric mean of Hup and Hdn [H ~> m or kg m-2].
897		real :: S2max ! An upper bound on the squared slopes [Z2 L-2 ~> nondim]
898		real :: wNE, wSE, wSW, wNW ! Weights of adjacent points [nondim]
899	24	real :: H_u(SZIB_(G)), H_v(SZI_(G)) ! Layer thicknesses at u- and v-points [H ~> m or kg m-2]
900
901		! Note that at some points in the code S2_u and S2_v hold the running depth
902		! integrals of the squared slope [H ~> m or kg m-2] before the average is taken.
903	24	real :: S2_u(SZIB_(G),SZJ_(G)) ! At first the thickness-weighted depth integral of the squared
904		! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
905		! squared slope [Z2 L-2 ~> nondim] at u points.
906	24	real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared
907		! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
908		! squared slope [Z2 L-2 ~> nondim] at v points.
909	24	integer :: OBC_dir_u(SZIB_(G),SZJ_(G)) ! An integer indicating where there are u OBCs: +1 for
910		! eastern OBCs, -1 for western OBCs and 0 at points with no OBCs.
911	24	integer :: OBC_dir_v(SZI_(G),SZJB_(G)) ! An integer indicating where there are v OBCs: +1 for
912		! northern OBCs, -1 for southern OBCs and 0 at points with no OBCs.
913	12	real :: h4_u(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a u-point
914		! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
915	12	real :: h4_v(SZI_(G),SZJB_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a v-point
916		! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
917		integer :: i, j, k, is, ie, js, je, nz
918
919	12	if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// &
920	0	"Module must be initialized before it is used.")
921
922	12	if (.not. CS%calculate_Eady_growth_rate) return
923	12	if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function: "// &
924	0	"%SN_u is not associated with use_variable_mixing.")
925	12	if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function: "// &
926	0	"%SN_v is not associated with use_variable_mixing.")
927
928	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
929
930	12	S2max = CS%Visbeck_S_max**2
931
932	106092	CS%SN_u(:,:) = 0.0
933	106824	CS%SN_v(:,:) = 0.0
934
935		! These settings apply where there are not open boundary conditions.
936	212904	OBC_dir_u(:,:) = 0 ; OBC_dir_v(:,:) = 0
937
938	12	if (associated(OBC) .and. CS%OBC_friendly) then
939		! Store the direction of any OBC faces.
940		!$OMP parallel do default(shared)
941	0	do j=js-1,je+1 ; do I=is-1,ie ; if (OBC%segnum_u(I,j) /= 0) then
942	0	if (OBC%segnum_u(I,j) > 0) OBC_dir_u(I,j) = 1 ! OBC_DIRECTION_E
943	0	if (OBC%segnum_u(I,j) < 0) OBC_dir_u(I,j) = -1 ! OBC_DIRECTION_W
944		endif ; enddo ; enddo
945		!$OMP parallel do default(shared)
946	0	do J=js-1,je ; do i=is-1,ie+1 ; if (OBC%segnum_v(i,J) /= 0) then
947	0	if (OBC%segnum_v(i,J) > 0) OBC_dir_v(i,J) = 1 ! OBC_DIRECTION_N
948	0	if (OBC%segnum_v(i,J) < 0) OBC_dir_v(i,J) = -1 ! OBC_DIRECTION_S
949		endif ; enddo ; enddo
950
951		! Use the masked product of the 4 (or 2) thicknesses around a velocity-point interface for weights.
952		!$OMP parallel do default(shared)
953	0	do K=2,nz
954	0	do j=js-1,je+1 ; do I=is-1,ie
955	0	if (OBC_dir_u(I,j) == 0) then
956	0	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)h(i+1,j,k)) (h(i,j,k-1)*h(i+1,j,k-1)) )
957	0	elseif (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E
958	0	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)*2) (h(i,j,k-1)**2) )
959	0	elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W
960	0	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i+1,j,k)*2) (h(i+1,j,k-1)**2) )
961		endif
962		enddo ; enddo
963	0	do J=js-1,je ; do i=is-1,ie+1
964	0	if (OBC_dir_v(i,J) == 0) then
965	0	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)h(i,j+1,k)) (h(i,j,k-1)*h(i,j+1,k-1)) )
966	0	elseif (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N
967	0	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)*2) (h(i,j,k-1)**2) )
968	0	elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S
969	0	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j+1,k)*2) (h(i,j+1,k-1)**2) )
970		endif
971		enddo ; enddo
972		enddo
973		else ! The land mask is sufficient and there are no special considerations taken at OBC points.
974		! Use the masked product of the 4 thicknesses around a velocity-point interface for weights.
975		!$OMP parallel do default(shared)
976	900	do K=2,nz
977	6717720	do j=js-1,je+1 ; do I=is-1,ie
978	6716832	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)h(i+1,j,k)) (h(i,j,k-1)*h(i+1,j,k-1)) )
979		enddo ; enddo
980	6663564	do J=js-1,je ; do i=is-1,ie+1
981	6662664	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)h(i,j+1,k)) (h(i,j,k-1)*h(i,j+1,k-1)) )
982		enddo ; enddo
983		enddo
984		endif
985
986		! To set the length scale based on the deformation radius, use wave_speed to
987		! calculate the first-mode gravity wave speed and then blend the equatorial
988		! and midlatitude deformation radii, using calc_resoln_function as a template.
989
990		!$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
991	732	do j=js,je
992	87840	do I=is-1,ie
993	87840	CS%SN_u(I,j) = 0. ; H_u(I) = 0. ; S2_u(I,j) = 0.
994		enddo
995	6500880	do K=2,nz ; do I=is-1,ie
996	6446880	Hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
997	6446880	Hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
998	6446880	H_geom = sqrt( Hdn * Hup )
999		!H_geom = H_geom * sqrt(N2) ! WKB-ish
1000		!H_geom = H_geom * N2 ! WKB-ish
1001	6446880	wSE = h4_v(i+1,J-1,K)
1002	6446880	wNW = h4_v(i,J,K)
1003	6446880	wNE = h4_v(i+1,J,K)
1004	6446880	wSW = h4_v(i,J-1,K)
1005	6446880	if (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E
1006	0	wSE = 0.0 ; wNE = 0.0
1007	0	H_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
1008	6446880	elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W
1009	0	wSW = 0.0 ; wNW = 0.0
1010	0	H_geom = sqrt( h(i+1,j,k) * h(i+1,j,k-1) )
1011		endif
1012		S2 = slope_x(I,j,K)**2 + &
1013		(((wNWslope_y(i,J,K)2) + (wSEslope_y(i+1,J-1,K)**2)) + &
1014		((wNEslope_y(i+1,J,K)2) + (wSWslope_y(i,J-1,K)**2)) ) / &
1015	6446880	( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 )
1016	6446880	if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2
1017
1018	6446880	N2 = max(0., N2_u(I,j,k))
1019	6446880	CS%SN_u(I,j) = CS%SN_u(I,j) + sqrt( S2N2 )H_geom
1020	6446880	S2_u(I,j) = S2_u(I,j) + S2*H_geom
1021	6500160	H_u(I) = H_u(I) + H_geom
1022		enddo ; enddo
1023	87852	do I=is-1,ie
1024	87840	if (H_u(I)>0.) then
1025	87120	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * CS%SN_u(I,j) / H_u(I)
1026	87120	S2_u(I,j) = G%OBCmaskCu(I,j) * S2_u(I,j) / H_u(I)
1027		else
1028	0	CS%SN_u(I,j) = 0.
1029		endif
1030		enddo
1031		enddo
1032
1033		!$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
1034	744	do J=js-1,je
1035	88572	do i=is,ie
1036	88572	CS%SN_v(i,J) = 0. ; H_v(i) = 0. ; S2_v(i,J) = 0.
1037		enddo
1038	6555060	do K=2,nz ; do i=is,ie
1039	6500160	Hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
1040	6500160	Hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
1041	6500160	H_geom = sqrt( Hdn * Hup )
1042		!H_geom = H_geom * sqrt(N2) ! WKB-ish
1043		!H_geom = H_geom * N2 ! WKB-ish
1044	6500160	wSE = h4_u(I,j,K)
1045	6500160	wNW = h4_u(I-1,j+1,K)
1046	6500160	wNE = h4_u(I,j+1,K)
1047	6500160	wSW = h4_u(I-1,j,K)
1048	6500160	if (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N
1049	0	wNW = 0.0 ; wNE = 0.0
1050	0	H_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
1051	6500160	elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S
1052	0	wSW = 0.0 ; wSE = 0.0
1053	0	H_geom = sqrt( h(i,j+1,k) * h(i,j+1,k-1) )
1054		endif
1055		S2 = slope_y(i,J,K)**2 + &
1056		(((wSEslope_x(I,j,K)2) + (wNWslope_x(I-1,j+1,K)**2)) + &
1057		((wNEslope_x(I,j+1,K)2) + (wSWslope_x(I-1,j,K)**2)) ) / &
1058	6500160	( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 )
1059	6500160	if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2
1060
1061	6500160	N2 = max(0., N2_v(i,J,K))
1062	6500160	CS%SN_v(i,J) = CS%SN_v(i,J) + sqrt( S2N2 )H_geom
1063	6500160	S2_v(i,J) = S2_v(i,J) + S2*H_geom
1064	6554328	H_v(i) = H_v(i) + H_geom
1065		enddo ; enddo
1066	88584	do i=is,ie
1067	88572	if (H_v(i)>0.) then
1068	87840	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * CS%SN_v(i,J) / H_v(i)
1069	87840	S2_v(i,J) = G%OBCmaskCv(i,J) * S2_v(i,J) / H_v(i)
1070		else
1071	0	CS%SN_v(i,J) = 0.
1072		endif
1073		enddo
1074		enddo
1075
1076		! Offer diagnostic fields for averaging.
1077	12	if (query_averaging_enabled(CS%diag)) then
1078	12	if (CS%id_S2_u > 0) call post_data(CS%id_S2_u, S2_u, CS%diag)
1079	12	if (CS%id_S2_v > 0) call post_data(CS%id_S2_v, S2_v, CS%diag)
1080		endif
1081
1082	12	if (CS%debug) then
1083		call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, &
1084	0	unscale=US%Z_to_L, haloshift=1)
1085		! call uvchksum("calc_Visbeck_coeffs_old S2_[uv]", S2_u, S2_v, G%HI, &
1086		! unscale=US%Z_to_L**2, scalar_pair=.true.)
1087		call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, &
1088	0	unscale=US%L_to_Z*2US%s_to_T**2, scalar_pair=.true.)
1089		call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
1090	0	unscale=US%s_to_T, scalar_pair=.true.)
1091		endif
1092
1093		end subroutine calc_Visbeck_coeffs_old
1094
1095		!> Calculates the Eady growth rate (2D fields) for use in MEKE and the Visbeck schemes
1096	0	subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, SN_u, SN_v)
1097		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1098		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1099		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1100		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1101		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m]
1102		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m]
1103		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m]
1104		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzv !< dz at v-points [Z ~> m]
1105		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzSxN !< dz Sx N at u-points [Z T-1 ~> m s-1]
1106		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzSyN !< dz Sy N at v-points [Z T-1 ~> m s-1]
1107		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: SN_u !< SN at u-points [T-1 ~> s-1]
1108		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: SN_v !< SN at v-points [T-1 ~> s-1]
1109		! Local variables
1110		real :: D_scale ! The depth over which to average SN [Z ~> m]
1111		real :: dnew ! Depth of bottom of layer [Z ~> m]
1112		real :: dz ! Limited thickness of this layer [Z ~> m]
1113		real :: weight ! Fraction of this layer that contributes to integral [nondim]
1114	0	real :: sum_dz(SZI_(G)) ! Cumulative sum of z-thicknesses [Z ~> m]
1115	0	real :: vint_SN(SZIB_(G)) ! Cumulative integral of SN [Z T-1 ~> m s-1]
1116	0	real, dimension(SZIB_(G),SZJ_(G)) :: SN_cpy !< SN at u-points [T-1 ~> s-1]
1117		real :: dz_neglect ! A negligibly small distance to avoid division by zero [Z ~> m]
1118		real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~> m-1]
1119		real :: dB, dT ! Elevation variables used when cropping [Z ~> m]
1120		integer :: i, j, k
1121		logical :: crop
1122
1123	0	dz_neglect = GV%dZ_subroundoff
1124	0	D_scale = CS%Eady_GR_D_scale
1125	0	if (D_scale<=0.) D_scale = 64.*GV%max_depth ! 0 means use full depth so choose something big
1126	0	r_crp_dist = 1. / max( dz_neglect, CS%cropping_distance )
1127	0	crop = CS%cropping_distance>=0. ! Only filter out in-/out-cropped interface is parameter if non-negative
1128
1129	0	if (CS%debug) then
1130	0	call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.)
1131		call uvchksum("calc_Eady_growth_rate_2D dzS2N2[uv]", dzSxN, dzSyN, G%HI, &
1132	0	unscale=US%Z_to_m*US%s_to_T, scalar_pair=.true.)
1133		endif
1134
1135		!$OMP parallel do default(shared)
1136	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1137	0	CS%SN_u(i,j) = 0.0
1138	0	CS%SN_v(i,j) = 0.0
1139		enddo ; enddo
1140
1141		!$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1142	0	do j=G%jsc-1,G%jec+1
1143	0	do I=G%isc-1,G%iec
1144	0	vint_SN(I) = 0.
1145	0	sum_dz(I) = dz_neglect
1146		enddo
1147	0	if (crop) then
1148	0	do K=2,GV%ke ; do I=G%isc-1,G%iec
1149	0	dnew = sum_dz(I) + dzu(I,j,K) ! This is where the bottom of the layer is
1150	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1151	0	dz = max(0., dnew - sum_dz(I)) ! This is the part of the layer to be included in the integral.
1152		! When D_scale>dnew, dz=dzu (+roundoff error).
1153		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1154		! When D_scale<sum_dz, dz=0.
1155	0	weight = dz / ( dzu(I,j,K) + dz_neglect ) ! Fraction of this layer to include
1156	0	dT = min( e(i,j,1), e(i+1,j,1) ) ! Deepest sea surface
1157	0	dB = max( e(i,j,K), e(i+1,j,K) ) ! Shallowest interface
1158	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1159	0	dT = min( e(i,j,K), e(i+1,j,K) ) ! Deepest interface
1160	0	dB = max( e(i,j,GV%ke+1), e(i+1,j,GV%ke+1) ) ! Shallowest topography
1161	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1162	0	vint_SN(I) = vint_SN(I) + weight * dzSxN(I,j,K)
1163	0	sum_dz(I) = sum_dz(I) + weight * dzu(I,j,K)
1164		enddo ; enddo
1165		else
1166	0	do K=2,GV%ke ; do I=G%isc-1,G%iec
1167	0	dnew = sum_dz(I) + dzu(I,j,K) ! This is where the bottom of the layer is
1168	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1169	0	dz = max(0., dnew - sum_dz(I)) ! This is the part of the layer to be included in the integral.
1170		! When D_scale>dnew, dz=dzu (+roundoff error).
1171		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1172		! When D_scale<sum_dz, dz=0.
1173	0	weight = dz / ( dzu(I,j,K) + dz_neglect ) ! Fraction of this layer to include
1174	0	vint_SN(I) = vint_SN(I) + weight * dzSxN(I,j,K)
1175	0	sum_dz(I) = sum_dz(I) + weight * dzu(I,j,K)
1176		enddo ; enddo
1177		endif
1178	0	do I=G%isc-1,G%iec
1179	0	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) )
1180	0	SN_cpy(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) )
1181		enddo
1182		enddo
1183
1184		!$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1185	0	do J=G%jsc-1,G%jec
1186	0	do i=G%isc-1,G%iec+1
1187	0	vint_SN(i) = 0.
1188	0	sum_dz(i) = dz_neglect
1189		enddo
1190	0	if (crop) then
1191	0	do K=2,GV%ke ; do i=G%isc-1,G%iec+1
1192	0	dnew = sum_dz(i) + dzv(i,J,K) ! This is where the bottom of the layer is
1193	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1194	0	dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1195		! When D_scale>dnew, dz=dzu (+roundoff error).
1196		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1197		! When D_scale<sum_dz, dz=0.
1198	0	weight = dz / ( dzv(i,J,K) + dz_neglect ) ! Fraction of this layer to include
1199	0	dT = min( e(i,j,1), e(i,j+1,1) ) ! Deepest sea surface
1200	0	dB = max( e(i,j,K), e(i,j+1,K) ) ! Shallowest interface
1201	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1202	0	dT = min( e(i,j,K), e(i,j+1,K) )! Deepest interface
1203	0	dB = max( e(i,j,GV%ke+1), e(i,j+1,GV%ke+1) ) ! Shallowest topography
1204	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1205	0	vint_SN(I) = vint_SN(I) + weight*2 dzSyN(i,J,K)
1206	0	sum_dz(i) = sum_dz(i) + weight * dzv(i,J,K)
1207		enddo ; enddo
1208		else
1209	0	do K=2,GV%ke ; do i=G%isc-1,G%iec+1
1210	0	dnew = sum_dz(i) + dzv(i,J,K) ! This is where the bottom of the layer is
1211	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1212	0	dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1213		! When D_scale>dnew, dz=dzu (+roundoff error).
1214		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1215		! When D_scale<sum_dz, dz=0.
1216	0	weight = dz / ( dzv(i,J,K) + dz_neglect ) ! Fraction of this layer to include
1217	0	vint_SN(I) = vint_SN(I) + weight*2 dzSyN(i,J,K)
1218	0	sum_dz(i) = sum_dz(i) + weight * dzv(i,J,K)
1219		enddo ; enddo
1220		endif
1221	0	do i=G%isc-1,G%iec+1
1222	0	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * ( vint_SN(i) / sum_dz(i) )
1223		enddo
1224		enddo
1225
1226	0	do j=G%jsc,G%jec
1227	0	do I=G%isc-1,G%iec
1228		CS%SN_u(I,j) = sqrt( SN_cpy(I,j)**2 &
1229		+ 0.25( ((CS%SN_v(i,J)2) + (CS%SN_v(i+1,J-1)*2)) &
1230	0	+ ((CS%SN_v(i+1,J)2) + (CS%SN_v(i,J-1)2)) ) )
1231		enddo
1232		enddo
1233	0	do J=G%jsc-1,G%jec
1234	0	do i=G%isc,G%iec
1235		CS%SN_v(i,J) = sqrt( CS%SN_v(i,J)**2 &
1236		+ 0.25( ((SN_cpy(I,j)2) + (SN_cpy(I-1,j+1)*2)) &
1237	0	+ ((SN_cpy(I,j+1)2) + (SN_cpy(I-1,j)2)) ) )
1238		enddo
1239		enddo
1240
1241	0	if (CS%debug) then
1242		call uvchksum("calc_Eady_growth_rate_2D SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
1243	0	unscale=US%s_to_T, scalar_pair=.true.)
1244		endif
1245
1246	0	end subroutine calc_Eady_growth_rate_2D
1247
1248		!> The original calc_slope_function() that calculated slopes using
1249		!! interface positions only, not accounting for density variations.
1250	0	subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
1251		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1252		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1253		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1254		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1255		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
1256		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m]
1257		! type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1258		! Local variables
1259	0	real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics)
1260	0	real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics)
1261	0	real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m]
1262		! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m]
1263		real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
1264		real :: dZ_cutoff ! A minimum water column depth for masking [H ~> m or kg m-2]
1265		real :: h1, h2 ! Temporary total thicknesses [Z ~> m]
1266		real :: h_neglect ! A thickness that is so small it is usually lost
1267		! in roundoff and can be neglected [H ~> m or kg m-2].
1268		real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
1269		real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2]
1270		real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
1271		real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
1272	0	real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times
1273		! the buoyancy frequency squared at u-points [Z T-2 ~> m s-2]
1274	0	real :: S2N2_v_local(SZI_(G),SZJB_(G),SZK_(GV)) ! The depth integral of the slope times
1275		! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2]
1276		logical :: use_dztot ! If true, use the total water column thickness rather than the
1277		! bathymetric depth for certain calculations.
1278		integer :: is, ie, js, je, nz
1279		integer :: i, j, k
1280
1281	0	if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_using_just_e: "// &
1282	0	"Module must be initialized before it is used.")
1283
1284	0	if (.not. CS%calculate_Eady_growth_rate) return
1285	0	if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function: "// &
1286	0	"%SN_u is not associated with use_variable_mixing.")
1287	0	if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function: "// &
1288	0	"%SN_v is not associated with use_variable_mixing.")
1289
1290	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1291
1292	0	h_neglect = GV%H_subroundoff
1293	0	H_cutoff = real(2nz) (GV%Angstrom_H + h_neglect)
1294	0	dZ_cutoff = real(2nz) (GV%Angstrom_Z + GV%dz_subroundoff)
1295
1296	0	use_dztot = CS%full_depth_Eady_growth_rate ! .or. .not.(GV%Boussinesq or GV%semi_Boussinesq)
1297
1298	0	if (use_dztot) then
1299		!$OMP parallel do default(shared)
1300	0	do j=js-1,je+1 ; do i=is-1,ie+1
1301	0	dz_tot(i,j) = e(i,j,1) - e(i,j,nz+1)
1302		enddo ; enddo
1303		! The following mathematically equivalent expression is more expensive but is less
1304		! sensitive to roundoff for large Z_ref:
1305		! call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
1306		! do j=js-1,je+1
1307		! do i=is-1,ie+1 ; dz_tot(i,j) = 0.0 ; enddo
1308		! do k=1,nz ; do i=is-1,ie+1
1309		! dz_tot(i,j) = dz_tot(i,j) + dz(i,j,k)
1310		! enddo ; enddo
1311		! enddo
1312		endif
1313
1314		! To set the length scale based on the deformation radius, use wave_speed to
1315		! calculate the first-mode gravity wave speed and then blend the equatorial
1316		! and midlatitude deformation radii, using calc_resoln_function as a template.
1317
1318		!$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
1319	0	do k=nz,CS%VarMix_Ktop,-1
1320
1321		! Calculate the interface slopes E_x and E_y and u- and v- points respectively
1322	0	do j=js-1,je+1 ; do I=is-1,ie
1323	0	E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
1324		! Mask slopes where interface intersects topography
1325	0	if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0.
1326		enddo ; enddo
1327	0	do J=js-1,je ; do i=is-1,ie+1
1328	0	E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
1329		! Mask slopes where interface intersects topography
1330	0	if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
1331		enddo ; enddo
1332
1333		! Calculate NSh from this layer and add to the sum
1334	0	do j=js,je ; do I=is-1,ie
1335		S2 = ( E_x(I,j)*2 + 0.25( &
1336	0	((E_y(i,J)2) + (E_y(i+1,J-1)2)) + ((E_y(i+1,J)2) + (E_y(i,J-1)2)) ) )
1337	0	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
1338
1339	0	Hdn = 2.h(i,j,k)h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1340	0	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)
1341	0	H_geom = sqrt(Hdn*Hup)
1342		! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1343	0	S2N2_u_local(I,j,k) = (H_geom * S2) * (GV%g_prime(k) / max(Hdn, Hup, CS%h_min_N2) )
1344		enddo ; enddo
1345	0	do J=js-1,je ; do i=is,ie
1346		S2 = ( E_y(i,J)*2 + 0.25( &
1347	0	((E_x(I,j)2) + (E_x(I-1,j+1)2)) + ((E_x(I,j+1)2) + (E_x(I-1,j)2)) ) )
1348	0	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
1349
1350	0	Hdn = 2.h(i,j,k)h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1351	0	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)
1352	0	H_geom = sqrt(Hdn*Hup)
1353		! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1354	0	S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2)))
1355		enddo ; enddo
1356
1357		enddo ! k
1358
1359		!$OMP parallel do default(shared)
1360	0	do j=js,je
1361	0	do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo
1362	0	do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie
1363	0	CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k)
1364		enddo ; enddo
1365		! SN above contains S^2N^2H, convert to vertical average of S*N
1366
1367	0	if (use_dztot) then
1368	0	do I=is-1,ie
1369		CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / &
1370	0	max(dz_tot(i,j), dz_tot(i+1,j), GV%dz_subroundoff) )
1371		enddo
1372		else
1373	0	do I=is-1,ie
1374	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
1375	0	h2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)
1376	0	if ( min(h1, h2) > dZ_cutoff ) then
1377	0	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / max(h1, h2) )
1378		else
1379	0	CS%SN_u(I,j) = 0.0
1380		endif
1381		enddo
1382		endif
1383		enddo
1384		!$OMP parallel do default(shared)
1385	0	do J=js-1,je
1386	0	do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo
1387	0	do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie
1388	0	CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k)
1389		enddo ; enddo
1390	0	if (use_dztot) then
1391	0	do i=is,ie
1392		CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / &
1393	0	max(dz_tot(i,j), dz_tot(i,j+1), GV%dz_subroundoff) )
1394		enddo
1395		else
1396	0	do i=is,ie
1397		! There is a primordial horizontal indexing bug on the following line from the previous
1398		! versions of the code. This comment should be deleted by the end of 2024.
1399		! if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
1400	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
1401	0	h2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)
1402	0	if ( min(h1, h2) > dZ_cutoff ) then
1403	0	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / max(h1, h2) )
1404		else
1405	0	CS%SN_v(i,J) = 0.0
1406		endif
1407		enddo
1408		endif
1409		enddo
1410
1411		end subroutine calc_slope_functions_using_just_e
1412
1413
1414		!> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity.
1415	0	subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC)
1416		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1417		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1418		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1419		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1420		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1421		real, intent(in) :: dt !< Time increment [T ~> s]
1422		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1423		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1424		type(VarMix_CS), intent(in) :: CS !< Variable mixing control structure
1425		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
1426		! Local variables
1427	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
1428		integer :: niblock, njblock, nkblock
1429
1430	0	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_QG_slopes: "//&
1431	0	"Module must be initialized before it is used.")
1432
1433	0	niblock = CS%niblock
1434	0	njblock = CS%njblock
1435	0	nkblock = CS%nkblock
1436
1437	0	if (niblock == 0) niblock = G%iec - G%isc + 1
1438	0	if (njblock == 0) njblock = G%jec - G%jsc + 1
1439	0	if (nkblock == 0) nkblock = GV%ke
1440
1441		!$omp target update to(h)
1442		!$omp target enter data map(alloc: e)
1443	0	call find_eta(h, tv, G, GV, US, e, halo_size=3)
1444		!$omp target enter data map(to: tv%T, tv%S)
1445		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
1446		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
1447		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
1448		slope_x, slope_y, niblock, njblock, nkblock, halo=2, OBC=OBC, &
1449	0	OBC_N2=CS%OBC_friendly)
1450		!$omp target exit data map(release: tv%T, tv%S)
1451		!$omp target exit data map(delete: e)
1452		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
1453		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
1454
1455	0	end subroutine calc_QG_slopes
1456
1457		!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
1458	0	subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, slope_x, slope_y, &
1459	0	vort_xy_dx, vort_xy_dy)
1460		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1461		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1462		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1463		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1464		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1465		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer vertical extents [Z ~> m]
1466		integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
1467		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence
1468		!! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1469		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence
1470		!! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1471		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1472		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1473		real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
1474		!! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1475		real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
1476		!! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1477		! Local variables
1478		real, dimension(SZI_(G),SZJB_(G)) :: &
1479	0	dslopey_dz, & ! z-derivative of y-slope at v-points [L-1 ~> m-1]
1480	0	h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
1481	0	beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
1482	0	grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
1483	0	grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]
1484
1485		real, dimension(SZIB_(G),SZJ_(G)) :: &
1486	0	dslopex_dz, & ! z-derivative of x-slope at u-points [L-1 ~> m-1]
1487	0	h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
1488	0	beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
1489	0	grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
1490	0	grad_div_mag_u ! Magnitude of divergence gradient at u-points [T-1 L-1 ~> s-1 m-1]
1491		real :: h_at_slope_above ! The thickness above [H ~> m or kg m-2]
1492		real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2]
1493		real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
1494		real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1]
1495		real :: Z_to_H ! A local copy of depth to thickness conversion factors or the inverse of the
1496		! mass-weighted average specific volumes around an interface [H Z-1 ~> nondim or kg m-3]
1497		real :: inv_PI3 ! The inverse of pi cubed [nondim]
1498		integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1499
1500	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1501	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1502	0	nz = GV%ke
1503
1504	0	inv_PI3 = 1.0 / ((4.0atan(1.0))*3)
1505	0	Z_to_H = GV%Z_to_H ! This will be replaced with a varying value in non-Boussinesq mode.
1506
1507	0	if ((k > 1) .and. (k < nz)) then
1508
1509	0	do j=js-2,je+2 ; do I=is-2,ie+1
1510		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) ) / &
1511		( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
1512	0	+ ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**3 )
1513		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) ) / &
1514		( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
1515	0	+ ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**3 )
1516	0	Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
1517	0	if (.not.GV%Boussinesq) &
1518		Z_to_H = ( (h(i,j,k-1) + h(i+1,j,k-1)) + (h(i,j,k) + h(i+1,j,k)) ) / &
1519	0	( (dz(i,j,k-1) + dz(i+1,j,k-1)) + (dz(i,j,k) + dz(i+1,j,k)) + GV%dZ_subroundoff)
1520	0	dslopex_dz(I,j) = 2. * ( slope_x(I,j,k) - slope_x(I,j,k+1) ) * (Z_to_H * Ih)
1521	0	h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
1522		enddo ; enddo
1523
1524	0	do J=js-2,je+1 ; do i=is-2,ie+2
1525		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) ) / &
1526		( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
1527	0	+ ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**3 )
1528		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) ) / &
1529		( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
1530	0	+ ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**3 )
1531	0	Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
1532	0	if (.not.GV%Boussinesq) &
1533		Z_to_H = ( (h(i,j,k-1) + h(i,j+1,k-1)) + (h(i,j,k) + h(i,j+1,k)) ) / &
1534	0	( (dz(i,j,k-1) + dz(i,j+1,k-1)) + (dz(i,j,k) + dz(i,j+1,k)) + GV%dZ_subroundoff)
1535	0	dslopey_dz(i,J) = 2. * ( slope_y(i,J,k) - slope_y(i,J,k+1) ) * (Z_to_H * Ih)
1536	0	h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
1537		enddo ; enddo
1538
1539	0	do J=js-2,je+1 ; do i=is-1,ie+1
1540	0	f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) )
1541		vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * &
1542		( ( (h_at_u(I,j) * dslopex_dz(I,j)) + (h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1)) ) &
1543		+ ( (h_at_u(I-1,j) * dslopex_dz(I-1,j)) + (h_at_u(I,j+1) * dslopex_dz(I,j+1)) ) ) / &
1544	0	( ( 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)
1545		enddo ; enddo
1546
1547	0	do j=js-1,je+1 ; do I=is-2,ie+1
1548	0	f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
1549		vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * &
1550		( ( (h_at_v(i,J) * dslopey_dz(i,J)) + (h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1)) ) &
1551		+ ( (h_at_v(i,J-1) * dslopey_dz(i,J-1)) + (h_at_v(i+1,J) * dslopey_dz(i+1,J)) ) ) / &
1552	0	( ( 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)
1553		enddo ; enddo
1554		endif ! k > 1
1555
1556	0	if (CS%use_QG_Leith_GM) then
1557
1558	0	do j=js,je ; do I=is-1,Ieq
1559		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)) &
1560	0	+ (vort_xy_dx(i+1,J) + vort_xy_dx(i,J-1))))**2)
1561		grad_div_mag_u(I,j) = SQRT(div_xx_dx(I,j)*2 + (0.25((div_xx_dy(i,J) + div_xx_dy(i+1,J-1)) &
1562	0	+ (div_xx_dy(i+1,J) + div_xx_dy(i,J-1))))**2)
1563	0	if (CS%use_beta_in_QG_Leith) then
1564		beta_u(I,j) = sqrt((0.5(G%dF_dx(i,j)+G%dF_dx(i+1,j))*2) + &
1565	0	(0.5(G%dF_dy(i,j)+G%dF_dy(i+1,j))*2))
1566		CS%KH_u_QG(I,j,k) = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), 3.0beta_u(I,j)) &
1567	0	CS%Laplac3_const_u(I,j) * inv_PI3
1568		else
1569		CS%KH_u_QG(I,j,k) = (grad_vort_mag_u(I,j) + grad_div_mag_u(I,j)) * &
1570	0	CS%Laplac3_const_u(I,j) * inv_PI3
1571		endif
1572		enddo ; enddo
1573
1574	0	do J=js-1,Jeq ; do i=is,ie
1575		grad_vort_mag_v(i,J) = SQRT(vort_xy_dx(i,J)*2 + (0.25((vort_xy_dy(I,j) + vort_xy_dy(I-1,j+1)) &
1576	0	+ (vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j))))**2)
1577		grad_div_mag_v(i,J) = SQRT(div_xx_dy(i,J)*2 + (0.25((div_xx_dx(I,j) + div_xx_dx(I-1,j+1)) &
1578	0	+ (div_xx_dx(I,j+1) + div_xx_dx(I-1,j))))**2)
1579	0	if (CS%use_beta_in_QG_Leith) then
1580		beta_v(i,J) = sqrt((0.5(G%dF_dx(i,j)+G%dF_dx(i,j+1))*2) + &
1581	0	(0.5(G%dF_dy(i,j)+G%dF_dy(i,j+1))*2))
1582		CS%KH_v_QG(i,J,k) = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), 3.0beta_v(i,J)) &
1583	0	CS%Laplac3_const_v(i,J) * inv_PI3
1584		else
1585		CS%KH_v_QG(i,J,k) = (grad_vort_mag_v(i,J) + grad_div_mag_v(i,J)) * &
1586	0	CS%Laplac3_const_v(i,J) * inv_PI3
1587		endif
1588		enddo ; enddo
1589		! post diagnostics
1590
1591	0	if (k==nz) then
1592	0	if (CS%id_KH_v_QG > 0) call post_data(CS%id_KH_v_QG, CS%KH_v_QG, CS%diag)
1593	0	if (CS%id_KH_u_QG > 0) call post_data(CS%id_KH_u_QG, CS%KH_u_QG, CS%diag)
1594		endif
1595		endif
1596
1597	0	end subroutine calc_QG_Leith_viscosity
1598
1599		!> Initializes the variables mixing coefficients container
1600	1	subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
1601		type(time_type), intent(in) :: Time !< Current model time
1602		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1603		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1604		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1605		type(param_file_type), intent(in) :: param_file !< Parameter file handles
1606		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
1607		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1608
1609		! Local variables
1610		real :: KhTr_Slope_Cff ! The nondimensional coefficient in the Visbeck formula
1611		! for the epipycnal tracer diffusivity [nondim]
1612		real :: KhTh_Slope_Cff ! The nondimensional coefficient in the Visbeck formula
1613		! for the interface depth diffusivity [nondim]
1614		real :: oneOrTwo ! A variable that may be 1 or 2, depending on which form
1615		! of the equatorial deformation radius us used [nondim]
1616		real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when
1617		! calculating the first-mode wave speed [H ~> m or kg m-2]
1618		real :: KhTr_passivity_coeff ! Coefficient setting the ratio between along-isopycnal tracer
1619		! mixing and interface height mixing [nondim]
1620		real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
1621		! default value is roughly (pi / (the age of the universe)).
1622		logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
1623		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1624		integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1625		! for remapping. Values below 20190101 recover the remapping
1626		! answers from 2018, while higher values use more robust
1627		! forms of the same remapping expressions.
1628		real :: MLE_front_length ! The frontal-length scale used to calculate the upscaling of
1629		! buoyancy gradients in boundary layer parameterizations [L ~> m]
1630		real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity [nondim]
1631		real :: grid_sp_u2, grid_sp_v2 ! Intermediate quantities for Leith metrics [L2 ~> m2]
1632		real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
1633		real :: wave_speed_min ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
1634		real :: wave_speed_tol ! The fractional tolerance for finding the wave speeds [nondim]
1635		logical :: Resoln_scaled_MEKE_visc ! If true, the viscosity contribution from MEKE is
1636		! scaled by the resolution function.
1637		logical :: better_speed_est ! If true, use a more robust estimate of the first
1638		! mode wave speed as the starting point for iterations.
1639		real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
1640		! temperature variance [nondim]
1641		logical :: use_SQG ! This is true if the SQG structure will be used for any parameterizations.
1642		logical :: om4_remap_via_sub_cells ! Use the OM4-era remap_via_sub_cells for calculating the EBT structure
1643		logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1644		! recreate the bugs, or if false bugs are only used if actively selected.
1645		logical :: mixing_coefs_OBC_bug ! If false, use only interior data for thickness weighting in
1646		! lateral mixing coefficient calculations and to calculate stratification
1647		! and other fields at open boundary condition faces.
1648		! This include declares and sets the variable "version".
1649		# include "version_variable.h"
1650		character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
1651		integer :: number_of_OBC_segments
1652		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j
1653		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
1654		#ifdef __NVCOMPILER_OPENMP_GPU
1655		integer, parameter :: default_niblock = 0
1656		integer, parameter :: default_njblock = 0
1657		integer, parameter :: default_nkblock = 0
1658		#else
1659		integer, parameter :: default_niblock = 0
1660		integer, parameter :: default_njblock = 1
1661		integer, parameter :: default_nkblock = 1
1662		#endif
1663
1664	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1665	1	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1666	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1667	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1668
1669	1	CS%initialized = .true.
1670	1	in_use = .false. ! Set to true to avoid deallocating
1671	1	CS%diag => diag ! Diagnostics pointer
1672	1	CS%calculate_cg1 = .false.
1673	1	CS%calculate_Rd_dx = .false.
1674	1	CS%calculate_res_fns = .false.
1675	1	CS%use_simpler_Eady_growth_rate = .false.
1676	1	CS%full_depth_Eady_growth_rate = .false.
1677	1	CS%calculate_depth_fns = .false.
1678		! Read all relevant parameters and write them to the model log.
1679	1	call log_version(param_file, mdl, version, "")
1680		call get_param(param_file, mdl, "USE_VARIABLE_MIXING", CS%use_variable_mixing,&
1681		"If true, the variable mixing code will be called. This "//&
1682		"allows diagnostics to be created even if the scheme is "//&
1683		"not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
1684		"this is set to true regardless of what is in the "//&
1685	1	"parameter file.", default=.false.)
1686		call get_param(param_file, mdl, "USE_VISBECK", CS%use_Visbeck,&
1687		"If true, use the Visbeck et al. (1997) formulation for \n"//&
1688	1	"thickness diffusivity.", default=.false.)
1689		call get_param(param_file, mdl, "RESOLN_SCALED_KH", CS%Resoln_scaled_Kh, &
1690		"If true, the Laplacian lateral viscosity is scaled away "//&
1691		"when the first baroclinic deformation radius is well "//&
1692	1	"resolved.", default=.false.)
1693		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", CS%Depth_scaled_KhTh, &
1694		"If true, KHTH is scaled away when the depth is shallower "//&
1695		"than a reference depth: KHTH = MIN(1,H/H0)*N KHTH, "//&
1696		"where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
1697		"and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
1698	1	default=.false.)
1699		call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", CS%Resoln_scaled_KhTh, &
1700		"If true, the interface depth diffusivity is scaled away "//&
1701		"when the first baroclinic deformation radius is well "//&
1702	1	"resolved.", default=.false.)
1703		call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", CS%Resoln_scaled_KhTr, &
1704		"If true, the epipycnal tracer diffusivity is scaled "//&
1705		"away when the first baroclinic deformation radius is "//&
1706	1	"well resolved.", default=.false.)
1707		call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
1708	1	default=.false., do_not_log=.true.)
1709		call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", Resoln_scaled_MEKE_visc, &
1710		"If true, the viscosity contribution from MEKE is scaled by "//&
1711	1	"the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
1712	1	if (.not.use_MEKE) Resoln_scaled_MEKE_visc = .false.
1713		call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, &
1714		"If true, uses the equivalent barotropic wave speed instead "//&
1715		"of first baroclinic wave for calculating the resolution function.",&
1716	1	default=.false.)
1717		call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, &
1718		"Power to raise EBT vertical structure to when backscatter "// &
1719	1	"has vertical structure.", units="nondim", default=0.0)
1720		call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", CS%BS_use_sqg_struct, &
1721		"If true, the SQG vertical structure is used for backscatter "//&
1722		"on the condition that BS_EBT_power=0", &
1723	1	default=.false.)
1724		call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, &
1725		"If true, uses the equivalent barotropic structure "//&
1726		"as the vertical structure of thickness diffusivity.",&
1727	1	default=.false.)
1728		call get_param(param_file, mdl, "KHTH_USE_SQG_STRUCT", CS%khth_use_sqg_struct, &
1729		"If true, uses the surface quasigeostrophic structure "//&
1730		"as the vertical structure of thickness diffusivity.",&
1731	1	default=.false.)
1732		call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%khtr_use_ebt_struct, &
1733		"If true, uses the equivalent barotropic structure "//&
1734		"as the vertical structure of tracer diffusivity.",&
1735	1	default=.false.)
1736		call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", CS%khtr_use_sqg_struct, &
1737		"If true, uses the surface quasigeostrophic structure "//&
1738		"as the vertical structure of tracer diffusivity.",&
1739	1	default=.false.)
1740		call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", CS%kdgl90_use_ebt_struct, &
1741		"If true, uses the equivalent barotropic structure "//&
1742		"as the vertical structure of diffusivity in the GL90 scheme.",&
1743	1	default=.false.)
1744		call get_param(param_file, mdl, "KD_GL90_USE_SQG_STRUCT", CS%kdgl90_use_sqg_struct, &
1745		"If true, uses the equivalent barotropic structure "//&
1746		"as the vertical structure of diffusivity in the GL90 scheme.",&
1747	1	default=.false.)
1748		call get_param(param_file, mdl, "KHTH_SLOPE_CFF", KhTh_Slope_Cff, &
1749		"The nondimensional coefficient in the Visbeck formula "//&
1750	1	"for the interface depth diffusivity", units="nondim", default=0.0)
1751		call get_param(param_file, mdl, "KHTR_SLOPE_CFF", KhTr_Slope_Cff, &
1752		"The nondimensional coefficient in the Visbeck formula "//&
1753	1	"for the epipycnal tracer diffusivity", units="nondim", default=0.0)
1754		call get_param(param_file, mdl, "USE_STORED_SLOPES", CS%use_stored_slopes,&
1755		"If true, the isopycnal slopes are calculated once and "//&
1756		"stored for re-use. This uses more memory but avoids calling "//&
1757		"the equation of state more times than should be necessary.", &
1758	1	default=.false.)
1759		call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1760		"A miniscule frequency that is used to avoid division by 0. The default "//&
1761		"value is roughly (pi / (the age of the universe)).", &
1762	1	default=1.0e-17, units="s-1", scale=US%T_to_s)
1763		call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, &
1764	1	default=.false., do_not_log=.true.)
1765		CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct &
1766	1	.or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.
1767	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE
1768		! Indicate whether to calculate the Eady growth rate
1769	1	CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.)
1770		call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, &
1771	1	units="nondim", default=0., do_not_log=.true.)
1772	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (KhTr_passivity_coeff>0.)
1773		call get_param(param_file, mdl, "MLE_FRONT_LENGTH", MLE_front_length, &
1774	1	units="m", default=0.0, scale=US%m_to_L, do_not_log=.true.)
1775	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (MLE_front_length>0.)
1776
1777	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
1778
1779		! Isopycnal blocking parameters
1780		call get_param(param_file, mdl, "ISOPYCNAL_NIBLOCK", CS%niblock, &
1781		"The i-direction block size used to calculate isopycnal slopes. "//&
1782		"If 0, or when running with OpenMP offload, "//&
1783		"the full computational domain width is used. "//&
1784		"If USE_STANLEY_ISO is true, ISOPYCNAL_NIBLOCK cannot equal 1.", &
1785	1	default=default_niblock, layoutParam=.true.)
1786		call get_param(param_file, mdl, "ISOPYCNAL_NJBLOCK", CS%njblock, &
1787		"The j-direction block size used to calculate isopycnal slopes. "//&
1788		"If 0, defaults to 1, except when running with OpenMP offload, "//&
1789		"in which case the full computational domain height is used. " //&
1790		"If USE_STANLEY_ISO is true, ISOPYCNAL_NJBLOCK cannot equal 1.", &
1791	1	default=default_njblock, layoutParam=.true.)
1792		call get_param(param_file, mdl, "ISOPYCNAL_NKBLOCK", CS%nkblock, &
1793		"The k-direction block size used to calculate isopycnal slopes. "//&
1794		"If 0, defaults to 1, except when "//&
1795		"running with OpenMP offload, in which case the full computational "//&
1796	1	"domain depth is used.", default=default_nkblock, layoutParam=.true.)
1797	1	if (CS%niblock < 0) &
1798		call MOM_error(FATAL, "ISOPYCNAL_NIBLOCK must be nonnegative; "//&
1799	0	"use 0 to select the default block size.")
1800	1	if (CS%njblock < 0) &
1801		call MOM_error(FATAL, "ISOPYCNAL_NJBLOCK must be nonnegative; "//&
1802	0	"use 0 to select the default block size.")
1803	1	if (CS%nkblock < 0) &
1804		call MOM_error(FATAL, "ISOPYCNAL_NKBLOCK must be nonnegative; "//&
1805	0	"use 0 to select the default block size.")
1806
1807		call get_param(param_file, mdl, "USE_STANLEY_ISO", CS%use_stanley_iso, &
1808		"If true, turn on Stanley SGS T variance parameterization "// &
1809	1	"in isopycnal slope code.", default=.false.)
1810	1	if (CS%use_stanley_iso) then
1811		call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
1812		"Coefficient correlating the temperature gradient and SGS T variance.", &
1813	0	units="nondim", default=-1.0, do_not_log=.true.)
1814	0	if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
1815	0	"STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.")
1816	0	if (CS%njblock == 1) then
1817		call MOM_error(WARNING, "ISOPYCNAL_NJBLOCK must be >= 2 or 0 if USE_STANLEY_ISO is true."//&
1818	0	" Changing block size from 1 to 2 for this run.")
1819	0	CS%njblock = 2
1820		endif
1821	0	if (CS%niblock == 1) then
1822		call MOM_error(WARNING, "ISOPYCNAL_NIBLOCK must be >= 2 or 0 if USE_STANLEY_ISO is true."//&
1823	0	" Changing block size from 1 to 2 for this run.")
1824	0	CS%niblock = 2
1825		endif
1826		endif
1827		call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_OBC_segments, &
1828	1	default=0, do_not_log=.true.)
1829		call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1830	1	default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1831		call get_param(param_file, mdl, "MIXING_COEFS_OBC_BUG", mixing_coefs_OBC_bug, &
1832		"If false, use only interior data for thickness weighting in lateral mixing "//&
1833		"coefficient calculations and to calculate stratification and other fields at "//&
1834		"open boundary condition faces.", &
1835	1	default=enable_bugs, do_not_log=(number_of_OBC_segments<=0))
1836	1	CS%OBC_friendly = .not. MIXING_COEFS_OBC_BUG
1837		call get_param(param_file, mdl, "RESOLN_FUNCTION_OBC_BUG", CS%res_fn_OBC_bug, &
1838		"If false, use only interior data for calculating the resolution functions at "//&
1839		"open boundary condition faces and vertices.", &
1840	1	default=enable_bugs, do_not_log=(number_of_OBC_segments<=0))
1841
1842		if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct &
1843	1	.or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) then
1844	0	in_use = .true.
1845		call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, &
1846		"The depth below which N2 is monotonized to avoid stratification "//&
1847		"artifacts from altering the equivalent barotropic mode structure. "//&
1848		"This monotonzization is disabled if this parameter is negative.", &
1849	0	units="m", default=-1.0, scale=GV%m_to_H)
1850	0	allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1851		endif
1852
1853		use_SQG = CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct .or. &
1854	1	CS%kdgl90_use_sqg_struct
1855		call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, &
1856		"Nondimensional exponent coeffecient of the SQG mode that is used for the "//&
1857		"vertical struture of diffusivities.", &
1858	1	units="nondim", default=1.0, do_not_log=.not.use_SQG)
1859		call get_param(param_file, mdl, "INTERPOLATED_SQG_STRUCTURE", CS%interpolated_sqg_struct, &
1860		"If true, interpolate properties to velocity points and then interpolate the "//&
1861		"buoyancy frequencies and layer thicknesses back to tracer points when "//&
1862		"calculating the SQG vertical structure.", &
1863	1	default=.true., do_not_log=.not.use_SQG)
1864		!### Consider changing the default for INTERPOLATED_SQG_STRUCTURE to false.
1865
1866	1	if ((CS%BS_EBT_power>0.) .and. CS%BS_use_sqg_struct) call MOM_error(FATAL, &
1867	0	"calc_resoln_function: BS_EBT_POWER>0. and BS_USE_SQG=True cannot be set together")
1868
1869	1	if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) call MOM_error(FATAL, &
1870	0	"calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT and KHTH_USE_SQG_STRUCT can be true")
1871
1872	1	if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) call MOM_error(FATAL, &
1873	0	"calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT and KHTR_USE_SQG_STRUCT can be true")
1874
1875	1	if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) call MOM_error(FATAL, &
1876	0	"calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT and KD_GL90_USE_SQG_STRUCT can be true")
1877
1878	1	if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
1879	0	allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1880		endif
1881
1882	1	if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
1883	0	allocate(CS%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1884		endif
1885
1886	1	if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then
1887	0	allocate(CS%khtr_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1888		endif
1889
1890	1	if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then
1891	0	allocate(CS%kdgl90_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1892		endif
1893
1894	1	if (CS%use_stored_slopes) then
1895	1	if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then
1896		call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", CS%Visbeck_S_max, &
1897		"If non-zero, is an upper bound on slopes used in the "//&
1898		"Visbeck formula for diffusivity. This does not affect the "//&
1899		"isopycnal slope calculation used within thickness diffusion.", &
1900	0	units="nondim", default=0.0, scale=US%L_to_Z)
1901		else
1902	1	CS%Visbeck_S_max = 0.
1903		endif
1904		endif
1905
1906	1	if (CS%use_stored_slopes .or. (CS%interpolated_sqg_struct .and. (CS%sqg_expo>0.0))) then
1907		! CS%calculate_Eady_growth_rate=.true.
1908	1	in_use = .true.
1909	671917	allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0)
1910	676553	allocate(CS%slope_y(isd:ied,JsdB:JedB,GV%ke+1), source=0.0)
1911		call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, &
1912		"A diapycnal diffusivity that is used to interpolate "//&
1913		"more sensible values of T & S into thin layers.", &
1914	1	units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T)
1915		endif
1916
1917	1	if (CS%calculate_Eady_growth_rate) then
1918	1	in_use = .true.
1919	8841	allocate(CS%SN_u(IsdB:IedB,jsd:jed), source=0.0)
1920	8902	allocate(CS%SN_v(isd:ied,JsdB:JedB), source=0.0)
1921		CS%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, Time, &
1922	1	'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=US%s_to_T)
1923		CS%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, Time, &
1924	1	'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=US%s_to_T)
1925		call get_param(param_file, mdl, "USE_SIMPLER_EADY_GROWTH_RATE", CS%use_simpler_Eady_growth_rate, &
1926		"If true, use a simpler method to calculate the Eady growth rate "//&
1927	1	"that avoids division by layer thickness. Recommended.", default=.false.)
1928	1	if (CS%use_simpler_Eady_growth_rate) then
1929	0	if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
1930		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1931	0	"When USE_SIMPLER_EADY_GROWTH_RATE=True, USE_STORED_SLOPES must also be True.")
1932		call get_param(param_file, mdl, "EADY_GROWTH_RATE_D_SCALE", CS%Eady_GR_D_scale, &
1933		"The depth from surface over which to average SN when calculating "//&
1934		"a 2D Eady growth rate. Zero mean use full depth.", &
1935	0	units="m", default=0., scale=US%m_to_Z)
1936		call get_param(param_file, mdl, "EADY_GROWTH_RATE_CROPPING_DISTANCE", CS%cropping_distance, &
1937		"Distance from surface or bottom to filter out outcropped or "//&
1938		"incropped interfaces for the Eady growth rate calc. "//&
1939	0	"Negative values disables cropping.", units="m", default=0., scale=US%m_to_Z)
1940		else
1941		call get_param(param_file, mdl, "VARMIX_KTOP", CS%VarMix_Ktop, &
1942		"The layer number at which to start vertical integration "//&
1943		"of S*N for purposes of finding the Eady growth rate.", &
1944	1	units="nondim", default=2)
1945		call get_param(param_file, mdl, "MIN_DZ_FOR_SLOPE_N2", CS%h_min_N2, &
1946		"The minimum vertical distance to use in the denominator of the "//&
1947		"bouyancy frequency used in the slope calculation.", &
1948	1	units="m", default=1.0, scale=GV%m_to_H, do_not_log=CS%use_stored_slopes)
1949
1950		call get_param(param_file, mdl, "FULL_DEPTH_EADY_GROWTH_RATE", CS%full_depth_Eady_growth_rate, &
1951		"If true, calculate the Eady growth rate based on average slope times "//&
1952		"stratification that includes contributions from sea-level changes "//&
1953		"in its denominator, rather than just the nominal depth of the bathymetry. "//&
1954		"This only applies when using the model interface heights as a proxy for "//&
1955		"isopycnal slopes.", default=.not.(GV%Boussinesq.or.GV%semi_Boussinesq), &
1956	1	do_not_log=CS%use_stored_slopes)
1957		endif
1958		endif
1959
1960	1	if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then
1961	0	in_use = .true.
1962		call get_param(param_file, mdl, "VISBECK_L_SCALE", CS%Visbeck_L_scale, &
1963		"The fixed length scale in the Visbeck formula, or if negative a nondimensional "//&
1964		"scaling factor relating this length scale squared to the cell areas.", &
1965	0	units="m or nondim", default=0.0, scale=US%m_to_L)
1966	0	allocate(CS%L2u(IsdB:IedB,jsd:jed), source=0.0)
1967	0	allocate(CS%L2v(isd:ied,JsdB:JedB), source=0.0)
1968	0	if (CS%Visbeck_L_scale<0) then
1969		! Undo the rescaling of CS%Visbeck_L_scale.
1970	0	do j=js,je ; do I=is-1,Ieq
1971	0	CS%L2u(I,j) = (US%L_to_mCS%Visbeck_L_scale)2 G%areaCu(I,j)
1972		enddo ; enddo
1973	0	do J=js-1,Jeq ; do i=is,ie
1974	0	CS%L2v(i,J) = (US%L_to_mCS%Visbeck_L_scale)2 G%areaCv(i,J)
1975		enddo ; enddo
1976		else
1977	0	CS%L2u(:,:) = CS%Visbeck_L_scale**2
1978	0	CS%L2v(:,:) = CS%Visbeck_L_scale**2
1979		endif
1980
1981		CS%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, Time, &
1982		'Length scale squared for mixing coefficient, at u-points', &
1983	0	'm2', conversion=US%L_to_m**2)
1984		CS%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, Time, &
1985		'Length scale squared for mixing coefficient, at v-points', &
1986	0	'm2', conversion=US%L_to_m**2)
1987		endif
1988
1989		CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, &
1990	1	'Vertical structure of SQG mode', 'nondim')
1991		if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
1992	1	.or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then
1993	0	allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1994		endif
1995
1996	1	if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
1997		CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
1998	0	'Vertical structure of backscatter', 'nondim')
1999		endif
2000	1	if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
2001		CS%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, Time, &
2002	0	'Vertical structure of thickness diffusivity', 'nondim')
2003		endif
2004	1	if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then
2005		CS%id_khtr_struct = register_diag_field('ocean_model', 'khtr_struct', diag%axesTl, Time, &
2006	0	'Vertical structure of tracer diffusivity', 'nondim')
2007		endif
2008	1	if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then
2009		CS%id_kdgl90_struct = register_diag_field('ocean_model', 'kdgl90_struct', diag%axesTl, Time, &
2010	0	'Vertical structure of GL90 diffusivity', 'nondim')
2011		endif
2012
2013	1	if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) ) then
2014		CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, &
2015		'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
2016	1	's-2', conversion=(US%L_to_ZUS%s_to_T)*2)
2017		CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, &
2018		'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', &
2019	1	's-2', conversion=(US%L_to_ZUS%s_to_T)*2)
2020		endif
2021	1	if (CS%use_simpler_Eady_growth_rate) then
2022		CS%id_dzu = register_diag_field('ocean_model', 'dzu_Visbeck', diag%axesCui, Time, &
2023		'dz at u-points, used in calculating Eady growth rate in Visbeck et al..', &
2024	0	'm', conversion=US%Z_to_m)
2025		CS%id_dzv = register_diag_field('ocean_model', 'dzv_Visbeck', diag%axesCvi, Time, &
2026		'dz at v-points, used in calculating Eady growth rate in Visbeck et al..', &
2027	0	'm', conversion=US%Z_to_m)
2028		CS%id_dzSxN = register_diag_field('ocean_model', 'dzSxN', diag%axesCui, Time, &
2029		'dz * \|slope_x\| * N, used in calculating Eady growth rate in '//&
2030	0	'Visbeck et al..', 'm s-1', conversion=US%Z_to_m*US%s_to_T)
2031		CS%id_dzSyN = register_diag_field('ocean_model', 'dzSyN', diag%axesCvi, Time, &
2032		'dz * \|slope_y\| * N, used in calculating Eady growth rate in '//&
2033	0	'Visbeck et al..', 'm s-1', conversion=US%Z_to_m*US%s_to_T)
2034		endif
2035	1	if (CS%use_stored_slopes) then
2036		CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, &
2037		'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', &
2038	1	'nondim', conversion=US%Z_to_L**2)
2039		CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, &
2040		'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', &
2041	1	'nondim', conversion=US%Z_to_L**2)
2042		endif
2043
2044	1	oneOrTwo = 1.0
2045		CS%Resoln_scaling_used = CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. &
2046	1	CS%Resoln_scaled_KhTr .or. Resoln_scaled_MEKE_visc
2047	1	if (CS%Resoln_scaling_used) then
2048	1	CS%calculate_Rd_dx = .true.
2049	1	CS%calculate_res_fns = .true.
2050	8773	allocate(CS%Res_fn_h(isd:ied,jsd:jed), source=0.0)
2051	8971	allocate(CS%Res_fn_q(IsdB:IedB,JsdB:JedB), source=0.0)
2052	8841	allocate(CS%Res_fn_u(IsdB:IedB,jsd:jed), source=0.0)
2053	8902	allocate(CS%Res_fn_v(isd:ied,JsdB:JedB), source=0.0)
2054	8971	allocate(CS%beta_dx2_q(IsdB:IedB,JsdB:JedB), source=0.0)
2055	8841	allocate(CS%beta_dx2_u(IsdB:IedB,jsd:jed), source=0.0)
2056	8902	allocate(CS%beta_dx2_v(isd:ied,JsdB:JedB), source=0.0)
2057	8971	allocate(CS%f2_dx2_q(IsdB:IedB,JsdB:JedB), source=0.0)
2058	8841	allocate(CS%f2_dx2_u(IsdB:IedB,jsd:jed), source=0.0)
2059	8902	allocate(CS%f2_dx2_v(isd:ied,JsdB:JedB), source=0.0)
2060
2061		CS%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, Time, &
2062	1	'Resolution function for scaling diffusivities', 'nondim')
2063
2064		call get_param(param_file, mdl, "KH_RES_SCALE_COEF", CS%Res_coef_khth, &
2065		"A coefficient that determines how KhTh is scaled away if "//&
2066		"RESOLN_SCALED_... is true, as "//&
2067		"F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER).", &
2068	1	units="nondim", default=1.0)
2069		call get_param(param_file, mdl, "KH_RES_FN_POWER", CS%Res_fn_power_khth, &
2070		"The power of dx/Ld in the Kh resolution function. Any "//&
2071		"positive integer may be used, although even integers "//&
2072		"are more efficient to calculate. Setting this greater "//&
2073		"than 100 results in a step-function being used.", &
2074	1	default=2)
2075		call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", CS%Res_coef_visc, &
2076		"A coefficient that determines how Kh is scaled away if "//&
2077		"RESOLN_SCALED_... is true, as "//&
2078		"F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
2079		"This function affects lateral viscosity, Kh, and not KhTh.", &
2080	1	units="nondim", default=CS%Res_coef_khth)
2081		call get_param(param_file, mdl, "VISC_RES_FN_POWER", CS%Res_fn_power_visc, &
2082		"The power of dx/Ld in the Kh resolution function. Any "//&
2083		"positive integer may be used, although even integers "//&
2084		"are more efficient to calculate. Setting this greater "//&
2085		"than 100 results in a step-function being used. "//&
2086		"This function affects lateral viscosity, Kh, and not KhTh.", &
2087	1	default=CS%Res_fn_power_khth)
2088		call get_param(param_file, mdl, "INTERPOLATE_RES_FN", CS%interpolate_Res_fn, &
2089		"If true, interpolate the resolution function to the "//&
2090		"velocity points from the thickness points; otherwise "//&
2091		"interpolate the wave speed and calculate the resolution "//&
2092	1	"function independently at each point.", default=.false.)
2093	1	if (CS%interpolate_Res_fn) then
2094	0	if (CS%Res_coef_visc /= CS%Res_coef_khth) call MOM_error(FATAL, &
2095		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2096	0	"When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
2097	0	if (CS%Res_fn_power_visc /= CS%Res_fn_power_khth) call MOM_error(FATAL, &
2098		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2099	0	"When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
2100		endif
2101		call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", Gill_equatorial_Ld, &
2102		"If true, uses Gill's definition of the baroclinic "//&
2103		"equatorial deformation radius, otherwise, if false, use "//&
2104		"Pedlosky's definition. These definitions differ by a factor "//&
2105		"of 2 in front of the beta term in the denominator. Gill's "//&
2106	1	"is the more appropriate definition.", default=.true.)
2107	1	if (Gill_equatorial_Ld) then
2108	1	oneOrTwo = 2.0
2109		endif
2110
2111	7443	do J=js-1,Jeq ; do I=is-1,Ieq
2112		CS%f2_dx2_q(I,J) = ((G%dxBu(I,J)2) + (G%dyBu(I,J)2)) * &
2113	7381	max(G%Coriolis2Bu(I,J), absurdly_small_freq**2)
2114		CS%beta_dx2_q(I,J) = oneOrTwo * ((G%dxBu(I,J)2) + (G%dyBu(I,J)2)) * (sqrt(0.5 * &
2115		( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
2116		(((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2)) + &
2117		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
2118	7442	(((G%CoriolisBu(I,J+1)-G%CoriolisBu(I,J)) * G%IdyCu(I,j+1))**2)) ) ))
2119		enddo ; enddo
2120
2121	7321	do j=js,je ; do I=is-1,Ieq
2122		CS%f2_dx2_u(I,j) = ((G%dxCu(I,j)2) + (G%dyCu(I,j)2)) * &
2123	7260	max(0.5* (G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I,J-1)), absurdly_small_freq**2)
2124		CS%beta_dx2_u(I,j) = oneOrTwo * ((G%dxCu(I,j)2) + (G%dyCu(I,j)2)) * (sqrt( &
2125		((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2 + &
2126		0.25( ((((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) G%IdxCv(i,J-1))**2) + &
2127		(((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2)) + &
2128		((((G%CoriolisBu(I+1,J-1)-G%CoriolisBu(I,J-1)) * G%IdxCv(i+1,J-1))**2) + &
2129	7320	(((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2)) ) ))
2130		enddo ; enddo
2131
2132	7382	do J=js-1,Jeq ; do i=is,ie
2133		CS%f2_dx2_v(i,J) = ((G%dxCv(i,J)2) + (G%dyCv(i,J)2)) * &
2134	7320	max(0.5(G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I-1,J)), absurdly_small_freq*2)
2135		CS%beta_dx2_v(i,J) = oneOrTwo * ((G%dxCv(i,J)2) + (G%dyCv(i,J)2)) * (sqrt( &
2136		((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
2137		0.25( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) G%IdyCu(I,j))**2) + &
2138		(((G%CoriolisBu(I-1,J+1)-G%CoriolisBu(I-1,J)) * G%IdyCu(I-1,j+1))**2)) + &
2139		((((G%CoriolisBu(I,J+1)-G%CoriolisBu(I,J)) * G%IdyCu(I,j+1))**2) + &
2140	7381	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ) ))
2141		enddo ; enddo
2142
2143		endif
2144
2145	1	if (CS%Depth_scaled_KhTh) then
2146	0	CS%calculate_depth_fns = .true.
2147	0	allocate(CS%Depth_fn_u(IsdB:IedB,jsd:jed), source=0.0)
2148	0	allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB), source=0.0)
2149		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, &
2150		"The depth above which KHTH is scaled away.", &
2151	0	units="m", scale=US%m_to_Z, default=1000.)
2152		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, &
2153		"The exponent used in the depth dependent scaling function for KHTH.", &
2154	0	units="nondim", default=3.0)
2155		endif
2156
2157		! Resolution %Rd_dx_h
2158		CS%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, Time, &
2159	1	'Ratio between deformation radius and grid spacing', 'm m-1')
2160	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (CS%id_Rd_dx>0)
2161
2162	1	if (CS%calculate_Rd_dx) then
2163	1	CS%calculate_cg1 = .true. ! We will need %cg1
2164	8773	allocate(CS%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
2165	8773	allocate(CS%beta_dx2_h(isd:ied,jsd:jed), source=0.0)
2166	8773	allocate(CS%f2_dx2_h(isd:ied,jsd:jed), source=0.0)
2167	7627	do j=js-1,je+1 ; do i=is-1,ie+1
2168		CS%f2_dx2_h(i,j) = ((G%dxT(i,j)2) + (G%dyT(i,j)2)) * &
2169		max(0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
2170		(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))), &
2171	7564	absurdly_small_freq**2)
2172		CS%beta_dx2_h(i,j) = oneOrTwo * ((G%dxT(i,j)2) + (G%dyT(i,j)2)) * (sqrt(0.5 * &
2173		( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
2174		(((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2)) + &
2175		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
2176	7626	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ) ))
2177		enddo ; enddo
2178		endif
2179
2180	1	if (CS%calculate_cg1) then
2181	1	in_use = .true.
2182	8773	allocate(CS%cg1(isd:ied,jsd:jed), source=0.0)
2183		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2184		"This sets the default value for the various _ANSWER_DATE parameters.", &
2185	1	default=99991231)
2186		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2187		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
2188		"Values below 20190101 result in the use of older, less accurate expressions "//&
2189		"that were in use at the end of 2018. Higher values result in the use of more "//&
2190		"robust and accurate forms of mathematically equivalent expressions.", &
2191	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
2192	1	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2193
2194		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
2195		"The fractional tolerance for finding the wave speeds.", &
2196	1	units="nondim", default=0.001)
2197		!### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
2198		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
2199		"A floor in the first mode speed below which 0 used instead.", &
2200	1	units="m s-1", default=0.0, scale=US%m_s_to_L_T)
2201		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
2202		"If true, use a more robust estimate of the first mode wave speed as the "//&
2203	1	"starting point for iterations.", default=.true.)
2204		call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2205	1	do_not_log=.true., default=.true.)
2206		call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2207		"If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
2208		"See REMAPPING_USE_OM4_SUBCELLS for details. "//&
2209	1	"We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2210		call wave_speed_init(CS%wave_speed, GV, use_ebt_mode=CS%Resoln_use_ebt, &
2211		mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, &
2212		better_speed_est=better_speed_est, min_speed=wave_speed_min, &
2213	1	om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol)
2214		endif
2215
2216		! Leith parameters
2217		call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%use_QG_Leith_GM, &
2218		"If true, use the QG Leith viscosity as the GM coefficient.", &
2219	1	default=.false.)
2220
2221	1	if (CS%Use_QG_Leith_GM) then
2222		call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
2223		"The nondimensional Laplacian Leith constant, \n"//&
2224	0	"often set to 1.0", units="nondim", default=0.0)
2225
2226		call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_QG_Leith, &
2227		"If true, include the beta term in the Leith nonlinear eddy viscosity.", &
2228	0	default=.true.)
2229
2230	0	allocate(CS%Laplac3_const_u(IsdB:IedB,jsd:jed), source=0.0)
2231	0	allocate(CS%Laplac3_const_v(isd:ied,JsdB:JedB), source=0.0)
2232	0	allocate(CS%KH_u_QG(IsdB:IedB,jsd:jed,GV%ke), source=0.0)
2233	0	allocate(CS%KH_v_QG(isd:ied,JsdB:JedB,GV%ke), source=0.0)
2234
2235		! register diagnostics
2236		CS%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, Time, &
2237	0	'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2238		CS%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, Time, &
2239	0	'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2240
2241	0	do j=Jsq,Jeq+1 ; do I=is-1,Ieq
2242		! Static factors in the Leith schemes
2243	0	grid_sp_u2 = G%dyCu(I,j)*G%dxCu(I,j)
2244	0	grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
2245	0	CS%Laplac3_const_u(I,j) = Leith_Lap_const * grid_sp_u3
2246		enddo ; enddo
2247	0	do j=js-1,Jeq ; do I=Isq,Ieq+1
2248		! Static factors in the Leith schemes
2249	0	grid_sp_v2 = G%dyCv(i,J)*G%dxCv(i,J)
2250	0	grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
2251	0	CS%Laplac3_const_v(i,J) = Leith_Lap_const * grid_sp_v3
2252		enddo ; enddo
2253
2254	0	if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
2255		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2256	0	"USE_STORED_SLOPES must be True when using QG Leith.")
2257		endif
2258
2259		! Re-enable variable mixing if one of the schemes was enabled
2260	1	CS%use_variable_mixing = in_use .or. CS%use_variable_mixing
2261
2262	1	end subroutine VarMix_init
2263
2264		!> Destructor for VarMix control structure
2265	1	subroutine VarMix_end(CS)
2266		type(VarMix_CS), intent(inout) :: CS
2267
2268	1	if (allocated(CS%ebt_struct)) deallocate(CS%ebt_struct)
2269	1	if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct)
2270	1	if (allocated(CS%BS_struct)) deallocate(CS%BS_struct)
2271	1	if (allocated(CS%khth_struct)) deallocate(CS%khth_struct)
2272	1	if (allocated(CS%khtr_struct)) deallocate(CS%khtr_struct)
2273	1	if (allocated(CS%kdgl90_struct)) deallocate(CS%kdgl90_struct)
2274
2275	1	if (allocated(CS%slope_x)) deallocate(CS%slope_x)
2276	1	if (allocated(CS%slope_y)) deallocate(CS%slope_y)
2277
2278	1	if (allocated(CS%SN_u)) deallocate(CS%SN_u)
2279	1	if (allocated(CS%SN_v)) deallocate(CS%SN_v)
2280
2281	1	if (allocated(CS%L2u)) deallocate(CS%L2u)
2282	1	if (allocated(CS%L2v)) deallocate(CS%L2v)
2283
2284	1	if (allocated(CS%Res_fn_h)) deallocate(CS%Res_fn_h)
2285	1	if (allocated(CS%Res_fn_q)) deallocate(CS%Res_fn_q)
2286	1	if (allocated(CS%Res_fn_u)) deallocate(CS%Res_fn_u)
2287	1	if (allocated(CS%Res_fn_v)) deallocate(CS%Res_fn_v)
2288	1	if (allocated(CS%beta_dx2_q)) deallocate(CS%beta_dx2_q)
2289	1	if (allocated(CS%beta_dx2_u)) deallocate(CS%beta_dx2_u)
2290	1	if (allocated(CS%beta_dx2_v)) deallocate(CS%beta_dx2_v)
2291	1	if (allocated(CS%f2_dx2_q)) deallocate(CS%f2_dx2_q)
2292	1	if (allocated(CS%f2_dx2_u)) deallocate(CS%f2_dx2_u)
2293	1	if (allocated(CS%f2_dx2_v)) deallocate(CS%f2_dx2_v)
2294
2295	1	if (allocated(CS%Depth_fn_u)) deallocate(CS%Depth_fn_u)
2296	1	if (allocated(CS%Depth_fn_v)) deallocate(CS%Depth_fn_v)
2297
2298	1	if (allocated(CS%Rd_dx_h)) deallocate(CS%Rd_dx_h)
2299	1	if (allocated(CS%beta_dx2_h)) deallocate(CS%beta_dx2_h)
2300	1	if (allocated(CS%f2_dx2_h)) deallocate(CS%f2_dx2_h)
2301
2302	1	if (allocated(CS%cg1)) deallocate(CS%cg1)
2303
2304	1	if (allocated(CS%Laplac3_const_u)) deallocate(CS%Laplac3_const_u)
2305	1	if (allocated(CS%Laplac3_const_v)) deallocate(CS%Laplac3_const_v)
2306	1	if (allocated(CS%KH_u_QG)) deallocate(CS%KH_u_QG)
2307	1	if (allocated(CS%KH_v_QG)) deallocate(CS%KH_v_QG)
2308
2309	1	end subroutine VarMix_end
2310
2311		!> \namespace mom_lateral_mixing_coeffs
2312		!!
2313		!! This module provides a container for various factors used in prescribing diffusivities, that are
2314		!! a function of the state (in particular the stratification and isoneutral slopes).
2315		!!
2316		!! \section section_Resolution_Function The resolution function
2317		!!
2318		!! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
2319		!! The square of the resolution parameter is
2320		!!
2321		!! \f[
2322		!! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
2323		!! \f]
2324		!!
2325		!! where the grid spacing is calculated as
2326		!!
2327		!! \f[
2328		!! \Delta^2 = \Delta x^2 + \Delta y^2 .
2329		!! \f]
2330		!!
2331		!! \todo Check this reference to Bob on/off paper.
2332		!! The resolution function used in scaling diffusivities (\cite hallberg2013) is
2333		!!
2334		!! \f[
2335		!! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
2336		!! \f]
2337		!!
2338		!! The resolution function can be applied independently to thickness diffusion $module mom_thickness_diffuse$,
2339		!! tracer diffusion $mom_tracer_hordiff$ lateral viscosity $mom_hor_visc$.
2340		!!
2341		!! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
2342		!! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007
2343		!!
2344		!! \| Symbol \| Module parameter \|
2345		!! \| ------ \| --------------- \|
2346		!! \| - \| <code>USE_VARIABLE_MIXING</code> \|
2347		!! \| - \| <code>RESOLN_SCALED_KH</code> \|
2348		!! \| - \| <code>RESOLN_SCALED_KHTH</code> \|
2349		!! \| - \| <code>RESOLN_SCALED_KHTR</code> \|
2350		!! \| \f$ \alpha \f$ \| <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) \|
2351		!! \| \f$ p \f$ \| <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) \|
2352		!! \| \f$ \alpha \f$ \| <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) \|
2353		!! \| \f$ p \f$ \| <code>VISC_RES_FN_POWER</code> (for lateral viscosity) \|
2354		!! \| - \| <code>GILL_EQUATORIAL_LD</code> \|
2355		!!
2356		!!
2357		!!
2358		!! \section section_Vicbeck Visbeck diffusivity
2359		!!
2360		!! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
2361		!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse but calculated in this module.
2362		!!
2363		!! \f[
2364		!! \kappa_h = \alpha_s L_s^2 S N
2365		!! \f]
2366		!!
2367		!! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
2368		!!
2369		!! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
2370		!! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
2371		!!
2372		!! \| Symbol \| Module parameter \|
2373		!! \| ------ \| --------------- \|
2374		!! \| - \| <code>USE_VARIABLE_MIXING</code> \|
2375		!! \| \f$ \alpha_s \f$ \| <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)\|
2376		!! \| \f$ \alpha_s \f$ \| <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)\|
2377		!! \| \f$ L_{s} \f$ \| <code>VISBECK_L_SCALE</code> \|
2378		!! \| \f$ S_{max} \f$ \| <code>VISBECK_MAX_SLOPE</code> \|
2379		!!
2380		!!
2381		!! \section section_vertical_structure_khth Vertical structure function for KhTh
2382		!!
2383		!! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
2384		!! velocity mode. The structure function is stored in the control structure for this module (varmix_cs) but is
2385		!! calculated using subroutines in mom_wave_speed.
2386		!!
2387		!! \| Symbol \| Module parameter \|
2388		!! \| ------ \| --------------- \|
2389		!! \| - \| <code>KHTH_USE_EBT_STRUCT</code> \|
2390
2391	0	end module MOM_lateral_mixing_coeffs