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		integer :: isoneutral_halo
664
665	0	niblock = CS%niblock
666	0	njblock = CS%njblock
667	0	nkblock = CS%nkblock
668
669	0	isoneutral_halo = 1
670
671	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
672	0	f_subround = 1.0e-40 * US%s_to_T
673
674	0	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
675	0	"Module must be initialized before it is used.")
676
677	0	if (CS%sqg_expo <= 0.) then
678	0	CS%sqg_struct(:,:,:) = 1.
679		else
680	0	if (allocated(MEKE%Le)) then
681	0	do j=js,je ; do i=is,ie
682	0	Le(i,j) = MEKE%Le(i,j)
683		enddo ; enddo
684		else
685	0	do j=js,je ; do i=is,ie
686	0	Le(i,j) = sqrt(G%areaT(i,j))
687		enddo ; enddo
688		endif
689
690	0	do j=js,je ; do i=is,ie
691		! Setting the structure averaged over the top layer to 1 is consistent with it being well mixed.
692	0	CS%sqg_struct(i,j,1) = 1.0
693		enddo ; enddo
694
695	0	if (CS%interpolated_sqg_struct) then
696	0	do j=js,je ; do i=is,ie
697		f(i,j) = max(0.25 * abs((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + &
698	0	(G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))), f_subround)
699		enddo ; enddo
700		!$omp target update to(h)
701		!$omp target enter data map(alloc: e)
702	0	call find_eta(h, tv, G, GV, US, e, halo_size=2) !### Could be halo_size=1?
703
704		! +2 needed because loops in calc_isoneutral slopes run from is-1 to ie and js-1 to je
705	0	if (niblock == 0) niblock = ie - is + 2 + 2*isoneutral_halo
706	0	if (njblock == 0) njblock = je - js + 2 + 2*isoneutral_halo
707	0	if (nkblock == 0) nkblock = nz
708
709		!$omp target enter data map(to: tv, tv%T, tv%S)
710		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
711		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
712		!$omp target enter data map(alloc: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
713		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
714		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
715		N2_v=N2_v, dzu=dzu, dzv=dzv, dzSxN=dzSxN, dzSyN=dzSyN, &
716	0	halo=isoneutral_halo, OBC=OBC, OBC_N2=CS%OBC_friendly)
717		!$omp target exit data map(release: tv%T, tv%S, tv)
718		!$omp target exit data map(delete: e)
719		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
720		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
721		!$omp target exit data map(from: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
722	0	do k=2,nz ; do j=js,je ; do i=is,ie
723	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)
724	0	dzc = 0.25 * ((dzu(I-1,j,K) + dzu(I,j,K)) + (dzv(i,J-1,K) + dzv(i,J,K)))
725		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
726	0	exp(-CS%sqg_expo * (dzc * sqrt(N2)/(f(i,j) * Le(i,j))))
727		enddo ; enddo ; enddo
728		else
729	0	do j=js,je ; do i=is,ie
730		I_f_Le(i,j) = 1.0 / &
731		(Le(i,j) * max(0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + &
732	0	(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1)))), f_subround))
733		enddo ; enddo
734
735	0	call thickness_to_dz(h, tv, dz, G, GV, US)
736
737	0	if (associated(tv%eqn_of_state)) then
738	0	EOSdom(:) = EOS_domain(G%HI)
739	0	H_to_pres = GV%H_to_RZ * GV%g_Earth
740		! Set the pressure at the topmost interior interface.
741	0	p_i(:,:) = 0.0
742	0	if (associated(tv%p_surf)) then
743	0	do j=js,je ; do i=is,ie ; p_i(i,j) = tv%p_surf(i,j) ; enddo ; enddo
744		endif
745	0	if (.not.allocated(tv%SpV_avg)) GxSpV = GV%g_Earth / GV%Rho0
746	0	do K=2,nz ; do j=js,je
747		! Find the derivatives of density with T and S at the interface.
748	0	do i=is,ie
749	0	p_i(i,j) = p_i(i,j) + H_to_pres * h(i,j,k-1)
750	0	T_i(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
751	0	S_i(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
752		enddo
753	0	call calculate_density_derivs(T_i, S_i, p_i(:,j), dRho_dT, dRho_dS, tv%eqn_of_state, EOSdom)
754
755	0	do i=is,ie
756	0	if (allocated(tv%SpV_avg)) & ! GxSpV is in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
757	0	GxSpV = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
758
759		drdk = max(dRho_dT(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
760	0	dRho_dS(i) * (tv%S(i,j,k)-tv%S(i,j,k-1)), 0.0) ! Density difference [R ~> kg m-3]
761	0	dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
762		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
763	0	exp(-CS%sqg_expo * (sqrt((GxSpV * drdk) * dz_int) * I_f_Le(i,j)) )
764		! To derive the expression above, note that
765		! N2 = GxSpV * drdk / dzh(i,j,K) ! Square of positive buoyancy freq. [L2 Z-2 T-2 ~> s-2]
766		! CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
767		! exp(-CS%sqg_expo * (dz_int(i,j,K) * sqrt(N2) * I_f_Le(i,j)) )
768		enddo
769		enddo ; enddo
770		else ! (GV%Boussinesq .and. .not.use_EOS) then
771	0	do K=2,nz ; do j=js,je ; do i=is,ie
772	0	dz_int = 0.5*(dz(i,j,k-1) + dz(i,j,k)) ! Thickness around interface [Z ~> m]
773		CS%sqg_struct(i,j,k) = CS%sqg_struct(i,j,k-1) * &
774	0	exp(-CS%sqg_expo * (sqrt(GV%g_prime(K) * dz_int) * I_f_Le(i,j)) )
775		enddo ; enddo ; enddo
776		endif
777		endif
778		endif
779
780	0	if (query_averaging_enabled(CS%diag)) then
781	0	if (CS%id_sqg_struct > 0) call post_data(CS%id_sqg_struct, CS%sqg_struct, CS%diag)
782	0	if (CS%interpolated_sqg_struct .and. (CS%sqg_expo > 0.)) then
783	0	if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
784	0	if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
785		endif
786		endif
787
788	0	end subroutine calc_sqg_struct
789
790		!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
791		!! style scaling of diffusivity
792	12	subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
793		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
794		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
795		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
796		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
797		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
798		real, intent(in) :: dt !< Time increment [T ~> s]
799		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
800		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
801
802		! Local variables
803	24	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
804	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]
805	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]
806	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m]
807	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m]
808	24	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! \|Sx\| N times dz at u-points [Z T-1 ~> m s-1]
809	24	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzSyN ! \|Sy\| N times dz at v-points [Z T-1 ~> m s-1]
810		integer :: niblock, njblock, nkblock
811		integer :: isoneutral_halo
812
813	12	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
814	0	"Module must be initialized before it is used.")
815
816	12	niblock = CS%niblock
817	12	njblock = CS%njblock
818	12	nkblock = CS%nkblock
819
820	12	isoneutral_halo = 1
821
822		! +2 needed because loops in calc_isoneutral slopes run from is-1 to ie and js-1 to je
823	12	if (niblock == 0) niblock = G%iec - G%isc + 2 + 2*isoneutral_halo
824	12	if (njblock == 0) njblock = G%jec - G%jsc + 2 + 2*isoneutral_halo
825	12	if (nkblock == 0) nkblock = GV%ke
826
827	12	if (CS%calculate_Eady_growth_rate) then
828		!$omp target update to(h)
829		!$omp target enter data map(alloc: e)
830	12	call find_eta(h, tv, G, GV, US, e, halo_size=2)
831
832	12	if (CS%use_simpler_Eady_growth_rate) then
833		!$omp target enter data map(to: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
834		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
835		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
836		!$omp target enter data map(alloc: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
837		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
838		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
839		N2_v=N2_v, dzu=dzu, dzv=dzv, dzSxN=dzSxN, dzSyN=dzSyN, &
840	0	halo=isoneutral_halo, OBC=OBC, OBC_N2=CS%OBC_friendly)
841		!$omp target update from(e, dzu, dzv, dzSxN, dzSyN)
842	0	call calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, CS%SN_u, CS%SN_v)
843		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
844		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
845		!$omp target exit data map(release: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
846		!$omp target exit data map(delete: N2_u, N2_v, dzu, dzv, dzSxN, dzSyN)
847	12	elseif (CS%use_stored_slopes) then
848		!$omp target enter data map(to: tv, tv%T, tv%S, CS%slope_x, CS%slope_y)
849		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
850		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
851		!$omp target enter data map(alloc: N2_u, N2_v)
852		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
853		CS%slope_x, CS%slope_y, niblock, njblock, nkblock, N2_u=N2_u, &
854	12	N2_v=N2_v, halo=isoneutral_halo, OBC=OBC, OBC_N2=CS%OBC_friendly)
855		!$omp target exit data map(from: CS%slope_x, CS%slope_y, N2_u, N2_v)
856	12	call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
857		!$omp target exit data map(release: tv%T, tv%S, tv, CS%slope_x, CS%slope_y)
858		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
859		!$omp target exit data map(delete: N2_u, N2_v )
860		else
861		!$omp target update from(e)
862	0	call calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
863		endif
864		!$omp target exit data map(delete: e)
865		endif
866
867	12	if (query_averaging_enabled(CS%diag)) then
868	12	if (CS%id_dzu > 0) call post_data(CS%id_dzu, dzu, CS%diag)
869	12	if (CS%id_dzv > 0) call post_data(CS%id_dzv, dzv, CS%diag)
870	12	if (CS%id_dzSxN > 0) call post_data(CS%id_dzSxN, dzSxN, CS%diag)
871	12	if (CS%id_dzSyN > 0) call post_data(CS%id_dzSyN, dzSyN, CS%diag)
872	12	if (CS%id_SN_u > 0) call post_data(CS%id_SN_u, CS%SN_u, CS%diag)
873	12	if (CS%id_SN_v > 0) call post_data(CS%id_SN_v, CS%SN_v, CS%diag)
874	12	if (CS%id_L2u > 0) call post_data(CS%id_L2u, CS%L2u, CS%diag)
875	12	if (CS%id_L2v > 0) call post_data(CS%id_L2v, CS%L2v, CS%diag)
876	12	if (CS%id_N2_u > 0) call post_data(CS%id_N2_u, N2_u, CS%diag)
877	12	if (CS%id_N2_v > 0) call post_data(CS%id_N2_v, N2_v, CS%diag)
878		endif
879
880	12	end subroutine calc_slope_functions
881
882		!> Calculates factors used when setting diffusivity coefficients similar to Visbeck et al., 1997.
883		!! This is on older implementation that is susceptible to large values of Eady growth rate
884		!! for incropping layers.
885	12	subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
886		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
887		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
888		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
889		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope [Z L-1 ~> nondim]
890		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency
891		!! at u-points [L2 Z-2 T-2 ~> s-2]
892		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope
893		!! [Z L-1 ~> nondim]
894		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency
895		!! at v-points [L2 Z-2 T-2 ~> s-2]
896		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
897		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
898		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
899
900		! Local variables
901		real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
902		real :: N2 ! Positive buoyancy frequency or zero [L2 Z-2 T-2 ~> s-2]
903		real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
904		real :: H_geom ! The geometric mean of Hup and Hdn [H ~> m or kg m-2].
905		real :: S2max ! An upper bound on the squared slopes [Z2 L-2 ~> nondim]
906		real :: wNE, wSE, wSW, wNW ! Weights of adjacent points [nondim]
907	24	real :: H_u(SZIB_(G)), H_v(SZI_(G)) ! Layer thicknesses at u- and v-points [H ~> m or kg m-2]
908
909		! Note that at some points in the code S2_u and S2_v hold the running depth
910		! integrals of the squared slope [H ~> m or kg m-2] before the average is taken.
911	24	real :: S2_u(SZIB_(G),SZJ_(G)) ! At first the thickness-weighted depth integral of the squared
912		! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
913		! squared slope [Z2 L-2 ~> nondim] at u points.
914	24	real :: S2_v(SZI_(G),SZJB_(G)) ! At first the thickness-weighted depth integral of the squared
915		! slope [H Z2 L-2 ~> m or kg m-2] and then the average of the
916		! squared slope [Z2 L-2 ~> nondim] at v points.
917	24	integer :: OBC_dir_u(SZIB_(G),SZJ_(G)) ! An integer indicating where there are u OBCs: +1 for
918		! eastern OBCs, -1 for western OBCs and 0 at points with no OBCs.
919	24	integer :: OBC_dir_v(SZI_(G),SZJB_(G)) ! An integer indicating where there are v OBCs: +1 for
920		! northern OBCs, -1 for southern OBCs and 0 at points with no OBCs.
921	12	real :: h4_u(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a u-point
922		! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
923	12	real :: h4_v(SZI_(G),SZJB_(G),SZK_(GV)+1) ! The product of the 4 thicknesses surrounding a v-point
924		! interface or the inward equivalent with OBCs [H4 ~> m4 or kg4 m-8]
925		integer :: i, j, k, is, ie, js, je, nz
926
927	12	if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// &
928	0	"Module must be initialized before it is used.")
929
930	12	if (.not. CS%calculate_Eady_growth_rate) return
931	12	if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function: "// &
932	0	"%SN_u is not associated with use_variable_mixing.")
933	12	if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function: "// &
934	0	"%SN_v is not associated with use_variable_mixing.")
935
936	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
937
938	12	S2max = CS%Visbeck_S_max**2
939
940	106092	CS%SN_u(:,:) = 0.0
941	106824	CS%SN_v(:,:) = 0.0
942
943		! These settings apply where there are not open boundary conditions.
944	212904	OBC_dir_u(:,:) = 0 ; OBC_dir_v(:,:) = 0
945
946	12	if (associated(OBC) .and. CS%OBC_friendly) then
947		! Store the direction of any OBC faces.
948		!$OMP parallel do default(shared)
949	0	do j=js-1,je+1 ; do I=is-1,ie ; if (OBC%segnum_u(I,j) /= 0) then
950	0	if (OBC%segnum_u(I,j) > 0) OBC_dir_u(I,j) = 1 ! OBC_DIRECTION_E
951	0	if (OBC%segnum_u(I,j) < 0) OBC_dir_u(I,j) = -1 ! OBC_DIRECTION_W
952		endif ; enddo ; enddo
953		!$OMP parallel do default(shared)
954	0	do J=js-1,je ; do i=is-1,ie+1 ; if (OBC%segnum_v(i,J) /= 0) then
955	0	if (OBC%segnum_v(i,J) > 0) OBC_dir_v(i,J) = 1 ! OBC_DIRECTION_N
956	0	if (OBC%segnum_v(i,J) < 0) OBC_dir_v(i,J) = -1 ! OBC_DIRECTION_S
957		endif ; enddo ; enddo
958
959		! Use the masked product of the 4 (or 2) thicknesses around a velocity-point interface for weights.
960		!$OMP parallel do default(shared)
961	0	do K=2,nz
962	0	do j=js-1,je+1 ; do I=is-1,ie
963	0	if (OBC_dir_u(I,j) == 0) then
964	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)) )
965	0	elseif (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E
966	0	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i,j,k)*2) (h(i,j,k-1)**2) )
967	0	elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W
968	0	h4_u(I,j,K) = G%mask2dCu(I,j) * ( (h(i+1,j,k)*2) (h(i+1,j,k-1)**2) )
969		endif
970		enddo ; enddo
971	0	do J=js-1,je ; do i=is-1,ie+1
972	0	if (OBC_dir_v(i,J) == 0) then
973	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)) )
974	0	elseif (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N
975	0	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j,k)*2) (h(i,j,k-1)**2) )
976	0	elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S
977	0	h4_v(i,J,K) = G%mask2dCv(i,J) * ( (h(i,j+1,k)*2) (h(i,j+1,k-1)**2) )
978		endif
979		enddo ; enddo
980		enddo
981		else ! The land mask is sufficient and there are no special considerations taken at OBC points.
982		! Use the masked product of the 4 thicknesses around a velocity-point interface for weights.
983		!$OMP parallel do default(shared)
984	900	do K=2,nz
985	6717720	do j=js-1,je+1 ; do I=is-1,ie
986	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)) )
987		enddo ; enddo
988	6663564	do J=js-1,je ; do i=is-1,ie+1
989	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)) )
990		enddo ; enddo
991		enddo
992		endif
993
994		! To set the length scale based on the deformation radius, use wave_speed to
995		! calculate the first-mode gravity wave speed and then blend the equatorial
996		! and midlatitude deformation radii, using calc_resoln_function as a template.
997
998		!$OMP parallel do default(shared) private(S2,H_u,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
999	732	do j=js,je
1000	87840	do I=is-1,ie
1001	87840	CS%SN_u(I,j) = 0. ; H_u(I) = 0. ; S2_u(I,j) = 0.
1002		enddo
1003	6500880	do K=2,nz ; do I=is-1,ie
1004	6446880	Hdn = sqrt( h(i,j,k) * h(i+1,j,k) )
1005	6446880	Hup = sqrt( h(i,j,k-1) * h(i+1,j,k-1) )
1006	6446880	H_geom = sqrt( Hdn * Hup )
1007		!H_geom = H_geom * sqrt(N2) ! WKB-ish
1008		!H_geom = H_geom * N2 ! WKB-ish
1009	6446880	wSE = h4_v(i+1,J-1,K)
1010	6446880	wNW = h4_v(i,J,K)
1011	6446880	wNE = h4_v(i+1,J,K)
1012	6446880	wSW = h4_v(i,J-1,K)
1013	6446880	if (OBC_dir_u(I,j) == 1) then ! OBC_DIRECTION_E
1014	0	wSE = 0.0 ; wNE = 0.0
1015	0	H_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
1016	6446880	elseif (OBC_dir_u(I,j) == -1) then ! OBC_DIRECTION_W
1017	0	wSW = 0.0 ; wNW = 0.0
1018	0	H_geom = sqrt( h(i+1,j,k) * h(i+1,j,k-1) )
1019		endif
1020		S2 = slope_x(I,j,K)**2 + &
1021		(((wNWslope_y(i,J,K)2) + (wSEslope_y(i+1,J-1,K)**2)) + &
1022		((wNEslope_y(i+1,J,K)2) + (wSWslope_y(i,J-1,K)**2)) ) / &
1023	6446880	( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 )
1024	6446880	if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2
1025
1026	6446880	N2 = max(0., N2_u(I,j,k))
1027	6446880	CS%SN_u(I,j) = CS%SN_u(I,j) + sqrt( S2N2 )H_geom
1028	6446880	S2_u(I,j) = S2_u(I,j) + S2*H_geom
1029	6500160	H_u(I) = H_u(I) + H_geom
1030		enddo ; enddo
1031	87852	do I=is-1,ie
1032	87840	if (H_u(I)>0.) then
1033	87120	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * CS%SN_u(I,j) / H_u(I)
1034	87120	S2_u(I,j) = G%OBCmaskCu(I,j) * S2_u(I,j) / H_u(I)
1035		else
1036	0	CS%SN_u(I,j) = 0.
1037		endif
1038		enddo
1039		enddo
1040
1041		!$OMP parallel do default(shared) private(S2,H_v,Hdn,Hup,H_geom,N2,wNE,wSE,wSW,wNW)
1042	744	do J=js-1,je
1043	88572	do i=is,ie
1044	88572	CS%SN_v(i,J) = 0. ; H_v(i) = 0. ; S2_v(i,J) = 0.
1045		enddo
1046	6555060	do K=2,nz ; do i=is,ie
1047	6500160	Hdn = sqrt( h(i,j,k) * h(i,j+1,k) )
1048	6500160	Hup = sqrt( h(i,j,k-1) * h(i,j+1,k-1) )
1049	6500160	H_geom = sqrt( Hdn * Hup )
1050		!H_geom = H_geom * sqrt(N2) ! WKB-ish
1051		!H_geom = H_geom * N2 ! WKB-ish
1052	6500160	wSE = h4_u(I,j,K)
1053	6500160	wNW = h4_u(I-1,j+1,K)
1054	6500160	wNE = h4_u(I,j+1,K)
1055	6500160	wSW = h4_u(I-1,j,K)
1056	6500160	if (OBC_dir_v(i,J) == 1) then ! OBC_DIRECTION_N
1057	0	wNW = 0.0 ; wNE = 0.0
1058	0	H_geom = sqrt( h(i,j,k) * h(i,j,k-1) )
1059	6500160	elseif (OBC_dir_v(i,J) == -1) then ! OBC_DIRECTION_S
1060	0	wSW = 0.0 ; wSE = 0.0
1061	0	H_geom = sqrt( h(i,j+1,k) * h(i,j+1,k-1) )
1062		endif
1063		S2 = slope_y(i,J,K)**2 + &
1064		(((wSEslope_x(I,j,K)2) + (wNWslope_x(I-1,j+1,K)**2)) + &
1065		((wNEslope_x(I,j+1,K)2) + (wSWslope_x(I-1,j,K)**2)) ) / &
1066	6500160	( ((wSE+wNW) + (wNE+wSW)) + GV%H_subroundoff**4 )
1067	6500160	if (S2max>0.) S2 = S2 * S2max / (S2 + S2max) ! Limit S2
1068
1069	6500160	N2 = max(0., N2_v(i,J,K))
1070	6500160	CS%SN_v(i,J) = CS%SN_v(i,J) + sqrt( S2N2 )H_geom
1071	6500160	S2_v(i,J) = S2_v(i,J) + S2*H_geom
1072	6554328	H_v(i) = H_v(i) + H_geom
1073		enddo ; enddo
1074	88584	do i=is,ie
1075	88572	if (H_v(i)>0.) then
1076	87840	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * CS%SN_v(i,J) / H_v(i)
1077	87840	S2_v(i,J) = G%OBCmaskCv(i,J) * S2_v(i,J) / H_v(i)
1078		else
1079	0	CS%SN_v(i,J) = 0.
1080		endif
1081		enddo
1082		enddo
1083
1084		! Offer diagnostic fields for averaging.
1085	12	if (query_averaging_enabled(CS%diag)) then
1086	12	if (CS%id_S2_u > 0) call post_data(CS%id_S2_u, S2_u, CS%diag)
1087	12	if (CS%id_S2_v > 0) call post_data(CS%id_S2_v, S2_v, CS%diag)
1088		endif
1089
1090	12	if (CS%debug) then
1091		call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, &
1092	0	unscale=US%Z_to_L, haloshift=1)
1093		! call uvchksum("calc_Visbeck_coeffs_old S2_[uv]", S2_u, S2_v, G%HI, &
1094		! unscale=US%Z_to_L**2, scalar_pair=.true.)
1095		call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, &
1096	0	unscale=US%L_to_Z*2US%s_to_T**2, scalar_pair=.true.)
1097		call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
1098	0	unscale=US%s_to_T, scalar_pair=.true.)
1099		endif
1100
1101		end subroutine calc_Visbeck_coeffs_old
1102
1103		!> Calculates the Eady growth rate (2D fields) for use in MEKE and the Visbeck schemes
1104	0	subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN, SN_u, SN_v)
1105		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1106		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1107		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1108		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1109		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m]
1110		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m]
1111		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m]
1112		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzv !< dz at v-points [Z ~> m]
1113		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzSxN !< dz Sx N at u-points [Z T-1 ~> m s-1]
1114		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzSyN !< dz Sy N at v-points [Z T-1 ~> m s-1]
1115		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: SN_u !< SN at u-points [T-1 ~> s-1]
1116		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: SN_v !< SN at v-points [T-1 ~> s-1]
1117		! Local variables
1118		real :: D_scale ! The depth over which to average SN [Z ~> m]
1119		real :: dnew ! Depth of bottom of layer [Z ~> m]
1120		real :: dz ! Limited thickness of this layer [Z ~> m]
1121		real :: weight ! Fraction of this layer that contributes to integral [nondim]
1122	0	real :: sum_dz(SZI_(G)) ! Cumulative sum of z-thicknesses [Z ~> m]
1123	0	real :: vint_SN(SZIB_(G)) ! Cumulative integral of SN [Z T-1 ~> m s-1]
1124	0	real, dimension(SZIB_(G),SZJ_(G)) :: SN_cpy !< SN at u-points [T-1 ~> s-1]
1125		real :: dz_neglect ! A negligibly small distance to avoid division by zero [Z ~> m]
1126		real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~> m-1]
1127		real :: dB, dT ! Elevation variables used when cropping [Z ~> m]
1128		integer :: i, j, k
1129		logical :: crop
1130
1131	0	dz_neglect = GV%dZ_subroundoff
1132	0	D_scale = CS%Eady_GR_D_scale
1133	0	if (D_scale<=0.) D_scale = 64.*GV%max_depth ! 0 means use full depth so choose something big
1134	0	r_crp_dist = 1. / max( dz_neglect, CS%cropping_distance )
1135	0	crop = CS%cropping_distance>=0. ! Only filter out in-/out-cropped interface is parameter if non-negative
1136
1137	0	if (CS%debug) then
1138	0	call uvchksum("calc_Eady_growth_rate_2D dz[uv]", dzu, dzv, G%HI, unscale=US%Z_to_m, scalar_pair=.true.)
1139		call uvchksum("calc_Eady_growth_rate_2D dzS2N2[uv]", dzSxN, dzSyN, G%HI, &
1140	0	unscale=US%Z_to_m*US%s_to_T, scalar_pair=.true.)
1141		endif
1142
1143		!$OMP parallel do default(shared)
1144	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1145	0	CS%SN_u(i,j) = 0.0
1146	0	CS%SN_v(i,j) = 0.0
1147		enddo ; enddo
1148
1149		!$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1150	0	do j=G%jsc-1,G%jec+1
1151	0	do I=G%isc-1,G%iec
1152	0	vint_SN(I) = 0.
1153	0	sum_dz(I) = dz_neglect
1154		enddo
1155	0	if (crop) then
1156	0	do K=2,GV%ke ; do I=G%isc-1,G%iec
1157	0	dnew = sum_dz(I) + dzu(I,j,K) ! This is where the bottom of the layer is
1158	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1159	0	dz = max(0., dnew - sum_dz(I)) ! This is the part of the layer to be included in the integral.
1160		! When D_scale>dnew, dz=dzu (+roundoff error).
1161		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1162		! When D_scale<sum_dz, dz=0.
1163	0	weight = dz / ( dzu(I,j,K) + dz_neglect ) ! Fraction of this layer to include
1164	0	dT = min( e(i,j,1), e(i+1,j,1) ) ! Deepest sea surface
1165	0	dB = max( e(i,j,K), e(i+1,j,K) ) ! Shallowest interface
1166	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1167	0	dT = min( e(i,j,K), e(i+1,j,K) ) ! Deepest interface
1168	0	dB = max( e(i,j,GV%ke+1), e(i+1,j,GV%ke+1) ) ! Shallowest topography
1169	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1170	0	vint_SN(I) = vint_SN(I) + weight * dzSxN(I,j,K)
1171	0	sum_dz(I) = sum_dz(I) + weight * dzu(I,j,K)
1172		enddo ; enddo
1173		else
1174	0	do K=2,GV%ke ; do I=G%isc-1,G%iec
1175	0	dnew = sum_dz(I) + dzu(I,j,K) ! This is where the bottom of the layer is
1176	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1177	0	dz = max(0., dnew - sum_dz(I)) ! This is the part of the layer to be included in the integral.
1178		! When D_scale>dnew, dz=dzu (+roundoff error).
1179		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1180		! When D_scale<sum_dz, dz=0.
1181	0	weight = dz / ( dzu(I,j,K) + dz_neglect ) ! Fraction of this layer to include
1182	0	vint_SN(I) = vint_SN(I) + weight * dzSxN(I,j,K)
1183	0	sum_dz(I) = sum_dz(I) + weight * dzu(I,j,K)
1184		enddo ; enddo
1185		endif
1186	0	do I=G%isc-1,G%iec
1187	0	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) )
1188	0	SN_cpy(I,j) = G%OBCmaskCu(I,j) * ( vint_SN(I) / sum_dz(I) )
1189		enddo
1190		enddo
1191
1192		!$OMP parallel do default(shared) private(dnew,dz,weight,vint_SN,sum_dz,dT,dB)
1193	0	do J=G%jsc-1,G%jec
1194	0	do i=G%isc-1,G%iec+1
1195	0	vint_SN(i) = 0.
1196	0	sum_dz(i) = dz_neglect
1197		enddo
1198	0	if (crop) then
1199	0	do K=2,GV%ke ; do i=G%isc-1,G%iec+1
1200	0	dnew = sum_dz(i) + dzv(i,J,K) ! This is where the bottom of the layer is
1201	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1202	0	dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1203		! When D_scale>dnew, dz=dzu (+roundoff error).
1204		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1205		! When D_scale<sum_dz, dz=0.
1206	0	weight = dz / ( dzv(i,J,K) + dz_neglect ) ! Fraction of this layer to include
1207	0	dT = min( e(i,j,1), e(i,j+1,1) ) ! Deepest sea surface
1208	0	dB = max( e(i,j,K), e(i,j+1,K) ) ! Shallowest interface
1209	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1210	0	dT = min( e(i,j,K), e(i,j+1,K) )! Deepest interface
1211	0	dB = max( e(i,j,GV%ke+1), e(i,j+1,GV%ke+1) ) ! Shallowest topography
1212	0	weight = weight * min( max( 0., (dT-dB)*r_crp_dist ), 1. )
1213	0	vint_SN(I) = vint_SN(I) + weight*2 dzSyN(i,J,K)
1214	0	sum_dz(i) = sum_dz(i) + weight * dzv(i,J,K)
1215		enddo ; enddo
1216		else
1217	0	do K=2,GV%ke ; do i=G%isc-1,G%iec+1
1218	0	dnew = sum_dz(i) + dzv(i,J,K) ! This is where the bottom of the layer is
1219	0	dnew = min(dnew, D_scale) ! This limits the depth to D_scale
1220	0	dz = max(0., dnew - sum_dz(i)) ! This is the part of the layer to be included in the integral.
1221		! When D_scale>dnew, dz=dzu (+roundoff error).
1222		! When sum_dz<D_scale<dnew, 0<dz<dzu.
1223		! When D_scale<sum_dz, dz=0.
1224	0	weight = dz / ( dzv(i,J,K) + dz_neglect ) ! Fraction of this layer to include
1225	0	vint_SN(I) = vint_SN(I) + weight*2 dzSyN(i,J,K)
1226	0	sum_dz(i) = sum_dz(i) + weight * dzv(i,J,K)
1227		enddo ; enddo
1228		endif
1229	0	do i=G%isc-1,G%iec+1
1230	0	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * ( vint_SN(i) / sum_dz(i) )
1231		enddo
1232		enddo
1233
1234	0	do j=G%jsc,G%jec
1235	0	do I=G%isc-1,G%iec
1236		CS%SN_u(I,j) = sqrt( SN_cpy(I,j)**2 &
1237		+ 0.25( ((CS%SN_v(i,J)2) + (CS%SN_v(i+1,J-1)*2)) &
1238	0	+ ((CS%SN_v(i+1,J)2) + (CS%SN_v(i,J-1)2)) ) )
1239		enddo
1240		enddo
1241	0	do J=G%jsc-1,G%jec
1242	0	do i=G%isc,G%iec
1243		CS%SN_v(i,J) = sqrt( CS%SN_v(i,J)**2 &
1244		+ 0.25( ((SN_cpy(I,j)2) + (SN_cpy(I-1,j+1)*2)) &
1245	0	+ ((SN_cpy(I,j+1)2) + (SN_cpy(I-1,j)2)) ) )
1246		enddo
1247		enddo
1248
1249	0	if (CS%debug) then
1250		call uvchksum("calc_Eady_growth_rate_2D SN_[uv]", CS%SN_u, CS%SN_v, G%HI, &
1251	0	unscale=US%s_to_T, scalar_pair=.true.)
1252		endif
1253
1254	0	end subroutine calc_Eady_growth_rate_2D
1255
1256		!> The original calc_slope_function() that calculated slopes using
1257		!! interface positions only, not accounting for density variations.
1258	0	subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
1259		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1260		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1261		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1262		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1263		type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
1264		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m]
1265		! type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1266		! Local variables
1267	0	real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics)
1268	0	real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics)
1269	0	real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m]
1270		! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m]
1271		real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
1272		real :: dZ_cutoff ! A minimum water column depth for masking [H ~> m or kg m-2]
1273		real :: h1, h2 ! Temporary total thicknesses [Z ~> m]
1274		real :: h_neglect ! A thickness that is so small it is usually lost
1275		! in roundoff and can be neglected [H ~> m or kg m-2].
1276		real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
1277		real :: N2 ! Brunt-Vaisala frequency squared [L2 Z-2 T-2 ~> s-2]
1278		real :: Hup, Hdn ! Thickness from above, below [H ~> m or kg m-2]
1279		real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2].
1280	0	real :: S2N2_u_local(SZIB_(G),SZJ_(G),SZK_(GV)) ! The depth integral of the slope times
1281		! the buoyancy frequency squared at u-points [Z T-2 ~> m s-2]
1282	0	real :: S2N2_v_local(SZI_(G),SZJB_(G),SZK_(GV)) ! The depth integral of the slope times
1283		! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2]
1284		logical :: use_dztot ! If true, use the total water column thickness rather than the
1285		! bathymetric depth for certain calculations.
1286		integer :: is, ie, js, je, nz
1287		integer :: i, j, k
1288
1289	0	if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_using_just_e: "// &
1290	0	"Module must be initialized before it is used.")
1291
1292	0	if (.not. CS%calculate_Eady_growth_rate) return
1293	0	if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function: "// &
1294	0	"%SN_u is not associated with use_variable_mixing.")
1295	0	if (.not. allocated(CS%SN_v)) call MOM_error(FATAL, "calc_slope_function: "// &
1296	0	"%SN_v is not associated with use_variable_mixing.")
1297
1298	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1299
1300	0	h_neglect = GV%H_subroundoff
1301	0	H_cutoff = real(2nz) (GV%Angstrom_H + h_neglect)
1302	0	dZ_cutoff = real(2nz) (GV%Angstrom_Z + GV%dz_subroundoff)
1303
1304	0	use_dztot = CS%full_depth_Eady_growth_rate ! .or. .not.(GV%Boussinesq or GV%semi_Boussinesq)
1305
1306	0	if (use_dztot) then
1307		!$OMP parallel do default(shared)
1308	0	do j=js-1,je+1 ; do i=is-1,ie+1
1309	0	dz_tot(i,j) = e(i,j,1) - e(i,j,nz+1)
1310		enddo ; enddo
1311		! The following mathematically equivalent expression is more expensive but is less
1312		! sensitive to roundoff for large Z_ref:
1313		! call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
1314		! do j=js-1,je+1
1315		! do i=is-1,ie+1 ; dz_tot(i,j) = 0.0 ; enddo
1316		! do k=1,nz ; do i=is-1,ie+1
1317		! dz_tot(i,j) = dz_tot(i,j) + dz(i,j,k)
1318		! enddo ; enddo
1319		! enddo
1320		endif
1321
1322		! To set the length scale based on the deformation radius, use wave_speed to
1323		! calculate the first-mode gravity wave speed and then blend the equatorial
1324		! and midlatitude deformation radii, using calc_resoln_function as a template.
1325
1326		!$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
1327	0	do k=nz,CS%VarMix_Ktop,-1
1328
1329		! Calculate the interface slopes E_x and E_y and u- and v- points respectively
1330	0	do j=js-1,je+1 ; do I=is-1,ie
1331	0	E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
1332		! Mask slopes where interface intersects topography
1333	0	if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0.
1334		enddo ; enddo
1335	0	do J=js-1,je ; do i=is-1,ie+1
1336	0	E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
1337		! Mask slopes where interface intersects topography
1338	0	if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
1339		enddo ; enddo
1340
1341		! Calculate NSh from this layer and add to the sum
1342	0	do j=js,je ; do I=is-1,ie
1343		S2 = ( E_x(I,j)*2 + 0.25( &
1344	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)) ) )
1345	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
1346
1347	0	Hdn = 2.h(i,j,k)h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1348	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)
1349	0	H_geom = sqrt(Hdn*Hup)
1350		! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1351	0	S2N2_u_local(I,j,k) = (H_geom * S2) * (GV%g_prime(k) / max(Hdn, Hup, CS%h_min_N2) )
1352		enddo ; enddo
1353	0	do J=js-1,je ; do i=is,ie
1354		S2 = ( E_y(i,J)*2 + 0.25( &
1355	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)) ) )
1356	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
1357
1358	0	Hdn = 2.h(i,j,k)h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
1359	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)
1360	0	H_geom = sqrt(Hdn*Hup)
1361		! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
1362	0	S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2)))
1363		enddo ; enddo
1364
1365		enddo ! k
1366
1367		!$OMP parallel do default(shared)
1368	0	do j=js,je
1369	0	do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo
1370	0	do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie
1371	0	CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k)
1372		enddo ; enddo
1373		! SN above contains S^2N^2H, convert to vertical average of S*N
1374
1375	0	if (use_dztot) then
1376	0	do I=is-1,ie
1377		CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / &
1378	0	max(dz_tot(i,j), dz_tot(i+1,j), GV%dz_subroundoff) )
1379		enddo
1380		else
1381	0	do I=is-1,ie
1382	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
1383	0	h2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)
1384	0	if ( min(h1, h2) > dZ_cutoff ) then
1385	0	CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / max(h1, h2) )
1386		else
1387	0	CS%SN_u(I,j) = 0.0
1388		endif
1389		enddo
1390		endif
1391		enddo
1392		!$OMP parallel do default(shared)
1393	0	do J=js-1,je
1394	0	do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo
1395	0	do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie
1396	0	CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k)
1397		enddo ; enddo
1398	0	if (use_dztot) then
1399	0	do i=is,ie
1400		CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / &
1401	0	max(dz_tot(i,j), dz_tot(i,j+1), GV%dz_subroundoff) )
1402		enddo
1403		else
1404	0	do i=is,ie
1405		! There is a primordial horizontal indexing bug on the following line from the previous
1406		! versions of the code. This comment should be deleted by the end of 2024.
1407		! if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
1408	0	h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
1409	0	h2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)
1410	0	if ( min(h1, h2) > dZ_cutoff ) then
1411	0	CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / max(h1, h2) )
1412		else
1413	0	CS%SN_v(i,J) = 0.0
1414		endif
1415		enddo
1416		endif
1417		enddo
1418
1419		end subroutine calc_slope_functions_using_just_e
1420
1421
1422		!> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity.
1423	0	subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC)
1424		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1425		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1426		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1427		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1428		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
1429		real, intent(in) :: dt !< Time increment [T ~> s]
1430		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1431		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1432		type(VarMix_CS), intent(in) :: CS !< Variable mixing control structure
1433		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure
1434		! Local variables
1435	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m]
1436		integer :: niblock, njblock, nkblock
1437		integer :: isoneutral_halo
1438
1439	0	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_QG_slopes: "//&
1440	0	"Module must be initialized before it is used.")
1441
1442	0	niblock = CS%niblock
1443	0	njblock = CS%njblock
1444	0	nkblock = CS%nkblock
1445	0	isoneutral_halo = 2
1446
1447		! +2 needed because loops in calc_isoneutral slopes run from is-1 to ie and js-1 to je
1448	0	if (niblock == 0) niblock = G%iec - G%isc + 2 + 2*isoneutral_halo
1449	0	if (njblock == 0) njblock = G%jec - G%jsc + 2 + 2*isoneutral_halo
1450	0	if (nkblock == 0) nkblock = GV%ke
1451
1452		!$omp target update to(h)
1453		!$omp target enter data map(alloc: e)
1454	0	call find_eta(h, tv, G, GV, US, e, halo_size=3)
1455		!$omp target enter data map(to: tv%T, tv%S)
1456		!$omp target enter data map(to: tv%SpV_avg) if (allocated(tv%SpV_avg))
1457		!$omp target enter data map(to: tv%p_surf) if (associated(tv%p_surf))
1458		call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, &
1459		slope_x, slope_y, niblock, njblock, nkblock, halo=isoneutral_halo, &
1460	0	OBC=OBC, OBC_N2=CS%OBC_friendly)
1461		!$omp target exit data map(release: tv%T, tv%S)
1462		!$omp target exit data map(delete: e)
1463		!$omp target exit data map(release: tv%SpV_avg) if (allocated(tv%SpV_avg))
1464		!$omp target exit data map(release: tv%p_surf) if (associated(tv%p_surf))
1465
1466	0	end subroutine calc_QG_slopes
1467
1468		!> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients
1469	0	subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, slope_x, slope_y, &
1470	0	vort_xy_dx, vort_xy_dy)
1471		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1472		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1473		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1474		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1475		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1476		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer vertical extents [Z ~> m]
1477		integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude
1478		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence
1479		!! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1480		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence
1481		!! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
1482		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim]
1483		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim]
1484		real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity
1485		!! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1486		real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity
1487		!! (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
1488		! Local variables
1489		real, dimension(SZI_(G),SZJB_(G)) :: &
1490	0	dslopey_dz, & ! z-derivative of y-slope at v-points [L-1 ~> m-1]
1491	0	h_at_v, & ! Thickness at v-points [H ~> m or kg m-2]
1492	0	beta_v, & ! Beta at v-points [T-1 L-1 ~> s-1 m-1]
1493	0	grad_vort_mag_v, & ! Magnitude of vorticity gradient at v-points [T-1 L-1 ~> s-1 m-1]
1494	0	grad_div_mag_v ! Magnitude of divergence gradient at v-points [T-1 L-1 ~> s-1 m-1]
1495
1496		real, dimension(SZIB_(G),SZJ_(G)) :: &
1497	0	dslopex_dz, & ! z-derivative of x-slope at u-points [L-1 ~> m-1]
1498	0	h_at_u, & ! Thickness at u-points [H ~> m or kg m-2]
1499	0	beta_u, & ! Beta at u-points [T-1 L-1 ~> s-1 m-1]
1500	0	grad_vort_mag_u, & ! Magnitude of vorticity gradient at u-points [T-1 L-1 ~> s-1 m-1]
1501	0	grad_div_mag_u ! Magnitude of divergence gradient at u-points [T-1 L-1 ~> s-1 m-1]
1502		real :: h_at_slope_above ! The thickness above [H ~> m or kg m-2]
1503		real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2]
1504		real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
1505		real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1]
1506		real :: Z_to_H ! A local copy of depth to thickness conversion factors or the inverse of the
1507		! mass-weighted average specific volumes around an interface [H Z-1 ~> nondim or kg m-3]
1508		real :: inv_PI3 ! The inverse of pi cubed [nondim]
1509		integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1510
1511	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1512	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1513	0	nz = GV%ke
1514
1515	0	inv_PI3 = 1.0 / ((4.0atan(1.0))*3)
1516	0	Z_to_H = GV%Z_to_H ! This will be replaced with a varying value in non-Boussinesq mode.
1517
1518	0	if ((k > 1) .and. (k < nz)) then
1519
1520	0	do j=js-2,je+2 ; do I=is-2,ie+1
1521		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) ) / &
1522		( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) &
1523	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 )
1524		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) ) / &
1525		( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
1526	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 )
1527	0	Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
1528	0	if (.not.GV%Boussinesq) &
1529		Z_to_H = ( (h(i,j,k-1) + h(i+1,j,k-1)) + (h(i,j,k) + h(i+1,j,k)) ) / &
1530	0	( (dz(i,j,k-1) + dz(i+1,j,k-1)) + (dz(i,j,k) + dz(i+1,j,k)) + GV%dZ_subroundoff)
1531	0	dslopex_dz(I,j) = 2. * ( slope_x(I,j,k) - slope_x(I,j,k+1) ) * (Z_to_H * Ih)
1532	0	h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
1533		enddo ; enddo
1534
1535	0	do J=js-2,je+1 ; do i=is-2,ie+2
1536		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) ) / &
1537		( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) &
1538	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 )
1539		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) ) / &
1540		( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
1541	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 )
1542	0	Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
1543	0	if (.not.GV%Boussinesq) &
1544		Z_to_H = ( (h(i,j,k-1) + h(i,j+1,k-1)) + (h(i,j,k) + h(i,j+1,k)) ) / &
1545	0	( (dz(i,j,k-1) + dz(i,j+1,k-1)) + (dz(i,j,k) + dz(i,j+1,k)) + GV%dZ_subroundoff)
1546	0	dslopey_dz(i,J) = 2. * ( slope_y(i,J,k) - slope_y(i,J,k+1) ) * (Z_to_H * Ih)
1547	0	h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
1548		enddo ; enddo
1549
1550	0	do J=js-2,je+1 ; do i=is-1,ie+1
1551	0	f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) )
1552		vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * &
1553		( ( (h_at_u(I,j) * dslopex_dz(I,j)) + (h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1)) ) &
1554		+ ( (h_at_u(I-1,j) * dslopex_dz(I-1,j)) + (h_at_u(I,j+1) * dslopex_dz(I,j+1)) ) ) / &
1555	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)
1556		enddo ; enddo
1557
1558	0	do j=js-1,je+1 ; do I=is-2,ie+1
1559	0	f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
1560		vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * &
1561		( ( (h_at_v(i,J) * dslopey_dz(i,J)) + (h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1)) ) &
1562		+ ( (h_at_v(i,J-1) * dslopey_dz(i,J-1)) + (h_at_v(i+1,J) * dslopey_dz(i+1,J)) ) ) / &
1563	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)
1564		enddo ; enddo
1565		endif ! k > 1
1566
1567	0	if (CS%use_QG_Leith_GM) then
1568
1569	0	do j=js,je ; do I=is-1,Ieq
1570		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)) &
1571	0	+ (vort_xy_dx(i+1,J) + vort_xy_dx(i,J-1))))**2)
1572		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)) &
1573	0	+ (div_xx_dy(i+1,J) + div_xx_dy(i,J-1))))**2)
1574	0	if (CS%use_beta_in_QG_Leith) then
1575		beta_u(I,j) = sqrt((0.5(G%dF_dx(i,j)+G%dF_dx(i+1,j))*2) + &
1576	0	(0.5(G%dF_dy(i,j)+G%dF_dy(i+1,j))*2))
1577		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)) &
1578	0	CS%Laplac3_const_u(I,j) * inv_PI3
1579		else
1580		CS%KH_u_QG(I,j,k) = (grad_vort_mag_u(I,j) + grad_div_mag_u(I,j)) * &
1581	0	CS%Laplac3_const_u(I,j) * inv_PI3
1582		endif
1583		enddo ; enddo
1584
1585	0	do J=js-1,Jeq ; do i=is,ie
1586		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)) &
1587	0	+ (vort_xy_dy(I,j+1) + vort_xy_dy(I-1,j))))**2)
1588		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)) &
1589	0	+ (div_xx_dx(I,j+1) + div_xx_dx(I-1,j))))**2)
1590	0	if (CS%use_beta_in_QG_Leith) then
1591		beta_v(i,J) = sqrt((0.5(G%dF_dx(i,j)+G%dF_dx(i,j+1))*2) + &
1592	0	(0.5(G%dF_dy(i,j)+G%dF_dy(i,j+1))*2))
1593		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)) &
1594	0	CS%Laplac3_const_v(i,J) * inv_PI3
1595		else
1596		CS%KH_v_QG(i,J,k) = (grad_vort_mag_v(i,J) + grad_div_mag_v(i,J)) * &
1597	0	CS%Laplac3_const_v(i,J) * inv_PI3
1598		endif
1599		enddo ; enddo
1600		! post diagnostics
1601
1602	0	if (k==nz) then
1603	0	if (CS%id_KH_v_QG > 0) call post_data(CS%id_KH_v_QG, CS%KH_v_QG, CS%diag)
1604	0	if (CS%id_KH_u_QG > 0) call post_data(CS%id_KH_u_QG, CS%KH_u_QG, CS%diag)
1605		endif
1606		endif
1607
1608	0	end subroutine calc_QG_Leith_viscosity
1609
1610		!> Initializes the variables mixing coefficients container
1611	1	subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
1612		type(time_type), intent(in) :: Time !< Current model time
1613		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1614		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1615		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1616		type(param_file_type), intent(in) :: param_file !< Parameter file handles
1617		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
1618		type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients
1619
1620		! Local variables
1621		real :: KhTr_Slope_Cff ! The nondimensional coefficient in the Visbeck formula
1622		! for the epipycnal tracer diffusivity [nondim]
1623		real :: KhTh_Slope_Cff ! The nondimensional coefficient in the Visbeck formula
1624		! for the interface depth diffusivity [nondim]
1625		real :: oneOrTwo ! A variable that may be 1 or 2, depending on which form
1626		! of the equatorial deformation radius us used [nondim]
1627		real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when
1628		! calculating the first-mode wave speed [H ~> m or kg m-2]
1629		real :: KhTr_passivity_coeff ! Coefficient setting the ratio between along-isopycnal tracer
1630		! mixing and interface height mixing [nondim]
1631		real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
1632		! default value is roughly (pi / (the age of the universe)).
1633		logical :: Gill_equatorial_Ld, use_FGNV_streamfn, use_MEKE, in_use
1634		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1635		integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1636		! for remapping. Values below 20190101 recover the remapping
1637		! answers from 2018, while higher values use more robust
1638		! forms of the same remapping expressions.
1639		real :: MLE_front_length ! The frontal-length scale used to calculate the upscaling of
1640		! buoyancy gradients in boundary layer parameterizations [L ~> m]
1641		real :: Leith_Lap_const ! The non-dimensional coefficient in the Leith viscosity [nondim]
1642		real :: grid_sp_u2, grid_sp_v2 ! Intermediate quantities for Leith metrics [L2 ~> m2]
1643		real :: grid_sp_u3, grid_sp_v3 ! Intermediate quantities for Leith metrics [L3 ~> m3]
1644		real :: wave_speed_min ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
1645		real :: wave_speed_tol ! The fractional tolerance for finding the wave speeds [nondim]
1646		logical :: Resoln_scaled_MEKE_visc ! If true, the viscosity contribution from MEKE is
1647		! scaled by the resolution function.
1648		logical :: better_speed_est ! If true, use a more robust estimate of the first
1649		! mode wave speed as the starting point for iterations.
1650		real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
1651		! temperature variance [nondim]
1652		logical :: use_SQG ! This is true if the SQG structure will be used for any parameterizations.
1653		logical :: om4_remap_via_sub_cells ! Use the OM4-era remap_via_sub_cells for calculating the EBT structure
1654		logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1655		! recreate the bugs, or if false bugs are only used if actively selected.
1656		logical :: mixing_coefs_OBC_bug ! If false, use only interior data for thickness weighting in
1657		! lateral mixing coefficient calculations and to calculate stratification
1658		! and other fields at open boundary condition faces.
1659		! This include declares and sets the variable "version".
1660		# include "version_variable.h"
1661		character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
1662		integer :: number_of_OBC_segments
1663		integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j
1664		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
1665		#ifdef __NVCOMPILER_OPENMP_GPU
1666		integer, parameter :: default_niblock = 0
1667		integer, parameter :: default_njblock = 0
1668		integer, parameter :: default_nkblock = 0
1669		#else
1670		integer, parameter :: default_niblock = 0
1671		integer, parameter :: default_njblock = 1
1672		integer, parameter :: default_nkblock = 1
1673		#endif
1674
1675	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1676	1	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1677	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1678	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1679
1680	1	CS%initialized = .true.
1681	1	in_use = .false. ! Set to true to avoid deallocating
1682	1	CS%diag => diag ! Diagnostics pointer
1683	1	CS%calculate_cg1 = .false.
1684	1	CS%calculate_Rd_dx = .false.
1685	1	CS%calculate_res_fns = .false.
1686	1	CS%use_simpler_Eady_growth_rate = .false.
1687	1	CS%full_depth_Eady_growth_rate = .false.
1688	1	CS%calculate_depth_fns = .false.
1689		! Read all relevant parameters and write them to the model log.
1690	1	call log_version(param_file, mdl, version, "")
1691		call get_param(param_file, mdl, "USE_VARIABLE_MIXING", CS%use_variable_mixing,&
1692		"If true, the variable mixing code will be called. This "//&
1693		"allows diagnostics to be created even if the scheme is "//&
1694		"not used. If KHTR_SLOPE_CFF>0 or KhTh_Slope_Cff>0, "//&
1695		"this is set to true regardless of what is in the "//&
1696	1	"parameter file.", default=.false.)
1697		call get_param(param_file, mdl, "USE_VISBECK", CS%use_Visbeck,&
1698		"If true, use the Visbeck et al. (1997) formulation for \n"//&
1699	1	"thickness diffusivity.", default=.false.)
1700		call get_param(param_file, mdl, "RESOLN_SCALED_KH", CS%Resoln_scaled_Kh, &
1701		"If true, the Laplacian lateral viscosity is scaled away "//&
1702		"when the first baroclinic deformation radius is well "//&
1703	1	"resolved.", default=.false.)
1704		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH", CS%Depth_scaled_KhTh, &
1705		"If true, KHTH is scaled away when the depth is shallower "//&
1706		"than a reference depth: KHTH = MIN(1,H/H0)*N KHTH, "//&
1707		"where H0 is a reference depth, controlled via DEPTH_SCALED_KHTH_H0, "//&
1708		"and the exponent (N) is controlled via DEPTH_SCALED_KHTH_EXP.",&
1709	1	default=.false.)
1710		call get_param(param_file, mdl, "RESOLN_SCALED_KHTH", CS%Resoln_scaled_KhTh, &
1711		"If true, the interface depth diffusivity is scaled away "//&
1712		"when the first baroclinic deformation radius is well "//&
1713	1	"resolved.", default=.false.)
1714		call get_param(param_file, mdl, "RESOLN_SCALED_KHTR", CS%Resoln_scaled_KhTr, &
1715		"If true, the epipycnal tracer diffusivity is scaled "//&
1716		"away when the first baroclinic deformation radius is "//&
1717	1	"well resolved.", default=.false.)
1718		call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
1719	1	default=.false., do_not_log=.true.)
1720		call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", Resoln_scaled_MEKE_visc, &
1721		"If true, the viscosity contribution from MEKE is scaled by "//&
1722	1	"the resolution function.", default=.false., do_not_log=.true.) ! Logged elsewhere.
1723	1	if (.not.use_MEKE) Resoln_scaled_MEKE_visc = .false.
1724		call get_param(param_file, mdl, "RESOLN_USE_EBT", CS%Resoln_use_ebt, &
1725		"If true, uses the equivalent barotropic wave speed instead "//&
1726		"of first baroclinic wave for calculating the resolution function.",&
1727	1	default=.false.)
1728		call get_param(param_file, mdl, "BACKSCAT_EBT_POWER", CS%BS_EBT_power, &
1729		"Power to raise EBT vertical structure to when backscatter "// &
1730	1	"has vertical structure.", units="nondim", default=0.0)
1731		call get_param(param_file, mdl, "BS_USE_SQG_STRUCT", CS%BS_use_sqg_struct, &
1732		"If true, the SQG vertical structure is used for backscatter "//&
1733		"on the condition that BS_EBT_power=0", &
1734	1	default=.false.)
1735		call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", CS%khth_use_ebt_struct, &
1736		"If true, uses the equivalent barotropic structure "//&
1737		"as the vertical structure of thickness diffusivity.",&
1738	1	default=.false.)
1739		call get_param(param_file, mdl, "KHTH_USE_SQG_STRUCT", CS%khth_use_sqg_struct, &
1740		"If true, uses the surface quasigeostrophic structure "//&
1741		"as the vertical structure of thickness diffusivity.",&
1742	1	default=.false.)
1743		call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%khtr_use_ebt_struct, &
1744		"If true, uses the equivalent barotropic structure "//&
1745		"as the vertical structure of tracer diffusivity.",&
1746	1	default=.false.)
1747		call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", CS%khtr_use_sqg_struct, &
1748		"If true, uses the surface quasigeostrophic structure "//&
1749		"as the vertical structure of tracer diffusivity.",&
1750	1	default=.false.)
1751		call get_param(param_file, mdl, "KD_GL90_USE_EBT_STRUCT", CS%kdgl90_use_ebt_struct, &
1752		"If true, uses the equivalent barotropic structure "//&
1753		"as the vertical structure of diffusivity in the GL90 scheme.",&
1754	1	default=.false.)
1755		call get_param(param_file, mdl, "KD_GL90_USE_SQG_STRUCT", CS%kdgl90_use_sqg_struct, &
1756		"If true, uses the equivalent barotropic structure "//&
1757		"as the vertical structure of diffusivity in the GL90 scheme.",&
1758	1	default=.false.)
1759		call get_param(param_file, mdl, "KHTH_SLOPE_CFF", KhTh_Slope_Cff, &
1760		"The nondimensional coefficient in the Visbeck formula "//&
1761	1	"for the interface depth diffusivity", units="nondim", default=0.0)
1762		call get_param(param_file, mdl, "KHTR_SLOPE_CFF", KhTr_Slope_Cff, &
1763		"The nondimensional coefficient in the Visbeck formula "//&
1764	1	"for the epipycnal tracer diffusivity", units="nondim", default=0.0)
1765		call get_param(param_file, mdl, "USE_STORED_SLOPES", CS%use_stored_slopes,&
1766		"If true, the isopycnal slopes are calculated once and "//&
1767		"stored for re-use. This uses more memory but avoids calling "//&
1768		"the equation of state more times than should be necessary.", &
1769	1	default=.false.)
1770		call get_param(param_file, mdl, "VERY_SMALL_FREQUENCY", absurdly_small_freq, &
1771		"A miniscule frequency that is used to avoid division by 0. The default "//&
1772		"value is roughly (pi / (the age of the universe)).", &
1773	1	default=1.0e-17, units="s-1", scale=US%T_to_s)
1774		call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", use_FGNV_streamfn, &
1775	1	default=.false., do_not_log=.true.)
1776		CS%calculate_cg1 = CS%calculate_cg1 .or. use_FGNV_streamfn .or. CS%khth_use_ebt_struct &
1777	1	.or. CS%kdgl90_use_ebt_struct .or. CS%BS_EBT_power>0.
1778	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. use_MEKE
1779		! Indicate whether to calculate the Eady growth rate
1780	1	CS%calculate_Eady_growth_rate = use_MEKE .or. (KhTr_Slope_Cff>0.) .or. (KhTh_Slope_Cff>0.)
1781		call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", KhTr_passivity_coeff, &
1782	1	units="nondim", default=0., do_not_log=.true.)
1783	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (KhTr_passivity_coeff>0.)
1784		call get_param(param_file, mdl, "MLE_FRONT_LENGTH", MLE_front_length, &
1785	1	units="m", default=0.0, scale=US%m_to_L, do_not_log=.true.)
1786	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (MLE_front_length>0.)
1787
1788	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
1789
1790		! Isopycnal blocking parameters
1791		call get_param(param_file, mdl, "ISOPYCNAL_NIBLOCK", CS%niblock, &
1792		"The i-direction block size used to calculate isopycnal slopes. "//&
1793		"If 0, or when running with OpenMP offload, "//&
1794		"the full computational domain width is used. "//&
1795		"If USE_STANLEY_ISO is true, ISOPYCNAL_NIBLOCK cannot equal 1.", &
1796	1	default=default_niblock, layoutParam=.true.)
1797		call get_param(param_file, mdl, "ISOPYCNAL_NJBLOCK", CS%njblock, &
1798		"The j-direction block size used to calculate isopycnal slopes. "//&
1799		"If 0, defaults to 1, except when running with OpenMP offload, "//&
1800		"in which case the full computational domain height is used. " //&
1801		"If USE_STANLEY_ISO is true, ISOPYCNAL_NJBLOCK cannot equal 1.", &
1802	1	default=default_njblock, layoutParam=.true.)
1803		call get_param(param_file, mdl, "ISOPYCNAL_NKBLOCK", CS%nkblock, &
1804		"The k-direction block size used to calculate isopycnal slopes. "//&
1805		"If 0, defaults to 1, except when "//&
1806		"running with OpenMP offload, in which case the full computational "//&
1807	1	"domain depth is used.", default=default_nkblock, layoutParam=.true.)
1808	1	if (CS%niblock < 0) &
1809		call MOM_error(FATAL, "ISOPYCNAL_NIBLOCK must be nonnegative; "//&
1810	0	"use 0 to select the default block size.")
1811	1	if (CS%njblock < 0) &
1812		call MOM_error(FATAL, "ISOPYCNAL_NJBLOCK must be nonnegative; "//&
1813	0	"use 0 to select the default block size.")
1814	1	if (CS%nkblock < 0) &
1815		call MOM_error(FATAL, "ISOPYCNAL_NKBLOCK must be nonnegative; "//&
1816	0	"use 0 to select the default block size.")
1817
1818		call get_param(param_file, mdl, "USE_STANLEY_ISO", CS%use_stanley_iso, &
1819		"If true, turn on Stanley SGS T variance parameterization "// &
1820	1	"in isopycnal slope code.", default=.false.)
1821	1	if (CS%use_stanley_iso) then
1822		call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
1823		"Coefficient correlating the temperature gradient and SGS T variance.", &
1824	0	units="nondim", default=-1.0, do_not_log=.true.)
1825	0	if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
1826	0	"STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.")
1827	0	if (CS%njblock == 1) then
1828		call MOM_error(WARNING, "ISOPYCNAL_NJBLOCK must be >= 2 or 0 if USE_STANLEY_ISO is true."//&
1829	0	" Changing block size from 1 to 2 for this run.")
1830	0	CS%njblock = 2
1831		endif
1832	0	if (CS%niblock == 1) then
1833		call MOM_error(WARNING, "ISOPYCNAL_NIBLOCK must be >= 2 or 0 if USE_STANLEY_ISO is true."//&
1834	0	" Changing block size from 1 to 2 for this run.")
1835	0	CS%niblock = 2
1836		endif
1837		endif
1838		call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_OBC_segments, &
1839	1	default=0, do_not_log=.true.)
1840		call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1841	1	default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1842		call get_param(param_file, mdl, "MIXING_COEFS_OBC_BUG", mixing_coefs_OBC_bug, &
1843		"If false, use only interior data for thickness weighting in lateral mixing "//&
1844		"coefficient calculations and to calculate stratification and other fields at "//&
1845		"open boundary condition faces.", &
1846	1	default=enable_bugs, do_not_log=(number_of_OBC_segments<=0))
1847	1	CS%OBC_friendly = .not. MIXING_COEFS_OBC_BUG
1848		call get_param(param_file, mdl, "RESOLN_FUNCTION_OBC_BUG", CS%res_fn_OBC_bug, &
1849		"If false, use only interior data for calculating the resolution functions at "//&
1850		"open boundary condition faces and vertices.", &
1851	1	default=enable_bugs, do_not_log=(number_of_OBC_segments<=0))
1852
1853		if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct &
1854	1	.or. CS%BS_EBT_power>0. .or. CS%khtr_use_ebt_struct) then
1855	0	in_use = .true.
1856		call get_param(param_file, mdl, "RESOLN_N2_FILTER_DEPTH", N2_filter_depth, &
1857		"The depth below which N2 is monotonized to avoid stratification "//&
1858		"artifacts from altering the equivalent barotropic mode structure. "//&
1859		"This monotonzization is disabled if this parameter is negative.", &
1860	0	units="m", default=-1.0, scale=GV%m_to_H)
1861	0	allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1862		endif
1863
1864		use_SQG = CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct .or. &
1865	1	CS%kdgl90_use_sqg_struct
1866		call get_param(param_file, mdl, "SQG_EXPO", CS%sqg_expo, &
1867		"Nondimensional exponent coeffecient of the SQG mode that is used for the "//&
1868		"vertical struture of diffusivities.", &
1869	1	units="nondim", default=1.0, do_not_log=.not.use_SQG)
1870		call get_param(param_file, mdl, "INTERPOLATED_SQG_STRUCTURE", CS%interpolated_sqg_struct, &
1871		"If true, interpolate properties to velocity points and then interpolate the "//&
1872		"buoyancy frequencies and layer thicknesses back to tracer points when "//&
1873		"calculating the SQG vertical structure.", &
1874	1	default=.true., do_not_log=.not.use_SQG)
1875		!### Consider changing the default for INTERPOLATED_SQG_STRUCTURE to false.
1876
1877	1	if ((CS%BS_EBT_power>0.) .and. CS%BS_use_sqg_struct) call MOM_error(FATAL, &
1878	0	"calc_resoln_function: BS_EBT_POWER>0. and BS_USE_SQG=True cannot be set together")
1879
1880	1	if (CS%khth_use_ebt_struct .and. CS%khth_use_sqg_struct) call MOM_error(FATAL, &
1881	0	"calc_resoln_function: Only one of KHTH_USE_EBT_STRUCT and KHTH_USE_SQG_STRUCT can be true")
1882
1883	1	if (CS%khtr_use_ebt_struct .and. CS%khtr_use_sqg_struct) call MOM_error(FATAL, &
1884	0	"calc_resoln_function: Only one of KHTR_USE_EBT_STRUCT and KHTR_USE_SQG_STRUCT can be true")
1885
1886	1	if (CS%kdgl90_use_ebt_struct .and. CS%kdgl90_use_sqg_struct) call MOM_error(FATAL, &
1887	0	"calc_resoln_function: Only one of KD_GL90_USE_EBT_STRUCT and KD_GL90_USE_SQG_STRUCT can be true")
1888
1889	1	if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
1890	0	allocate(CS%BS_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
1891		endif
1892
1893	1	if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
1894	0	allocate(CS%khth_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1895		endif
1896
1897	1	if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then
1898	0	allocate(CS%khtr_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1899		endif
1900
1901	1	if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then
1902	0	allocate(CS%kdgl90_struct(isd:ied, jsd:jed, gv%ke), source=0.0)
1903		endif
1904
1905	1	if (CS%use_stored_slopes) then
1906	1	if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then
1907		call get_param(param_file, mdl, "VISBECK_MAX_SLOPE", CS%Visbeck_S_max, &
1908		"If non-zero, is an upper bound on slopes used in the "//&
1909		"Visbeck formula for diffusivity. This does not affect the "//&
1910		"isopycnal slope calculation used within thickness diffusion.", &
1911	0	units="nondim", default=0.0, scale=US%L_to_Z)
1912		else
1913	1	CS%Visbeck_S_max = 0.
1914		endif
1915		endif
1916
1917	1	if (CS%use_stored_slopes .or. (CS%interpolated_sqg_struct .and. (CS%sqg_expo>0.0))) then
1918		! CS%calculate_Eady_growth_rate=.true.
1919	1	in_use = .true.
1920	671917	allocate(CS%slope_x(IsdB:IedB,jsd:jed,GV%ke+1), source=0.0)
1921	676553	allocate(CS%slope_y(isd:ied,JsdB:JedB,GV%ke+1), source=0.0)
1922		call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, &
1923		"A diapycnal diffusivity that is used to interpolate "//&
1924		"more sensible values of T & S into thin layers.", &
1925	1	units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T)
1926		endif
1927
1928	1	if (CS%calculate_Eady_growth_rate) then
1929	1	in_use = .true.
1930	8841	allocate(CS%SN_u(IsdB:IedB,jsd:jed), source=0.0)
1931	8902	allocate(CS%SN_v(isd:ied,JsdB:JedB), source=0.0)
1932		CS%id_SN_u = register_diag_field('ocean_model', 'SN_u', diag%axesCu1, Time, &
1933	1	'Inverse eddy time-scale, S*N, at u-points', 's-1', conversion=US%s_to_T)
1934		CS%id_SN_v = register_diag_field('ocean_model', 'SN_v', diag%axesCv1, Time, &
1935	1	'Inverse eddy time-scale, S*N, at v-points', 's-1', conversion=US%s_to_T)
1936		call get_param(param_file, mdl, "USE_SIMPLER_EADY_GROWTH_RATE", CS%use_simpler_Eady_growth_rate, &
1937		"If true, use a simpler method to calculate the Eady growth rate "//&
1938	1	"that avoids division by layer thickness. Recommended.", default=.false.)
1939	1	if (CS%use_simpler_Eady_growth_rate) then
1940	0	if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
1941		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
1942	0	"When USE_SIMPLER_EADY_GROWTH_RATE=True, USE_STORED_SLOPES must also be True.")
1943		call get_param(param_file, mdl, "EADY_GROWTH_RATE_D_SCALE", CS%Eady_GR_D_scale, &
1944		"The depth from surface over which to average SN when calculating "//&
1945		"a 2D Eady growth rate. Zero mean use full depth.", &
1946	0	units="m", default=0., scale=US%m_to_Z)
1947		call get_param(param_file, mdl, "EADY_GROWTH_RATE_CROPPING_DISTANCE", CS%cropping_distance, &
1948		"Distance from surface or bottom to filter out outcropped or "//&
1949		"incropped interfaces for the Eady growth rate calc. "//&
1950	0	"Negative values disables cropping.", units="m", default=0., scale=US%m_to_Z)
1951		else
1952		call get_param(param_file, mdl, "VARMIX_KTOP", CS%VarMix_Ktop, &
1953		"The layer number at which to start vertical integration "//&
1954		"of S*N for purposes of finding the Eady growth rate.", &
1955	1	units="nondim", default=2)
1956		call get_param(param_file, mdl, "MIN_DZ_FOR_SLOPE_N2", CS%h_min_N2, &
1957		"The minimum vertical distance to use in the denominator of the "//&
1958		"bouyancy frequency used in the slope calculation.", &
1959	1	units="m", default=1.0, scale=GV%m_to_H, do_not_log=CS%use_stored_slopes)
1960
1961		call get_param(param_file, mdl, "FULL_DEPTH_EADY_GROWTH_RATE", CS%full_depth_Eady_growth_rate, &
1962		"If true, calculate the Eady growth rate based on average slope times "//&
1963		"stratification that includes contributions from sea-level changes "//&
1964		"in its denominator, rather than just the nominal depth of the bathymetry. "//&
1965		"This only applies when using the model interface heights as a proxy for "//&
1966		"isopycnal slopes.", default=.not.(GV%Boussinesq.or.GV%semi_Boussinesq), &
1967	1	do_not_log=CS%use_stored_slopes)
1968		endif
1969		endif
1970
1971	1	if (KhTr_Slope_Cff>0. .or. KhTh_Slope_Cff>0.) then
1972	0	in_use = .true.
1973		call get_param(param_file, mdl, "VISBECK_L_SCALE", CS%Visbeck_L_scale, &
1974		"The fixed length scale in the Visbeck formula, or if negative a nondimensional "//&
1975		"scaling factor relating this length scale squared to the cell areas.", &
1976	0	units="m or nondim", default=0.0, scale=US%m_to_L)
1977	0	allocate(CS%L2u(IsdB:IedB,jsd:jed), source=0.0)
1978	0	allocate(CS%L2v(isd:ied,JsdB:JedB), source=0.0)
1979	0	if (CS%Visbeck_L_scale<0) then
1980		! Undo the rescaling of CS%Visbeck_L_scale.
1981	0	do j=js,je ; do I=is-1,Ieq
1982	0	CS%L2u(I,j) = (US%L_to_mCS%Visbeck_L_scale)2 G%areaCu(I,j)
1983		enddo ; enddo
1984	0	do J=js-1,Jeq ; do i=is,ie
1985	0	CS%L2v(i,J) = (US%L_to_mCS%Visbeck_L_scale)2 G%areaCv(i,J)
1986		enddo ; enddo
1987		else
1988	0	CS%L2u(:,:) = CS%Visbeck_L_scale**2
1989	0	CS%L2v(:,:) = CS%Visbeck_L_scale**2
1990		endif
1991
1992		CS%id_L2u = register_diag_field('ocean_model', 'L2u', diag%axesCu1, Time, &
1993		'Length scale squared for mixing coefficient, at u-points', &
1994	0	'm2', conversion=US%L_to_m**2)
1995		CS%id_L2v = register_diag_field('ocean_model', 'L2v', diag%axesCv1, Time, &
1996		'Length scale squared for mixing coefficient, at v-points', &
1997	0	'm2', conversion=US%L_to_m**2)
1998		endif
1999
2000		CS%id_sqg_struct = register_diag_field('ocean_model', 'sqg_struct', diag%axesTl, Time, &
2001	1	'Vertical structure of SQG mode', 'nondim')
2002		if (CS%BS_use_sqg_struct .or. CS%khth_use_sqg_struct .or. CS%khtr_use_sqg_struct &
2003	1	.or. CS%kdgl90_use_sqg_struct .or. CS%id_sqg_struct>0) then
2004	0	allocate(CS%sqg_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
2005		endif
2006
2007	1	if (CS%BS_EBT_power>0. .or. CS%BS_use_sqg_struct) then
2008		CS%id_BS_struct = register_diag_field('ocean_model', 'BS_struct', diag%axesTl, Time, &
2009	0	'Vertical structure of backscatter', 'nondim')
2010		endif
2011	1	if (CS%khth_use_ebt_struct .or. CS%khth_use_sqg_struct) then
2012		CS%id_khth_struct = register_diag_field('ocean_model', 'khth_struct', diag%axesTl, Time, &
2013	0	'Vertical structure of thickness diffusivity', 'nondim')
2014		endif
2015	1	if (CS%khtr_use_ebt_struct .or. CS%khtr_use_sqg_struct) then
2016		CS%id_khtr_struct = register_diag_field('ocean_model', 'khtr_struct', diag%axesTl, Time, &
2017	0	'Vertical structure of tracer diffusivity', 'nondim')
2018		endif
2019	1	if (CS%kdgl90_use_ebt_struct .or. CS%kdgl90_use_sqg_struct) then
2020		CS%id_kdgl90_struct = register_diag_field('ocean_model', 'kdgl90_struct', diag%axesTl, Time, &
2021	0	'Vertical structure of GL90 diffusivity', 'nondim')
2022		endif
2023
2024	1	if ((CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) ) then
2025		CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, &
2026		'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', &
2027	1	's-2', conversion=(US%L_to_ZUS%s_to_T)*2)
2028		CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, &
2029		'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', &
2030	1	's-2', conversion=(US%L_to_ZUS%s_to_T)*2)
2031		endif
2032	1	if (CS%use_simpler_Eady_growth_rate) then
2033		CS%id_dzu = register_diag_field('ocean_model', 'dzu_Visbeck', diag%axesCui, Time, &
2034		'dz at u-points, used in calculating Eady growth rate in Visbeck et al..', &
2035	0	'm', conversion=US%Z_to_m)
2036		CS%id_dzv = register_diag_field('ocean_model', 'dzv_Visbeck', diag%axesCvi, Time, &
2037		'dz at v-points, used in calculating Eady growth rate in Visbeck et al..', &
2038	0	'm', conversion=US%Z_to_m)
2039		CS%id_dzSxN = register_diag_field('ocean_model', 'dzSxN', diag%axesCui, Time, &
2040		'dz * \|slope_x\| * N, used in calculating Eady growth rate in '//&
2041	0	'Visbeck et al..', 'm s-1', conversion=US%Z_to_m*US%s_to_T)
2042		CS%id_dzSyN = register_diag_field('ocean_model', 'dzSyN', diag%axesCvi, Time, &
2043		'dz * \|slope_y\| * N, used in calculating Eady growth rate in '//&
2044	0	'Visbeck et al..', 'm s-1', conversion=US%Z_to_m*US%s_to_T)
2045		endif
2046	1	if (CS%use_stored_slopes) then
2047		CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, &
2048		'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', &
2049	1	'nondim', conversion=US%Z_to_L**2)
2050		CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, &
2051		'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', &
2052	1	'nondim', conversion=US%Z_to_L**2)
2053		endif
2054
2055	1	oneOrTwo = 1.0
2056		CS%Resoln_scaling_used = CS%Resoln_scaled_Kh .or. CS%Resoln_scaled_KhTh .or. &
2057	1	CS%Resoln_scaled_KhTr .or. Resoln_scaled_MEKE_visc
2058	1	if (CS%Resoln_scaling_used) then
2059	1	CS%calculate_Rd_dx = .true.
2060	1	CS%calculate_res_fns = .true.
2061	8773	allocate(CS%Res_fn_h(isd:ied,jsd:jed), source=0.0)
2062	8971	allocate(CS%Res_fn_q(IsdB:IedB,JsdB:JedB), source=0.0)
2063	8841	allocate(CS%Res_fn_u(IsdB:IedB,jsd:jed), source=0.0)
2064	8902	allocate(CS%Res_fn_v(isd:ied,JsdB:JedB), source=0.0)
2065	8971	allocate(CS%beta_dx2_q(IsdB:IedB,JsdB:JedB), source=0.0)
2066	8841	allocate(CS%beta_dx2_u(IsdB:IedB,jsd:jed), source=0.0)
2067	8902	allocate(CS%beta_dx2_v(isd:ied,JsdB:JedB), source=0.0)
2068	8971	allocate(CS%f2_dx2_q(IsdB:IedB,JsdB:JedB), source=0.0)
2069	8841	allocate(CS%f2_dx2_u(IsdB:IedB,jsd:jed), source=0.0)
2070	8902	allocate(CS%f2_dx2_v(isd:ied,JsdB:JedB), source=0.0)
2071
2072		CS%id_Res_fn = register_diag_field('ocean_model', 'Res_fn', diag%axesT1, Time, &
2073	1	'Resolution function for scaling diffusivities', 'nondim')
2074
2075		call get_param(param_file, mdl, "KH_RES_SCALE_COEF", CS%Res_coef_khth, &
2076		"A coefficient that determines how KhTh 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	1	units="nondim", default=1.0)
2080		call get_param(param_file, mdl, "KH_RES_FN_POWER", CS%Res_fn_power_khth, &
2081		"The power of dx/Ld in the Kh resolution function. Any "//&
2082		"positive integer may be used, although even integers "//&
2083		"are more efficient to calculate. Setting this greater "//&
2084		"than 100 results in a step-function being used.", &
2085	1	default=2)
2086		call get_param(param_file, mdl, "VISC_RES_SCALE_COEF", CS%Res_coef_visc, &
2087		"A coefficient that determines how Kh is scaled away if "//&
2088		"RESOLN_SCALED_... is true, as "//&
2089		"F = 1 / (1 + (KH_RES_SCALE_COEF*Rd/dx)^KH_RES_FN_POWER). "//&
2090		"This function affects lateral viscosity, Kh, and not KhTh.", &
2091	1	units="nondim", default=CS%Res_coef_khth)
2092		call get_param(param_file, mdl, "VISC_RES_FN_POWER", CS%Res_fn_power_visc, &
2093		"The power of dx/Ld in the Kh resolution function. Any "//&
2094		"positive integer may be used, although even integers "//&
2095		"are more efficient to calculate. Setting this greater "//&
2096		"than 100 results in a step-function being used. "//&
2097		"This function affects lateral viscosity, Kh, and not KhTh.", &
2098	1	default=CS%Res_fn_power_khth)
2099		call get_param(param_file, mdl, "INTERPOLATE_RES_FN", CS%interpolate_Res_fn, &
2100		"If true, interpolate the resolution function to the "//&
2101		"velocity points from the thickness points; otherwise "//&
2102		"interpolate the wave speed and calculate the resolution "//&
2103	1	"function independently at each point.", default=.false.)
2104	1	if (CS%interpolate_Res_fn) then
2105	0	if (CS%Res_coef_visc /= CS%Res_coef_khth) call MOM_error(FATAL, &
2106		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2107	0	"When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_SCALE_COEF.")
2108	0	if (CS%Res_fn_power_visc /= CS%Res_fn_power_khth) call MOM_error(FATAL, &
2109		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2110	0	"When INTERPOLATE_RES_FN=True, VISC_RES_FN_POWER must equal KH_RES_FN_POWER.")
2111		endif
2112		call get_param(param_file, mdl, "GILL_EQUATORIAL_LD", Gill_equatorial_Ld, &
2113		"If true, uses Gill's definition of the baroclinic "//&
2114		"equatorial deformation radius, otherwise, if false, use "//&
2115		"Pedlosky's definition. These definitions differ by a factor "//&
2116		"of 2 in front of the beta term in the denominator. Gill's "//&
2117	1	"is the more appropriate definition.", default=.true.)
2118	1	if (Gill_equatorial_Ld) then
2119	1	oneOrTwo = 2.0
2120		endif
2121
2122	7443	do J=js-1,Jeq ; do I=is-1,Ieq
2123		CS%f2_dx2_q(I,J) = ((G%dxBu(I,J)2) + (G%dyBu(I,J)2)) * &
2124	7381	max(G%Coriolis2Bu(I,J), absurdly_small_freq**2)
2125		CS%beta_dx2_q(I,J) = oneOrTwo * ((G%dxBu(I,J)2) + (G%dyBu(I,J)2)) * (sqrt(0.5 * &
2126		( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
2127		(((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2)) + &
2128		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
2129	7442	(((G%CoriolisBu(I,J+1)-G%CoriolisBu(I,J)) * G%IdyCu(I,j+1))**2)) ) ))
2130		enddo ; enddo
2131
2132	7321	do j=js,je ; do I=is-1,Ieq
2133		CS%f2_dx2_u(I,j) = ((G%dxCu(I,j)2) + (G%dyCu(I,j)2)) * &
2134	7260	max(0.5* (G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I,J-1)), absurdly_small_freq**2)
2135		CS%beta_dx2_u(I,j) = oneOrTwo * ((G%dxCu(I,j)2) + (G%dyCu(I,j)2)) * (sqrt( &
2136		((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2 + &
2137		0.25( ((((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) G%IdxCv(i,J-1))**2) + &
2138		(((G%CoriolisBu(I+1,J)-G%CoriolisBu(I,J)) * G%IdxCv(i+1,J))**2)) + &
2139		((((G%CoriolisBu(I+1,J-1)-G%CoriolisBu(I,J-1)) * G%IdxCv(i+1,J-1))**2) + &
2140	7320	(((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2)) ) ))
2141		enddo ; enddo
2142
2143	7382	do J=js-1,Jeq ; do i=is,ie
2144		CS%f2_dx2_v(i,J) = ((G%dxCv(i,J)2) + (G%dyCv(i,J)2)) * &
2145	7320	max(0.5(G%Coriolis2Bu(I,J)+G%Coriolis2Bu(I-1,J)), absurdly_small_freq*2)
2146		CS%beta_dx2_v(i,J) = oneOrTwo * ((G%dxCv(i,J)2) + (G%dyCv(i,J)2)) * (sqrt( &
2147		((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2 + &
2148		0.25( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) G%IdyCu(I,j))**2) + &
2149		(((G%CoriolisBu(I-1,J+1)-G%CoriolisBu(I-1,J)) * G%IdyCu(I-1,j+1))**2)) + &
2150		((((G%CoriolisBu(I,J+1)-G%CoriolisBu(I,J)) * G%IdyCu(I,j+1))**2) + &
2151	7381	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ) ))
2152		enddo ; enddo
2153
2154		endif
2155
2156	1	if (CS%Depth_scaled_KhTh) then
2157	0	CS%calculate_depth_fns = .true.
2158	0	allocate(CS%Depth_fn_u(IsdB:IedB,jsd:jed), source=0.0)
2159	0	allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB), source=0.0)
2160		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, &
2161		"The depth above which KHTH is scaled away.", &
2162	0	units="m", scale=US%m_to_Z, default=1000.)
2163		call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, &
2164		"The exponent used in the depth dependent scaling function for KHTH.", &
2165	0	units="nondim", default=3.0)
2166		endif
2167
2168		! Resolution %Rd_dx_h
2169		CS%id_Rd_dx = register_diag_field('ocean_model', 'Rd_dx', diag%axesT1, Time, &
2170	1	'Ratio between deformation radius and grid spacing', 'm m-1')
2171	1	CS%calculate_Rd_dx = CS%calculate_Rd_dx .or. (CS%id_Rd_dx>0)
2172
2173	1	if (CS%calculate_Rd_dx) then
2174	1	CS%calculate_cg1 = .true. ! We will need %cg1
2175	8773	allocate(CS%Rd_dx_h(isd:ied,jsd:jed), source=0.0)
2176	8773	allocate(CS%beta_dx2_h(isd:ied,jsd:jed), source=0.0)
2177	8773	allocate(CS%f2_dx2_h(isd:ied,jsd:jed), source=0.0)
2178	7627	do j=js-1,je+1 ; do i=is-1,ie+1
2179		CS%f2_dx2_h(i,j) = ((G%dxT(i,j)2) + (G%dyT(i,j)2)) * &
2180		max(0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
2181		(G%Coriolis2Bu(I-1,J) + G%Coriolis2Bu(I,J-1))), &
2182	7564	absurdly_small_freq**2)
2183		CS%beta_dx2_h(i,j) = oneOrTwo * ((G%dxT(i,j)2) + (G%dyT(i,j)2)) * (sqrt(0.5 * &
2184		( ((((G%CoriolisBu(I,J)-G%CoriolisBu(I-1,J)) * G%IdxCv(i,J))**2) + &
2185		(((G%CoriolisBu(I,J-1)-G%CoriolisBu(I-1,J-1)) * G%IdxCv(i,J-1))**2)) + &
2186		((((G%CoriolisBu(I,J)-G%CoriolisBu(I,J-1)) * G%IdyCu(I,j))**2) + &
2187	7626	(((G%CoriolisBu(I-1,J)-G%CoriolisBu(I-1,J-1)) * G%IdyCu(I-1,j))**2)) ) ))
2188		enddo ; enddo
2189		endif
2190
2191	1	if (CS%calculate_cg1) then
2192	1	in_use = .true.
2193	8773	allocate(CS%cg1(isd:ied,jsd:jed), source=0.0)
2194		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2195		"This sets the default value for the various _ANSWER_DATE parameters.", &
2196	1	default=99991231)
2197		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2198		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
2199		"Values below 20190101 result in the use of older, less accurate expressions "//&
2200		"that were in use at the end of 2018. Higher values result in the use of more "//&
2201		"robust and accurate forms of mathematically equivalent expressions.", &
2202	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
2203	1	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2204
2205		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
2206		"The fractional tolerance for finding the wave speeds.", &
2207	1	units="nondim", default=0.001)
2208		!### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
2209		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
2210		"A floor in the first mode speed below which 0 used instead.", &
2211	1	units="m s-1", default=0.0, scale=US%m_s_to_L_T)
2212		call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
2213		"If true, use a more robust estimate of the first mode wave speed as the "//&
2214	1	"starting point for iterations.", default=.true.)
2215		call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2216	1	do_not_log=.true., default=.true.)
2217		call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2218		"If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
2219		"See REMAPPING_USE_OM4_SUBCELLS for details. "//&
2220	1	"We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2221		call wave_speed_init(CS%wave_speed, GV, use_ebt_mode=CS%Resoln_use_ebt, &
2222		mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, &
2223		better_speed_est=better_speed_est, min_speed=wave_speed_min, &
2224	1	om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol)
2225		endif
2226
2227		! Leith parameters
2228		call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%use_QG_Leith_GM, &
2229		"If true, use the QG Leith viscosity as the GM coefficient.", &
2230	1	default=.false.)
2231
2232	1	if (CS%Use_QG_Leith_GM) then
2233		call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
2234		"The nondimensional Laplacian Leith constant, \n"//&
2235	0	"often set to 1.0", units="nondim", default=0.0)
2236
2237		call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_QG_Leith, &
2238		"If true, include the beta term in the Leith nonlinear eddy viscosity.", &
2239	0	default=.true.)
2240
2241	0	allocate(CS%Laplac3_const_u(IsdB:IedB,jsd:jed), source=0.0)
2242	0	allocate(CS%Laplac3_const_v(isd:ied,JsdB:JedB), source=0.0)
2243	0	allocate(CS%KH_u_QG(IsdB:IedB,jsd:jed,GV%ke), source=0.0)
2244	0	allocate(CS%KH_v_QG(isd:ied,JsdB:JedB,GV%ke), source=0.0)
2245
2246		! register diagnostics
2247		CS%id_KH_u_QG = register_diag_field('ocean_model', 'KH_u_QG', diag%axesCuL, Time, &
2248	0	'Horizontal viscosity from Leith QG, at u-points', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2249		CS%id_KH_v_QG = register_diag_field('ocean_model', 'KH_v_QG', diag%axesCvL, Time, &
2250	0	'Horizontal viscosity from Leith QG, at v-points', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2251
2252	0	do j=Jsq,Jeq+1 ; do I=is-1,Ieq
2253		! Static factors in the Leith schemes
2254	0	grid_sp_u2 = G%dyCu(I,j)*G%dxCu(I,j)
2255	0	grid_sp_u3 = grid_sp_u2*sqrt(grid_sp_u2)
2256	0	CS%Laplac3_const_u(I,j) = Leith_Lap_const * grid_sp_u3
2257		enddo ; enddo
2258	0	do j=js-1,Jeq ; do I=Isq,Ieq+1
2259		! Static factors in the Leith schemes
2260	0	grid_sp_v2 = G%dyCv(i,J)*G%dxCv(i,J)
2261	0	grid_sp_v3 = grid_sp_v2*sqrt(grid_sp_v2)
2262	0	CS%Laplac3_const_v(i,J) = Leith_Lap_const * grid_sp_v3
2263		enddo ; enddo
2264
2265	0	if (.not. CS%use_stored_slopes) call MOM_error(FATAL, &
2266		"MOM_lateral_mixing_coeffs.F90, VarMix_init: "//&
2267	0	"USE_STORED_SLOPES must be True when using QG Leith.")
2268		endif
2269
2270		! Re-enable variable mixing if one of the schemes was enabled
2271	1	CS%use_variable_mixing = in_use .or. CS%use_variable_mixing
2272
2273	1	end subroutine VarMix_init
2274
2275		!> Destructor for VarMix control structure
2276	1	subroutine VarMix_end(CS)
2277		type(VarMix_CS), intent(inout) :: CS
2278
2279	1	if (allocated(CS%ebt_struct)) deallocate(CS%ebt_struct)
2280	1	if (allocated(CS%sqg_struct)) deallocate(CS%sqg_struct)
2281	1	if (allocated(CS%BS_struct)) deallocate(CS%BS_struct)
2282	1	if (allocated(CS%khth_struct)) deallocate(CS%khth_struct)
2283	1	if (allocated(CS%khtr_struct)) deallocate(CS%khtr_struct)
2284	1	if (allocated(CS%kdgl90_struct)) deallocate(CS%kdgl90_struct)
2285
2286	1	if (allocated(CS%slope_x)) deallocate(CS%slope_x)
2287	1	if (allocated(CS%slope_y)) deallocate(CS%slope_y)
2288
2289	1	if (allocated(CS%SN_u)) deallocate(CS%SN_u)
2290	1	if (allocated(CS%SN_v)) deallocate(CS%SN_v)
2291
2292	1	if (allocated(CS%L2u)) deallocate(CS%L2u)
2293	1	if (allocated(CS%L2v)) deallocate(CS%L2v)
2294
2295	1	if (allocated(CS%Res_fn_h)) deallocate(CS%Res_fn_h)
2296	1	if (allocated(CS%Res_fn_q)) deallocate(CS%Res_fn_q)
2297	1	if (allocated(CS%Res_fn_u)) deallocate(CS%Res_fn_u)
2298	1	if (allocated(CS%Res_fn_v)) deallocate(CS%Res_fn_v)
2299	1	if (allocated(CS%beta_dx2_q)) deallocate(CS%beta_dx2_q)
2300	1	if (allocated(CS%beta_dx2_u)) deallocate(CS%beta_dx2_u)
2301	1	if (allocated(CS%beta_dx2_v)) deallocate(CS%beta_dx2_v)
2302	1	if (allocated(CS%f2_dx2_q)) deallocate(CS%f2_dx2_q)
2303	1	if (allocated(CS%f2_dx2_u)) deallocate(CS%f2_dx2_u)
2304	1	if (allocated(CS%f2_dx2_v)) deallocate(CS%f2_dx2_v)
2305
2306	1	if (allocated(CS%Depth_fn_u)) deallocate(CS%Depth_fn_u)
2307	1	if (allocated(CS%Depth_fn_v)) deallocate(CS%Depth_fn_v)
2308
2309	1	if (allocated(CS%Rd_dx_h)) deallocate(CS%Rd_dx_h)
2310	1	if (allocated(CS%beta_dx2_h)) deallocate(CS%beta_dx2_h)
2311	1	if (allocated(CS%f2_dx2_h)) deallocate(CS%f2_dx2_h)
2312
2313	1	if (allocated(CS%cg1)) deallocate(CS%cg1)
2314
2315	1	if (allocated(CS%Laplac3_const_u)) deallocate(CS%Laplac3_const_u)
2316	1	if (allocated(CS%Laplac3_const_v)) deallocate(CS%Laplac3_const_v)
2317	1	if (allocated(CS%KH_u_QG)) deallocate(CS%KH_u_QG)
2318	1	if (allocated(CS%KH_v_QG)) deallocate(CS%KH_v_QG)
2319
2320	1	end subroutine VarMix_end
2321
2322		!> \namespace mom_lateral_mixing_coeffs
2323		!!
2324		!! This module provides a container for various factors used in prescribing diffusivities, that are
2325		!! a function of the state (in particular the stratification and isoneutral slopes).
2326		!!
2327		!! \section section_Resolution_Function The resolution function
2328		!!
2329		!! The resolution function is expressed in terms of the ratio of grid-spacing to deformation radius.
2330		!! The square of the resolution parameter is
2331		!!
2332		!! \f[
2333		!! R^2 = \frac{L_d^2}{\Delta^2} = \frac{ c_g^2 }{ f^2 \Delta^2 + c_g \beta \Delta^2 }
2334		!! \f]
2335		!!
2336		!! where the grid spacing is calculated as
2337		!!
2338		!! \f[
2339		!! \Delta^2 = \Delta x^2 + \Delta y^2 .
2340		!! \f]
2341		!!
2342		!! \todo Check this reference to Bob on/off paper.
2343		!! The resolution function used in scaling diffusivities (\cite hallberg2013) is
2344		!!
2345		!! \f[
2346		!! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
2347		!! \f]
2348		!!
2349		!! The resolution function can be applied independently to thickness diffusion $module mom_thickness_diffuse$,
2350		!! tracer diffusion $mom_tracer_hordiff$ lateral viscosity $mom_hor_visc$.
2351		!!
2352		!! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
2353		!! Ocean Modelling, 71, pp 92-103. http://dx.doi.org/10.1016/j.ocemod.2013.08.007
2354		!!
2355		!! \| Symbol \| Module parameter \|
2356		!! \| ------ \| --------------- \|
2357		!! \| - \| <code>USE_VARIABLE_MIXING</code> \|
2358		!! \| - \| <code>RESOLN_SCALED_KH</code> \|
2359		!! \| - \| <code>RESOLN_SCALED_KHTH</code> \|
2360		!! \| - \| <code>RESOLN_SCALED_KHTR</code> \|
2361		!! \| \f$ \alpha \f$ \| <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) \|
2362		!! \| \f$ p \f$ \| <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) \|
2363		!! \| \f$ \alpha \f$ \| <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) \|
2364		!! \| \f$ p \f$ \| <code>VISC_RES_FN_POWER</code> (for lateral viscosity) \|
2365		!! \| - \| <code>GILL_EQUATORIAL_LD</code> \|
2366		!!
2367		!!
2368		!!
2369		!! \section section_Vicbeck Visbeck diffusivity
2370		!!
2371		!! This module also calculates factors used in setting the thickness diffusivity similar to a Visbeck et al., 1997,
2372		!! scheme. The factors are combined in mom_thickness_diffuse::thickness_diffuse but calculated in this module.
2373		!!
2374		!! \f[
2375		!! \kappa_h = \alpha_s L_s^2 S N
2376		!! \f]
2377		!!
2378		!! where \f$S\f$ is the magnitude of the isoneutral slope and \f$N\f$ is the Brunt-Vaisala frequency.
2379		!!
2380		!! Visbeck, Marshall, Haine and Spall, 1997: Specification of Eddy Transfer Coefficients in Coarse-Resolution
2381		!! Ocean Circulation Models. J. Phys. Oceanogr. http://dx.doi.org/10.1175/1520-0485(1997)027%3C0381:SOETCI%3E2.0.CO;2
2382		!!
2383		!! \| Symbol \| Module parameter \|
2384		!! \| ------ \| --------------- \|
2385		!! \| - \| <code>USE_VARIABLE_MIXING</code> \|
2386		!! \| \f$ \alpha_s \f$ \| <code>KHTH_SLOPE_CFF</code> (for mom_thickness_diffuse module)\|
2387		!! \| \f$ \alpha_s \f$ \| <code>KHTR_SLOPE_CFF</code> (for mom_tracer_hordiff module)\|
2388		!! \| \f$ L_{s} \f$ \| <code>VISBECK_L_SCALE</code> \|
2389		!! \| \f$ S_{max} \f$ \| <code>VISBECK_MAX_SLOPE</code> \|
2390		!!
2391		!!
2392		!! \section section_vertical_structure_khth Vertical structure function for KhTh
2393		!!
2394		!! The thickness diffusivity can be prescribed a vertical distribution with the shape of the equivalent barotropic
2395		!! velocity mode. The structure function is stored in the control structure for this module (varmix_cs) but is
2396		!! calculated using subroutines in mom_wave_speed.
2397		!!
2398		!! \| Symbol \| Module parameter \|
2399		!! \| ------ \| --------------- \|
2400		!! \| - \| <code>KHTH_USE_EBT_STRUCT</code> \|
2401
2402	0	end module MOM_lateral_mixing_coeffs