src/parameterizations/lateral/MOM_thickness_diffuse.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		!> Isopycnal height diffusion (or Gent McWilliams diffusion)
6		module MOM_thickness_diffuse
7
8		use MOM_debugging, only : hchksum, uvchksum
9		use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
10		use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
11		use MOM_diag_mediator, only : diag_update_remap_grids
12		use MOM_domains, only : pass_var, CORNER, pass_vector
13		use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
14		use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain
15		use MOM_EOS, only : calculate_density_second_derivs
16		use MOM_file_parser, only : get_param, log_version, param_file_type
17		use MOM_grid, only : ocean_grid_type
18		use MOM_io, only : MOM_read_data, slasher
19		use MOM_interface_heights, only : find_eta, thickness_to_dz
20		use MOM_isopycnal_slopes, only : vert_fill_TS
21		use MOM_lateral_mixing_coeffs, only : VarMix_CS
22		use MOM_MEKE_types, only : MEKE_type
23		use MOM_stochastics, only : stochastic_CS
24		use MOM_unit_scaling, only : unit_scale_type
25		use MOM_variables, only : thermo_var_ptrs, cont_diag_ptrs
26		use MOM_verticalGrid, only : verticalGrid_type
27		implicit none ; private
28
29		#include <MOM_memory.h>
30
31		public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
32		public thickness_diffuse_get_KH
33
34		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
35		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
36		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
37		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
38
39		!> Control structure for thickness_diffuse
40		type, public :: thickness_diffuse_CS ; private
41		logical :: initialized = .false. !< True if this control structure has been initialized.
42		real :: Khth !< Background isopycnal depth diffusivity [L2 T-1 ~> m2 s-1]
43		real :: Khth_Slope_Cff !< Slope dependence coefficient of Khth [nondim]
44		real :: max_Khth_CFL !< Maximum value of the diffusive CFL for isopycnal height diffusion [nondim]
45		real :: Khth_Min !< Minimum value of Khth [L2 T-1 ~> m2 s-1]
46		real :: Khth_Max !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max
47		real :: Kh_eta_bg !< Background isopycnal height diffusivity [L2 T-1 ~> m2 s-1]
48		real :: Kh_eta_vel !< Velocity scale that is multiplied by the grid spacing to give
49		!! the isopycnal height diffusivity [L T-1 ~> m s-1]
50		real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim]
51		real :: kappa_smooth !< Vertical diffusivity used to interpolate more sensible values
52		!! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
53		logical :: thickness_diffuse !< If true, interfaces heights are diffused.
54		logical :: full_depth_khth_min !< If true, KHTH_MIN is enforced throughout the whole water column.
55		!! Otherwise, KHTH_MIN is only enforced at the surface. This parameter
56		!! is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0.
57		logical :: use_FGNV_streamfn !< If true, use the streamfunction formulation of
58		!! Ferrari et al., 2010, which effectively emphasizes
59		!! graver vertical modes by smoothing in the vertical.
60		real :: FGNV_scale !< A coefficient scaling the vertical smoothing term in the
61		!! Ferrari et al., 2010, streamfunction formulation [nondim].
62		real :: FGNV_c_min !< A minimum wave speed used in the Ferrari et al., 2010,
63		!! streamfunction formulation [L T-1 ~> m s-1].
64		real :: N2_floor !< A floor for squared buoyancy frequency in the Ferrari et al., 2010,
65		!! streamfunction formulation divided by aspect ratio rescaling factors
66		!! [L2 Z-2 T-2 ~> s-2].
67		logical :: detangle_interfaces !< If true, add 3-d structured interface height
68		!! diffusivities to horizontally smooth jagged layers.
69		real :: detangle_time !< If detangle_interfaces is true, this is the
70		!! timescale over which maximally jagged grid-scale
71		!! thickness variations are suppressed [T ~> s]. This must be
72		!! longer than DT, or 0 (the default) to use DT.
73		integer :: nkml !< number of layers within mixed layer
74		logical :: debug !< write verbose checksums for debugging purposes
75		logical :: use_GME_thickness_diffuse !< If true, passes GM coefficients to MOM_hor_visc for use
76		!! with GME closure.
77		logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC
78		!! framework (Marshall et al., 2012)
79		real :: MEKE_GEOMETRIC_alpha!< The nondimensional coefficient governing the efficiency of
80		!! the GEOMETRIC isopycnal height diffusion [nondim]
81		real :: MEKE_GEOMETRIC_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness
82		!! diffusivity [T-1 ~> s-1].
83		integer :: MEKE_GEOM_answer_date !< The vintage of the expressions in the MEKE_GEOMETRIC
84		!! calculation. Values below 20190101 recover the answers from the
85		!! original implementation, while higher values use expressions that
86		!! satisfy rotational symmetry.
87		logical :: Use_KH_in_MEKE !< If true, uses the isopycnal height diffusivity calculated here to diffuse MEKE.
88		real :: MEKE_min_depth_diff !< The minimum total depth over which to average the diffusivity
89		!! used for MEKE [H ~> m or kg m-2]. When the total depth is less
90		!! than this, the diffusivity is scaled away.
91		logical :: GM_src_alt !< If true, use the GM energy conversion form S^2N^2kappa rather
92		!! than the streamfunction for the GM source term for MEKE.
93		integer :: MEKE_src_answer_date !< The vintage of the expressions in the GM energy conversion
94		!! calculation when MEKE_GM_SRC_ALT is true. Values below 20240601
95		!! recover the answers from the original implementation, while higher
96		!! values use expressions that satisfy rotational symmetry.
97		logical :: MEKE_src_slope_bug !< If true, use a bug that limits the positive values, but not the
98		!! negative values, of the slopes used when MEKE_GM_SRC_ALT is true.
99		!! When this is true, it breaks rotational symmetry.
100		logical :: use_GM_work_bug !< If true, use the incorrect sign for the
101		!! top-level work tendency on the top layer.
102		real :: Stanley_det_coeff !< The coefficient correlating SGS temperature variance with the mean
103		!! temperature gradient in the deterministic part of the Stanley parameterization.
104		!! Negative values disable the scheme. [nondim]
105		logical :: read_khth !< If true, read a file containing the spatially varying horizontal
106		!! isopycnal height diffusivity
107		logical :: use_stanley_gm !< If true, also use the Stanley parameterization in MOM_thickness_diffuse
108
109		type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics
110		real, allocatable :: GMwork(:,:) !< Work by isopycnal height diffusion [R Z L2 T-3 ~> W m-2]
111		real, allocatable :: diagSlopeX(:,:,:) !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim]
112		real, allocatable :: diagSlopeY(:,:,:) !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim]
113
114		real, allocatable :: Kh_eta_u(:,:) !< Isopycnal height diffusivities at u points [L2 T-1 ~> m2 s-1]
115		real, allocatable :: Kh_eta_v(:,:) !< Isopycnal height diffusivities in v points [L2 T-1 ~> m2 s-1]
116
117		real, allocatable :: KH_u_GME(:,:,:) !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
118		real, allocatable :: KH_v_GME(:,:,:) !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
119		real, allocatable :: khth2d(:,:) !< 2D isopycnal height diffusivity at h-points [L2 T-1 ~> m2 s-1]
120
121		!>@{
122		!! Diagnostic identifier
123		integer :: id_uhGM = -1, id_vhGM = -1, id_GMwork = -1
124		integer :: id_KH_u = -1, id_KH_v = -1, id_KH_t = -1
125		integer :: id_KH_u1 = -1, id_KH_v1 = -1, id_KH_t1 = -1
126		integer :: id_slope_x = -1, id_slope_y = -1
127		integer :: id_sfn_unlim_x = -1, id_sfn_unlim_y = -1, id_sfn_x = -1, id_sfn_y = -1
128		!>@}
129		end type thickness_diffuse_CS
130
131		contains
132
133		!> Calculates isopycnal height diffusion coefficients and applies isopycnal height diffusion
134		!! by modifying to the layer thicknesses, h. Diffusivities are limited to ensure stability.
135		!! Also returns along-layer mass fluxes used in the continuity equation.
136	12	subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS, STOCH)
137		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
138		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
139		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
140		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
141		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Accumulated zonal mass flux
142		!! [L2 H ~> m3 or kg]
143		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Accumulated meridional mass flux
144		!! [L2 H ~> m3 or kg]
145		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
146		real, intent(in) :: dt !< Time increment [T ~> s]
147		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
148		type(VarMix_CS), target, intent(in) :: VarMix !< Variable mixing coefficients
149		type(cont_diag_ptrs), intent(inout) :: CDp !< Diagnostics for the continuity equation
150		type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse
151		type(stochastic_CS), intent(inout) :: STOCH !< Stochastic control structure
152		! Local variables
153	12	real :: e(SZI_(G),SZJ_(G),SZK_(GV)+1) ! heights of interfaces, relative to mean
154		! sea level [Z ~> m], positive up.
155	24	real :: uhD(SZIB_(G),SZJ_(G),SZK_(GV)) ! Diffusive u*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]
156	24	real :: vhD(SZI_(G),SZJB_(G),SZK_(GV)) ! Diffusive v*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]
157
158		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
159	24	KH_u, & ! Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
160	24	int_slope_u ! A nondimensional ratio from 0 to 1 that gives the relative
161		! weighting of the interface slopes to that calculated also
162		! using density gradients at u points. The physically correct
163		! slopes occur at 0, while 1 is used for numerical closures [nondim].
164		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
165	24	KH_v, & ! Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
166	24	int_slope_v ! A nondimensional ratio from 0 to 1 that gives the relative
167		! weighting of the interface slopes to that calculated also
168		! using density gradients at v points. The physically correct
169		! slopes occur at 0, while 1 is used for numerical closures [nondim].
170		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
171	24	KH_t ! diagnosed diffusivity at tracer points [L2 T-1 ~> m2 s-1]
172
173		real, dimension(SZIB_(G),SZJ_(G)) :: &
174	24	KH_u_CFL ! The maximum stable isopycnal height diffusivity at u grid points [L2 T-1 ~> m2 s-1]
175		real, dimension(SZI_(G),SZJB_(G)) :: &
176	24	KH_v_CFL ! The maximum stable isopycnal height diffusivity at v grid points [L2 T-1 ~> m2 s-1]
177		real, dimension(SZI_(G),SZJ_(G)) :: &
178	24	htot ! The sum of the total layer thicknesses [H ~> m or kg m-2]
179	24	real :: Khth_Loc_u(SZIB_(G),SZJ_(G)) ! The isopycnal height diffusivity at u points [L2 T-1 ~> m2 s-1]
180	24	real :: Khth_Loc_v(SZI_(G),SZJB_(G)) ! The isopycnal height diffusivity at v points [L2 T-1 ~> m2 s-1]
181		real :: h_neglect ! A thickness that is so small it is usually lost
182		! in roundoff and can be neglected [H ~> m or kg m-2].
183		real, dimension(:,:), pointer :: cg1 => null() !< Wave speed [L T-1 ~> m s-1]
184	24	real :: hu(SZI_(G),SZJ_(G)) ! A thickness-based mask at u points, used for diagnostics [nondim]
185	24	real :: hv(SZI_(G),SZJ_(G)) ! A thickness-based mask at v points, used for diagnostics [nondim]
186	24	real :: KH_u_lay(SZI_(G),SZJ_(G)) ! Diagnostic of isopycnal height diffusivities at u-points averaged
187		! to layer centers [L2 T-1 ~> m2 s-1]
188	12	real :: KH_v_lay(SZI_(G),SZJ_(G)) ! Diagnostic of isopycnal height diffusivities at v-points averaged
189		! to layer centers [L2 T-1 ~> m2 s-1]
190		logical :: use_VarMix, Resoln_scaled, Depth_scaled, use_stored_slopes, khth_use_vert_struct, use_Visbeck
191		logical :: use_QG_Leith
192		integer :: i, j, k, is, ie, js, je, nz
193
194	12	if (.not. CS%initialized) call MOM_error(FATAL, "MOM_thickness_diffuse: "//&
195	0	"Module must be initialized before it is used.")
196
197		if ((.not.CS%thickness_diffuse) &
198	12	.or. .not. (CS%Khth > 0.0 .or. CS%read_khth &
199	0	.or. VarMix%use_variable_mixing)) return
200
201	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
202	12	h_neglect = GV%H_subroundoff
203
204	12	if (allocated(MEKE%GM_src)) then
205	87132	do j=js,je ; do i=is,ie ; MEKE%GM_src(i,j) = 0. ; enddo ; enddo
206		endif
207
208	12	use_VarMix = .false. ; Resoln_scaled = .false. ; use_stored_slopes = .false.
209	12	khth_use_vert_struct = .false. ; use_Visbeck = .false. ; use_QG_Leith = .false.
210	12	Depth_scaled = .false.
211
212	12	if (VarMix%use_variable_mixing) then
213	12	use_VarMix = VarMix%use_variable_mixing .and. (CS%KHTH_Slope_Cff > 0.)
214	12	Resoln_scaled = VarMix%Resoln_scaled_KhTh
215	12	Depth_scaled = VarMix%Depth_scaled_KhTh
216	12	use_stored_slopes = VarMix%use_stored_slopes
217	12	khth_use_vert_struct = allocated(VarMix%khth_struct)
218	12	use_Visbeck = VarMix%use_Visbeck
219	12	use_QG_Leith = VarMix%use_QG_Leith_GM
220	12	if (allocated(VarMix%cg1)) cg1 => VarMix%cg1
221		else
222	0	cg1 => null()
223		endif
224
225
226		!$OMP parallel do default(shared)
227	87852	do j=js,je ; do I=is-1,ie
228		KH_u_CFL(I,j) = (0.25*CS%max_Khth_CFL) / &
229	87840	(dt * ((G%IdxCu(I,j)G%IdxCu(I,j)) + (G%IdyCu(I,j)G%IdyCu(I,j))))
230		enddo ; enddo
231		!$OMP parallel do default(shared)
232	88584	do J=js-1,je ; do i=is,ie
233		KH_v_CFL(i,J) = (0.25*CS%max_Khth_CFL) / &
234	88572	(dt * ((G%IdxCv(i,J)G%IdxCv(i,J)) + (G%IdyCv(i,J)G%IdyCv(i,J))))
235		enddo ; enddo
236
237		! Calculates interface heights, e, in [Z ~> m].
238		!$omp target update to(h)
239		!$omp target enter data map(alloc: e)
240	12	call find_eta(h, tv, G, GV, US, e, halo_size=1)
241		!$omp target exit data map(from: e)
242
243		! Set the diffusivities.
244		!$OMP parallel default(shared)
245	12	if (.not. CS%read_khth) then
246		!$OMP do
247	87852	do j=js,je ; do I=is-1,ie
248	87840	Khth_loc_u(I,j) = CS%Khth
249		enddo ; enddo
250		else ! use 2d KHTH that was read in from file
251		!$OMP do
252	0	do j=js,je ; do I=is-1,ie
253	0	Khth_loc_u(I,j) = 0.5 * (CS%khth2d(i,j) + CS%khth2d(i+1,j))
254		enddo ; enddo
255		endif
256
257	12	if (use_VarMix) then
258	0	if (use_Visbeck) then
259		!$OMP do
260	0	do j=js,je ; do I=is-1,ie
261		Khth_loc_u(I,j) = Khth_loc_u(I,j) + &
262	0	CS%KHTH_Slope_CffVarMix%L2u(I,j) VarMix%SN_u(I,j)
263		enddo ; enddo
264		endif
265		endif
266
267	12	if (allocated(MEKE%Kh)) then
268	12	if (CS%MEKE_GEOMETRIC) then
269		!$OMP do
270	0	do j=js,je ; do I=is-1,ie
271		Khth_loc_u(I,j) = Khth_loc_u(I,j) + G%OBCmaskCu(I,j) * CS%MEKE_GEOMETRIC_alpha * &
272		0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)) / &
273	0	(VarMix%SN_u(I,j) + CS%MEKE_GEOMETRIC_epsilon)
274		enddo ; enddo
275		else
276	87852	do j=js,je ; do I=is-1,ie
277	87840	Khth_loc_u(I,j) = Khth_loc_u(I,j) + MEKE%KhTh_facsqrt(MEKE%Kh(i,j)MEKE%Kh(i+1,j))
278		enddo ; enddo
279		endif
280		endif
281
282	12	if (Resoln_scaled) then
283		!$OMP do
284	87852	do j=js,je ; do I=is-1,ie
285	87840	Khth_loc_u(I,j) = Khth_loc_u(I,j) * VarMix%Res_fn_u(I,j)
286		enddo ; enddo
287		endif
288
289	12	if (Depth_scaled) then
290		!$OMP do
291	0	do j=js,je ; do I=is-1,ie
292	0	Khth_loc_u(I,j) = Khth_loc_u(I,j) * VarMix%Depth_fn_u(I,j)
293		enddo ; enddo
294		endif
295
296	12	if (CS%Khth_Max > 0) then
297		!$OMP do
298	0	do j=js,je ; do I=is-1,ie
299	0	Khth_loc_u(I,j) = max(CS%Khth_Min, min(Khth_loc_u(I,j), CS%Khth_Max))
300		enddo ; enddo
301		else
302		!$OMP do
303	87852	do j=js,je ; do I=is-1,ie
304	87840	Khth_loc_u(I,j) = max(CS%Khth_Min, Khth_loc_u(I,j))
305		enddo ; enddo
306		endif
307		!$OMP do
308	87852	do j=js,je ; do I=is-1,ie
309	87840	KH_u(I,j,1) = min(KH_u_CFL(I,j), Khth_loc_u(I,j))
310		enddo ; enddo
311
312	12	if (khth_use_vert_struct) then
313	0	if (CS%full_depth_khth_min) then
314		!$OMP do
315	0	do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
316	0	KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i+1,j,k-1) )
317	0	KH_u(I,j,K) = max(KH_u(I,j,K), CS%Khth_Min)
318		enddo ; enddo ; enddo
319		else
320		!$OMP do
321	0	do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
322	0	KH_u(I,j,K) = KH_u(I,j,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i+1,j,k-1) )
323		enddo ; enddo ; enddo
324		endif
325		else
326		!$OMP do
327	6588912	do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
328	6588000	KH_u(I,j,K) = KH_u(I,j,1)
329		enddo ; enddo ; enddo
330		endif
331
332	12	if (use_VarMix) then
333	0	if (use_QG_Leith) then
334		!$OMP do
335	0	do k=1,nz ; do j=js,je ; do I=is-1,ie
336	0	KH_u(I,j,k) = VarMix%KH_u_QG(I,j,k)
337		enddo ; enddo ; enddo
338		endif
339		endif
340
341	12	if (CS%use_GME_thickness_diffuse) then
342		!$OMP do
343	0	do k=1,nz+1 ; do j=js,je ; do I=is-1,ie
344	0	CS%KH_u_GME(I,j,k) = KH_u(I,j,k)
345		enddo ; enddo ; enddo
346		endif
347
348	12	if (.not. CS%read_khth) then
349		!$OMP do
350	88584	do J=js-1,je ; do i=is,ie
351	88572	Khth_loc_v(i,J) = CS%Khth
352		enddo ; enddo
353		else ! read KHTH from file
354		!$OMP do
355	0	do J=js-1,je ; do i=is,ie
356	0	Khth_loc_v(i,J) = 0.5 * (CS%khth2d(i,j) + CS%khth2d(i,j+1))
357		enddo ; enddo
358		endif
359
360	12	if (use_VarMix) then
361	0	if (use_Visbeck) then
362		!$OMP do
363	0	do J=js-1,je ; do i=is,ie
364	0	Khth_loc_v(i,J) = Khth_loc_v(i,J) + CS%KHTH_Slope_CffVarMix%L2v(i,J)VarMix%SN_v(i,J)
365		enddo ; enddo
366		endif
367		endif
368	12	if (allocated(MEKE%Kh)) then
369	12	if (CS%MEKE_GEOMETRIC) then
370		!$OMP do
371	0	do J=js-1,je ; do i=is,ie
372		Khth_loc_v(i,J) = Khth_loc_v(i,J) + G%OBCmaskCv(i,J) * CS%MEKE_GEOMETRIC_alpha * &
373		0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / &
374	0	(VarMix%SN_v(i,J) + CS%MEKE_GEOMETRIC_epsilon)
375		enddo ; enddo
376		else
377	88584	do J=js-1,je ; do i=is,ie
378	88572	Khth_loc_v(i,J) = Khth_loc_v(i,J) + MEKE%KhTh_facsqrt(MEKE%Kh(i,j)MEKE%Kh(i,j+1))
379		enddo ; enddo
380		endif
381		endif
382
383	12	if (Resoln_scaled) then
384		!$OMP do
385	88584	do J=js-1,je ; do i=is,ie
386	88572	Khth_loc_v(i,J) = Khth_loc_v(i,J) * VarMix%Res_fn_v(i,J)
387		enddo ; enddo
388		endif
389
390	12	if (Depth_scaled) then
391		!$OMP do
392	0	do J=js-1,je ; do i=is,ie
393	0	Khth_loc_v(i,J) = Khth_loc_v(i,J) * VarMix%Depth_fn_v(i,J)
394		enddo ; enddo
395		endif
396
397	12	if (CS%Khth_Max > 0) then
398		!$OMP do
399	0	do J=js-1,je ; do i=is,ie
400	0	Khth_loc_v(i,J) = max(CS%Khth_Min, min(Khth_loc_v(i,J), CS%Khth_Max))
401		enddo ; enddo
402		else
403		!$OMP do
404	88584	do J=js-1,je ; do i=is,ie
405	88572	Khth_loc_v(i,J) = max(CS%Khth_Min, Khth_loc_v(i,J))
406		enddo ; enddo
407		endif
408
409	12	if (CS%max_Khth_CFL > 0.0) then
410		!$OMP do
411	88584	do J=js-1,je ; do i=is,ie
412	88572	KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_loc_v(i,J))
413		enddo ; enddo
414		endif
415
416	12	if (khth_use_vert_struct) then
417	0	if (CS%full_depth_khth_min) then
418		!$OMP do
419	0	do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
420	0	KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i,j+1,k-1) )
421	0	KH_v(i,J,K) = max(KH_v(i,J,K), CS%Khth_Min)
422		enddo ; enddo ; enddo
423		else
424		!$OMP do
425	0	do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
426	0	KH_v(i,J,K) = KH_v(i,J,1) * 0.5 * ( VarMix%khth_struct(i,j,k-1) + VarMix%khth_struct(i,j+1,k-1) )
427		enddo ; enddo ; enddo
428		endif
429		else
430		!$OMP do
431	6643812	do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
432	6642900	KH_v(i,J,K) = KH_v(i,J,1)
433		enddo ; enddo ; enddo
434		endif
435
436	12	if (use_VarMix) then
437	0	if (use_QG_Leith) then
438		!$OMP do
439	0	do k=1,nz ; do J=js-1,je ; do i=is,ie
440	0	KH_v(i,J,k) = VarMix%KH_v_QG(i,J,k)
441		enddo ; enddo ; enddo
442		endif
443		endif
444
445	12	if (CS%use_GME_thickness_diffuse) then
446		!$OMP do
447	0	do k=1,nz+1 ; do J=js-1,je ; do i=is,ie
448	0	CS%KH_v_GME(i,J,k) = KH_v(i,J,k)
449		enddo ; enddo ; enddo
450		endif
451
452	12	if (allocated(MEKE%Kh)) then
453	12	if (CS%MEKE_GEOMETRIC) then
454	0	if (CS%MEKE_GEOM_answer_date < 20190101) then
455		!$OMP do
456	0	do j=js,je ; do i=is,ie
457		! This does not give bitwise rotational symmetry.
458		MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / &
459		(0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j) + &
460		VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + &
461	0	CS%MEKE_GEOMETRIC_epsilon)
462		enddo ; enddo
463		else
464		!$OMP do
465	0	do j=js,je ; do i=is,ie
466		! With the additional parentheses this gives bitwise rotational symmetry.
467		MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / &
468		(0.25*((VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)) + &
469		(VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1))) + &
470	0	CS%MEKE_GEOMETRIC_epsilon)
471		enddo ; enddo
472		endif
473		endif
474		endif
475
476		!$OMP do
477	6676764	do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo
478		!$OMP do
479	6732396	do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo
480		!$OMP end parallel
481
482	12	if (CS%detangle_interfaces) then
483		call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, &
484	0	CS, int_slope_u, int_slope_v)
485		endif
486
487	12	if ((CS%Kh_eta_bg > 0.0) .or. (CS%Kh_eta_vel > 0.0)) then
488	0	call add_interface_Kh(G, GV, US, CS, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, int_slope_u, int_slope_v)
489		endif
490
491	12	if (CS%debug) then
492		call uvchksum("Kh_[uv]", Kh_u, Kh_v, G%HI, haloshift=0, &
493	0	unscale=(US%L_to_m*2)US%s_to_T, scalar_pair=.true.)
494		call uvchksum("Kh_[uv]_CFL", Kh_u_CFL, Kh_v_CFL, G%HI, haloshift=0, &
495	0	unscale=(US%L_to_m*2)US%s_to_T, scalar_pair=.true.)
496	0	if (Resoln_scaled) then
497		call uvchksum("Res_fn_[uv]", VarMix%Res_fn_u, VarMix%Res_fn_v, G%HI, haloshift=0, &
498	0	unscale=1.0, scalar_pair=.true.)
499		endif
500	0	call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, G%HI, haloshift=0)
501	0	call hchksum(h, "thickness_diffuse_1 h", G%HI, haloshift=1, unscale=GV%H_to_m)
502	0	call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, unscale=US%Z_to_m)
503	0	if (use_stored_slopes) then
504		call uvchksum("VarMix%slope_[xy]", VarMix%slope_x, VarMix%slope_y, &
505	0	G%HI, haloshift=0, unscale=US%Z_to_L)
506		endif
507	0	if (associated(tv%eqn_of_state)) then
508	0	call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1, unscale=US%C_to_degC)
509	0	call hchksum(tv%S, "thickness_diffuse S", G%HI, haloshift=1, unscale=US%S_to_ppt)
510		endif
511		endif
512
513		! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S
514	12	if (STOCH%skeb_use_gm) then
515	0	if (use_stored_slopes) then
516		call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, &
517		int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y, &
518	0	STOCH=STOCH, VarMix=VarMix)
519		else
520		call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, &
521	0	int_slope_u, int_slope_v, STOCH=STOCH, VarMix=VarMix)
522		endif
523		else
524	12	if (use_stored_slopes) then
525		call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, &
526	12	int_slope_u, int_slope_v, VarMix%slope_x, VarMix%slope_y)
527		else
528		call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, CS, &
529	0	int_slope_u, int_slope_v)
530		endif
531		endif
532
533	12	if (VarMix%use_variable_mixing) then
534	12	if (allocated(MEKE%Rd_dx_h) .and. allocated(VarMix%Rd_dx_h)) then
535		!$OMP parallel do default(shared)
536	87132	do j=js,je ; do i=is,ie
537	87120	MEKE%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j)
538		enddo ; enddo
539		endif
540		endif
541
542		! offer diagnostic fields for averaging
543	12	if (query_averaging_enabled(CS%diag)) then
544	12	if (CS%id_uhGM > 0) call post_data(CS%id_uhGM, uhD, CS%diag)
545	12	if (CS%id_vhGM > 0) call post_data(CS%id_vhGM, vhD, CS%diag)
546	12	if (CS%id_GMwork > 0) call post_data(CS%id_GMwork, CS%GMwork, CS%diag)
547	12	if (CS%id_KH_u > 0) call post_data(CS%id_KH_u, KH_u, CS%diag)
548	12	if (CS%id_KH_v > 0) call post_data(CS%id_KH_v, KH_v, CS%diag)
549	12	if (CS%id_KH_u1 > 0) call post_data(CS%id_KH_u1, KH_u(:,:,1), CS%diag)
550	12	if (CS%id_KH_v1 > 0) call post_data(CS%id_KH_v1, KH_v(:,:,1), CS%diag)
551
552		! Diagnose diffusivity at T-cell point. Do a simple average, rather than a
553		! thickness-weighted average, so that KH_t is depth-independent when KH_u and KH_v
554		! are depth independent. If a thickness-weighted average were used, the variations
555		! of thickness could give a spurious depth dependence to the diagnosed KH_t.
556	12	if (CS%id_KH_t > 0 .or. CS%id_KH_t1 > 0 .or. CS%Use_KH_in_MEKE) then
557	0	do k=1,nz
558		! thicknesses across u and v faces, converted to 0/1 mask
559		! layer average of the interface diffusivities KH_u and KH_v
560	0	do j=js,je ; do I=is-1,ie
561		! This expression uses harmonic mean thicknesses:
562		! hu(I,j) = 2.0h(i,j,k)h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
563		! This expression is a 0/1 mask based on depths where there are thick layers:
564	0	hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0
565	0	KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1))
566		enddo ; enddo
567	0	do J=js-1,je ; do i=is,ie
568		! This expression uses harmonic mean thicknesses:
569		! hv(i,J) = 2.0h(i,j,k)h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
570		! This expression is a 0/1 mask based on depths where there are thick layers:
571	0	hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0
572	0	KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1))
573		enddo ; enddo
574		! diagnose diffusivity at T-points
575	0	do j=js,je ; do i=is,ie
576		Kh_t(i,j,k) = (((hu(I-1,j)KH_u_lay(i-1,j)) + (hu(I,j)KH_u_lay(I,j))) + &
577		((hv(i,J-1)KH_v_lay(i,J-1)) + (hv(i,J)KH_v_lay(i,J)))) / &
578	0	((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + 1.0e-20)
579		! Use this denominator instead if hu and hv are actual thicknesses rather than a 0/1 mask:
580		! ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + h_neglect)
581		enddo ; enddo
582		enddo
583
584	0	if (CS%Use_KH_in_MEKE) then
585	0	MEKE%Kh_diff(:,:) = 0.0
586	0	htot(:,:) = 0.0
587	0	do k=1,nz
588	0	do j=js,je ; do i=is,ie
589	0	MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + Kh_t(i,j,k) * h(i,j,k)
590	0	htot(i,j) = htot(i,j) + h(i,j,k)
591		enddo ; enddo
592		enddo
593
594	0	do j=js,je ; do i=is,ie
595	0	MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(CS%MEKE_min_depth_diff, htot(i,j))
596		enddo ; enddo
597		endif
598
599	0	if (CS%id_KH_t > 0) call post_data(CS%id_KH_t, KH_t, CS%diag)
600	0	if (CS%id_KH_t1 > 0) call post_data(CS%id_KH_t1, KH_t(:,:,1), CS%diag)
601		endif
602
603		endif
604
605		!$OMP parallel do default(shared)
606	912	do k=1,nz
607	6588900	do j=js,je ; do I=is-1,ie
608	6534000	uhtr(I,j,k) = uhtr(I,j,k) + uhD(I,j,k) * dt
609	6588000	if (associated(CDp%uhGM)) CDp%uhGM(I,j,k) = uhD(I,j,k)
610		enddo ; enddo
611	6643800	do J=js-1,je ; do i=is,ie
612	6588000	vhtr(i,J,k) = vhtr(i,J,k) + vhD(i,J,k) * dt
613	6642900	if (associated(CDp%vhGM)) CDp%vhGM(i,J,k) = vhD(i,J,k)
614		enddo ; enddo
615	6534912	do j=js,je ; do i=is,ie
616		h(i,j,k) = h(i,j,k) - dt * G%IareaT(i,j) * &
617	6480000	((uhD(I,j,k) - uhD(I-1,j,k)) + (vhD(i,J,k) - vhD(i,J-1,k)))
618	6534000	if (h(i,j,k) < GV%Angstrom_H) h(i,j,k) = GV%Angstrom_H
619		enddo ; enddo
620		enddo
621
622		! Whenever thickness changes let the diag manager know, target grids
623		! for vertical remapping may need to be regenerated.
624		! This needs to happen after the H update and before the next post_data.
625	12	call diag_update_remap_grids(CS%diag)
626
627	12	if (CS%debug) then
628		call uvchksum("thickness_diffuse [uv]hD", uhD, vhD, &
629	0	G%HI, haloshift=0, unscale=GV%H_to_mUS%L_to_m2US%s_to_T)
630		call uvchksum("thickness_diffuse [uv]htr", uhtr, vhtr, &
631	0	G%HI, haloshift=0, unscale=US%L_to_m*2GV%H_to_m)
632	0	call hchksum(h, "thickness_diffuse h", G%HI, haloshift=0, unscale=GV%H_to_m)
633		endif
634
635		end subroutine thickness_diffuse
636
637		!> Calculates parameterized layer transports for use in the continuity equation.
638		!! Fluxes are limited to give positive definite thicknesses.
639		!! Called by thickness_diffuse().
640	24	subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, &
641	24	CS, int_slope_u, int_slope_v, slope_x, slope_y, STOCH, VarMix)
642		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
643		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
644		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
645		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
646		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface positions [Z ~> m]
647		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: Kh_u !< Isopycnal height diffusivity
648		!! at u points [L2 T-1 ~> m2 s-1]
649		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: Kh_v !< Isopycnal height diffusivity
650		!! at v points [L2 T-1 ~> m2 s-1]
651		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
652		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: uhD !< Zonal mass fluxes
653		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
654		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: vhD !< Meridional mass fluxes
655		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
656		real, dimension(:,:), pointer :: cg1 !< Wave speed [L T-1 ~> m s-1]
657		real, intent(in) :: dt !< Time increment [T ~> s]
658		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
659		type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse
660		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: int_slope_u !< Ratio that determine how much of
661		!! the isopycnal slopes are taken directly from
662		!! the interface slopes without consideration of
663		!! density gradients [nondim].
664		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: int_slope_v !< Ratio that determine how much of
665		!! the isopycnal slopes are taken directly from
666		!! the interface slopes without consideration of
667		!! density gradients [nondim].
668		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopyc. slope at u [Z L-1 ~> nondim]
669		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopyc. slope at v [Z L-1 ~> nondim]
670		type(stochastic_CS), optional, intent(inout) :: STOCH !< Stochastic control structure
671		type(VarMix_CS), target, optional, intent(in) :: VarMix !< Variable mixing coefficents
672
673		! Local variables
674		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
675	24	T, & ! The temperature [C ~> degC], with the values in
676		! in massless layers filled vertically by diffusion.
677	24	S, & ! The filled salinity [S ~> ppt], with the values in
678		! in massless layers filled vertically by diffusion.
679	24	h_avail, & ! The mass available for diffusion out of each face, divided
680		! by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
681	24	h_frac ! The fraction of the mass in the column above the bottom
682		! interface of a layer that is within a layer [nondim]. 0<h_frac<=1
683	24	real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
684		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
685	24	Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below) [Z L-1 ~> nondim]
686	24	hN2_y_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency
687		! at v-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2],
688		! used for calculating the potential energy release
689		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
690	24	Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below) [Z L-1 ~> nondim]
691	24	hN2_x_PE ! Harmonic mean of thicknesses around the interfaces times the buoyancy frequency
692		! at u-points with unit conversion factors [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2],
693		! used for calculating the potential energy release
694		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
695	24	pres, & ! The pressure at an interface [R L2 T-2 ~> Pa].
696	24	h_avail_rsum ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1].
697		real, dimension(SZIB_(G)) :: &
698	24	drho_dT_u, & ! The derivative of density with temperature at u points [R C-1 ~> kg m-3 degC-1]
699	24	drho_dS_u ! The derivative of density with salinity at u points [R S-1 ~> kg m-3 ppt-1].
700	24	real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs()
701		! with various units that will be ignored [various]
702		real, dimension(SZI_(G)) :: &
703	24	drho_dT_v, & ! The derivative of density with temperature at v points [R C-1 ~> kg m-3 degC-1]
704	24	drho_dS_v, & ! The derivative of density with salinity at v points [R S-1 ~> kg m-3 ppt-1].
705	24	drho_dT_dT_h, & ! The second derivative of density with temperature at h points [R C-2 ~> kg m-3 degC-2]
706	24	drho_dT_dT_hr ! The second derivative of density with temperature at h (+1) points [R C-2 ~> kg m-3 degC-2]
707	24	real :: uhtot(SZIB_(G),SZJ_(G)) ! The vertical sum of uhD [H L2 T-1 ~> m3 s-1 or kg s-1].
708	24	real :: vhtot(SZI_(G),SZJB_(G)) ! The vertical sum of vhD [H L2 T-1 ~> m3 s-1 or kg s-1].
709		real, dimension(SZIB_(G)) :: &
710	24	T_u, & ! Temperature on the interface at the u-point [C ~> degC].
711	24	S_u, & ! Salinity on the interface at the u-point [S ~> ppt].
712	24	pres_u ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
713		real, dimension(SZI_(G)) :: &
714	24	T_v, & ! Temperature on the interface at the v-point [C ~> degC].
715	24	S_v, & ! Salinity on the interface at the v-point [S ~> ppt].
716	24	pres_v, & ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
717	24	T_h, & ! Temperature on the interface at the h-point [C ~> degC].
718	24	S_h, & ! Salinity on the interface at the h-point [S ~> ppt].
719	24	pres_h, & ! Pressure on the interface at the h-point [R L2 T-2 ~> Pa].
720	24	T_hr, & ! Temperature on the interface at the h (+1) point [C ~> degC].
721	24	S_hr, & ! Salinity on the interface at the h (+1) point [S ~> ppt].
722	24	pres_hr ! Pressure on the interface at the h (+1) point [R L2 T-2 ~> Pa].
723	24	real :: Work_u(SZIB_(G),SZJ_(G)) ! The work done by the isopycnal height diffusion
724		! integrated over u-point water columns [R Z L4 T-3 ~> W]
725	24	real :: Work_v(SZI_(G),SZJB_(G)) ! The work done by the isopycnal height diffusion
726		! integrated over v-point water columns [R Z L4 T-3 ~> W]
727		real :: Work_h ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2].
728		real :: PE_release_h ! The amount of potential energy released by GM averaged over an h-cell
729		! [R Z L2 T-3 ~> W m-2]. The calculation equals rho0 * h * S^2 * N^2 * kappa_GM.
730		real :: I4dt ! 1 / 4 dt [T-1 ~> s-1].
731		real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A)
732		! and below (B) the interface times the grid spacing [R ~> kg m-3].
733		real :: drdjA, drdjB ! Along layer meridional potential density gradients in the layers above (A)
734		! and below (B) the interface times the grid spacing [R ~> kg m-3].
735		real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3].
736	24	real :: drdi_u(SZIB_(G),SZK_(GV)) ! Copy of drdi at u-points [R ~> kg m-3].
737	24	real :: drdj_v(SZI_(G),SZK_(GV)) ! Copy of drdj at v-points [R ~> kg m-3].
738	24	real :: drdkDe_u(SZIB_(G),SZK_(GV)+1) ! Lateral difference of product of drdk and e at u-points
739		! [Z R ~> kg m-2].
740	24	real :: drdkDe_v(SZI_(G),SZK_(GV)+1) ! Lateral difference of product of drdk and e at v-points
741		! [Z R ~> kg m-2].
742		real :: hg2A, hg2B, hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
743		real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2].
744		real :: dzg2A, dzg2B ! Squares of geometric mean vertical layer extents [Z2 ~> m2].
745		real :: dzaA, dzaB ! Arithmetic mean vertical layer extents [Z ~> m].
746		real :: dzaL, dzaR ! Temporary vertical layer extents [Z ~> m]
747		real :: wtA, wtB ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6]
748		real :: wtL, wtR ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5]
749		real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
750		real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4].
751		real :: dz_harm ! Harmonic mean layer vertical extent [Z ~> m].
752	24	real :: c2_dz_u(SZIB_(G),SZK_(GV)+1) ! Wave speed squared divided by dz at u-points [L2 Z-1 T-2 ~> m s-2]
753	24	real :: c2_dz_v(SZI_(G),SZK_(GV)+1) ! Wave speed squared divided by dz at v-points [L2 Z-1 T-2 ~> m s-2]
754	24	real :: dzN2_u(SZIB_(G),SZK_(GV)+1) ! Vertical extent times N2 at interfaces above u-points times
755		! rescaling factors from vertical to horizontal distances [L2 Z-1 T-2 ~> m s-2]
756	24	real :: dzN2_v(SZI_(G),SZK_(GV)+1) ! Vertical extent times N2 at interfaces above v-points times
757		! rescaling factors from vertical to horizontal distances [L2 Z-1 T-2 ~> m s-2]
758		real :: Sfn_est ! A preliminary estimate (before limiting) of the overturning
759		! streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1].
760	24	real :: Sfn_unlim_u(SZIB_(G),SZK_(GV)+1) ! Volume streamfunction for u-points [Z L2 T-1 ~> m3 s-1]
761	24	real :: Sfn_unlim_v(SZI_(G),SZK_(GV)+1) ! Volume streamfunction for v-points [Z L2 T-1 ~> m3 s-1]
762	24	real :: slope2_Ratio_u(SZIB_(G),SZK_(GV)+1) ! The ratio of the slope squared to slope_max squared [nondim]
763	24	real :: slope2_Ratio_v(SZI_(G),SZK_(GV)+1) ! The ratio of the slope squared to slope_max squared [nondim]
764		real :: Sfn_in_h ! The overturning streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] (note that
765		! the units are different from other Sfn vars).
766		real :: Sfn_safe ! The streamfunction that goes linearly back to 0 at the surface
767		! [H L2 T-1 ~> m3 s-1 or kg s-1]. This is a good value to use when the
768		! slope is so large as to be meaningless, usually due to weak stratification.
769		real :: Slope ! The slope of density surfaces, calculated in a way that is always
770		! between -1 and 1 after undoing dimensional scaling, [Z L-1 ~> nondim]
771		real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
772		real :: I_slope_max2 ! The inverse of slope_max squared [L2 Z-2 ~> nondim].
773		real :: h_neglect ! A thickness that is so small it is usually lost
774		! in roundoff and can be neglected [H ~> m or kg m-2].
775		real :: hn_2 ! Half of h_neglect [H ~> m or kg m-2].
776		real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
777		real :: dz_neglect ! A thickness [Z ~> m], that is so small it is usually lost
778		! in roundoff and can be neglected [Z ~> m].
779		real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2]
780		real :: G_scale ! The gravitational acceleration times a unit conversion
781		! factor [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
782		logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
783		logical :: find_work ! If true, find the change in energy due to the fluxes.
784		integer :: nk_linear ! The number of layers over which the streamfunction goes to 0.
785		real :: G_rho0 ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2].
786		real :: Rho_avg ! The in situ density averaged to an interface [R ~> kg m-3]
787		real :: N2_floor ! A floor for N2 to avoid degeneracy in the elliptic solver
788		! times unit conversion factors [L2 Z-2 T-2 ~> s-2]
789		real :: N2_unlim ! An unlimited estimate of the buoyancy frequency
790		! times unit conversion factors [L2 Z-2 T-2 ~> s-2]
791		real :: Z_to_H ! A conversion factor from heights to thicknesses, perhaps based on
792		! a spatially variable local density [H Z-1 ~> nondim or kg m-3]
793	24	real :: diag_sfn_x(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! Diagnostic of the x-face streamfunction
794		! [H L2 T-1 ~> m3 s-1 or kg s-1]
795	24	real :: diag_sfn_unlim_x(SZIB_(G),SZJ_(G),SZK_(GV)+1) ! Diagnostic of the x-face streamfunction before
796		! applying limiters [Z L2 T-1 ~> m3 s-1]
797	24	real :: diag_sfn_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction
798		! [H L2 T-1 ~> m3 s-1 or kg s-1]
799	24	real :: diag_sfn_unlim_y(SZI_(G),SZJB_(G),SZK_(GV)+1) ! Diagnostic of the y-face streamfunction before
800		! applying limiters [Z L2 T-1 ~> m3 s-1]
801		! applying limiters [H L2 T-1 ~> m3 s-1 or kg s-1]
802	12	real, allocatable :: skeb_gm_work(:,:) ! Temp array to hold GM work for SKEB
803	12	real, allocatable :: skeb_ebt_norm2(:,:) ! Used to normalize EBT for SKEB
804
805		logical :: present_slope_x, present_slope_y, calc_derivatives
806		integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of
807		! state calculations at u-points.
808		integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of
809		! state calculations at v-points.
810		integer, dimension(2) :: EOSdom_h1 ! The shifted i-computational domain to use for equation of
811		! state calculations at h points with 1 extra halo point
812		logical :: use_stanley, skeb_use_gm
813		integer :: is, ie, js, je, nz, IsdB, halo
814		integer :: i, j, k
815	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke ; IsdB = G%IsdB
816
817	12	I4dt = 0.25 / dt
818	12	I_slope_max2 = 1.0 / (CS%slope_max**2)
819
820	12	h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*2 ; hn_2 = 0.5h_neglect
821	12	dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect**2
822	12	if (GV%Boussinesq) G_rho0 = GV%g_Earth / GV%Rho0
823	12	N2_floor = CS%N2_floor
824
825	12	use_EOS = associated(tv%eqn_of_state)
826	12	present_slope_x = PRESENT(slope_x)
827	12	present_slope_y = PRESENT(slope_y)
828
829	12	use_stanley = CS%use_stanley_gm
830
831	12	skeb_use_gm = .false.
832	12	if (present(STOCH)) skeb_use_gm = STOCH%skeb_use_gm
833	12	if (skeb_use_gm) then
834	0	allocate(skeb_gm_work(is:ie,js:je), source=0.)
835	0	allocate(skeb_ebt_norm2(is:ie,js:je), source=0.)
836		endif
837
838	12	nk_linear = max(GV%nkml, 1)
839
840	8063004	Slope_x_PE(:,:,:) = 0.0
841	8118636	Slope_y_PE(:,:,:) = 0.0
842	8063004	hN2_x_PE(:,:,:) = 0.0
843	8118636	hN2_y_PE(:,:,:) = 0.0
844
845	12	find_work = allocated(MEKE%GM_src)
846	12	find_work = (allocated(CS%GMwork) .or. find_work)
847	12	find_work = (skeb_use_gm .or. find_work)
848
849	12	if (use_EOS) then
850	12	halo = 1 ! Default halo to fill is 1
851		!$omp target enter data map(to: h, tv%T, tv%S)
852		!$omp target enter data map(alloc: T, S)
853	12	call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo, larger_h_denom=.true.)
854		!$omp target exit data map(from: T, S)
855		!$omp target exit data map(release: h, tv%T, tv%S)
856		endif
857
858		! Rescale the thicknesses, perhaps using the specific volume.
859	12	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
860
861	12	if (CS%use_FGNV_streamfn .and. .not. associated(cg1)) call MOM_error(FATAL, &
862	0	"cg1 must be associated when using FGNV streamfunction.")
863
864		!$OMP parallel default(shared)
865		! Find the maximum and minimum permitted streamfunction.
866		!$OMP do
867	91524	do j=js-1,je+1 ; do i=is-1,ie+1
868	90768	h_avail_rsum(i,j,1) = 0.0
869	90768	pres(i,j,1) = 0.0
870	90768	if (associated(tv%p_surf)) then ; pres(i,j,1) = tv%p_surf(i,j) ; endif
871
872	90768	h_avail(i,j,1) = max(I4dtG%areaT(i,j)(h(i,j,1)-GV%Angstrom_H),0.0)
873	90768	h_avail_rsum(i,j,2) = h_avail(i,j,1)
874	90768	h_frac(i,j,1) = 1.0
875	91512	pres(i,j,2) = pres(i,j,1) + (GV%g_EarthGV%H_to_RZ) h(i,j,1)
876		enddo ; enddo
877	756	do j=js-1,je+1
878	6772644	do k=2,nz ; do i=is-1,ie+1
879	6716832	h_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H),0.0)
880	6716832	h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
881	6716832	h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) &
882	6716832	h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
883	6771888	pres(i,j,K+1) = pres(i,j,K) + (GV%g_EarthGV%H_to_RZ) h(i,j,k)
884		enddo ; enddo
885		enddo
886		!$OMP do
887	87852	do j=js,je ; do I=is-1,ie
888	87840	uhtot(I,j) = 0.0 ; Work_u(I,j) = 0.0
889		enddo ; enddo
890		!$OMP do
891	88584	do J=js-1,je ; do i=is,ie
892	88572	vhtot(i,J) = 0.0 ; Work_v(i,J) = 0.0
893		enddo ; enddo
894		!$OMP end parallel
895
896	12	if (CS%id_sfn_x > 0) then ; diag_sfn_x(:,:,1) = 0.0 ; diag_sfn_x(:,:,nz+1) = 0.0 ; endif
897	12	if (CS%id_sfn_y > 0) then ; diag_sfn_y(:,:,1) = 0.0 ; diag_sfn_y(:,:,nz+1) = 0.0 ; endif
898	12	if (CS%id_sfn_unlim_x > 0) then ; diag_sfn_unlim_x(:,:,1) = 0.0 ; diag_sfn_unlim_x(:,:,nz+1) = 0.0 ; endif
899	12	if (CS%id_sfn_unlim_y > 0) then ; diag_sfn_unlim_y(:,:,1) = 0.0 ; diag_sfn_unlim_y(:,:,nz+1) = 0.0 ; endif
900
901	12	EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1)
902	12	EOSdom_v(:) = EOS_domain(G%HI)
903	12	EOSdom_h1(:) = EOS_domain(G%HI, halo=1)
904
905		!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
906		!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz,dz_neglect,dz_neglect2, &
907		!$OMP h_neglect2,hn_2,I_slope_max2,int_slope_u,KH_u,uhtot, &
908		!$OMP h_frac,h_avail_rsum,uhD,h_avail,Work_u,CS,slope_x,cg1, &
909		!$OMP diag_sfn_x,diag_sfn_unlim_x,N2_floor,EOSdom_u,EOSdom_h1, &
910		!$OMP use_stanley,present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) &
911		!$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u,G_scale, &
912		!$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
913		!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h,N2_unlim, &
914		!$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
915		!$OMP dzg2A,dzg2B,dzaA,dzaB,dz_harm,Z_to_H, &
916		!$OMP drdx,mag_grad2,Slope,slope2_Ratio_u,dzN2_u, &
917		!$OMP Sfn_unlim_u,Rho_avg,drdi_u,drdkDe_u,c2_dz_u, &
918		!$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
919	732	do j=js,je
920	87840	do I=is-1,ie ; dzN2_u(I,1) = 0. ; dzN2_u(I,nz+1) = 0. ; enddo
921	54000	do K=nz,2,-1
922	53280	if (find_work .and. .not.(use_EOS)) then
923	0	drdiA = 0.0 ; drdiB = 0.0
924	0	drdkL = GV%Rlay(k) - GV%Rlay(k-1) ; drdkR = drdkL
925		endif
926
927		calc_derivatives = use_EOS .and. (k >= nk_linear) .and. &
928	53280	(find_work .or. .not. present_slope_x .or. CS%use_FGNV_streamfn .or. use_stanley)
929
930		! Calculate the zonal fluxes and gradients.
931	53280	if (calc_derivatives) then
932	6500160	do I=is-1,ie
933	6446880	pres_u(I) = 0.5*(pres(i,j,K) + pres(i+1,j,K))
934	6446880	T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1)))
935	6500160	S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
936		enddo
937		call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, &
938	53280	tv%eqn_of_state, EOSdom_u)
939		endif
940	53280	if (use_stanley) then
941	0	do i=is-1,ie+1
942	0	pres_h(i) = pres(i,j,K)
943	0	T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1))
944	0	S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1))
945		enddo
946
947		! The second line below would correspond to arguments
948		! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
949		call calculate_density_second_derivs(T_h, S_h, pres_h, &
950		scrap, scrap, drho_dT_dT_h, scrap, scrap, &
951	0	tv%eqn_of_state, EOSdom_h1)
952		endif
953
954	6500880	do I=is-1,ie
955	6446880	if (calc_derivatives) then
956		! Estimate the horizontal density gradients along layers.
957		drdiA = drho_dT_u(I) * (T(i+1,j,k-1)-T(i,j,k-1)) + &
958	6446880	drho_dS_u(I) * (S(i+1,j,k-1)-S(i,j,k-1))
959		drdiB = drho_dT_u(I) * (T(i+1,j,k)-T(i,j,k)) + &
960	6446880	drho_dS_u(I) * (S(i+1,j,k)-S(i,j,k))
961
962		! Estimate the vertical density gradients times the grid spacing.
963		drdkL = (drho_dT_u(I) * (T(i,j,k)-T(i,j,k-1)) + &
964	6446880	drho_dS_u(I) * (S(i,j,k)-S(i,j,k-1)))
965		drdkR = (drho_dT_u(I) * (T(i+1,j,k)-T(i+1,j,k-1)) + &
966	6446880	drho_dS_u(I) * (S(i+1,j,k)-S(i+1,j,k-1)))
967	6446880	drdkDe_u(I,K) = (drdkR * e(i+1,j,K)) - (drdkL * e(i,j,K))
968	0	elseif (find_work) then ! This is used in pure stacked SW mode
969	0	drdkDe_u(I,K) = (drdkR * e(i+1,j,K)) - (drdkL * e(i,j,K))
970		endif
971	6446880	if (use_stanley) then
972		! Correction to the horizontal density gradient due to nonlinearity in
973		! the EOS rectifying SGS temperature anomalies
974		drdiA = drdiA + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k-1)) - &
975	0	(drho_dT_dT_h(i) * tv%varT(i,j,k-1)) )
976		drdiB = drdiB + 0.5 * ((drho_dT_dT_h(i+1) * tv%varT(i+1,j,k)) - &
977	0	(drho_dT_dT_h(i) * tv%varT(i,j,k)) )
978		endif
979	6446880	if (find_work) drdi_u(I,k) = drdiB
980
981	6446880	if (k > nk_linear) then
982	6446880	if (use_EOS) then
983	6446880	if (CS%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x) then
984	6446880	hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
985	6446880	hg2R = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
986	6446880	haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
987	6446880	haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
988	6446880	if (GV%Boussinesq) then
989	6446880	dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z
990	0	elseif (GV%semi_Boussinesq) then
991	0	dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
992	0	dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect
993		else
994	0	dzaL = 0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect
995	0	dzaR = 0.5*(dz(i+1,j,k-1) + dz(i+1,j,k)) + dz_neglect
996		endif
997		! Use the harmonic mean thicknesses to weight the horizontal gradients.
998		! These unnormalized weights have been rearranged to minimize divisions.
999	6446880	wtL = hg2L(haRdzaR) ; wtR = hg2R(haLdzaL)
1000
1001	6446880	drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaLwtL) + (dzaRwtR))
1002		! The expression for drdz above is mathematically equivalent to:
1003		! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
1004		! ((hg2L/haL) + (hg2R/haR))
1005	6446880	hg2A = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
1006	6446880	hg2B = h(i,j,k)*h(i+1,j,k) + h_neglect2
1007	6446880	haA = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
1008	6446880	haB = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
1009
1010	6446880	if (GV%Boussinesq) then
1011	6446880	N2_unlim = drdz*G_rho0
1012		else
1013		N2_unlim = (GV%g_EarthGV%RZ_to_H) &
1014	0	(((wtL * drdkL) + (wtR * drdkR)) / ((haLwtL) + (haRwtR)))
1015		endif
1016
1017	6446880	dzg2A = dz(i,j,k-1)*dz(i+1,j,k-1) + dz_neglect2
1018	6446880	dzg2B = dz(i,j,k)*dz(i+1,j,k) + dz_neglect2
1019	6446880	dzaA = 0.5*(dz(i,j,k-1) + dz(i+1,j,k-1)) + dz_neglect
1020	6446880	dzaB = 0.5*(dz(i,j,k) + dz(i+1,j,k)) + dz_neglect
1021		! dzN2_u is used with the FGNV streamfunction formulation
1022	6446880	dzN2_u(I,K) = (0.5 * ( dzg2A / dzaA + dzg2B / dzaB )) * max(N2_unlim, N2_floor)
1023	6446880	if (find_work .and. CS%GM_src_alt) &
1024	0	hN2_x_PE(I,j,k) = (0.5 * ( hg2A / haA + hg2B / haB )) * max(N2_unlim, N2_floor)
1025		endif
1026
1027	6446880	if (present_slope_x) then
1028	6446880	Slope = slope_x(I,j,k)
1029	6446880	slope2_Ratio_u(I,K) = Slope*2 I_slope_max2
1030		else
1031		! Use the harmonic mean thicknesses to weight the horizontal gradients.
1032		! These unnormalized weights have been rearranged to minimize divisions.
1033	0	wtA = hg2AhaB ; wtB = hg2BhaA
1034		! This is the gradient of density along geopotentials.
1035		drdx = ((wtA * drdiA + wtB * drdiB) / (wtA + wtB) - &
1036	0	drdz * (e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j)
1037
1038		! This estimate of slope is accurate for small slopes, but bounded
1039		! to be between -1 and 1.
1040	0	mag_grad2 = (US%Z_to_Ldrdx)2 + drdz*2
1041	0	if (mag_grad2 > 0.0) then
1042	0	Slope = drdx / sqrt(mag_grad2)
1043	0	slope2_Ratio_u(I,K) = Slope*2 I_slope_max2
1044		else ! Just in case mag_grad2 = 0 ever.
1045	0	Slope = 0.0
1046	0	slope2_Ratio_u(I,K) = 1.0e20 ! Force the use of the safe streamfunction.
1047		endif
1048		endif
1049
1050		! Adjust real slope by weights that bias towards slope of interfaces
1051		! that ignore density gradients along layers.
1052		Slope = (1.0 - int_slope_u(I,j,K)) * Slope + &
1053	6446880	int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j))
1054	6446880	slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K)
1055
1056	6446880	if (CS%MEKE_src_slope_bug) then
1057	0	Slope_x_PE(I,j,k) = MIN(Slope, CS%slope_max)
1058		else
1059	6446880	Slope_x_PE(I,j,k) = Slope
1060	6446880	if (Slope > CS%slope_max) Slope_x_PE(I,j,k) = CS%slope_max
1061	6446880	if (Slope < -CS%slope_max) Slope_x_PE(I,j,k) = -CS%slope_max
1062		endif
1063	6446880	if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope
1064
1065		! Estimate the streamfunction at each interface [H L2 T-1 ~> m3 s-1 or kg s-1].
1066	6446880	Sfn_unlim_u(I,K) = -(KH_u(I,j,K)G%dy_Cu(I,j))Slope
1067
1068		! Avoid moving dense water upslope from below the level of
1069		! the bottom on the receiving side.
1070	6446880	if (Sfn_unlim_u(I,K) > 0.0) then ! The flow below this interface is positive.
1071	0	if (e(i,j,K) < e(i+1,j,nz+1)) then
1072	0	Sfn_unlim_u(I,K) = 0.0 ! This is not uhtot, because it may compensate for
1073		! deeper flow in very unusual cases.
1074	0	elseif (e(i+1,j,nz+1) > e(i,j,K+1)) then
1075		! Scale the transport with the fraction of the donor layer above
1076		! the bottom on the receiving side.
1077		Sfn_unlim_u(I,K) = Sfn_unlim_u(I,K) * ((e(i,j,K) - e(i+1,j,nz+1)) / &
1078	0	((e(i,j,K) - e(i,j,K+1)) + dz_neglect))
1079		endif
1080		else
1081	6446880	if (e(i+1,j,K) < e(i,j,nz+1)) then ; Sfn_unlim_u(I,K) = 0.0
1082	5679480	elseif (e(i,j,nz+1) > e(i+1,j,K+1)) then
1083		Sfn_unlim_u(I,K) = Sfn_unlim_u(I,K) * ((e(i+1,j,K) - e(i,j,nz+1)) / &
1084	29316	((e(i+1,j,K) - e(i+1,j,K+1)) + dz_neglect))
1085		endif
1086		endif
1087
1088		else ! .not. use_EOS
1089	0	if (present_slope_x) then
1090	0	Slope = slope_x(I,j,k)
1091		else
1092	0	Slope = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu_OBCmask(I,j)
1093		endif
1094	0	if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope
1095	0	Sfn_unlim_u(I,K) = -(KH_u(I,j,K)G%dy_Cu(I,j))Slope
1096	0	dzN2_u(I,K) = GV%g_prime(K)
1097		endif ! if (use_EOS)
1098		else ! if (k > nk_linear)
1099	0	dzN2_u(I,K) = N2_floor * dz_neglect
1100	0	Sfn_unlim_u(I,K) = 0.
1101		endif ! if (k > nk_linear)
1102	6500160	if (CS%id_sfn_unlim_x>0) diag_sfn_unlim_x(I,j,K) = Sfn_unlim_u(I,K)
1103		enddo ! i-loop
1104		enddo ! k-loop
1105
1106	720	if (CS%use_FGNV_streamfn) then
1107	0	do k=1,nz ; do I=is-1,ie ; if (G%OBCmaskCu(I,j)>0.) then
1108		dz_harm = max( dz_neglect, &
1109	0	2. * dz(i,j,k) * dz(i+1,j,k) / ( ( dz(i,j,k) + dz(i+1,j,k) ) + dz_neglect ) )
1110	0	c2_dz_u(I,k) = CS%FGNV_scale * ( 0.5( cg1(i,j) + cg1(i+1,j) ) )*2 / dz_harm
1111		endif ; enddo ; enddo
1112
1113		! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
1114	0	do I=is-1,ie
1115	0	if (G%OBCmaskCu(I,j)>0.) then
1116	0	do K=2,nz
1117	0	Sfn_unlim_u(I,K) = (1. + CS%FGNV_scale) * Sfn_unlim_u(I,K)
1118		enddo
1119	0	call streamfn_solver(nz, c2_dz_u(I,:), dzN2_u(I,:), Sfn_unlim_u(I,:))
1120		else
1121	0	do K=2,nz
1122	0	Sfn_unlim_u(I,K) = 0.
1123		enddo
1124		endif
1125		enddo
1126		endif
1127
1128	54012	do K=nz,2,-1
1129	6500880	do I=is-1,ie
1130
1131	6446880	if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then
1132		Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i+1,j,k) + h(i+1,j,k-1))) + 4.0*hn_2 ) / &
1133		( (((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k)) + ((h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1))) + &
1134	0	(((h(i+1,j,k)+hn_2)tv%SpV_avg(i+1,j,k)) + ((h(i+1,j,k-1)+hn_2)tv%SpV_avg(i+1,j,k-1))) )
1135		! Use an average density to convert the volume streamfunction estimate into a mass streamfunction.
1136	0	Z_to_H = GV%RZ_to_H*Rho_avg
1137		else
1138	6446880	Z_to_H = GV%Z_to_H
1139		endif
1140
1141	6446880	if (k > nk_linear) then
1142	6446880	if (use_EOS) then
1143
1144	6446880	if (uhtot(I,j) <= 0.0) then
1145		! The transport that must balance the transport below is positive.
1146	6446880	Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i,j,k))
1147		else ! (uhtot(I,j) > 0.0)
1148	0	Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i+1,j,k))
1149		endif
1150
1151		! Determine the actual streamfunction at each interface.
1152	6446880	Sfn_est = (Z_to_HSfn_unlim_u(I,K) + slope2_Ratio_u(I,K)Sfn_safe) / (1.0 + slope2_Ratio_u(I,K))
1153		else ! When use_EOS is false, the layers are constant density.
1154	0	Sfn_est = Z_to_H*Sfn_unlim_u(I,K)
1155		endif
1156
1157		! Make sure that there is enough mass above to allow the streamfunction
1158		! to satisfy the boundary condition of 0 at the surface.
1159	6446880	Sfn_in_H = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i+1,j,K))
1160
1161		! The actual transport is limited by the mass available in the two
1162		! neighboring grid cells.
1163		uhD(I,j,k) = max(min((Sfn_in_H - uhtot(I,j)), h_avail(i,j,k)), &
1164	6446880	-h_avail(i+1,j,k))
1165
1166	6446880	if (CS%id_sfn_x>0) diag_sfn_x(I,j,K) = diag_sfn_x(I,j,K+1) + uhD(I,j,k)
1167		! sfn_x(I,j,K) = max(min(Sfn_in_h, uhtot(I,j)+h_avail(i,j,k)), &
1168		! uhtot(I,j)-h_avail(i+1,j,K))
1169		! sfn_slope_x(I,j,K) = max(uhtot(I,j)-h_avail(i+1,j,k), &
1170		! min(uhtot(I,j)+h_avail(i,j,k), &
1171		! min(h_avail_rsum(i+1,j,K), max(-h_avail_rsum(i,j,K), &
1172		! (KH_u(I,j,K)G%dy_Cu(I,j)) ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) )) ))
1173		else ! k <= nk_linear
1174		! Balance the deeper flow with a return flow uniformly distributed
1175		! though the remaining near-surface layers. This is the same as
1176		! using Sfn_safe above. There is no need to apply the limiters in
1177		! this case.
1178	0	if (uhtot(I,j) <= 0.0) then
1179	0	uhD(I,j,k) = -uhtot(I,j) * h_frac(i,j,k)
1180		else ! (uhtot(I,j) > 0.0)
1181	0	uhD(I,j,k) = -uhtot(I,j) * h_frac(i+1,j,k)
1182		endif
1183
1184	0	if (CS%id_sfn_x>0) diag_sfn_x(I,j,K) = diag_sfn_x(I,j,K+1) + uhD(I,j,k)
1185		! if (sfn_slope_x(I,j,K+1) <= 0.0) then
1186		! sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i,j,k))
1187		! else
1188		! sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i+1,j,k))
1189		! endif
1190
1191		endif
1192
1193	6446880	uhtot(I,j) = uhtot(I,j) + uhD(I,j,k)
1194
1195	6500160	if (find_work) then
1196		! This is the energy tendency based on the original profiles, and does
1197		! not include any nonlinear terms due to a finite time step (which would
1198		! involve interactions between the fluxes through the different faces.
1199		! A second order centered estimate is used for the density transferred
1200		! between water columns.
1201
1202	6446880	if (allocated(tv%SpV_avg)) then
1203	0	G_scale = GV%H_to_RZ * GV%g_Earth / Rho_avg
1204		else
1205	6446880	G_scale = GV%g_Earth * GV%H_to_Z
1206		endif
1207
1208		Work_u(I,j) = Work_u(I,j) + G_scale * &
1209		( uhtot(I,j) * drdkDe_u(I,K) - &
1210		(uhD(I,j,k) * drdi_u(I,k)) * 0.25 * &
1211	6446880	((e(i,j,K) + e(i,j,K+1)) + (e(i+1,j,K) + e(i+1,j,K+1))) )
1212		endif
1213
1214		enddo
1215		enddo ! end of k-loop
1216		enddo ! end of j-loop
1217
1218		! Calculate the meridional fluxes and gradients.
1219
1220		!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S,dz, &
1221		!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,dz_neglect2, &
1222		!$OMP h_neglect2,int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, &
1223		!$OMP I_slope_max2,vhD,h_avail,Work_v,CS,slope_y,cg1,hn_2,&
1224		!$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,&
1225		!$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) &
1226		!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v,S_h,S_hr, &
1227		!$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA,G_scale, &
1228		!$OMP drho_dT_dT_h,drho_dT_dT_hr,scrap,pres_h,T_h,T_hr, &
1229		!$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,pres_hr, &
1230		!$OMP dzg2A,dzg2B,dzaA,dzaB,dz_harm,Z_to_H, &
1231		!$OMP drdy,mag_grad2,Slope,slope2_Ratio_v,dzN2_v,N2_unlim, &
1232		!$OMP Sfn_unlim_v,Rho_avg,drdj_v,drdkDe_v,c2_dz_v, &
1233		!$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
1234	744	do J=js-1,je
1235	54900	do K=nz,2,-1
1236	54168	if (find_work .and. .not.(use_EOS)) then
1237	0	drdjA = 0.0 ; drdjB = 0.0
1238	0	drdkL = GV%Rlay(k) - GV%Rlay(k-1) ; drdkR = drdkL
1239		endif
1240
1241		calc_derivatives = use_EOS .and. (k >= nk_linear) .and. &
1242	54168	(find_work .or. .not. present_slope_y .or. CS%use_FGNV_streamfn .or. use_stanley)
1243
1244	54168	if (calc_derivatives) then
1245	6554328	do i=is,ie
1246	6500160	pres_v(i) = 0.5*(pres(i,j,K) + pres(i,j+1,K))
1247	6500160	T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
1248	6554328	S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
1249		enddo
1250		call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, &
1251	54168	tv%eqn_of_state, EOSdom_v)
1252		endif
1253	54168	if (use_stanley) then
1254	0	do i=is,ie
1255	0	pres_h(i) = pres(i,j,K)
1256	0	T_h(i) = 0.5*(T(i,j,k) + T(i,j,k-1))
1257	0	S_h(i) = 0.5*(S(i,j,k) + S(i,j,k-1))
1258
1259	0	pres_hr(i) = pres(i,j+1,K)
1260	0	T_hr(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1))
1261	0	S_hr(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1))
1262		enddo
1263
1264		! The second line below would correspond to arguments
1265		! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
1266		call calculate_density_second_derivs(T_h, S_h, pres_h, &
1267		scrap, scrap, drho_dT_dT_h, scrap, scrap, &
1268	0	tv%eqn_of_state, EOSdom_v)
1269		call calculate_density_second_derivs(T_hr, S_hr, pres_hr, &
1270		scrap, scrap, drho_dT_dT_hr, scrap, scrap, &
1271	0	tv%eqn_of_state, EOSdom_v)
1272		endif
1273	6555060	do i=is,ie
1274	6500160	if (calc_derivatives) then
1275		! Estimate the horizontal density gradients along layers.
1276		drdjA = drho_dT_v(i) * (T(i,j+1,k-1)-T(i,j,k-1)) + &
1277	6500160	drho_dS_v(i) * (S(i,j+1,k-1)-S(i,j,k-1))
1278		drdjB = drho_dT_v(i) * (T(i,j+1,k)-T(i,j,k)) + &
1279	6500160	drho_dS_v(i) * (S(i,j+1,k)-S(i,j,k))
1280
1281		! Estimate the vertical density gradients times the grid spacing.
1282		drdkL = (drho_dT_v(i) * (T(i,j,k)-T(i,j,k-1)) + &
1283	6500160	drho_dS_v(i) * (S(i,j,k)-S(i,j,k-1)))
1284		drdkR = (drho_dT_v(i) * (T(i,j+1,k)-T(i,j+1,k-1)) + &
1285	6500160	drho_dS_v(i) * (S(i,j+1,k)-S(i,j+1,k-1)))
1286	6500160	drdkDe_v(i,K) = (drdkR * e(i,j+1,K)) - (drdkL * e(i,j,K))
1287	0	elseif (find_work) then ! This is used in pure stacked SW mode
1288	0	drdkDe_v(i,K) = (drdkR * e(i,j+1,K)) - (drdkL * e(i,j,K))
1289		endif
1290	6500160	if (use_stanley) then
1291		! Correction to the horizontal density gradient due to nonlinearity in
1292		! the EOS rectifying SGS temperature anomalies
1293		drdjA = drdjA + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k-1)) - &
1294	0	(drho_dT_dT_h(i) * tv%varT(i,j,k-1)) )
1295		drdjB = drdjB + 0.5 * ((drho_dT_dT_hr(i) * tv%varT(i,j+1,k)) - &
1296	0	(drho_dT_dT_h(i) * tv%varT(i,j,k)) )
1297		endif
1298
1299	6500160	if (find_work) drdj_v(i,k) = drdjB
1300
1301	6500160	if (k > nk_linear) then
1302	6500160	if (use_EOS) then
1303	6500160	if (CS%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y) then
1304	6500160	hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
1305	6500160	hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
1306	6500160	haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
1307	6500160	haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
1308
1309	6500160	if (GV%Boussinesq) then
1310	6500160	dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z
1311	0	elseif (GV%semi_Boussinesq) then
1312	0	dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
1313	0	dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect
1314		else
1315	0	dzaL = 0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect
1316	0	dzaR = 0.5*(dz(i,j+1,k-1) + dz(i,j+1,k)) + dz_neglect
1317		endif
1318		! Use the harmonic mean thicknesses to weight the horizontal gradients.
1319		! These unnormalized weights have been rearranged to minimize divisions.
1320	6500160	wtL = hg2L(haRdzaR) ; wtR = hg2R(haLdzaL)
1321
1322	6500160	drdz = ((wtL * drdkL) + (wtR * drdkR)) / ((dzaLwtL) + (dzaRwtR))
1323		! The expression for drdz above is mathematically equivalent to:
1324		! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
1325		! ((hg2L/haL) + (hg2R/haR))
1326	6500160	hg2A = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
1327	6500160	hg2B = h(i,j,k)*h(i,j+1,k) + h_neglect2
1328	6500160	haA = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
1329	6500160	haB = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
1330
1331	6500160	if (GV%Boussinesq) then
1332	6500160	N2_unlim = drdz*G_rho0
1333		else
1334		N2_unlim = (GV%g_EarthGV%RZ_to_H) &
1335	0	(((wtL * drdkL) + (wtR * drdkR)) / ((haLwtL) + (haRwtR)))
1336		endif
1337
1338	6500160	dzg2A = dz(i,j,k-1)*dz(i,j+1,k-1) + dz_neglect2
1339	6500160	dzg2B = dz(i,j,k)*dz(i,j+1,k) + dz_neglect2
1340	6500160	dzaA = 0.5*(dz(i,j,k-1) + dz(i,j+1,k-1)) + dz_neglect
1341	6500160	dzaB = 0.5*(dz(i,j,k) + dz(i,j+1,k)) + dz_neglect
1342
1343		! dzN2_v is used with the FGNV streamfunction formulation
1344	6500160	dzN2_v(i,K) = (0.5( dzg2A / dzaA + dzg2B / dzaB )) max(N2_unlim, N2_floor)
1345	6500160	if (find_work .and. CS%GM_src_alt) &
1346	0	hN2_y_PE(i,J,k) = (0.5( hg2A / haA + hg2B / haB )) max(N2_unlim, N2_floor)
1347		endif
1348	6500160	if (present_slope_y) then
1349	6500160	Slope = slope_y(i,J,k)
1350	6500160	slope2_Ratio_v(i,K) = Slope*2 I_slope_max2
1351		else
1352		! Use the harmonic mean thicknesses to weight the horizontal gradients.
1353		! These unnormalized weights have been rearranged to minimize divisions.
1354	0	wtA = hg2AhaB ; wtB = hg2BhaA
1355		! This is the gradient of density along geopotentials.
1356		drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - &
1357	0	drdz * (e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J)
1358
1359		! This estimate of slope is accurate for small slopes, but bounded
1360		! to be between -1 and 1.
1361	0	mag_grad2 = (US%Z_to_Ldrdy)2 + drdz*2
1362	0	if (mag_grad2 > 0.0) then
1363	0	Slope = drdy / sqrt(mag_grad2)
1364	0	slope2_Ratio_v(i,K) = Slope*2 I_slope_max2
1365		else ! Just in case mag_grad2 = 0 ever.
1366	0	Slope = 0.0
1367	0	slope2_Ratio_v(i,K) = 1.0e20 ! Force the use of the safe streamfunction.
1368		endif
1369		endif
1370
1371		! Adjust real slope by weights that bias towards slope of interfaces
1372		! that ignore density gradients along layers.
1373		Slope = (1.0 - int_slope_v(i,J,K)) * Slope + &
1374	6500160	int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J))
1375	6500160	slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K)
1376
1377	6500160	if (CS%MEKE_src_slope_bug) then
1378	0	Slope_y_PE(i,J,k) = MIN(Slope, CS%slope_max)
1379		else
1380	6500160	Slope_y_PE(i,J,k) = Slope
1381	6500160	if (Slope > CS%slope_max) Slope_y_PE(i,J,k) = CS%slope_max
1382	6500160	if (Slope < -CS%slope_max) Slope_y_PE(i,J,k) = -CS%slope_max
1383		endif
1384	6500160	if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope
1385
1386	6500160	Sfn_unlim_v(i,K) = -((KH_v(i,J,K)G%dx_Cv(i,J))Slope)
1387
1388		! Avoid moving dense water upslope from below the level of
1389		! the bottom on the receiving side.
1390	6500160	if (Sfn_unlim_v(i,K) > 0.0) then ! The flow below this interface is positive.
1391	0	if (e(i,j,K) < e(i,j+1,nz+1)) then
1392	0	Sfn_unlim_v(i,K) = 0.0 ! This is not vhtot, because it may compensate for
1393		! deeper flow in very unusual cases.
1394	0	elseif (e(i,j+1,nz+1) > e(i,j,K+1)) then
1395		! Scale the transport with the fraction of the donor layer above
1396		! the bottom on the receiving side.
1397		Sfn_unlim_v(i,K) = Sfn_unlim_v(i,K) * ((e(i,j,K) - e(i,j+1,nz+1)) / &
1398	0	((e(i,j,K) - e(i,j,K+1)) + dz_neglect))
1399		endif
1400		else
1401	6500160	if (e(i,j+1,K) < e(i,j,nz+1)) then ; Sfn_unlim_v(i,K) = 0.0
1402	5904816	elseif (e(i,j,nz+1) > e(i,j+1,K+1)) then
1403		Sfn_unlim_v(i,K) = Sfn_unlim_v(i,K) * ((e(i,j+1,K) - e(i,j,nz+1)) / &
1404	21144	((e(i,j+1,K) - e(i,j+1,K+1)) + dz_neglect))
1405		endif
1406		endif
1407
1408		else ! .not. use_EOS
1409	0	if (present_slope_y) then
1410	0	Slope = slope_y(i,J,k)
1411		else
1412	0	Slope = (e(i,j+1,K)-e(i,j,K)) * G%IdyCv_OBCmask(i,J)
1413		endif
1414	0	if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope
1415	0	Sfn_unlim_v(i,K) = -((KH_v(i,J,K)G%dx_Cv(i,J))Slope)
1416	0	dzN2_v(i,K) = GV%g_prime(K)
1417		endif ! if (use_EOS)
1418		else ! if (k > nk_linear)
1419	0	dzN2_v(i,K) = N2_floor * dz_neglect
1420	0	Sfn_unlim_v(i,K) = 0.
1421		endif ! if (k > nk_linear)
1422	6554328	if (CS%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,J,K) = Sfn_unlim_v(i,K)
1423		enddo ! i-loop
1424		enddo ! k-loop
1425
1426	732	if (CS%use_FGNV_streamfn) then
1427	0	do k=1,nz ; do i=is,ie ; if (G%OBCmaskCv(i,J)>0.) then
1428		dz_harm = max( dz_neglect, &
1429	0	2. * dz(i,j,k) * dz(i,j+1,k) / ( ( dz(i,j,k) + dz(i,j+1,k) ) + dz_neglect ) )
1430	0	c2_dz_v(i,k) = CS%FGNV_scale * ( 0.5( cg1(i,j) + cg1(i,j+1) ) )*2 / dz_harm
1431		endif ; enddo ; enddo
1432
1433		! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
1434	0	do i=is,ie
1435	0	if (G%OBCmaskCv(i,J)>0.) then
1436	0	do K=2,nz
1437	0	Sfn_unlim_v(i,K) = (1. + CS%FGNV_scale) * Sfn_unlim_v(i,K)
1438		enddo
1439	0	call streamfn_solver(nz, c2_dz_v(i,:), dzN2_v(i,:), Sfn_unlim_v(i,:))
1440		else
1441	0	do K=2,nz
1442	0	Sfn_unlim_v(i,K) = 0.
1443		enddo
1444		endif
1445		enddo
1446		endif
1447
1448	54912	do K=nz,2,-1
1449	6555060	do i=is,ie
1450	6500160	if (allocated(tv%SpV_avg) .and. (find_work .or. (k > nk_linear)) ) then
1451		Rho_avg = ( ((h(i,j,k) + h(i,j,k-1)) + (h(i,j+1,k) + h(i,j+1,k-1))) + 4.0*hn_2 ) / &
1452		( (((h(i,j,k)+hn_2) * tv%SpV_avg(i,j,k)) + ((h(i,j,k-1)+hn_2) * tv%SpV_avg(i,j,k-1))) + &
1453	0	(((h(i,j+1,k)+hn_2)tv%SpV_avg(i,j+1,k)) + ((h(i,j+1,k-1)+hn_2)tv%SpV_avg(i,j+1,k-1))) )
1454		! Use an average density to convert the volume streamfunction estimate into a mass streamfunction.
1455	0	Z_to_H = GV%RZ_to_H*Rho_avg
1456		else
1457	6500160	Z_to_H = GV%Z_to_H
1458		endif
1459
1460	6500160	if (k > nk_linear) then
1461	6500160	if (use_EOS) then
1462
1463	6500160	if (vhtot(i,J) <= 0.0) then
1464		! The transport that must balance the transport below is positive.
1465	6500160	Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j,k))
1466		else ! (vhtot(I,j) > 0.0)
1467	0	Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j+1,k))
1468		endif
1469
1470		! Find the actual streamfunction at each interface.
1471	6500160	Sfn_est = (Z_to_HSfn_unlim_v(i,K) + slope2_Ratio_v(i,K)Sfn_safe) / (1.0 + slope2_Ratio_v(i,K))
1472		else ! When use_EOS is false, the layers are constant density.
1473	0	Sfn_est = Z_to_H*Sfn_unlim_v(i,K)
1474		endif
1475
1476		! Make sure that there is enough mass above to allow the streamfunction
1477		! to satisfy the boundary condition of 0 at the surface.
1478	6500160	Sfn_in_H = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i,j+1,K))
1479
1480		! The actual transport is limited by the mass available in the two
1481		! neighboring grid cells.
1482	6500160	vhD(i,J,k) = max(min((Sfn_in_H - vhtot(i,J)), h_avail(i,j,k)), -h_avail(i,j+1,k))
1483
1484	6500160	if (CS%id_sfn_y>0) diag_sfn_y(i,J,K) = diag_sfn_y(i,J,K+1) + vhD(i,J,k)
1485		! sfn_y(i,J,K) = max(min(Sfn_in_h, vhtot(i,J)+h_avail(i,j,k)), &
1486		! vhtot(i,J)-h_avail(i,j+1,k))
1487		! sfn_slope_y(i,J,K) = max(vhtot(i,J)-h_avail(i,j+1,k), &
1488		! min(vhtot(i,J)+h_avail(i,j,k), &
1489		! min(h_avail_rsum(i,j+1,K), max(-h_avail_rsum(i,j,K), &
1490		! (KH_v(i,J,K)G%dx_Cv(i,J)) ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) )) ))
1491		else ! k <= nk_linear
1492		! Balance the deeper flow with a return flow uniformly distributed
1493		! though the remaining near-surface layers. This is the same as
1494		! using Sfn_safe above. There is no need to apply the limiters in
1495		! this case.
1496	0	if (vhtot(i,J) <= 0.0) then
1497	0	vhD(i,J,k) = -vhtot(i,J) * h_frac(i,j,k)
1498		else ! (vhtot(i,J) > 0.0)
1499	0	vhD(i,J,k) = -vhtot(i,J) * h_frac(i,j+1,k)
1500		endif
1501
1502	0	if (CS%id_sfn_y>0) diag_sfn_y(i,J,K) = diag_sfn_y(i,J,K+1) + vhD(i,J,k)
1503		! if (sfn_slope_y(i,J,K+1) <= 0.0) then
1504		! sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j,k))
1505		! else
1506		! sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j+1,k))
1507		! endif
1508		endif
1509
1510	6500160	vhtot(i,J) = vhtot(i,J) + vhD(i,J,k)
1511
1512	6554328	if (find_work) then
1513		! This is the energy tendency based on the original profiles, and does
1514		! not include any nonlinear terms due to a finite time step (which would
1515		! involve interactions between the fluxes through the different faces.
1516		! A second order centered estimate is used for the density transferred
1517		! between water columns.
1518
1519	6500160	if (allocated(tv%SpV_avg)) then
1520	0	G_scale = GV%H_to_RZ * GV%g_Earth / Rho_avg
1521		else
1522	6500160	G_scale = GV%g_Earth * GV%H_to_Z
1523		endif
1524
1525		Work_v(i,J) = Work_v(i,J) + G_scale * &
1526		( vhtot(i,J) * drdkDe_v(i,K) - &
1527		(vhD(i,J,k) * drdj_v(i,k)) * 0.25 * &
1528	6500160	((e(i,j,K) + e(i,j,K+1)) + (e(i,j+1,K) + e(i,j+1,K+1))) )
1529		endif
1530
1531		enddo
1532		enddo ! end of k-loop
1533		enddo ! end of j-loop
1534
1535		! In layer 1, enforce the boundary conditions that Sfn(z=0) = 0.0
1536	12	if (.not.find_work .or. .not.(use_EOS)) then
1537	0	do j=js,je ; do I=is-1,ie ; uhD(I,j,1) = -uhtot(I,j) ; enddo ; enddo
1538	0	do J=js-1,je ; do i=is,ie ; vhD(i,J,1) = -vhtot(i,J) ; enddo ; enddo
1539		else
1540	12	EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1)
1541		!$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB,G_scale)
1542	732	do j=js,je
1543	720	if (use_EOS) then
1544	87840	do I=is-1,ie
1545	87120	pres_u(I) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
1546	87120	T_u(I) = 0.5*(T(i,j,1) + T(i+1,j,1))
1547	87840	S_u(I) = 0.5*(S(i,j,1) + S(i+1,j,1))
1548		enddo
1549		call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, &
1550	720	tv%eqn_of_state, EOSdom_u )
1551		endif
1552	87852	do I=is-1,ie
1553	87120	uhD(I,j,1) = -uhtot(I,j)
1554
1555	87120	G_scale = GV%g_Earth * GV%H_to_Z
1556	87120	if (use_EOS) then
1557		drdiB = drho_dT_u(I) * (T(i+1,j,1)-T(i,j,1)) + &
1558	87120	drho_dS_u(I) * (S(i+1,j,1)-S(i,j,1))
1559	87120	if (allocated(tv%SpV_avg)) then
1560		G_scale = GV%H_to_RZ * GV%g_Earth * &
1561		( ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1)) + ((h(i+1,j,1)+hn_2) * tv%SpV_avg(i+1,j,1)) ) / &
1562	0	( (h(i,j,1) + h(i+1,j,1)) + 2.0*hn_2 ) )
1563		endif
1564		endif
1565	87840	if (CS%use_GM_work_bug) then
1566		Work_u(I,j) = Work_u(I,j) + G_scale * &
1567		( (uhD(I,j,1) * drdiB) * 0.25 * &
1568	0	((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1569		else
1570		Work_u(I,j) = Work_u(I,j) - G_scale * &
1571		( (uhD(I,j,1) * drdiB) * 0.25 * &
1572	87120	((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1573		endif
1574		enddo
1575		enddo
1576
1577	12	EOSdom_v(:) = EOS_domain(G%HI)
1578		!$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB,G_scale)
1579	744	do J=js-1,je
1580	732	if (use_EOS) then
1581	88572	do i=is,ie
1582	87840	pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
1583	87840	T_v(i) = 0.5*(T(i,j,1) + T(i,j+1,1))
1584	88572	S_v(i) = 0.5*(S(i,j,1) + S(i,j+1,1))
1585		enddo
1586		call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, &
1587	732	tv%eqn_of_state, EOSdom_v)
1588		endif
1589	88584	do i=is,ie
1590	87840	vhD(i,J,1) = -vhtot(i,J)
1591
1592	87840	G_scale = GV%g_Earth * GV%H_to_Z
1593	87840	if (use_EOS) then
1594		drdjB = drho_dT_v(i) * (T(i,j+1,1)-T(i,j,1)) + &
1595	87840	drho_dS_v(i) * (S(i,j+1,1)-S(i,j,1))
1596	87840	if (allocated(tv%SpV_avg)) then
1597		G_scale = GV%H_to_RZ * GV%g_Earth * &
1598		( ( ((h(i,j,1)+hn_2) * tv%SpV_avg(i,j,1)) + ((h(i,j+1,1)+hn_2) * tv%SpV_avg(i,j+1,1)) ) / &
1599	0	( (h(i,j,1) + h(i,j+1,1)) + 2.0*hn_2 ) )
1600		endif
1601		endif
1602		Work_v(i,J) = Work_v(i,J) - G_scale * &
1603		( (vhD(i,J,1) * drdjB) * 0.25 * &
1604	88572	((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
1605		enddo
1606		enddo
1607		endif
1608
1609	87132	if (find_work) then ; do j=js,je ; do i=is,ie
1610		! Note that the units of Work_v and Work_u are [R Z L4 T-3 ~> W], while Work_h is in [R Z L2 T-3 ~> W m-2].
1611		Work_h = 0.5 * G%IareaT(i,j) * &
1612	86400	((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J)))
1613	86400	if (allocated(CS%GMwork)) CS%GMwork(i,j) = Work_h
1614	86400	if (.not. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then
1615	86400	MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
1616		endif ; endif
1617	87120	if (skeb_use_gm) then
1618	0	skeb_gm_work(i,j) = STOCH%skeb_gm_coef * Work_h
1619	0	skeb_ebt_norm2(i,j) = 0.0
1620	0	do k=1,nz
1621	0	skeb_ebt_norm2(i,j) = skeb_ebt_norm2(i,j) + h(i,j,k) * VarMix%ebt_struct(i,j,k)**2
1622		enddo
1623	0	skeb_ebt_norm2(i,j) = GV%H_to_RZ * (skeb_ebt_norm2(i,j) + h_neglect)
1624		endif
1625		enddo ; enddo ; endif
1626
1627	12	if (skeb_use_gm) then
1628		! This block spreads the GM work down through the column using the ebt vertical structure, squared.
1629		! Note the sign convention.
1630	0	do k=1,nz ; do j=js,je ; do i=is,ie
1631		STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - skeb_gm_work(i,j) * &
1632	0	VarMix%ebt_struct(i,j,k)**2 / skeb_ebt_norm2(i,j)
1633		enddo ; enddo ; enddo
1634		endif
1635
1636	12	if (find_work .and. CS%GM_src_alt) then ; if (allocated(MEKE%GM_src)) then
1637	0	if (CS%MEKE_src_answer_date >= 20240601) then
1638	0	do j=js,je ; do i=is,ie ; do k=nz,1,-1
1639		PE_release_h = -0.25 * GV%H_to_RZ * &
1640		( ((KH_u(I,j,k)(Slope_x_PE(I,j,k)2) hN2_x_PE(I,j,k)) + &
1641		(Kh_u(I-1,j,k)(Slope_x_PE(I-1,j,k)2) hN2_x_PE(I-1,j,k))) + &
1642		((Kh_v(i,J,k)(Slope_y_PE(i,J,k)2) hN2_y_PE(i,J,k)) + &
1643	0	(Kh_v(i,J-1,k)(Slope_y_PE(i,J-1,k)2) hN2_y_PE(i,J-1,k))) )
1644	0	MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
1645		enddo ; enddo ; enddo
1646		else
1647	0	do j=js,je ; do i=is,ie ; do k=nz,1,-1
1648		PE_release_h = -0.25 * GV%H_to_RZ * &
1649		((KH_u(I,j,k)(Slope_x_PE(I,j,k)2) hN2_x_PE(I,j,k)) + &
1650		(Kh_u(I-1,j,k)(Slope_x_PE(I-1,j,k)2) hN2_x_PE(I-1,j,k)) + &
1651		(Kh_v(i,J,k)(Slope_y_PE(i,J,k)2) hN2_y_PE(i,J,k)) + &
1652	0	(Kh_v(i,J-1,k)(Slope_y_PE(i,J-1,k)2) hN2_y_PE(i,J-1,k)))
1653	0	MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + PE_release_h
1654		enddo ; enddo ; enddo
1655		endif
1656
1657	0	if (CS%debug) then
1658	0	call hchksum(MEKE%GM_src, 'MEKE%GM_src', G%HI, unscale=US%RZ3_T3_to_W_m2US%L_to_Z*2)
1659		call uvchksum("KH_[uv]", Kh_u, Kh_v, G%HI, unscale=US%L_to_m*2US%s_to_T, &
1660	0	scalar_pair=.true.)
1661	0	call uvchksum("Slope_[xy]_PE", Slope_x_PE, Slope_y_PE, G%HI, unscale=US%Z_to_L)
1662		call uvchksum("hN2_[xy]_PE", hN2_x_PE, hN2_y_PE, G%HI, unscale=GV%H_to_mksUS%L_to_Z2US%s_to_T**2, &
1663	0	scalar_pair=.true.)
1664		endif
1665		endif ; endif
1666
1667	12	if (CS%id_slope_x > 0) call post_data(CS%id_slope_x, CS%diagSlopeX, CS%diag)
1668	12	if (CS%id_slope_y > 0) call post_data(CS%id_slope_y, CS%diagSlopeY, CS%diag)
1669	12	if (CS%id_sfn_x > 0) call post_data(CS%id_sfn_x, diag_sfn_x, CS%diag)
1670	12	if (CS%id_sfn_y > 0) call post_data(CS%id_sfn_y, diag_sfn_y, CS%diag)
1671	12	if (CS%id_sfn_unlim_x > 0) call post_data(CS%id_sfn_unlim_x, diag_sfn_unlim_x, CS%diag)
1672	12	if (CS%id_sfn_unlim_y > 0) call post_data(CS%id_sfn_unlim_y, diag_sfn_unlim_y, CS%diag)
1673
1674	36	end subroutine thickness_diffuse_full
1675
1676		!> Tridiagonal solver for streamfunction at interfaces
1677	0	subroutine streamfn_solver(nk, c2_h, hN2, sfn)
1678		integer, intent(in) :: nk !< Number of layers
1679		real, dimension(nk), intent(in) :: c2_h !< Wave speed squared over thickness in layers, rescaled to
1680		!! [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2]
1681		real, dimension(nk+1), intent(in) :: hN2 !< Thickness times N2 at interfaces times rescaling factors
1682		!! [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2]
1683		real, dimension(nk+1), intent(inout) :: sfn !< Streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] or arbitrary units
1684		!! On entry, equals diffusivity times slope.
1685		!! On exit, equals the streamfunction.
1686		! Local variables
1687	0	real :: c1(nk) ! The dependence of the final streamfunction on the values below [nondim]
1688		real :: d1 ! The complement of c1(k) (i.e., 1 - c1(k)) [nondim]
1689		real :: b_denom ! A term in the denominator of beta [H L2 Z-2 T-2 ~> m s-2 or kg m-2 s-2]
1690		real :: beta ! The normalization for the pivot [Z2 T2 H-1 L-2 ~> s2 m-1 or m2 s2 kg-1]
1691		integer :: k
1692
1693	0	sfn(1) = 0.
1694	0	b_denom = hN2(2) + c2_h(1)
1695	0	beta = 1.0 / ( b_denom + c2_h(2) )
1696	0	d1 = beta * b_denom
1697	0	sfn(2) = ( beta * hN2(2) )*sfn(2)
1698	0	do K=3,nk
1699	0	c1(k-1) = beta * c2_h(k-1)
1700	0	b_denom = hN2(K) + d1*c2_h(k-1)
1701	0	beta = 1.0 / (b_denom + c2_h(k))
1702	0	d1 = beta * b_denom
1703	0	sfn(K) = beta * (hN2(K)sfn(K) + c2_h(k-1)sfn(K-1))
1704		enddo
1705	0	c1(nk) = beta * c2_h(nk)
1706	0	sfn(nk+1) = 0.
1707	0	do K=nk,2,-1
1708	0	sfn(K) = sfn(K) + c1(k)*sfn(K+1)
1709		enddo
1710
1711	0	end subroutine streamfn_solver
1712
1713		!> Add a diffusivity that acts on the isopycnal heights, regardless of the densities
1714	0	subroutine add_interface_Kh(G, GV, US, CS, Kh_u, Kh_v, Kh_u_CFL, Kh_v_CFL, int_slope_u, int_slope_v)
1715		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1716		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1717		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1718		type(thickness_diffuse_CS), intent(in) :: CS !< Control structure for thickness_diffuse
1719		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kh_u !< Isopycnal height diffusivity
1720		!! at u points [L2 T-1 ~> m2 s-1]
1721		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: Kh_v !< Isopycnal height diffusivity
1722		!! at v points [L2 T-1 ~> m2 s-1]
1723		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Kh_u_CFL !< Maximum stable isopycnal height
1724		!! diffusivity at u points [L2 T-1 ~> m2 s-1]
1725		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Kh_v_CFL !< Maximum stable isopycnal height
1726		!! diffusivity at v points [L2 T-1 ~> m2 s-1]
1727		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: int_slope_u !< Ratio that determine how much of
1728		!! the isopycnal slopes are taken directly from
1729		!! the interface slopes without consideration
1730		!! of density gradients [nondim].
1731		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: int_slope_v !< Ratio that determine how much of
1732		!! the isopycnal slopes are taken directly from
1733		!! the interface slopes without consideration
1734		!! of density gradients [nondim].
1735
1736		! Local variables
1737		integer :: i, j, k, is, ie, js, je, nz
1738
1739	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1740
1741	0	do k=1,nz+1 ; do j=js,je ; do I=is-1,ie ; if (CS%Kh_eta_u(I,j) > 0.0) then
1742		int_slope_u(I,j,K) = (int_slope_u(I,j,K)*Kh_u(I,j,K) + CS%Kh_eta_u(I,j)) / &
1743	0	(Kh_u(I,j,K) + CS%Kh_eta_u(I,j))
1744	0	Kh_u(I,j,K) = min(Kh_u(I,j,K) + CS%Kh_eta_u(I,j), Kh_u_CFL(I,j))
1745		endif ; enddo ; enddo ; enddo
1746
1747	0	do k=1,nz+1 ; do J=js-1,je ; do i=is,ie ; if (CS%Kh_eta_v(i,J) > 0.0) then
1748		int_slope_v(i,J,K) = (int_slope_v(i,J,K)*Kh_v(i,J,K) + CS%Kh_eta_v(i,J)) / &
1749	0	(Kh_v(i,J,K) + CS%Kh_eta_v(i,J))
1750	0	Kh_v(i,J,K) = min(Kh_v(i,J,K) + CS%Kh_eta_v(i,J), Kh_v_CFL(i,J))
1751		endif ; enddo ; enddo ; enddo
1752
1753	0	end subroutine add_interface_Kh
1754
1755		!> Modifies isopycnal height diffusivities to untangle layer structures
1756	0	subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, CS, &
1757	0	int_slope_u, int_slope_v)
1758		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1759		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1760		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1761		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1762		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface positions [Z ~> m]
1763		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kh_u !< Isopycnal height diffusivity
1764		!! at u points [L2 T-1 ~> m2 s-1]
1765		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: Kh_v !< Isopycnal height diffusivity
1766		!! at v points [L2 T-1 ~> m2 s-1]
1767		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Kh_u_CFL !< Maximum stable isopycnal height
1768		!! diffusivity at u points [L2 T-1 ~> m2 s-1]
1769		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Kh_v_CFL !< Maximum stable isopycnal height
1770		!! diffusivity at v points [L2 T-1 ~> m2 s-1]
1771		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1772		real, intent(in) :: dt !< Time increment [T ~> s]
1773		type(thickness_diffuse_CS), intent(in) :: CS !< Control structure for thickness_diffuse
1774		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: int_slope_u !< Ratio that determine how much of
1775		!! the isopycnal slopes are taken directly from
1776		!! the interface slopes without consideration
1777		!! of density gradients [nondim].
1778		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: int_slope_v !< Ratio that determine how much of
1779		!! the isopycnal slopes are taken directly from
1780		!! the interface slopes without consideration
1781		!! of density gradients [nondim].
1782		! Local variables
1783		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1784	0	de_top ! The distances between the top of a layer and the top of the
1785		! region where the detangling is applied [H ~> m or kg m-2].
1786		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
1787	0	Kh_lay_u ! The tentative isopycnal height diffusivity for each layer at
1788		! u points [L2 T-1 ~> m2 s-1].
1789		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
1790	0	Kh_lay_v ! The tentative isopycnal height diffusivity for each layer at
1791		! v points [L2 T-1 ~> m2 s-1].
1792		real, dimension(SZI_(G),SZJ_(G)) :: &
1793	0	de_bot ! The distances from the bottom of the region where the
1794		! detangling is applied [H ~> m or kg m-2].
1795		real :: h1, h2 ! The thinner and thicker surrounding thicknesses [H ~> m or kg m-2],
1796		! with the thinner modified near the boundaries to mask out
1797		! thickness variations due to topography, etc.
1798		real :: jag_Rat ! The nondimensional jaggedness ratio for a layer, going
1799		! from 0 (smooth) to 1 (jagged) [nondim]. This is the difference
1800		! between the arithmetic and harmonic mean thicknesses
1801		! normalized by the arithmetic mean thickness.
1802		real :: Kh_scale ! A ratio by which Kh_u_CFL is scaled for maximally jagged
1803		! layers [nondim].
1804		real :: h_neglect ! A thickness that is so small it is usually lost
1805		! in roundoff and can be neglected [H ~> m or kg m-2].
1806
1807		real :: I_sl ! The absolute value of the larger in magnitude of the slopes
1808		! above and below [L Z-1 ~> nondim].
1809		real :: Rsl ! The ratio of the smaller magnitude slope to the larger
1810		! magnitude one [nondim]. 0 <= Rsl <1.
1811		real :: IRsl ! The (limited) inverse of Rsl [nondim]. 1 < IRsl <= 1e9.
1812		real :: dH ! The thickness gradient divided by the damping timescale
1813		! and the ratio of the face length to the adjacent cell
1814		! areas for comparability with the diffusivities [L Z T-1 ~> m2 s-1].
1815		real :: adH ! The absolute value of dH [L Z T-1 ~> m2 s-1].
1816		real :: sign ! 1 or -1, with the same sign as the layer thickness gradient [nondim].
1817		real :: sl_K ! The sign-corrected slope of the interface above [Z L-1 ~> nondim].
1818		real :: sl_Kp1 ! The sign-corrected slope of the interface below [Z L-1 ~> nondim].
1819		real :: I_sl_K ! The (limited) inverse of sl_K [L Z-1 ~> nondim].
1820		real :: I_sl_Kp1 ! The (limited) inverse of sl_Kp1 [L Z-1 ~> nondim].
1821		real :: I_4t ! A quarter of a flux scaling factor divided by
1822		! the damping timescale [T-1 ~> s-1].
1823		real :: Fn_R ! A function of Rsl, such that Rsl < Fn_R < 1 [nondim]
1824		real :: Idx_eff ! The effective inverse x-grid spacing at a u-point [L-1 ~> m-1]
1825		real :: Idy_eff ! The effective inverse y-grid spacing at a v-point [L-1 ~> m-1]
1826		real :: slope_sq ! The sum of the squared slopes above and below a layer [Z2 L-2 ~> nondim]
1827		real :: Kh_max ! A local ceiling on the diffusivity [L2 T-1 ~> m2 s-1].
1828		real :: wt1, wt2 ! Nondimensional weights [nondim].
1829		! Variables used only in testing code.
1830		! real, dimension(SZK_(GV)) :: uh_here ! The transport in a layer [Z L2 T-1 ~> m3 s-1]
1831		! real, dimension(SZK_(GV)+1) :: Sfn ! The streamfunction at an interface [Z L T-1 ~> m2 s-1]
1832		real :: dKh ! An increment in the diffusivity [L2 T-1 ~> m2 s-1].
1833
1834		real, dimension(SZIB_(G),SZK_(GV)+1) :: &
1835	0	Kh_bg, & ! The background (floor) value of Kh [L2 T-1 ~> m2 s-1].
1836	0	Kh, & ! The tentative value of Kh [L2 T-1 ~> m2 s-1].
1837	0	Kh_detangle, & ! The detangling diffusivity that could be used [L2 T-1 ~> m2 s-1].
1838	0	Kh_min_max_p, & ! The smallest ceiling that can be placed on Kh(I,K)
1839		! based on the value of Kh(I,K+1) [L2 T-1 ~> m2 s-1].
1840	0	Kh_min_max_m, & ! The smallest ceiling that can be placed on Kh(I,K)
1841		! based on the value of Kh(I,K-1) [L2 T-1 ~> m2 s-1].
1842		! The following are variables that define the relationships between
1843		! successive values of Kh.
1844		! Search for Kh that satisfy...
1845		! Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
1846		! Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
1847		! Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
1848		! Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
1849	0	Kh_min_m , & ! See above [nondim].
1850	0	Kh0_min_m , & ! See above [L2 T-1 ~> m2 s-1].
1851	0	Kh_max_m , & ! See above [nondim].
1852	0	Kh0_max_m, & ! See above [L2 T-1 ~> m2 s-1].
1853	0	Kh_min_p , & ! See above [nondim].
1854	0	Kh0_min_p , & ! See above [L2 T-1 ~> m2 s-1].
1855	0	Kh_max_p , & ! See above [nondim].
1856	0	Kh0_max_p ! See above [L2 T-1 ~> m2 s-1].
1857		real, dimension(SZIB_(G)) :: &
1858	0	Kh_max_max ! The maximum diffusivity permitted in a column [L2 T-1 ~> m2 s-1]
1859		logical, dimension(SZIB_(G)) :: &
1860	0	do_i ! If true, work on a column.
1861		integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
1862	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1863
1864	0	k_top = GV%nk_rho_varies + 1
1865	0	h_neglect = GV%H_subroundoff
1866		! The 0.5 is because we are not using uniform weightings, but are
1867		! distributing the diffusivities more effectively (with wt1 & wt2), but this
1868		! means that the additions to a single interface can be up to twice as large.
1869	0	Kh_scale = 0.5
1870	0	if (CS%detangle_time > dt) Kh_scale = 0.5 * dt / CS%detangle_time
1871
1872	0	do j=js-1,je+1 ; do i=is-1,ie+1
1873	0	de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
1874		enddo ; enddo
1875	0	do k=k_top+1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
1876	0	de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
1877		enddo ; enddo ; enddo
1878
1879	0	do j=js,je ; do I=is-1,ie
1880	0	Kh_lay_u(I,j,nz) = 0.0 ; Kh_lay_u(I,j,k_top) = 0.0
1881		enddo ; enddo
1882	0	do J=js-1,je ; do i=is,ie
1883	0	Kh_lay_v(i,J,nz) = 0.0 ; Kh_lay_v(i,J,k_top) = 0.0
1884		enddo ; enddo
1885
1886	0	do k=nz-1,k_top+1,-1
1887		! Find the diffusivities associated with each layer.
1888	0	do j=js-1,je+1 ; do i=is-1,ie+1
1889	0	de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
1890		enddo ; enddo
1891
1892	0	do j=js,je ; do I=is-1,ie ; if (G%OBCmaskCu(I,j) > 0.0) then
1893	0	if (h(i,j,k) > h(i+1,j,k)) then
1894	0	h2 = h(i,j,k)
1895	0	h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
1896		else
1897	0	h2 = h(i+1,j,k)
1898	0	h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1899		endif
1900	0	jag_Rat = (h2 - h1)2 / (h2 + h1 + h_neglect)2
1901	0	KH_lay_u(I,j,k) = (Kh_scale * KH_u_CFL(I,j)) * jag_Rat**2
1902		endif ; enddo ; enddo
1903
1904	0	do J=js-1,je ; do i=is,ie ; if (G%OBCmaskCv(i,J) > 0.0) then
1905	0	if (h(i,j,k) > h(i,j+1,k)) then
1906	0	h2 = h(i,j,k)
1907	0	h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
1908		else
1909	0	h2 = h(i,j+1,k)
1910	0	h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1911		endif
1912	0	jag_Rat = (h2 - h1)2 / (h2 + h1 + h_neglect)2
1913	0	KH_lay_v(i,J,k) = (Kh_scale * KH_v_CFL(i,J)) * jag_Rat**2
1914		endif ; enddo ; enddo
1915		enddo
1916
1917		! Limit the diffusivities
1918
1919	0	I_4t = Kh_scale / (4.0 * dt)
1920
1921	0	do n=1,2
1922	0	if (n==1) then ; jsh = js ; ish = is-1
1923	0	else ; jsh = js-1 ; ish = is ; endif
1924
1925	0	do j=jsh,je
1926
1927		! First, populate the diffusivities
1928	0	if (n==1) then ! This is a u-column.
1929	0	do i=ish,ie
1930	0	do_i(I) = (G%OBCmaskCu(I,j) > 0.0)
1931	0	Kh_Max_max(I) = KH_u_CFL(I,j)
1932		enddo
1933	0	do K=1,nz+1 ; do i=ish,ie
1934	0	Kh_bg(I,K) = KH_u(I,j,K) ; Kh(I,K) = Kh_bg(I,K)
1935	0	Kh_min_max_p(I,K) = Kh_bg(I,K) ; Kh_min_max_m(I,K) = Kh_bg(I,K)
1936	0	Kh_detangle(I,K) = 0.0
1937		enddo ; enddo
1938		else ! This is a v-column.
1939	0	do i=ish,ie
1940	0	do_i(i) = (G%OBCmaskCv(i,J) > 0.0) ; Kh_Max_max(I) = KH_v_CFL(i,J)
1941		enddo
1942	0	do K=1,nz+1 ; do i=ish,ie
1943	0	Kh_bg(I,K) = KH_v(I,j,K) ; Kh(I,K) = Kh_bg(I,K)
1944	0	Kh_min_max_p(I,K) = Kh_bg(I,K) ; Kh_min_max_m(I,K) = Kh_bg(I,K)
1945	0	Kh_detangle(I,K) = 0.0
1946		enddo ; enddo
1947		endif
1948
1949		! Determine the limits on the diffusivities.
1950	0	do k=k_top,nz ; do i=ish,ie ; if (do_i(i)) then
1951	0	if (n==1) then ! This is a u-column.
1952	0	dH = 0.0
1953	0	Idx_eff = ((G%IareaT(i+1,j) + G%IareaT(i,j)) * G%dy_Cu(I,j))
1954		! This expression uses differences in e in place of h for better
1955		! consistency with the slopes.
1956	0	if (Idx_eff > 0.0) &
1957		dH = I_4t * ((e(i+1,j,K) - e(i+1,j,K+1)) - &
1958	0	(e(i,j,K) - e(i,j,K+1))) / Idx_eff
1959		! dH = I_4t * (h(i+1,j,k) - h(i,j,k)) / Idx_eff
1960
1961	0	adH = abs(dH)
1962	0	sign = 1.0 ; if (dH < 0) sign = -1.0
1963	0	sl_K = sign * (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
1964	0	sl_Kp1 = sign * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IdxCu(I,j)
1965
1966		! Add the incremental diffusivities to the surrounding interfaces.
1967		! Adding more to the more steeply sloping layers (as below) makes
1968		! the diffusivities more than twice as effective.
1969	0	slope_sq = (sl_K2 + sl_Kp12)
1970	0	wt1 = 0.5 ; wt2 = 0.5
1971	0	if (slope_sq > 0.0) then
1972	0	wt1 = sl_K2 / slope_sq ; wt2 = sl_Kp12 / slope_sq
1973		endif
1974	0	Kh_detangle(I,K) = Kh_detangle(I,K) + wt1*KH_lay_u(I,j,k)
1975	0	Kh_detangle(I,K+1) = Kh_detangle(I,K+1) + wt2*KH_lay_u(I,j,k)
1976		else ! This is a v-column.
1977	0	dH = 0.0
1978	0	Idy_eff = ((G%IareaT(i,j+1) + G%IareaT(i,j)) * G%dx_Cv(I,j))
1979	0	if (Idy_eff > 0.0) &
1980		dH = I_4t * ((e(i,j+1,K) - e(i,j+1,K+1)) - &
1981	0	(e(i,j,K) - e(i,j,K+1))) / Idy_eff
1982		! dH = I_4t * (h(i,j+1,k) - h(i,j,k)) / Idy_eff
1983
1984	0	adH = abs(dH)
1985	0	sign = 1.0 ; if (dH < 0) sign = -1.0
1986	0	sl_K = sign * (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)
1987	0	sl_Kp1 = sign * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IdyCv(i,J)
1988
1989		! Add the incremental diffusivities to the surrounding interfaces.
1990		! Adding more to the more steeply sloping layers (as below) makes
1991		! the diffusivities more than twice as effective.
1992	0	slope_sq = (sl_K2 + sl_Kp12)
1993	0	wt1 = 0.5 ; wt2 = 0.5
1994	0	if (slope_sq > 0.0) then
1995	0	wt1 = sl_K2 / slope_sq ; wt2 = sl_Kp12 / slope_sq
1996		endif
1997	0	Kh_detangle(I,K) = Kh_detangle(I,K) + wt1*KH_lay_v(i,J,k)
1998	0	Kh_detangle(I,K+1) = Kh_detangle(I,K+1) + wt2*KH_lay_v(i,J,k)
1999		endif
2000
2001	0	if (adH == 0.0) then
2002	0	Kh_min_m(I,K+1) = 1.0 ; Kh0_min_m(I,K+1) = 0.0
2003	0	Kh_max_m(I,K+1) = 1.0 ; Kh0_max_m(I,K+1) = 0.0
2004	0	Kh_min_p(I,K) = 1.0 ; Kh0_min_p(I,K) = 0.0
2005	0	Kh_max_p(I,K) = 1.0 ; Kh0_max_p(I,K) = 0.0
2006	0	elseif (adH > 0.0) then
2007	0	if (sl_K <= sl_Kp1) then
2008		! This case should only arise from nonlinearities in the equation of state.
2009		! Treat it as though dedx(K) = dedx(K+1) & dH = 0.
2010	0	Kh_min_m(I,K+1) = 1.0 ; Kh0_min_m(I,K+1) = 0.0
2011	0	Kh_max_m(I,K+1) = 1.0 ; Kh0_max_m(I,K+1) = 0.0
2012	0	Kh_min_p(I,K) = 1.0 ; Kh0_min_p(I,K) = 0.0
2013	0	Kh_max_p(I,K) = 1.0 ; Kh0_max_p(I,K) = 0.0
2014	0	elseif (sl_K <= 0.0) then ! Both slopes are opposite to dH
2015	0	I_sl = -1.0 / sl_Kp1
2016	0	Rsl = -sl_K * I_sl ! 0 <= Rsl < 1
2017	0	IRsl = 1e9 ; if (Rsl > 1e-9) IRsl = 1.0/Rsl ! 1 < IRsl <= 1e9
2018
2019	0	Fn_R = Rsl
2020	0	if (Kh_max_max(I) > 0) &
2021	0	Fn_R = min(sqrt(Rsl), Rsl + (adH * I_sl) / (Kh_Max_max(I)))
2022
2023	0	Kh_min_m(I,K+1) = Fn_R ; Kh0_min_m(I,K+1) = 0.0
2024	0	Kh_max_m(I,K+1) = Rsl ; Kh0_max_m(I,K+1) = adH * I_sl
2025	0	Kh_min_p(I,K) = IRsl ; Kh0_min_p(I,K) = -adH * (I_sl*IRsl)
2026	0	Kh_max_p(I,K) = 1.0/(Fn_R + 1.0e-30) ; Kh0_max_p(I,K) = 0.0
2027	0	elseif (sl_Kp1 < 0.0) then ! Opposite (nonzero) signs of slopes.
2028	0	I_sl_K = 1e18US%Z_to_L ; if (sl_K > 1e-18US%L_to_Z) I_sl_K = 1.0 / sl_K
2029	0	I_sl_Kp1 = 1e18US%Z_to_L ; if (-sl_Kp1 > 1e-18US%L_to_Z) I_sl_Kp1 = -1.0 / sl_Kp1
2030
2031	0	Kh_min_m(I,K+1) = 0.0 ; Kh0_min_m(I,K+1) = 0.0
2032	0	Kh_max_m(I,K+1) = - sl_KI_sl_Kp1 ; Kh0_max_m(I,K+1) = adHI_sl_Kp1
2033	0	Kh_min_p(I,K) = 0.0 ; Kh0_min_p(I,K) = 0.0
2034	0	Kh_max_p(I,K) = sl_Kp1I_sl_K ; Kh0_max_p(I,K) = adHI_sl_K
2035
2036		! This limit does not use the slope weighting so that potentially
2037		! sharp gradients in diffusivities are not forced to occur.
2038	0	Kh_Max = adH / (sl_K - sl_Kp1)
2039	0	Kh_min_max_p(I,K) = max(Kh_min_max_p(I,K), Kh_Max)
2040	0	Kh_min_max_m(I,K+1) = max(Kh_min_max_m(I,K+1), Kh_Max)
2041		else ! Both slopes are of the same sign as dH.
2042	0	I_sl = 1.0 / sl_K
2043	0	Rsl = sl_Kp1 * I_sl ! 0 <= Rsl < 1
2044	0	IRsl = 1e9 ; if (Rsl > 1e-9) IRsl = 1.0/Rsl ! 1 < IRsl <= 1e9
2045
2046		! Rsl <= Fn_R <= 1
2047	0	Fn_R = Rsl
2048	0	if (Kh_max_max(I) > 0) &
2049	0	Fn_R = min(sqrt(Rsl), Rsl + (adH * I_sl) / Kh_Max_max(I))
2050
2051	0	Kh_min_m(I,K+1) = IRsl ; Kh0_min_m(I,K+1) = -adH * (I_sl*IRsl)
2052	0	Kh_max_m(I,K+1) = 1.0/(Fn_R + 1.0e-30) ; Kh0_max_m(I,K+1) = 0.0
2053	0	Kh_min_p(I,K) = Fn_R ; Kh0_min_p(I,K) = 0.0
2054	0	Kh_max_p(I,K) = Rsl ; Kh0_max_p(I,K) = adH * I_sl
2055		endif
2056		endif
2057		endif ; enddo ; enddo ! I-loop & k-loop
2058
2059	0	do k=k_top,nz+1,nz+1-k_top ; do i=ish,ie ; if (do_i(i)) then
2060		! The diffusivities at k_top and nz+1 are both fixed.
2061	0	Kh_min_m(I,k) = 0.0 ; Kh0_min_m(I,k) = 0.0
2062	0	Kh_max_m(I,k) = 0.0 ; Kh0_max_m(I,k) = 0.0
2063	0	Kh_min_p(I,k) = 0.0 ; Kh0_min_p(I,k) = 0.0
2064	0	Kh_max_p(I,k) = 0.0 ; Kh0_max_p(I,k) = 0.0
2065	0	Kh_min_max_p(I,K) = Kh_bg(I,K)
2066	0	Kh_min_max_m(I,K) = Kh_bg(I,K)
2067		endif ; enddo ; enddo ! I-loop and k_top/nz+1 loop
2068
2069		! Search for Kh that satisfy...
2070		! Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
2071		! Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
2072		! Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
2073		! Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
2074
2075		! Increase the diffusivities to satisfy the min constraints.
2076		! All non-zero min constraints on one diffusivity are max constraints on another.
2077	0	do K=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
2078		Kh(I,K) = max(Kh_bg(I,K), Kh_detangle(I,K), &
2079	0	min(Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K), Kh(I,K-1)))
2080
2081	0	if (Kh0_max_m(I,K) > Kh_bg(I,K)) Kh(I,K) = min(Kh(I,K), Kh0_max_m(I,K))
2082	0	if (Kh0_max_p(I,K) > Kh_bg(I,K)) Kh(I,K) = min(Kh(I,K), Kh0_max_p(I,K))
2083		endif ; enddo ; enddo ! I-loop & k-loop
2084		! This is still true... do i=ish,ie ; Kh(I,nz+1) = Kh_bg(I,nz+1) ; enddo
2085	0	do K=nz,k_top+1,-1 ; do i=ish,ie ; if (do_i(i)) then
2086	0	Kh(I,k) = max(Kh(I,K), min(Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K), Kh(I,K+1)))
2087
2088	0	Kh_Max = max(Kh_min_max_p(I,K), Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K))
2089	0	Kh(I,k) = min(Kh(I,k), Kh_Max)
2090		endif ; enddo ; enddo ! I-loop & k-loop
2091		! All non-zero min constraints on one diffusivity are max constraints on
2092		! another layer, so the min constraints can now be discounted.
2093
2094		! Decrease the diffusivities to satisfy the max constraints.
2095	0	do K=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
2096	0	Kh_Max = max(Kh_min_max_m(I,K), Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K))
2097	0	if (Kh(I,k) > Kh_Max) Kh(I,k) = Kh_Max
2098		endif ; enddo ; enddo ! i- and K-loops
2099
2100		! This code tests the solutions...
2101		! do i=ish,ie
2102		! Sfn(:) = 0.0 ; uh_here(:) = 0.0
2103		! do K=k_top,nz
2104		! if ((Kh(i,K) > Kh_bg(i,K)) .or. (Kh(i,K+1) > Kh_bg(i,K+1))) then
2105		! if (n==1) then ! u-point.
2106		! if ((h(i+1,j,k) - h(i,j,k)) * &
2107		! ((e(i+1,j,K)-e(i+1,j,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
2108		! Sfn(K) = -Kh(i,K) * (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
2109		! Sfn(K+1) = -Kh(i,K+1) * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IdxCu(I,j)
2110		! uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dy_Cu(I,j)
2111		! if (abs(uh_here(k)) * min(G%IareaT(i,j), G%IareaT(i+1,j)) > &
2112		! (1e-10*GV%m_to_H)) then
2113		! if (uh_here(k) * (h(i+1,j,k) - h(i,j,k)) > 0.0) then
2114		! call MOM_error(WARNING, "Corrective u-transport is up the thickness gradient.", .true.)
2115		! endif
2116		! if (((h(i,j,k) - 4.0dtG%IareaT(i,j)*uh_here(k)) - &
2117		! (h(i+1,j,k) + 4.0dtG%IareaT(i+1,j)uh_here(k))) &
2118		! (h(i,j,k) - h(i+1,j,k)) < 0.0) then
2119		! call MOM_error(WARNING, "Corrective u-transport is too large.", .true.)
2120		! endif
2121		! endif
2122		! endif
2123		! else ! v-point
2124		! if ((h(i,j+1,k) - h(i,j,k)) * &
2125		! ((e(i,j+1,K)-e(i,j+1,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
2126		! Sfn(K) = -Kh(i,K) * (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)
2127		! Sfn(K+1) = -Kh(i,K+1) * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IdyCv(i,J)
2128		! uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dx_Cv(i,J)
2129		! if (abs(uh_here(K)) * min(G%IareaT(i,j), G%IareaT(i,j+1)) > &
2130		! (1e-10*GV%m_to_H)) then
2131		! if (uh_here(K) * (h(i,j+1,k) - h(i,j,k)) > 0.0) then
2132		! call MOM_error(WARNING, &
2133		! "Corrective v-transport is up the thickness gradient.", .true.)
2134		! endif
2135		! if (((h(i,j,k) - 4.0dtG%IareaT(i,j)*uh_here(K)) - &
2136		! (h(i,j+1,k) + 4.0dtG%IareaT(i,j+1)uh_here(K))) &
2137		! (h(i,j,k) - h(i,j+1,k)) < 0.0) then
2138		! call MOM_error(WARNING, &
2139		! "Corrective v-transport is too large.", .true.)
2140		! endif
2141		! endif
2142		! endif
2143		! endif ! u- or v- selection.
2144		! ! de_dx(I,K) = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
2145		! endif
2146		! enddo
2147		! enddo
2148
2149	0	if (n==1) then ! This is a u-column.
2150	0	do K=k_top+1,nz ; do i=ish,ie
2151	0	if (Kh(I,K) > KH_u(I,j,K)) then
2152	0	dKh = (Kh(I,K) - KH_u(I,j,K))
2153	0	int_slope_u(I,j,K) = dKh / Kh(I,K)
2154	0	KH_u(I,j,K) = Kh(I,K)
2155		endif
2156		enddo ; enddo
2157		else ! This is a v-column.
2158	0	do K=k_top+1,nz ; do i=ish,ie
2159	0	if (Kh(i,K) > KH_v(i,J,K)) then
2160	0	dKh = Kh(i,K) - KH_v(i,J,K)
2161	0	int_slope_v(i,J,K) = dKh / Kh(i,K)
2162	0	KH_v(i,J,K) = Kh(i,K)
2163		endif
2164		enddo ; enddo
2165		endif
2166
2167		enddo ! j-loop
2168		enddo ! n-loop over u- and v- directions.
2169
2170	0	end subroutine add_detangling_Kh
2171
2172		!> Initialize the isopycnal height diffusion module and its control structure
2173	1	subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
2174		type(time_type), intent(in) :: Time !< Current model time
2175		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2176		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2177		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2178		type(param_file_type), intent(in) :: param_file !< Parameter file handles
2179		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
2180		type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity equation diagnostics
2181		type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse
2182
2183		! Local variables
2184		character(len=40) :: mdl = "MOM_thickness_diffuse" ! This module's name.
2185		character(len=200) :: khth_file, inputdir, khth_varname
2186		! This include declares and sets the variable "version".
2187		# include "version_variable.h"
2188		real :: grid_sp ! The local grid spacing [L ~> m]
2189		real :: omega ! The Earth's rotation rate [T-1 ~> s-1]
2190		real :: strat_floor ! A floor for buoyancy frequency in the Ferrari et al. 2010,
2191		! streamfunction formulation, expressed as a fraction of planetary
2192		! rotation divided by an aspect ratio rescaling factor [L Z-1 ~> nondim]
2193		real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
2194		! temperature variance [nondim]
2195		logical :: khth_use_ebt_struct ! If true, uses the equivalent barotropic structure
2196		! as the vertical structure of thickness diffusivity.
2197		! Used to determine if FULL_DEPTH_KHTH_MIN should be
2198		! available.
2199		logical :: use_meke = .false. ! If true, use the MEKE formulation for the thickness diffusivity.
2200		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2201		integer :: i, j
2202
2203	1	CS%initialized = .true.
2204	1	CS%diag => diag
2205
2206		! Read all relevant parameters and write them to the model log.
2207	1	call log_version(param_file, mdl, version, "")
2208		call get_param(param_file, mdl, "THICKNESSDIFFUSE", CS%thickness_diffuse, &
2209		"If true, interface heights are diffused with a "//&
2210	1	"coefficient of KHTH.", default=.false.)
2211		call get_param(param_file, mdl, "KHTH", CS%Khth, &
2212		"The background horizontal thickness diffusivity.", &
2213	1	default=0.0, units="m2 s-1", scale=US%m_to_L*2US%T_to_s)
2214		call get_param(param_file, mdl, "READ_KHTH", CS%read_khth, &
2215		"If true, read a file (given by KHTH_FILE) containing the "//&
2216		"spatially varying horizontal isopycnal height diffusivity.", &
2217	1	default=.false.)
2218	1	if (CS%read_khth) then
2219	0	if (CS%Khth > 0) then
2220		call MOM_error(FATAL, "thickness_diffuse_init: KHTH > 0 is not "// &
2221	0	"compatible with READ_KHTH = TRUE. ")
2222		endif
2223		call get_param(param_file, mdl, "INPUTDIR", inputdir, &
2224		"The directory in which all input files are found.", &
2225	0	default=".", do_not_log=.true.)
2226	0	inputdir = slasher(inputdir)
2227		call get_param(param_file, mdl, "KHTH_FILE", khth_file, &
2228		"The file containing the spatially varying horizontal "//&
2229	0	"isopycnal height diffusivity.", default="khth.nc")
2230		call get_param(param_file, mdl, "KHTH_VARIABLE", khth_varname, &
2231		"The name of the isopycnal height diffusivity variable to read "//&
2232		"from KHTH_FILE.", &
2233	0	default="khth")
2234	0	khth_file = trim(inputdir) // trim(khth_file)
2235
2236	0	allocate(CS%khth2d(G%isd:G%ied, G%jsd:G%jed), source=0.0)
2237	0	call MOM_read_data(khth_file, khth_varname, CS%khth2d(:,:), G%domain, scale=US%m_to_L*2US%T_to_s)
2238	0	call pass_var(CS%khth2d, G%domain)
2239		endif
2240		call get_param(param_file, mdl, "KHTH_SLOPE_CFF", CS%KHTH_Slope_Cff, &
2241		"The nondimensional coefficient in the Visbeck formula for "//&
2242	1	"the interface depth diffusivity", units="nondim", default=0.0)
2243		call get_param(param_file, mdl, "KHTH_MIN", CS%KHTH_Min, &
2244		"The minimum horizontal thickness diffusivity.", &
2245	1	default=0.0, units="m2 s-1", scale=US%m_to_L*2US%T_to_s)
2246		call get_param(param_file, mdl, "KHTH_USE_EBT_STRUCT", khth_use_ebt_struct, &
2247		"If true, uses the equivalent barotropic structure "//&
2248		"as the vertical structure of thickness diffusivity.",&
2249	1	default=.false., do_not_log=.true.)
2250	1	if (khth_use_ebt_struct .and. CS%KHTH_Min>0.0) then
2251		call get_param(param_file, mdl, "FULL_DEPTH_KHTH_MIN", CS%full_depth_khth_min, &
2252		"If true, KHTH_MIN is enforced throughout the whole water column. "//&
2253		"Otherwise, KHTH_MIN is only enforced at the surface. This parameter "//&
2254		"is only available when KHTH_USE_EBT_STRUCT=True and KHTH_MIN>0.", &
2255	0	default=.false.)
2256		endif
2257		call get_param(param_file, mdl, "KHTH_MAX", CS%KHTH_Max, &
2258		"The maximum horizontal thickness diffusivity.", &
2259	1	default=0.0, units="m2 s-1", scale=US%m_to_L*2US%T_to_s)
2260		call get_param(param_file, mdl, "KHTH_MAX_CFL", CS%max_Khth_CFL, &
2261		"The maximum value of the local diffusive CFL ratio that "//&
2262		"is permitted for the thickness diffusivity. 1.0 is the "//&
2263		"marginally unstable value in a pure layered model, but "//&
2264		"much smaller numbers (e.g. 0.1) seem to work better for "//&
2265	1	"ALE-based models.", units="nondimensional", default=0.8)
2266
2267		call get_param(param_file, mdl, "KH_ETA_CONST", CS%Kh_eta_bg, &
2268		"The background horizontal diffusivity of the interface heights (without "//&
2269		"considering the layer density structure). If diffusive CFL limits are "//&
2270		"encountered, the diffusivities of the isopycnals and the interfaces heights "//&
2271		"are scaled back proportionately.", &
2272	1	default=0.0, units="m2 s-1", scale=US%m_to_L*2US%T_to_s)
2273		call get_param(param_file, mdl, "KH_ETA_VEL_SCALE", CS%Kh_eta_vel, &
2274		"A velocity scale that is multiplied by the grid spacing to give a contribution "//&
2275		"to the horizontal diffusivity of the interface heights (without considering "//&
2276		"the layer density structure).", &
2277	1	default=0.0, units="m s-1", scale=US%m_to_L*US%T_to_s)
2278
2279	1	if ((CS%Kh_eta_bg > 0.0) .or. (CS%Kh_eta_vel > 0.0)) then
2280	0	allocate(CS%Kh_eta_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.)
2281	0	allocate(CS%Kh_eta_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.)
2282	0	do j=G%jsc,G%jec ; do I=G%isc-1,G%iec
2283	0	grid_sp = sqrt((2.0G%dxCu(I,j)2 G%dyCu(I,j)2) / ((G%dxCu(I,j)2) + (G%dyCu(I,j)**2)))
2284	0	CS%Kh_eta_u(I,j) = G%OBCmaskCu(I,j) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp)
2285		enddo ; enddo
2286	0	do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
2287	0	grid_sp = sqrt((2.0G%dxCv(i,J)2 G%dyCv(i,J)2) / ((G%dxCv(i,J)2) + (G%dyCv(i,J)**2)))
2288	0	CS%Kh_eta_v(i,J) = G%OBCmaskCv(i,J) * MAX(0.0, CS%Kh_eta_bg + CS%Kh_eta_vel * grid_sp)
2289		enddo ; enddo
2290		endif
2291
2292	1	if (CS%max_Khth_CFL < 0.0) CS%max_Khth_CFL = 0.0
2293		call get_param(param_file, mdl, "DETANGLE_INTERFACES", CS%detangle_interfaces, &
2294		"If defined add 3-d structured enhanced interface height "//&
2295		"diffusivities to horizontally smooth jagged layers.", &
2296	1	default=.false.)
2297	1	CS%detangle_time = 0.0
2298	1	if (CS%detangle_interfaces) &
2299		call get_param(param_file, mdl, "DETANGLE_TIMESCALE", CS%detangle_time, &
2300		"A timescale over which maximally jagged grid-scale "//&
2301		"thickness variations are suppressed. This must be "//&
2302	0	"longer than DT, or 0 to use DT.", units="s", default=0.0, scale=US%s_to_T)
2303		call get_param(param_file, mdl, "KHTH_SLOPE_MAX", CS%slope_max, &
2304		"A slope beyond which the calculated isopycnal slope is "//&
2305	1	"not reliable and is scaled away.", units="nondim", default=0.01, scale=US%L_to_Z)
2306		call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, &
2307		"A diapycnal diffusivity that is used to interpolate "//&
2308		"more sensible values of T & S into thin layers.", &
2309	1	units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T)
2310		call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", CS%use_FGNV_streamfn, &
2311		"If true, use the streamfunction formulation of "//&
2312		"Ferrari et al., 2010, which effectively emphasizes "//&
2313		"graver vertical modes by smoothing in the vertical.", &
2314	1	default=.false.)
2315		call get_param(param_file, mdl, "FGNV_FILTER_SCALE", CS%FGNV_scale, &
2316		"A coefficient scaling the vertical smoothing term in the "//&
2317		"Ferrari et al., 2010, streamfunction formulation.", &
2318	1	units="nondim", default=1., do_not_log=.not.CS%use_FGNV_streamfn)
2319		call get_param(param_file, mdl, "FGNV_C_MIN", CS%FGNV_c_min, &
2320		"A minium wave speed used in the Ferrari et al., 2010, "//&
2321		"streamfunction formulation.", &
2322	1	default=0., units="m s-1", scale=US%m_s_to_L_T, do_not_log=.not.CS%use_FGNV_streamfn)
2323		call get_param(param_file, mdl, "FGNV_STRAT_FLOOR", strat_floor, &
2324		"A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "//&
2325		"streamfunction formulation, expressed as a fraction of planetary "//&
2326		"rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
2327	1	default=1.e-15, units="nondim", scale=US%Z_to_L, do_not_log=.not.CS%use_FGNV_streamfn)
2328		call get_param(param_file, mdl, "USE_STANLEY_GM", CS%use_stanley_gm, &
2329		"If true, turn on Stanley SGS T variance parameterization "// &
2330	1	"in GM code.", default=.false.)
2331	1	if (CS%use_stanley_gm) then
2332		call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
2333		"Coefficient correlating the temperature gradient and SGS T variance.", &
2334	0	units="nondim", default=-1.0, do_not_log=.true.)
2335	0	if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
2336	0	"STANLEY_COEFF must be set >= 0 if USE_STANLEY_GM is true.")
2337		endif
2338		call get_param(param_file, mdl, "OMEGA", omega, &
2339		"The rotation rate of the earth.", &
2340	1	default=7.2921e-5, units="s-1", scale=US%T_to_s, do_not_log=.not.CS%use_FGNV_streamfn)
2341	1	CS%N2_floor = 0.
2342	1	if (CS%use_FGNV_streamfn) CS%N2_floor = (strat_flooromega)*2
2343		call get_param(param_file, mdl, "DEBUG", CS%debug, &
2344		"If true, write out verbose debugging data.", &
2345	1	default=.false., debuggingParam=.true.)
2346
2347		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2348		"This sets the default value for the various _ANSWER_DATE parameters.", &
2349	1	default=99991231, do_not_log=.true.)
2350
2351		call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", CS%GM_src_alt, &
2352		"If true, use the GM energy conversion form S^2N^2kappa rather "//&
2353	1	"than the streamfunction for the GM source term.", default=.false.)
2354		call get_param(param_file, mdl, "MEKE_GM_SRC_ANSWER_DATE", CS%MEKE_src_answer_date, &
2355		"The vintage of the expressions in the GM energy conversion calculation when "//&
2356		"MEKE_GM_SRC_ALT is true. Values below 20240601 recover the answers from the "//&
2357		"original implementation, while higher values use expressions that satisfy "//&
2358		"rotational symmetry.", &
2359	1	default=default_answer_date, do_not_log=.not.CS%GM_src_alt)
2360		call get_param(param_file, mdl, "MEKE_GM_SRC_ALT_SLOPE_BUG", CS%MEKE_src_slope_bug, &
2361		"If true, use a bug that limits the positive values, but not the negative values, "//&
2362		"of the slopes used when MEKE_GM_SRC_ALT is true. When this is true, it breaks "//&
2363		"all of the symmetry rules that MOM6 is supposed to obey.", &
2364	1	default=.false., do_not_log=.not.CS%GM_src_alt)
2365
2366		call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, &
2367		"If true, uses the GM coefficient formulation from the GEOMETRIC "//&
2368	1	"framework (Marshall et al., 2012).", default=.false.)
2369	1	if (CS%MEKE_GEOMETRIC) then
2370		call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", CS%MEKE_GEOMETRIC_epsilon, &
2371		"Minimum Eady growth rate used in the calculation of GEOMETRIC "//&
2372	0	"thickness diffusivity.", units="s-1", default=1.0e-7, scale=US%T_to_s)
2373		call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, &
2374		"The nondimensional coefficient governing the efficiency of the GEOMETRIC "//&
2375	0	"thickness diffusion.", units="nondim", default=0.05)
2376
2377		call get_param(param_file, mdl, "MEKE_GEOMETRIC_ANSWER_DATE", CS%MEKE_GEOM_answer_date, &
2378		"The vintage of the expressions in the MEKE_GEOMETRIC calculation. "//&
2379		"Values below 20190101 recover the answers from the original implementation, "//&
2380		"while higher values use expressions that satisfy rotational symmetry.", &
2381	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
2382	0	if (.not.GV%Boussinesq) CS%MEKE_GEOM_answer_date = max(CS%MEKE_GEOM_answer_date, 20230701)
2383		endif
2384
2385	1	call get_param(param_file, mdl, "USE_MEKE", use_meke, default=.false., do_not_log=.true.)
2386	1	if (use_meke) then
2387		call get_param(param_file, mdl, "USE_KH_IN_MEKE", CS%Use_KH_in_MEKE, &
2388		"If true, uses the thickness diffusivity calculated here to diffuse MEKE.", &
2389	1	default=.false.)
2390		call get_param(param_file, mdl, "MEKE_MIN_DEPTH_DIFF", CS%MEKE_min_depth_diff, &
2391		"The minimum total depth over which to average the diffusivity used for MEKE. "//&
2392		"When the total depth is less than this, the diffusivity is scaled away.", &
2393	1	units="m", default=1.0, scale=GV%m_to_H, do_not_log=.not.CS%Use_KH_in_MEKE)
2394		endif
2395
2396		call get_param(param_file, mdl, "USE_GME", CS%use_GME_thickness_diffuse, &
2397		"If true, use the GM+E backscatter scheme in association "//&
2398	1	"with the Gent and McWilliams parameterization.", default=.false.)
2399
2400		call get_param(param_file, mdl, "USE_GM_WORK_BUG", CS%use_GM_work_bug, &
2401		"If true, compute the top-layer work tendency on the u-grid "//&
2402		"with the incorrect sign, for legacy reproducibility.", &
2403	1	default=.false.)
2404
2405	1	if (CS%use_GME_thickness_diffuse) then
2406	0	allocate(CS%KH_u_GME(G%IsdB:G%IedB, G%jsd:G%jed, GV%ke+1), source=0.)
2407	0	allocate(CS%KH_v_GME(G%isd:G%ied, G%JsdB:G%JedB, GV%ke+1), source=0.)
2408		endif
2409
2410		CS%id_uhGM = register_diag_field('ocean_model', 'uhGM', diag%axesCuL, Time, &
2411		'Time Mean Diffusive Zonal Thickness Flux', &
2412		'kg s-1', conversion=GV%H_to_kg_m2US%L_to_m2US%s_to_T, &
2413	1	y_cell_method='sum', v_extensive=.true.)
2414	1	if (CS%id_uhGM > 0) call safe_alloc_ptr(CDp%uhGM,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
2415		CS%id_vhGM = register_diag_field('ocean_model', 'vhGM', diag%axesCvL, Time, &
2416		'Time Mean Diffusive Meridional Thickness Flux', &
2417		'kg s-1', conversion=GV%H_to_kg_m2US%L_to_m2US%s_to_T, &
2418	1	x_cell_method='sum', v_extensive=.true.)
2419	1	if (CS%id_vhGM > 0) call safe_alloc_ptr(CDp%vhGM,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
2420
2421		CS%id_GMwork = register_diag_field('ocean_model', 'GMwork', diag%axesT1, Time, &
2422		'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
2423		'W m-2', conversion=US%RZ3_T3_to_W_m2US%L_to_Z*2, cmor_field_name='tnkebto', &
2424		cmor_long_name='Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
2425	1	cmor_standard_name='tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection')
2426	1	if (CS%id_GMwork > 0) &
2427	0	allocate(CS%GMwork(G%isd:G%ied,G%jsd:G%jed), source=0.)
2428
2429		CS%id_KH_u = register_diag_field('ocean_model', 'KHTH_u', diag%axesCui, Time, &
2430		'Parameterized mesoscale eddy advection diffusivity at U-point', &
2431	1	'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2432		CS%id_KH_v = register_diag_field('ocean_model', 'KHTH_v', diag%axesCvi, Time, &
2433		'Parameterized mesoscale eddy advection diffusivity at V-point', &
2434	1	'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2435		CS%id_KH_t = register_diag_field('ocean_model', 'KHTH_t', diag%axesTL, Time, &
2436		'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
2437		'm2 s-1', conversion=US%L_to_m*2US%s_to_T, &
2438		cmor_field_name='diftrblo', &
2439		cmor_long_name='Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
2440	1	cmor_standard_name='ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection')
2441
2442		CS%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', diag%axesCu1, Time, &
2443		'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)', &
2444	1	'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2445		CS%id_KH_v1 = register_diag_field('ocean_model', 'KHTH_v1', diag%axesCv1, Time, &
2446		'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)', &
2447	1	'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2448		CS%id_KH_t1 = register_diag_field('ocean_model', 'KHTH_t1', diag%axesT1, Time, &
2449		'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)', &
2450	1	'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
2451
2452		CS%id_slope_x = register_diag_field('ocean_model', 'neutral_slope_x', diag%axesCui, Time, &
2453	1	'Zonal slope of neutral surface', 'nondim', conversion=US%Z_to_L)
2454	1	if (CS%id_slope_x > 0) &
2455	0	allocate(CS%diagSlopeX(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke+1), source=0.)
2456
2457		CS%id_slope_y = register_diag_field('ocean_model', 'neutral_slope_y', diag%axesCvi, Time, &
2458	1	'Meridional slope of neutral surface', 'nondim', conversion=US%Z_to_L)
2459	1	if (CS%id_slope_y > 0) &
2460	0	allocate(CS%diagSlopeY(G%isd:G%ied,G%JsdB:G%JedB,GV%ke+1), source=0.)
2461
2462		CS%id_sfn_x = register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, Time, &
2463		'Parameterized Zonal Overturning Streamfunction', &
2464	1	'm3 s-1', conversion=GV%H_to_mUS%L_to_m2US%s_to_T)
2465		CS%id_sfn_y = register_diag_field('ocean_model', 'GM_sfn_y', diag%axesCvi, Time, &
2466		'Parameterized Meridional Overturning Streamfunction', &
2467	1	'm3 s-1', conversion=GV%H_to_mUS%L_to_m2US%s_to_T)
2468		CS%id_sfn_unlim_x = register_diag_field('ocean_model', 'GM_sfn_unlim_x', diag%axesCui, Time, &
2469		'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
2470	1	'm3 s-1', conversion=US%Z_to_mUS%L_to_m2US%s_to_T)
2471		CS%id_sfn_unlim_y = register_diag_field('ocean_model', 'GM_sfn_unlim_y', diag%axesCvi, Time, &
2472		'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
2473	1	'm3 s-1', conversion=US%Z_to_mUS%L_to_m2US%s_to_T)
2474
2475	1	end subroutine thickness_diffuse_init
2476
2477		!> Copies KH_u_GME and KH_v_GME from private type into arrays provided as arguments
2478	0	subroutine thickness_diffuse_get_KH(CS, KH_u_GME, KH_v_GME, G, GV)
2479		type(thickness_diffuse_CS), intent(in) :: CS !< Control structure for this module
2480		type(ocean_grid_type), intent(in) :: G !< Grid structure
2481		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2482		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: KH_u_GME !< Isopycnal height
2483		!! diffusivities at u-faces [L2 T-1 ~> m2 s-1]
2484		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: KH_v_GME !< Isopycnal height
2485		!! diffusivities at v-faces [L2 T-1 ~> m2 s-1]
2486		! Local variables
2487		integer :: i,j,k
2488
2489	0	do k=1,GV%ke+1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
2490	0	KH_u_GME(I,j,k) = CS%KH_u_GME(I,j,k)
2491		enddo ; enddo ; enddo
2492
2493	0	do k=1,GV%ke+1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
2494	0	KH_v_GME(i,J,k) = CS%KH_v_GME(i,J,k)
2495		enddo ; enddo ; enddo
2496
2497	0	end subroutine thickness_diffuse_get_KH
2498
2499		!> Deallocate the thickness_diffus3 control structure
2500	1	subroutine thickness_diffuse_end(CS, CDp)
2501		type(thickness_diffuse_CS), intent(inout) :: CS !< Control structure for thickness_diffuse
2502		type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity diagnostic control structure
2503
2504	1	if (CS%id_slope_x > 0) deallocate(CS%diagSlopeX)
2505	1	if (CS%id_slope_y > 0) deallocate(CS%diagSlopeY)
2506
2507	1	if (CS%id_GMwork > 0) deallocate(CS%GMwork)
2508
2509		! NOTE: [uv]hGM may be allocated either here or the diagnostic module
2510	1	if (associated(CDp%uhGM)) deallocate(CDp%uhGM)
2511	1	if (associated(CDp%vhGM)) deallocate(CDp%vhGM)
2512
2513	1	if (CS%use_GME_thickness_diffuse) then
2514	0	deallocate(CS%KH_u_GME)
2515	0	deallocate(CS%KH_v_GME)
2516		endif
2517
2518	1	if (allocated(CS%khth2d)) deallocate(CS%khth2d)
2519	1	end subroutine thickness_diffuse_end
2520
2521		!> \namespace mom_thickness_diffuse
2522		!!
2523		!! \section section_gm Isopycnal height diffusion (aka Gent-McWilliams)
2524		!!
2525		!! Isopycnal height diffusion is implemented via along-layer mass fluxes
2526		!! \f[
2527		!! h^\dagger \leftarrow h^n - \Delta t \nabla \cdot ( \vec{uh}^* )
2528		!! \f]
2529		!! where the mass fluxes are cast as the difference in vector streamfunction
2530		!!
2531		!! \f[
2532		!! \vec{uh}^* = \delta_k \vec{\psi} .
2533		!! \f]
2534		!!
2535		!! The GM implementation of isopycnal height diffusion made the streamfunction proportional
2536		!! to the potential density slope
2537		!! \f[
2538		!! \vec{\psi} = - \kappa_h \frac{\nabla_z \rho}{\partial_z \rho}
2539		!! = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = -\kappa_h \frac{M^2}{N^2}
2540		!! \f]
2541		!! but for robustness the scheme is implemented as
2542		!! \f[
2543		!! \vec{\psi} = -\kappa_h \frac{M^2}{\sqrt{N^4 + M^4}}
2544		!! \f]
2545		!! since the quantity \f$\frac{M^2}{\sqrt{N^4 + M^4}}\f$ is bounded between $-1$ and $1$ and does not change sign
2546		!! if \f$N^2<0\f$.
2547		!!
2548		!! Optionally, the method of \cite ferrari2010, can be used to obtain the streamfunction which solves the
2549		!! vertically elliptic equation:
2550		!! \f[
2551		!! \gamma_F \partial_z c^2 \partial_z \psi - N_^2 \psi = -( 1 + \gamma_F ) \kappa_h N_^2 \frac{M^2}{\sqrt{N^4+M^4}}
2552		!! \f]
2553		!! which recovers the previous streamfunction relation in the limit that \f$ c \rightarrow 0 \f$.
2554		!! Here, \f$c=\max(c_{min},c_g)\f$ is the maximum of either \f$c_{min}\f$ and either the first baroclinic mode
2555		!! wave-speed or the equivalent barotropic mode wave-speed.
2556		!! \f$N_*^2 = \max(N^2,0)\f$ is a non-negative form of the square of the buoyancy frequency.
2557		!! The parameter \f$\gamma_F\f$ is used to reduce the vertical smoothing length scale.
2558		!! \f[
2559		!! \kappa_h = \left( \kappa_o + \alpha_{s} L_{s}^2 < S N > + \alpha_{M} \kappa_{M} \right) r(\Delta x,L_d)
2560		!! \f]
2561		!! where \f$ S \f$ is the isoneutral slope magnitude, \f$ N \f$ is the buoyancy frequency,
2562		!! \f$\kappa_{M}\f$ is the diffusivity calculated by the MEKE parameterization (mom_meke module) and
2563		!! \f$ r(\Delta x,L_d) \f$ is a function of the local resolution (ratio of grid-spacing, \f$\Delta x\f$,
2564		!! to deformation radius, \f$L_d\f$). The length \f$L_s\f$ is provided by the mom_lateral_mixing_coeffs module
2565		!! (enabled with <code>USE_VARIABLE_MIXING=True</code> and the term \f$<SN>\f$ is the vertical average slope
2566		!! times the buoyancy frequency prescribed by \cite visbeck1996.
2567		!!
2568		!! The result of the above expression is subsequently bounded by minimum and maximum values, including an upper
2569		!! diffusivity consistent with numerical stability (\f$ \kappa_{cfl} \f$ is calculated internally).
2570		!! \f[
2571		!! \kappa_h \leftarrow \min{\left( \kappa_{max}, \kappa_{cfl}, \max{\left( \kappa_{min}, \kappa_h \right)} \right)}
2572		!! f(c_g,z)
2573		!! \f]
2574		!!
2575		!! where \f$f(c_g,z)\f$ is a vertical structure function.
2576		!! \f$f(c_g,z)\f$ is calculated in module mom_lateral_mixing_coeffs.
2577		!! If <code>KHTH_USE_EBT_STRUCT=True</code> then \f$f(c_g,z)\f$ is set to look like the equivalent barotropic
2578		!! modal velocity structure. Otherwise \f$f(c_g,z)=1\f$ and the diffusivity is independent of depth.
2579		!!
2580		!! In order to calculate meaningful slopes in vanished layers, temporary copies of the thermodynamic variables
2581		!! are passed through a vertical smoother, function vert_fill_ts():
2582		!! \f{eqnarray*}{
2583		!! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] \theta & \leftarrow & \theta \\
2584		!! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] s & \leftarrow & s
2585		!! \f}
2586		!!
2587		!! \subsection section_khth_module_parameters Module mom_thickness_diffuse parameters
2588		!!
2589		!! \| Symbol \| Module parameter \|
2590		!! \| ------ \| --------------- \|
2591		!! \| - \| <code>THICKNESSDIFFUSE</code> \|
2592		!! \| \f$ \kappa_o \f$ \| <code>KHTH</code> \|
2593		!! \| \f$ \alpha_{s} \f$ \| <code>KHTH_SLOPE_CFF</code> \|
2594		!! \| \f$ \kappa_{min} \f$ \| <code>KHTH_MIN</code> \|
2595		!! \| \f$ \kappa_{max} \f$ \| <code>KHTH_MAX</code> \|
2596		!! \| - \| <code>KHTH_MAX_CFL</code> \|
2597		!! \| \f$ \kappa_{smth} \f$ \| <code>KD_SMOOTH</code> \|
2598		!! \| \f$ \alpha_{M} \f$ \| <code>MEKE_KHTH_FAC</code> (from mom_meke module) \|
2599		!! \| - \| <code>KHTH_USE_EBT_STRUCT</code> (from mom_lateral_mixing_coeffs module) \|
2600		!! \| - \| <code>KHTH_USE_FGNV_STREAMFUNCTION</code> \|
2601		!! \| \f$ \gamma_F \f$ \| <code>FGNV_FILTER_SCALE</code> \|
2602		!! \| \f$ c_{min} \f$ \| <code>FGNV_C_MIN</code> \|
2603		!!
2604		!! \subsection section_khth_module_reference References
2605		!!
2606		!! Ferrari, R., S.M. Griffies, A.J.G. Nurser and G.K. Vallis, 2010:
2607		!! A boundary-value problem for the parameterized mesoscale eddy transport.
2608		!! Ocean Modelling, 32, 143-156. http://doi.org/10.1016/j.ocemod.2010.01.004
2609		!!
2610		!! Visbeck, M., J.C. Marshall, H. Jones, 1996:
2611		!! Dynamics of isolated convective regions in the ocean. J. Phys. Oceangr., 26, 1721-1734.
2612		!! http://dx.doi.org/10.1175/1520-0485(1996)026%3C1721:DOICRI%3E2.0.CO;2
2613
2614	0	end module MOM_thickness_diffuse