src/tracer/MOM_tracer_hor_diff.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		#include "do_concurrent_compat.h"
5
6		!> Main routine for lateral (along surface or neutral) diffusion of tracers
7		module MOM_tracer_hor_diff
8
9		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10		use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE
11		use MOM_diag_mediator, only : post_data, diag_ctrl
12		use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
13		use MOM_domains, only : sum_across_PEs, max_across_PEs
14		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
15		use MOM_domains, only : pass_vector
16		use MOM_debugging, only : hchksum, uvchksum
17		use MOM_diabatic_driver, only : diabatic_CS
18		use MOM_EOS, only : calculate_density, EOS_type, EOS_domain
19		use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe
20		use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery
21		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
22		use MOM_file_parser, only : get_param, log_version, param_file_type
23		use MOM_grid, only : ocean_grid_type
24		use MOM_lateral_mixing_coeffs, only : VarMix_CS
25		use MOM_MEKE_types, only : MEKE_type
26		use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end
27		use MOM_neutral_diffusion, only : neutral_diffusion_CS
28		use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion
29		use MOM_hor_bnd_diffusion, only : hbd_CS, hor_bnd_diffusion_init
30		use MOM_hor_bnd_diffusion, only : hor_bnd_diffusion, hor_bnd_diffusion_end
31		use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum
32		use MOM_unit_scaling, only : unit_scale_type
33		use MOM_variables, only : thermo_var_ptrs, vertvisc_type
34		use MOM_verticalGrid, only : verticalGrid_type
35
36		implicit none ; private
37
38		#include <MOM_memory.h>
39
40		public tracer_hordiff, tracer_hor_diff_init, tracer_hor_diff_end
41
42		!> The control structure for along-layer and epineutral tracer diffusion
43		type, public :: tracer_hor_diff_CS ; private
44		real :: KhTr !< The along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
45		real :: KhTr_Slope_Cff !< The non-dimensional coefficient in KhTr formula [nondim]
46		real :: KhTr_min !< Minimum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
47		real :: KhTr_max !< Maximum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
48		real :: KhTr_passivity_coeff !< Passivity coefficient that scales Rd/dx (default = 0)
49		!! where passivity is the ratio between along-isopycnal
50		!! tracer mixing and thickness mixing [nondim]
51		real :: KhTr_passivity_min !< Passivity minimum (default = 1/2) [nondim]
52		real :: ML_KhTR_scale !< With Diffuse_ML_interior, the ratio of the
53		!! truly horizontal diffusivity in the mixed
54		!! layer to the epipycnal diffusivity [nondim].
55		real :: max_diff_CFL !< If positive, locally limit the along-isopycnal
56		!! tracer diffusivity to keep the diffusive CFL
57		!! locally at or below this value [nondim].
58		logical :: KhTr_use_vert_struct !< If true, uses the equivalent barotropic structure
59		!! as the vertical structure of tracer diffusivity.
60		logical :: full_depth_khtr_min !< If true, KHTR_MIN is enforced throughout the whole water column.
61		!! Otherwise, KHTR_MIN is only enforced at the surface. This parameter
62		!! is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.
63		logical :: Diffuse_ML_interior !< If true, diffuse along isopycnals between
64		!! the mixed layer and the interior.
65		logical :: check_diffusive_CFL !< If true, automatically iterate the diffusion
66		!! to ensure that the diffusive equivalent of
67		!! the CFL limit is not violated.
68		logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within
69		!! tracer_hor_diff.
70		logical :: use_hor_bnd_diffusion !< If true, use the hor_bnd_diffusion module from within
71		!! tracer_hor_diff.
72		logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been
73		!! exceeded
74		logical :: limit_bug !< If true and the answer date is 20240330 or below, use a
75		!! rotational symmetry breaking bug when limiting the tracer
76		!! properties in tracer_epipycnal_ML_diff.
77		integer :: answer_date !< The vintage of the order of arithmetic to use for the tracer
78		!! diffusion. Values of 20240330 or below recover the answers
79		!! from the original form of this code, while higher values use
80		!! mathematically equivalent expressions that recover rotational symmetry
81		!! when DIFFUSE_ML_TO_INTERIOR is true.
82		type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion.
83		type(hbd_CS), pointer :: hor_bnd_diffusion_CSp => NULL() !< Control structure for
84		!! horizontal boundary diffusion.
85		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
86		!! regulate the timing of diagnostic output.
87		logical :: debug !< If true, write verbose checksums for debugging purposes.
88		logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level.
89		logical :: first_call = .true. !< This is true until after the first call
90		!>@{ Diagnostic IDs
91		integer :: id_KhTr_u = -1
92		integer :: id_KhTr_v = -1
93		integer :: id_KhTr_h = -1
94		integer :: id_CFL = -1
95		integer :: id_khdt_x = -1
96		integer :: id_khdt_y = -1
97		!>@}
98
99		type(group_pass_type) :: pass_t !< For group halo pass, used in both
100		!! tracer_hordiff and tracer_epipycnal_ML_diff
101		end type tracer_hor_diff_CS
102
103		!> A type that can be used to create arrays of pointers to 2D arrays
104		type p2d
105		real, dimension(:,:), pointer :: p => NULL() !< A pointer to a 2D array of reals [various]
106		end type p2d
107		!> A type that can be used to create arrays of pointers to 2D integer arrays
108		type p2di
109		integer, dimension(:,:), pointer :: p => NULL() !< A pointer to a 2D array of integers
110		end type p2di
111
112		!>@{ CPU time clocks
113		integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync
114		!>@}
115
116		contains
117
118		!> Compute along-coordinate diffusion of all tracers
119		!! using the diffusivity in CS%KhTr, or using space-dependent diffusivity.
120		!! Multiple iterations are used (if necessary) so that there is no limit
121		!! on the acceptable time increment.
122	0	subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y)
123		type(ocean_grid_type), intent(inout) :: G !< Grid type
124		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
125		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
126		intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
127		real, intent(in) :: dt !< time step [T ~> s]
128		type(MEKE_type), intent(in) :: MEKE !< MEKE fields
129		type(VarMix_CS), intent(in) :: VarMix !< Variable mixing type
130		type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities,
131		!! boundary layer properties and related fields
132		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
133		type(tracer_hor_diff_CS), pointer :: CS !< module control structure
134		type(tracer_registry_type), pointer :: Reg !< registered tracers
135		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
136		!! thermodynamic fields, including potential temp and
137		!! salinity or mixed layer density. Absent fields have
138		!! NULL ptrs, and these may (probably will) point to
139		!! some of the same arrays as Tr does. tv is required
140		!! for epipycnal mixing between mixed layer and the interior.
141		! Optional inputs for offline tracer transport
142		logical, optional, intent(in) :: do_online_flag !< If present and true, do online
143		!! tracer transport with stored velocities.
144		! The next two arguments do not appear to be used anywhere.
145		real, dimension(SZIB_(G),SZJ_(G)), &
146		optional, intent(in) :: read_khdt_x !< If present, these are the zonal diffusivities
147		!! times a timestep from a previous run [L2 ~> m2]
148		real, dimension(SZI_(G),SZJB_(G)), &
149		optional, intent(in) :: read_khdt_y !< If present, these are the meridional diffusivities
150		!! times a timestep from a previous run [L2 ~> m2]
151
152
153		real, dimension(SZI_(G),SZJ_(G)) :: &
154	24	Ihdxdy, & ! The inverse of the volume or mass of fluid in a layer in a
155		! grid cell [H-1 L-2 ~> m-3 or kg-1].
156	12	CFL, & ! A diffusive CFL number for each cell [nondim].
157	24	dTr ! The change in a tracer's concentration, in units of concentration [Conc].
158
159	24	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: Kh_h
160		! The tracer diffusivity averaged to tracer points [L2 T-1 ~> m2 s-1].
161		real, dimension(SZIB_(G),SZJ_(G)) :: &
162	24	khdt_x ! The value of Khtr*dt times the open face width divided by
163		! the distance between adjacent tracer points [L2 ~> m2].
164		real, dimension(SZI_(G),SZJB_(G)) :: &
165	24	khdt_y ! The value of Khtr*dt times the open face width divided by
166		! the distance between adjacent tracer points [L2 ~> m2].
167		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
168	24	Coef_x, & ! The coefficients relating zonal tracer differences to time-integrated
169		! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others.
170	24	Kh_u ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
171		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
172	24	Coef_y, & ! The coefficients relating meridional tracer differences to time-integrated
173		! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others.
174	12	Kh_v ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
175
176		real :: khdt_max ! The local limiting value of khdt_x or khdt_y [L2 ~> m2].
177		real :: Coef_min ! The local limiting value of Coef_x or Coef_y, in [L2 ~> m2] for some
178		! schemes and [H L2 ~> m3 or kg] for others.
179		real :: max_CFL ! The global maximum of the diffusive CFL number [nondim]
180		logical :: use_VarMix, Resoln_scaled, do_online, use_Eady
181		integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts
182		real :: I_numitts ! The inverse of the number of iterations, num_itts [nondim]
183		real :: scale ! The fraction of khdt_x or khdt_y that is applied in this
184		! layer for this iteration [nondim].
185		real :: Idt ! The inverse of the time step [T-1 ~> s-1].
186		real :: h_neglect ! A thickness that is so small it is usually lost
187		! in roundoff and can be neglected [H ~> m or kg m-2].
188		real :: Kh_loc ! The local value of Kh [L2 T-1 ~> m2 s-1].
189		real :: Res_Fn ! The local value of the resolution function [nondim].
190		real :: Rd_dx ! The local value of deformation radius over grid-spacing [nondim].
191		real :: normalize ! normalization used for diagnostic Kh_h [nondim]; diffusivity averaged to h-points.
192		real :: MEKE_KhTr_fac
193
194	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
195
196	12	do_online = .true.
197	12	if (present(do_online_flag)) do_online = do_online_flag
198
199	12	if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// &
200	0	"register_tracer must be called before tracer_hordiff.")
201	12	if (.not. associated(Reg)) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// &
202	0	"register_tracer must be called before tracer_hordiff.")
203	12	if (Reg%ntr == 0 .or. (CS%KhTr <= 0.0 .and. .not. VarMix%use_variable_mixing)) return
204
205	12	if (CS%show_call_tree) call callTree_enter("tracer_hordiff(), MOM_tracer_hor_diff.F90")
206
207	12	call cpu_clock_begin(id_clock_diffuse)
208
209		!$omp target enter data map(to: Reg, Reg%Tr, CS) map(alloc: khdt_x, khdt_y, kh_u, kh_v)
210		!$ do m = 1, Reg%ntr
211		!$omp target enter data map(to: Reg%Tr(m)%t)
212		!$omp target enter data map(to: Reg%Tr(m)%df_x) if(associated(Reg%Tr(m)%df_x))
213		!$omp target enter data map(to: Reg%Tr(m)%df_y) if(associated(Reg%Tr(m)%df_y))
214		!$omp target enter data map(to: Reg%Tr(m)%df2d_x) if(associated(Reg%Tr(m)%df2d_x))
215		!$omp target enter data map(to: Reg%Tr(m)%df2d_y) if(associated(Reg%Tr(m)%df2d_y))
216		!$ enddo
217
218	12	ntr = Reg%ntr
219	12	Idt = 1.0 / dt
220	12	h_neglect = GV%H_subroundoff
221
222	12	if (CS%Diffuse_ML_interior .and. CS%first_call) then ; if (is_root_pe()) then
223	0	do m=1,ntr ; if (associated(Reg%Tr(m)%df_x) .or. associated(Reg%Tr(m)%df_y)) &
224		call MOM_error(WARNING, "tracer_hordiff: Tracer "//trim(Reg%Tr(m)%name)// &
225		" has associated 3-d diffusive flux diagnostics. These are not "//&
226		"valid when DIFFUSE_ML_TO_INTERIOR is defined. Use 2-d tracer "//&
227	0	"diffusion diagnostics instead to get accurate total fluxes.")
228		enddo
229		endif ; endif
230	12	CS%first_call = .false.
231
232	12	if (CS%debug) call MOM_tracer_chksum("Before tracer diffusion ", Reg, G)
233
234	12	use_VarMix = .false. ; Resoln_scaled = .false. ; use_Eady = .false.
235	12	if (VarMix%use_variable_mixing) then
236	12	use_VarMix = VarMix%use_variable_mixing
237	12	Resoln_scaled = VarMix%Resoln_scaled_KhTr
238	12	use_Eady = CS%KhTr_Slope_Cff > 0.
239	12	CS%KhTr_use_vert_struct = allocated(VarMix%khtr_struct)
240		endif
241
242	12	call cpu_clock_begin(id_clock_pass)
243	48	do m=1,ntr
244	48	call create_group_pass(CS%pass_t, Reg%Tr(m)%t(:,:,:), G%Domain)
245		enddo
246	12	call cpu_clock_end(id_clock_pass)
247
248	12	if (CS%show_call_tree) call callTree_waypoint("Calculating diffusivity (tracer_hordiff)")
249
250	12	if (do_online) then
251	12	if (use_VarMix) then
252	12	MEKE_KhTr_fac = MEKE%KhTr_fac
253		!$omp target enter data map(to: VarMix, VarMix%SN_u, VarMix%L2u, VarMix%Res_fn_h, &
254		!$omp VarMix%Rd_dx_h, MEKE, MEKE%Kh)
255	12	do concurrent (j=js:je, I=is-1:ie)
256	87120	Kh_loc = CS%KhTr
257	87120	if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_CffVarMix%L2u(I,j)VarMix%SN_u(I,j)
258	87120	if (allocated(MEKE%Kh)) &
259	87120	Kh_loc = Kh_loc + MEKE_KhTr_facsqrt(MEKE%Kh(i,j)MEKE%Kh(i+1,j))
260	87120	if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max)
261	87120	if (Resoln_scaled) &
262	0	Kh_loc = Kh_loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j))
263	87120	Kh_u(I,j,1) = max(Kh_loc, CS%KhTr_min)
264	88584	if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity
265	0	Rd_dx=0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points
266	0	Kh_loc = Kh_u(I,j,1)max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeffRd_dx )
267	0	if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max
268	0	Kh_u(I,j,1) = max(Kh_loc, CS%KhTr_min) ! Re-apply min
269		endif
270		enddo
271	12	do concurrent (J=js-1:je, i=is:ie)
272	87840	Kh_loc = CS%KhTr
273	87840	if (use_Eady) Kh_loc = Kh_loc + CS%KhTr_Slope_CffVarMix%L2v(i,J)VarMix%SN_v(i,J)
274	87840	if (allocated(MEKE%Kh)) &
275	87840	Kh_loc = Kh_loc + MEKE_KhTr_facsqrt(MEKE%Kh(i,j)MEKE%Kh(i,j+1))
276	87840	if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max)
277	87840	if (Resoln_scaled) &
278	0	Kh_loc = Kh_loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1))
279	87840	Kh_v(i,J,1) = max(Kh_loc, CS%KhTr_min)
280	89292	if (CS%KhTr_passivity_coeff>0.) then ! Apply passivity
281	0	Rd_dx = 0.5*( VarMix%Rd_dx_h(i,j)+VarMix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points
282	0	Kh_loc = Kh_v(i,J,1)max( CS%KhTr_passivity_min, CS%KhTr_passivity_coeffRd_dx )
283	0	if (CS%KhTr_max > 0.) Kh_loc = min(Kh_loc, CS%KhTr_max) ! Re-apply max
284	0	Kh_v(i,J,1) = max(Kh_loc, CS%KhTr_min) ! Re-apply min
285		endif
286		enddo
287		!$omp target exit data map(release: VarMix, VarMix%SN_u, VarMix%L2u, VarMix%SN_v, &
288		!$omp VarMix%L2v, VarMix%Res_fn_h, VarMix%Rd_dx_h, MEKE, MEKE%Kh)
289
290	12	do concurrent (j=js:je, I=is-1:ie)
291	88584	khdt_x(I,j) = dt(Kh_u(I,j,1)(G%dy_Cu(I,j)*G%IdxCu(I,j)))
292		enddo
293	12	do concurrent (J=js-1:je, i=is:ie)
294	89292	khdt_y(i,J) = dt(Kh_v(i,J,1)(G%dx_Cv(i,J)*G%IdyCv(i,J)))
295		enddo
296	0	elseif (Resoln_scaled) then
297		!$OMP parallel do default(shared) private(Res_fn)
298	0	do j=js,je ; do I=is-1,ie
299	0	Res_fn = 0.5 * (VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j))
300	0	Kh_u(I,j,1) = max(CS%KhTr * Res_fn, CS%KhTr_min)
301	0	khdt_x(I,j) = dt(CS%KhTr(G%dy_Cu(I,j)G%IdxCu(I,j))) Res_fn
302		enddo ; enddo
303		!$OMP parallel do default(shared) private(Res_fn)
304	0	do J=js-1,je ; do i=is,ie
305	0	Res_fn = 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1))
306	0	Kh_v(i,J,1) = max(CS%KhTr * Res_fn, CS%KhTr_min)
307	0	khdt_y(i,J) = dt(CS%KhTr(G%dx_Cv(i,J)G%IdyCv(i,J))) Res_fn
308		enddo ; enddo
309		!$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v)
310		else ! Use a simple constant diffusivity.
311	0	if (CS%id_KhTr_u > 0) then
312		!$OMP parallel do default(shared)
313	0	do j=js,je ; do I=is-1,ie
314	0	Kh_u(I,j,1) = CS%KhTr
315	0	khdt_x(I,j) = dt(CS%KhTr(G%dy_Cu(I,j)*G%IdxCu(I,j)))
316		enddo ; enddo
317		else
318		!$OMP parallel do default(shared)
319	0	do j=js,je ; do I=is-1,ie
320	0	khdt_x(I,j) = dt(CS%KhTr(G%dy_Cu(I,j)*G%IdxCu(I,j)))
321		enddo ; enddo
322		endif
323	0	if (CS%id_KhTr_v > 0) then
324		!$OMP parallel do default(shared)
325	0	do J=js-1,je ; do i=is,ie
326	0	Kh_v(i,J,1) = CS%KhTr
327	0	khdt_y(i,J) = dt(CS%KhTr(G%dx_Cv(i,J)*G%IdyCv(i,J)))
328		enddo ; enddo
329		else
330		!$OMP parallel do default(shared)
331	0	do J=js-1,je ; do i=is,ie
332	0	khdt_y(i,J) = dt(CS%KhTr(G%dx_Cv(i,J)*G%IdyCv(i,J)))
333		enddo ; enddo
334		endif
335		!$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v)
336		endif ! VarMix
337
338	12	if (CS%max_diff_CFL > 0.0) then
339		!$omp target update from(khdt_x, khdt_y, Kh_u, Kh_v)
340	0	if ((CS%id_KhTr_u > 0) .or. (CS%id_KhTr_h > 0)) then
341		!$OMP parallel do default(shared) private(khdt_max)
342	0	do j=js,je ; do I=is-1,ie
343	0	khdt_max = 0.125CS%max_diff_CFL min(G%areaT(i,j), G%areaT(i+1,j))
344	0	if (khdt_x(I,j) > khdt_max) then
345	0	khdt_x(I,j) = khdt_max
346	0	if (dt(G%dy_Cu(I,j)G%IdxCu(I,j)) > 0.0) &
347	0	Kh_u(I,j,1) = khdt_x(I,j) / (dt(G%dy_Cu(I,j)G%IdxCu(I,j)))
348		endif
349		enddo ; enddo
350		else
351		!$OMP parallel do default(shared) private(khdt_max)
352	0	do j=js,je ; do I=is-1,ie
353	0	khdt_max = 0.125CS%max_diff_CFL min(G%areaT(i,j), G%areaT(i+1,j))
354	0	khdt_x(I,j) = min(khdt_x(I,j), khdt_max)
355		enddo ; enddo
356		endif
357	0	if ((CS%id_KhTr_v > 0) .or. (CS%id_KhTr_h > 0)) then
358		!$OMP parallel do default(shared) private(khdt_max)
359	0	do J=js-1,je ; do i=is,ie
360	0	khdt_max = 0.125CS%max_diff_CFL min(G%areaT(i,j), G%areaT(i,j+1))
361	0	if (khdt_y(i,J) > khdt_max) then
362	0	khdt_y(i,J) = khdt_max
363	0	if (dt(G%dx_Cv(i,J)G%IdyCv(i,J)) > 0.0) &
364	0	Kh_v(i,J,1) = khdt_y(i,J) / (dt(G%dx_Cv(i,J)G%IdyCv(i,J)))
365		endif
366		enddo ; enddo
367		else
368		!$OMP parallel do default(shared) private(khdt_max)
369	0	do J=js-1,je ; do i=is,ie
370	0	khdt_max = 0.125CS%max_diff_CFL min(G%areaT(i,j), G%areaT(i,j+1))
371	0	khdt_y(i,J) = min(khdt_y(i,J), khdt_max)
372		enddo ; enddo
373		endif
374		!$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v)
375		endif
376
377		else ! .not. do_online
378		!$OMP parallel do default(shared)
379	0	do j=js,je ; do I=is-1,ie
380	0	khdt_x(I,j) = read_khdt_x(I,j)
381		enddo ; enddo
382		!$OMP parallel do default(shared)
383	0	do J=js-1,je ; do i=is,ie
384	0	khdt_y(i,J) = read_khdt_y(i,J)
385		enddo ; enddo
386	0	call pass_vector(khdt_x, khdt_y, G%Domain)
387		!$omp target update to(khdt_x, khdt_y, Kh_u, Kh_v)
388		endif ! do_online
389
390	12	if (CS%check_diffusive_CFL) then
391	12	if (CS%show_call_tree) call callTree_waypoint("Checking diffusive CFL (tracer_hordiff)")
392	12	max_CFL = 0.0
393	87132	do j=js,je ; do i=is,ie
394		CFL(i,j) = 2.0*((khdt_x(I-1,j) + khdt_x(I,j)) + &
395	86400	(khdt_y(i,J-1) + khdt_y(i,J))) * G%IareaT(i,j)
396	87120	if (max_CFL < CFL(i,j)) max_CFL = CFL(i,j)
397		enddo ; enddo
398	12	call cpu_clock_begin(id_clock_sync)
399	12	call max_across_PEs(max_CFL)
400	12	call cpu_clock_end(id_clock_sync)
401	12	num_itts = max(1, ceiling(max_CFL - 4.0*EPSILON(max_CFL)))
402	12	I_numitts = 1.0 / (real(num_itts))
403	12	if (CS%id_CFL > 0) call post_data(CS%id_CFL, CFL, CS%diag)
404	0	elseif (CS%max_diff_CFL > 0.0) then
405	0	num_itts = max(1, ceiling(CS%max_diff_CFL - 4.0*EPSILON(CS%max_diff_CFL)))
406	0	I_numitts = 1.0 / (real(num_itts))
407		else
408	0	num_itts = 1 ; I_numitts = 1.0
409		endif
410
411	48	do m=1,ntr
412	36	if (associated(Reg%Tr(m)%df_x)) then
413	0	do k=1,nz ; do j=js,je ; do I=is-1,ie
414	0	Reg%Tr(m)%df_x(I,j,k) = 0.0
415		enddo ; enddo ; enddo
416		!$omp target update to(Reg%Tr(m)%df_x)
417		endif
418	36	if (associated(Reg%Tr(m)%df_y)) then
419	0	do k=1,nz ; do J=js-1,je ; do i=is,ie
420	0	Reg%Tr(m)%df_y(i,J,k) = 0.0
421		enddo ; enddo ; enddo
422		!$omp target update to(Reg%Tr(m)%df_y)
423		endif
424	36	if (associated(Reg%Tr(m)%df2d_x)) then
425	0	do j=js,je ; do I=is-1,ie ; Reg%Tr(m)%df2d_x(I,j) = 0.0 ; enddo ; enddo
426		!$omp target update to(Reg%Tr(m)%df2d_x)
427		endif
428	48	if (associated(Reg%Tr(m)%df2d_y)) then
429	0	do J=js-1,je ; do i=is,ie ; Reg%Tr(m)%df2d_y(i,J) = 0.0 ; enddo ; enddo
430		!$omp target update to(Reg%Tr(m)%df2d_y)
431		endif
432		enddo
433
434	12	if (CS%use_hor_bnd_diffusion) then
435		!$omp target update from(h)
436
437	0	if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)")
438		!$omp target update from(khdt_x, khdt_y)
439	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass)
440
441	0	do k=1,nz+1
442	0	do J=js-1,je
443	0	do i=is,ie
444	0	Coef_y(i,J,K) = I_numitts * khdt_y(i,J)
445		enddo
446		enddo
447		enddo
448	0	do k=1,nz+1
449	0	do j=js,je
450	0	do I=is-1,ie
451	0	Coef_x(I,j,K) = I_numitts * khdt_x(I,j)
452		enddo
453		enddo
454		enddo
455	0	if (CS%KhTr_use_vert_struct) then
456	0	if (CS%full_depth_khtr_min) then
457	0	do K=2,nz+1
458	0	do J=js-1,je
459	0	do i=is,ie
460	0	Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) )
461	0	Coef_min = I_numitts * dt * (CS%KhTr_min(G%dx_Cv(i,J)G%IdyCv(i,J)))
462	0	Coef_y(i,J,K) = max(Coef_y(i,J,K), Coef_min)
463		enddo
464		enddo
465		enddo
466	0	do k=2,nz+1
467	0	do j=js,je
468	0	do I=is-1,ie
469	0	Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) )
470	0	Coef_min = I_numitts * dt * (CS%KhTr_min(G%dy_Cu(I,j)G%IdxCu(I,j)))
471	0	Coef_x(I,j,K) = max(Coef_x(I,j,K), Coef_min)
472		enddo
473		enddo
474		enddo
475		else
476	0	do K=2,nz+1
477	0	do J=js-1,je
478	0	do i=is,ie
479	0	Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i,j+1,k-1) )
480		enddo
481		enddo
482		enddo
483	0	do k=2,nz+1
484	0	do j=js,je
485	0	do I=is-1,ie
486	0	Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%ebt_struct(i,j,k-1) + VarMix%ebt_struct(i+1,j,k-1) )
487		enddo
488		enddo
489		enddo
490		endif
491		endif
492
493	0	do itt=1,num_itts
494	0	if (CS%show_call_tree) call callTree_waypoint("Calling horizontal boundary diffusion (tracer_hordiff)",itt)
495	0	if (itt>1) then ! Update halos for subsequent iterations
496	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass)
497		endif
498		call hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, visc, &
499	0	CS%hor_bnd_diffusion_CSp)
500		enddo ! itt
501		endif
502
503	12	if (CS%use_neutral_diffusion) then
504		!$omp target update from(h)
505
506	0	if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)")
507
508	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass)
509		! We are assuming that neutral surfaces do not evolve (much) as a result of multiple
510		!horizontal diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs()
511		! would be inside the itt-loop. -AJA
512
513	0	if (associated(tv%p_surf)) then
514	0	call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, p_surf=tv%p_surf)
515		else
516	0	call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp)
517		endif
518
519	0	do k=1,nz+1
520	0	do J=js-1,je
521	0	do i=is,ie
522	0	Coef_y(i,J,K) = I_numitts * khdt_y(i,J)
523		enddo
524		enddo
525		enddo
526	0	do k=1,nz+1
527	0	do j=js,je
528	0	do I=is-1,ie
529	0	Coef_x(I,j,K) = I_numitts * khdt_x(I,j)
530		enddo
531		enddo
532		enddo
533	0	if (CS%KhTr_use_vert_struct) then
534	0	do K=2,nz+1
535	0	do J=js-1,je
536	0	do i=is,ie
537	0	Coef_y(i,J,K) = Coef_y(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) )
538		enddo
539		enddo
540		enddo
541	0	do k=2,nz+1
542	0	do j=js,je
543	0	do I=is-1,ie
544	0	Coef_x(I,j,K) = Coef_x(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) )
545		enddo
546		enddo
547		enddo
548		endif
549
550	0	do itt=1,num_itts
551	0	if (CS%show_call_tree) call callTree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt)
552	0	if (itt>1) then ! Update halos for subsequent iterations
553	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass)
554	0	if (CS%recalc_neutral_surf) then
555	0	if (associated(tv%p_surf)) then
556		call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, &
557	0	p_surf=tv%p_surf)
558		else
559	0	call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp)
560		endif
561		endif
562		endif
563	0	call neutral_diffusion(G, GV, h, Coef_x, Coef_y, I_numitts*dt, Reg, US, CS%neutral_diffusion_CSp)
564		enddo ! itt
565
566		else ! following if not using neutral diffusion, but instead along-surface diffusion
567
568	12	if (CS%show_call_tree) call callTree_waypoint("Calculating horizontal diffusion (tracer_hordiff)")
569		!$omp target enter data map(alloc: dTr, Ihdxdy)
570		!$omp target enter data map(to: Coef_x, Coef_y) if(CS%use_hor_bnd_diffusion) ! copy from non-ported loops above
571		!$omp target enter data map(alloc: Coef_x, Coef_y) if(.not.CS%use_hor_bnd_diffusion)
572	24	do itt=1,num_itts
573	12	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.)
574		! loop should probably be reordered
575	912	do k=1,nz
576	900	scale = I_numitts
577	900	if (CS%Diffuse_ML_interior) then
578	0	if (k<=GV%nkml) then
579	0	if (CS%ML_KhTr_scale <= 0.0) cycle
580	0	scale = I_numitts * CS%ML_KhTr_scale
581		endif
582	0	if ((k>GV%nkml) .and. (k<=GV%nk_rho_varies)) cycle
583		endif
584
585	900	do concurrent (J=js-1:je, i=is:ie)
586		Coef_y(i,J,1) = ((scale * khdt_y(i,J))2.0(h(i,j,k)*h(i,j+1,k))) / &
587	6696900	(h(i,j,k)+h(i,j+1,k)+h_neglect)
588		enddo
589
590	900	do concurrent (j=js:je)
591	54000	do concurrent (I=is-1:ie)
592		Coef_x(I,j,1) = ((scale * khdt_x(I,j))2.0(h(i,j,k)*h(i+1,j,k))) / &
593	6588000	(h(i,j,k)+h(i+1,j,k)+h_neglect)
594		enddo
595
596	54900	do concurrent (i=is:ie)
597	6534000	Ihdxdy(i,j) = G%IareaT(i,j) / (h(i,j,k)+h_neglect)
598		enddo
599		enddo
600
601	3612	do m=1,ntr
602	2700	do concurrent (j=js:je, i=is:ie)
603		dTr(i,j) = Ihdxdy(i,j) * &
604		( ((Coef_x(I-1,j,1) * (Reg%Tr(m)%t(i-1,j,k) - Reg%Tr(m)%t(i,j,k))) - &
605		(Coef_x(I,j,1) * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)))) + &
606		((Coef_y(i,J-1,1) * (Reg%Tr(m)%t(i,j-1,k) - Reg%Tr(m)%t(i,j,k))) - &
607	19766700	(Coef_y(i,J,1) * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)))) )
608		enddo
609	2700	if (associated(Reg%Tr(m)%df_x)) then ; do j=js,je ; do I=G%IscB,G%IecB
610		Reg%Tr(m)%df_x(I,j,k) = Reg%Tr(m)%df_x(I,j,k) + Coef_x(I,j,1) &
611	0	* (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt
612		enddo ; enddo ; endif
613	2700	if (associated(Reg%Tr(m)%df_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie
614		Reg%Tr(m)%df_y(i,J,k) = Reg%Tr(m)%df_y(i,J,k) + Coef_y(i,J,1) &
615	0	* (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt
616		enddo ; enddo ; endif
617	2700	if (associated(Reg%Tr(m)%df2d_x)) then ; do j=js,je ; do I=G%IscB,G%IecB
618		Reg%Tr(m)%df2d_x(I,j) = Reg%Tr(m)%df2d_x(I,j) + Coef_x(I,j,1) &
619	0	* (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt
620		enddo ; enddo ; endif
621	2700	if (associated(Reg%Tr(m)%df2d_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie
622		Reg%Tr(m)%df2d_y(i,J) = Reg%Tr(m)%df2d_y(i,J) + Coef_y(i,J,1) &
623	0	* (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt
624		enddo ; enddo ; endif
625	900	do concurrent (j=js:je, i=is:ie)
626	19766700	Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j)
627		enddo
628		enddo
629
630		enddo ! End of k loop.
631
632		! Do user controlled underflow of the tracer concentrations.
633	60	do m=1,ntr ; if (Reg%Tr(m)%conc_underflow > 0.0) then
634	0	do concurrent (k=1:nz, j=js:je, i=is:ie, abs(Reg%Tr(m)%t(i,j,k)) < Reg%Tr(m)%conc_underflow)
635	0	Reg%Tr(m)%t(i,j,k) = 0.0
636		enddo
637		endif ; enddo
638
639		enddo ! End of "while" loop.
640		!$omp target exit data map(release: dTr, Ihdxdy, Coef_x, Coef_y)
641		endif ! endif for CS%use_neutral_diffusion
642	12	call cpu_clock_end(id_clock_diffuse)
643
644	12	if (CS%Diffuse_ML_interior) then
645	0	if (CS%show_call_tree) call callTree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)")
646	0	if (CS%debug) call MOM_tracer_chksum("Before epipycnal diff ", Reg, G)
647
648	0	call cpu_clock_begin(id_clock_epimix)
649		call tracer_epipycnal_ML_diff(h, dt, Reg%Tr, ntr, khdt_x, khdt_y, G, GV, US, &
650	0	CS, tv, num_itts)
651	0	call cpu_clock_end(id_clock_epimix)
652		endif
653		!$ do m = 1, Reg%ntr
654		!$omp target exit data map(from: Reg%Tr(m)%t)
655		!$omp target exit data map(from: Reg%Tr(m)%df_x) if(associated(Reg%Tr(m)%df_x))
656		!$omp target exit data map(from: Reg%Tr(m)%df_y) if(associated(Reg%Tr(m)%df_y))
657		!$omp target exit data map(from: Reg%Tr(m)%df2d_x) if(associated(Reg%Tr(m)%df2d_x))
658		!$omp target exit data map(from: Reg%Tr(m)%df2d_y) if(associated(Reg%Tr(m)%df2d_y))
659		!$ enddo
660		!$omp target exit data map(release: Reg%Tr, Reg)
661
662	12	if (CS%debug) call MOM_tracer_chksum("After tracer diffusion ", Reg, G)
663
664		! post diagnostics for 2d tracer diffusivity
665	12	if (CS%id_KhTr_u > 0) then
666		!$omp target exit data map(from: Kh_u)
667	0	do j=js,je ; do I=is-1,ie
668	0	Kh_u(I,j,:) = G%mask2dCu(I,j)*Kh_u(I,j,1)
669		enddo ; enddo
670	0	if (CS%KhTr_use_vert_struct) then
671	0	do K=2,nz+1
672	0	do j=js,je
673	0	do I=is-1,ie
674	0	Kh_u(I,j,K) = Kh_u(I,j,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i+1,j,k-1) )
675		enddo
676		enddo
677		enddo
678		endif
679	0	call post_data(CS%id_KhTr_u, Kh_u, CS%diag)
680		endif
681	12	if (CS%id_KhTr_v > 0) then
682		!$omp target exit data map(from: Kh_v)
683	0	do J=js-1,je ; do i=is,ie
684	0	Kh_v(i,J,:) = G%mask2dCv(i,J)*Kh_v(i,J,1)
685		enddo ; enddo
686	0	if (CS%KhTr_use_vert_struct) then
687	0	do K=2,nz+1
688	0	do J=js-1,je
689	0	do i=is,ie
690	0	Kh_v(i,J,K) = Kh_v(i,J,1) * 0.5 * ( VarMix%khtr_struct(i,j,k-1) + VarMix%khtr_struct(i,j+1,k-1) )
691		enddo
692		enddo
693		enddo
694		endif
695	0	call post_data(CS%id_KhTr_v, Kh_v, CS%diag)
696		endif
697	12	if (CS%id_KhTr_h > 0) then
698		!$omp target exit data map(from: Kh_u, Kh_v)
699	0	Kh_h(:,:,:) = 0.0
700	0	do j=js,je ; do I=is-1,ie
701	0	Kh_u(I,j,1) = G%mask2dCu(I,j)*Kh_u(I,j,1)
702		enddo ; enddo
703	0	do J=js-1,je ; do i=is,ie
704	0	Kh_v(i,J,1) = G%mask2dCv(i,J)*Kh_v(i,J,1)
705		enddo ; enddo
706
707	0	do j=js,je ; do i=is,ie
708		normalize = 1.0 / ((G%mask2dCu(I-1,j)+G%mask2dCu(I,j)) + &
709	0	(G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37)
710		Kh_h(i,j,:) = normalizeG%mask2dT(i,j)((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + &
711	0	(Kh_v(i,J-1,1)+Kh_v(i,J,1)))
712	0	if (CS%KhTr_use_vert_struct) then
713	0	do K=2,nz+1
714		Kh_h(i,j,K) = normalizeG%mask2dT(i,j)VarMix%khtr_struct(i,j,k-1)*((Kh_u(I-1,j,1)+Kh_u(I,j,1)) + &
715	0	(Kh_v(i,J-1,1)+Kh_v(i,J,1)))
716		enddo
717		endif
718		enddo ; enddo
719	0	call post_data(CS%id_KhTr_h, Kh_h, CS%diag)
720		endif
721
722		!$omp target update from(khdt_x, khdt_y) if(CS%debug .or. CS%id_khdt_x>0 .or. CS%id_khdt_y>0)
723		!$omp target exit data map(release: khdt_x, khdt_y, Kh_u, Kh_v) map(release: CS)
724
725	12	if (CS%debug) then
726		call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, &
727		G%HI, haloshift=0, symmetric=.true., unscale=US%L_to_m**2, &
728	0	scalar_pair=.true.)
729		endif
730
731	12	if (CS%id_khdt_x > 0) call post_data(CS%id_khdt_x, khdt_x, CS%diag)
732	12	if (CS%id_khdt_y > 0) call post_data(CS%id_khdt_y, khdt_y, CS%diag)
733
734	12	if (CS%show_call_tree) call callTree_leave("tracer_hordiff()")
735
736	12	end subroutine tracer_hordiff
737
738		!> This subroutine does epipycnal diffusion of all tracers between the mixed
739		!! and buffer layers and the interior, using the diffusivity in CS%KhTr.
740		!! Multiple iterations are used (if necessary) so that there is no limit on the
741		!! acceptable time increment.
742	0	subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
743		GV, US, CS, tv, num_itts)
744		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
745		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
746		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
747		real, intent(in) :: dt !< time step [T ~> s]
748		type(tracer_type), intent(inout) :: Tr(:) !< tracer array
749		integer, intent(in) :: ntr !< number of tracers
750		real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: khdt_epi_x !< Zonal epipycnal diffusivity times
751		!! a time step and the ratio of the open face width over
752		!! the distance between adjacent tracer points [L2 ~> m2]
753		real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times
754		!! a time step and the ratio of the open face width over
755		!! the distance between adjacent tracer points [L2 ~> m2]
756		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
757		type(tracer_hor_diff_CS), intent(inout) :: CS !< module control structure
758		type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure
759		integer, intent(in) :: num_itts !< number of iterations (usually=1)
760
761
762		real, dimension(SZI_(G), SZJ_(G)) :: &
763	0	Rml_max ! The maximum coordinate density within the mixed layer [R ~> kg m-3].
764		real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: &
765	0	rho_coord ! The coordinate density that is used to mix along [R ~> kg m-3].
766
767		! The naming mnemonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point.
768		! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage.
769		real, dimension(:,:,:), allocatable :: &
770	0	deep_wt_Lu, deep_wt_Ru, & ! The relative weighting of the deeper of a pair [nondim].
771	0	hP_Lu, hP_Ru, & ! The total thickness on each side for each pair [H ~> m or kg m-2].
772	0	deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair [nondim].
773	0	hP_Lv, hP_Rv ! The total thickness on each side for each pair [H ~> m or kg m-2].
774
775		integer, dimension(:,:,:), allocatable :: &
776	0	k0b_Lu, k0a_Lu, & ! The original k-indices of the layers that participate
777	0	k0b_Ru, k0a_Ru, & ! in each pair of mixing at u-faces.
778	0	k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate
779	0	k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces.
780
781		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
782	0	tr_flux_N, & ! The tracer flux through the northern face [conc H L2 ~> conc m3 or conc kg]
783	0	tr_flux_S, & ! The tracer flux through the southern face [conc H L2 ~> conc m3 or conc kg]
784	0	tr_flux_E, & ! The tracer flux through the eastern face [conc H L2 ~> conc m3 or conc kg]
785	0	tr_flux_W, & ! The tracer flux through the western face [conc H L2 ~> conc m3 or conc kg]
786	0	tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
787
788		! The following 3-d arrays were created in 2014 in MOM6 PR#12 to facilitate openMP threading
789		! on an i-loop, which might have been ill advised.
790		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)*2) :: &
791	0	Tr_flux_3d, & ! The tracer flux through pairings at meridional faces [conc H L2 ~> conc m3 or conc kg]
792	0	Tr_adj_vert_L, & ! Vertical adjustments to which layer the fluxes go into in the southern
793		! columns at meridional face [conc H L2 ~> conc m3 or conc kg]
794	0	Tr_adj_vert_R ! Vertical adjustments to which layer the fluxes go into in the northern
795		! columns at meridional face [conc H L2 ~> conc m3 or conc kg]
796
797		real, dimension(SZI_(G),SZK_(GV), SZJ_(G)) :: &
798	0	rho_srt, & ! The density of each layer of the sorted columns [R ~> kg m-3].
799	0	h_srt ! The thickness of each layer of the sorted columns [H ~> m or kg m-2].
800		integer, dimension(SZI_(G),SZK_(GV), SZJ_(G)) :: &
801	0	k0_srt ! The original k-index that each layer of the sorted column corresponds to.
802
803		real, dimension(SZK_(GV)) :: &
804	0	h_demand_L, & ! The thickness in the left column that is demanded to match the thickness
805		! in the counterpart [H ~> m or kg m-2].
806	0	h_demand_R, & ! The thickness in the right column that is demanded to match the thickness
807		! in the counterpart [H ~> m or kg m-2].
808	0	h_used_L, & ! The summed thickness from the left column that has actually been used [H ~> m or kg m-2]
809	0	h_used_R, & ! The summed thickness from the right columns that has actually been used [H ~> m or kg m-2]
810	0	h_supply_frac_L, & ! The fraction of the demanded thickness that can actually be supplied
811		! from a layer on the left [nondim].
812	0	h_supply_frac_R ! The fraction of the demanded thickness that can actually be supplied
813		! from a layer on the right [nondim].
814		integer, dimension(SZI_(G), SZJ_(G)) :: &
815	0	num_srt, & ! The number of layers that are sorted in each column.
816	0	k_end_srt, & ! The maximum index in each column that might need to be
817		! sorted, based on neighboring values of max_kRho
818	0	max_kRho ! The index of the layer whose target density is just denser
819		! than the densest part of the mixed layer.
820		integer, dimension(SZJ_(G)) :: &
821	0	max_srt ! The maximum value of num_srt in a k-row.
822		integer, dimension(SZIB_(G), SZJ_(G)) :: &
823	0	nPu ! The number of epipycnal pairings at each u-point.
824		integer, dimension(SZI_(G), SZJB_(G)) :: &
825	0	nPv ! The number of epipycnal pairings at each v-point.
826		real :: h_exclude ! A thickness that layers must attain to be considered
827		! for inclusion in mixing [H ~> m or kg m-2].
828		real :: Idt ! The inverse of the time step [T-1 ~> s-1].
829		real :: I_maxitt ! The inverse of the maximum number of iterations [nondim]
830		real :: rho_pair, rho_a, rho_b ! Temporary densities [R ~> kg m-3].
831		real :: Tr_min_face ! The minimum tracer concentration associated with a pairing [Conc]
832		real :: Tr_max_face ! The maximum tracer concentration associated with a pairing [Conc]
833		real :: Tr_La, Tr_Lb ! The 2 left-side tracer concentrations that might be associated with a pairing [Conc]
834		real :: Tr_Ra, Tr_Rb ! The 2 right-side tracer concentrations that might be associated with a pairing [Conc]
835		real :: Tr_av_L ! The average tracer concentrations on the left side of a pairing [Conc].
836		real :: Tr_av_R ! The average tracer concentrations on the right side of a pairing [Conc].
837		real :: Tr_flux ! The tracer flux from left to right in a pair [conc H L2 ~> conc m3 or conc kg].
838		real :: Tr_adj_vert ! A downward vertical adjustment to Tr_flux between the two cells that
839		! make up one side of the pairing [conc H L2 ~> conc m3 or conc kg].
840		real :: h_L, h_R ! Thicknesses to the left and right [H ~> m or kg m-2].
841		real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim].
842		real :: vol ! A cell volume or mass [H L2 ~> m3 or kg].
843
844		! The total number of pairings is usually much less than twice the number of layers, but
845		! the memory in these 1-d columns of pairings can be allocated generously for safety.
846		integer, dimension(SZK_(GV)*2) :: &
847	0	kbs_Lp, & ! The sorted indices of the Left and Right columns for
848	0	kbs_Rp ! each pairing.
849		logical, dimension(SZK_(GV)*2) :: &
850	0	left_set, & ! If true, the left or right point determines the density of
851	0	right_set ! of the trio. If densities are exactly equal, both are true.
852
853		real :: tmp ! A temporary variable used in swaps [various]
854	0	real :: p_ref_cv(SZI_(G),SZJ_(G)) ! The reference pressure for the coordinate density [R L2 T-2 ~> Pa]
855
856		integer, dimension(2,2) :: EOSdom ! The i-computational domain for the equation of state
857		integer :: k_max, k_min, k_test, itmp
858		integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb
859		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, k_size
860		integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt
861		integer :: PEmax_kRho
862
863	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
864	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
865	0	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
866	0	Idt = 1.0 / dt
867	0	nkmb = GV%nk_rho_varies
868
869	0	if (num_itts <= 1) then
870	0	max_itt = 1 ; I_maxitt = 1.0
871		else
872	0	max_itt = num_itts ; I_maxitt = 1.0 / (real(max_itt))
873		endif
874
875		!$omp target enter data map(alloc: rho_coord, Rml_max, max_kRho, rho_srt, k0_srt, num_srt, &
876		!$omp h_srt, p_ref_cv, k_end_srt, max_srt)
877
878	0	do concurrent (j=jsd:jed, i=isd:ied)
879	0	p_ref_cv(i,j) = tv%P_Ref
880		enddo
881
882	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.)
883		! Determine which layers the mixed- and buffer-layers map into...
884		!$OMP parallel do default(shared)
885	0	do k=1,nkmb
886		call calculate_density(tv%T(:,:,k), tv%S(:,:,k), p_ref_cv, rho_coord(:,:,k), &
887	0	tv%eqn_of_state)
888		enddo
889		!$omp target exit data map(release: p_ref_cv)
890
891	0	do concurrent (j=js-2:je+2, i=is-2:ie+2)
892	0	Rml_max(i,j) = rho_coord(i,j,1)
893	0	num_srt(i,j) = 0 ; max_kRho(i,j) = 0
894		enddo
895	0	do k=2,nkmb ; do concurrent (j=js-2:je+2, i=is-2:ie+2, Rml_max(i,j) < rho_coord(i,j,k))
896	0	Rml_max(i,j) = rho_coord(i,j,k)
897		enddo ; enddo
898
899		! Use bracketing and bisection to find the k-level that the densest of the
900		! mixed and buffer layer corresponds to, such that:
901		! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
902	0	do concurrent (j=js-2:je+2, i=is-2:ie+2, G%mask2dT(i,j) > 0.0)
903	0	if ((Rml_max(i,j) > GV%Rlay(nz)) .or. (nkmb+1 > nz)) then ; max_kRho(i,j) = nz+1
904	0	elseif ((Rml_max(i,j) <= GV%Rlay(nkmb+1)) .or. (nkmb+2 > nz)) then ; max_kRho(i,j) = nkmb+1
905		else
906	0	k_min = nkmb+2 ; k_max = nz
907	0	do
908	0	k_test = (k_min + k_max) / 2
909	0	if (Rml_max(i,j) <= GV%Rlay(k_test-1)) then ; k_max = k_test-1
910	0	elseif (GV%Rlay(k_test) < Rml_max(i,j)) then ; k_min = k_test+1
911	0	else ; max_kRho(i,j) = k_test ; exit ; endif
912
913	0	if (k_min == k_max) then ; max_kRho(i,j) = k_max ; exit ; endif
914		enddo
915		endif
916		enddo
917
918	0	PEmax_kRho = 0
919	0	do concurrent (j=js-1:je+1, i=is-1:ie+1) DO_LOCALITY(reduce(max:PEmax_kRho))
920		k_end_srt(i,j) = max(max_kRho(i,j), max_kRho(i-1,j), max_kRho(i+1,j), &
921	0	max_kRho(i,j-1), max_kRho(i,j+1))
922	0	if (PEmax_kRho < k_end_srt(i,j)) PEmax_kRho = k_end_srt(i,j)
923		enddo
924	0	if (PEmax_kRho > nz) PEmax_kRho = nz ! PEmax_kRho could have been nz+1.
925
926	0	h_exclude = 10.0*(GV%Angstrom_H + GV%H_subroundoff)
927
928	0	do concurrent (j=js-1:je+1) DO_LOCALITY(local(k, ns))
929	0	do k=1,nkmb ; do concurrent (i=is-1:ie+1, G%mask2dT(i,j) > 0.0) DO_LOCALITY(local(ns))
930	0	if (h(i,j,k) > h_exclude) then
931	0	num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
932	0	k0_srt(i,ns,j) = k
933	0	rho_srt(i,ns,j) = rho_coord(i,j,k)
934	0	h_srt(i,ns,j) = h(i,j,k)
935		endif
936		enddo ; enddo
937	0	do k=nkmb+1,PEmax_kRho ; do concurrent (i=is-1:ie+1, G%mask2dT(i,j) > 0.0) DO_LOCALITY(local(ns))
938	0	if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
939	0	num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
940	0	k0_srt(i,ns,j) = k
941	0	rho_srt(i,ns,j) = GV%Rlay(k)
942	0	h_srt(i,ns,j) = h(i,j,k)
943		endif
944		enddo ; enddo
945		enddo
946		!$omp target exit data map(release: rho_coord, Rml_max, k_end_srt)
947		! Sort each column by increasing density. This should already be close,
948		! and the size of the arrays are small, so straight insertion is used.
949	0	do concurrent (j=js-1:je+1, i=is-1:ie+1)
950	0	do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then
951		! The last segment needs to be shuffled earlier in the list.
952	0	do k2 = k,2,-1 ; if (rho_srt(i,k2,j) < rho_srt(i,k2-1,j)) then
953	0	itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp
954	0	tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp
955	0	tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp
956		endif ; enddo
957		endif ; enddo
958		enddo
959	0	do concurrent (j=js-1:je+1) DO_LOCALITY(local(itmp))
960		! max_srt(j) = 0
961	0	itmp = 0
962		! nvfortran do concurrent cannot reduce array elements
963	0	do concurrent (i=is-1:ie+1) DO_LOCALITY(reduce(max:itmp))
964	0	itmp = max(itmp, num_srt(i,j))
965		enddo
966	0	max_srt(j) = itmp
967		enddo
968	0	k_size = 1
969	0	do concurrent (j=js-1:je+1) DO_LOCALITY(reduce(max:k_size))
970	0	k_size = max(k_size, 2*max_srt(j))
971		enddo
972	0	allocate(k0a_Lu(IsdB:iedB,k_size,jsd:jed))
973	0	allocate(k0a_Ru(IsdB:iedB,k_size,jsd:jed))
974	0	allocate(k0b_Lu(IsdB:iedB,k_size,jsd:jed))
975	0	allocate(k0b_Ru(IsdB:iedB,k_size,jsd:jed))
976	0	allocate(deep_wt_Lu(IsdB:iedB,k_size,jsd:jed))
977	0	allocate(deep_wt_Ru(IsdB:iedB,k_size,jsd:jed))
978	0	allocate(hP_Lu(IsdB:iedB,k_size,jsd:jed))
979	0	allocate(hP_Ru(IsdB:iedB,k_size,jsd:jed))
980		!$omp target enter data map(alloc: deep_wt_Lu, deep_wt_Ru, hP_Lu, hP_Ru, k0a_Lu, k0a_Ru, k0b_Lu, &
981		!$omp k0b_Ru)
982		!$omp target enter data map(alloc: nPu)
983		!$OMP target teams loop collapse(2) &
984		!$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
985		!$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
986		!$OMP wt_b,left_set,right_set,h_supply_frac_R, &
987		!$OMP h_supply_frac_L, k)
988	0	do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
989		! Set up the pairings for fluxes through the zonal faces.
990
991	0	do k=1,num_srt(i,j) ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo
992	0	do k=1,num_srt(i+1,j) ; h_demand_R(k) = 0.0 ; h_used_R(k) = 0.0 ; enddo
993
994		! First merge the left and right lists into a single, sorted list.
995
996		! Discard any layers that are lighter than the lightest in the other
997		! column. They can only participate in mixing as the lighter part of a
998		! pair of points.
999	0	if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then
1000	0	kR = 1
1001	0	do kL=2,num_srt(i,j) ; if (rho_srt(i,kL,j) >= rho_srt(i+1,1,j)) exit ; enddo
1002	0	elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then
1003	0	kL = 1
1004	0	do kR=2,num_srt(i+1,j) ; if (rho_srt(i+1,kR,j) >= rho_srt(i,1,j)) exit ; enddo
1005		else
1006	0	kL = 1 ; kR = 1
1007		endif
1008	0	nP = 0
1009	0	do ! Loop to accumulate pairs of columns.
1010	0	if ((kL > num_srt(i,j)) .or. (kR > num_srt(i+1,j))) exit
1011
1012	0	if (rho_srt(i,kL,j) > rho_srt(i+1,kR,j)) then
1013		! The right point is lighter and defines the density for this trio.
1014	0	nP = nP+1 ; k = nP
1015	0	rho_pair = rho_srt(i+1,kR,j)
1016
1017	0	k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j)
1018	0	k0a_Lu(I,k,j) = k0_srt(i,kL-1,j) ; k0a_Ru(I,k,j) = k0b_Ru(I,k,j)
1019	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1020
1021	0	rho_a = rho_srt(i,kL-1,j) ; rho_b = rho_srt(i,kL,j)
1022	0	wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1023	0	wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1024	0	deep_wt_Lu(I,k,j) = wt_b ; deep_wt_Ru(I,k,j) = 1.0
1025
1026	0	h_demand_L(kL) = h_demand_L(kL) + 0.5h_srt(i+1,kR,j) wt_b
1027	0	h_demand_L(kL-1) = h_demand_L(kL-1) + 0.5h_srt(i+1,kR,j) (1.0-wt_b)
1028
1029	0	kR = kR+1 ; left_set(k) = .false. ; right_set(k) = .true.
1030	0	elseif (rho_srt(i,kL,j) < rho_srt(i+1,kR,j)) then
1031		! The left point is lighter and defines the density for this trio.
1032	0	nP = nP+1 ; k = nP
1033	0	rho_pair = rho_srt(i,kL,j)
1034	0	k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j)
1035	0	k0a_Lu(I,k,j) = k0b_Lu(I,k,j) ; k0a_Ru(I,k,j) = k0_srt(i+1,kR-1,j)
1036
1037	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1038
1039	0	rho_a = rho_srt(i+1,kR-1,j) ; rho_b = rho_srt(i+1,kR,j)
1040	0	wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1041	0	wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1042	0	deep_wt_Lu(I,k,j) = 1.0 ; deep_wt_Ru(I,k,j) = wt_b
1043
1044	0	h_demand_R(kR) = h_demand_R(kR) + 0.5h_srt(i,kL,j) wt_b
1045	0	h_demand_R(kR-1) = h_demand_R(kR-1) + 0.5h_srt(i,kL,j) (1.0-wt_b)
1046
1047	0	kL = kL+1 ; left_set(k) = .true. ; right_set(k) = .false.
1048	0	elseif ((k0_srt(i,kL,j) <= nkmb) .or. (k0_srt(i+1,kR,j) <= nkmb)) then
1049		! The densities are exactly equal and one layer is above the interior.
1050	0	nP = nP+1 ; k = nP
1051	0	k0b_Lu(I,k,j) = k0_srt(i,kL,j) ; k0b_Ru(I,k,j) = k0_srt(i+1,kR,j)
1052	0	k0a_Lu(I,k,j) = k0b_Lu(I,k,j) ; k0a_Ru(I,k,j) = k0b_Ru(I,k,j)
1053	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1054	0	deep_wt_Lu(I,k,j) = 1.0 ; deep_wt_Ru(I,k,j) = 1.0
1055
1056	0	h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i+1,kR,j)
1057	0	h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j)
1058
1059	0	kL = kL+1 ; kR = kR+1 ; left_set(k) = .true. ; right_set(k) = .true.
1060		else ! The densities are exactly equal and in the interior.
1061		! Mixing in this case has already occurred, so accumulate the thickness
1062		! demanded for that mixing and skip onward.
1063	0	h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i+1,kR,j)
1064	0	h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j)
1065
1066	0	kL = kL+1 ; kR = kR+1
1067		endif
1068		enddo ! Loop to accumulate pairs of columns.
1069	0	nPu(I,j) = nP ! This is the number of active pairings.
1070
1071		! Determine what fraction of the thickness "demand" can be supplied.
1072	0	do k=1,num_srt(i+1,j)
1073	0	h_supply_frac_R(k) = 1.0
1074	0	if (h_demand_R(k) > 0.5*h_srt(i+1,k,j)) &
1075	0	h_supply_frac_R(k) = 0.5*h_srt(i+1,k,j) / h_demand_R(k)
1076		enddo
1077	0	do k=1,num_srt(i,j)
1078	0	h_supply_frac_L(k) = 1.0
1079	0	if (h_demand_L(k) > 0.5*h_srt(i,k,j)) &
1080	0	h_supply_frac_L(k) = 0.5*h_srt(i,k,j) / h_demand_L(k)
1081		enddo
1082
1083		! Distribute the "exported" thicknesses proportionately.
1084	0	do k=1,nPu(I,j)
1085	0	kL = kbs_Lp(k) ; kR = kbs_Rp(k)
1086	0	hP_Lu(I,k,j) = 0.0 ; hP_Ru(I,k,j) = 0.0
1087	0	if (left_set(k)) then ! Add the contributing thicknesses on the right.
1088	0	if (deep_wt_Ru(I,k,j) < 1.0) then
1089	0	hP_Ru(I,k,j) = 0.5h_srt(i,kL,j) min(h_supply_frac_R(kR), h_supply_frac_R(kR-1))
1090	0	wt_b = deep_wt_Ru(I,k,j)
1091	0	h_used_R(kR-1) = h_used_R(kR-1) + (1.0 - wt_b)*hP_Ru(I,k,j)
1092	0	h_used_R(kR) = h_used_R(kR) + wt_b*hP_Ru(I,k,j)
1093		else
1094	0	hP_Ru(I,k,j) = 0.5h_srt(i,kL,j) h_supply_frac_R(kR)
1095	0	h_used_R(kR) = h_used_R(kR) + hP_Ru(I,k,j)
1096		endif
1097		endif
1098	0	if (right_set(k)) then ! Add the contributing thicknesses on the left.
1099	0	if (deep_wt_Lu(I,k,j) < 1.0) then
1100	0	hP_Lu(I,k,j) = 0.5h_srt(i+1,kR,j) min(h_supply_frac_L(kL), h_supply_frac_L(kL-1))
1101	0	wt_b = deep_wt_Lu(I,k,j)
1102	0	h_used_L(kL-1) = h_used_L(kL-1) + (1.0 - wt_b)*hP_Lu(I,k,j)
1103	0	h_used_L(kL) = h_used_L(kL) + wt_b*hP_Lu(I,k,j)
1104		else
1105	0	hP_Lu(I,k,j) = 0.5h_srt(i+1,kR,j) h_supply_frac_L(kL)
1106	0	h_used_L(kL) = h_used_L(kL) + hP_Lu(I,k,j)
1107		endif
1108		endif
1109		enddo
1110
1111		! The left-over thickness (at least half the layer thickness) is now
1112		! added to the thicknesses of the importing columns.
1113	0	do k=1,nPu(I,j)
1114	0	if (left_set(k)) hP_Lu(I,k,j) = hP_Lu(I,k,j) + &
1115	0	(h_srt(i,kbs_Lp(k),j) - h_used_L(kbs_Lp(k)))
1116	0	if (right_set(k)) hP_Ru(I,k,j) = hP_Ru(I,k,j) + &
1117	0	(h_srt(i+1,kbs_Rp(k),j) - h_used_R(kbs_Rp(k)))
1118		enddo
1119
1120		endif ; enddo ; enddo ! i- & j- loops over zonal faces.
1121
1122	0	allocate(deep_wt_Lv(isd:ied,k_size,JsdB:JedB))
1123	0	allocate(deep_wt_Rv(isd:ied,k_size,JsdB:JedB))
1124	0	allocate(hP_Lv(isd:ied,k_size,JsdB:JedB))
1125	0	allocate(hP_Rv(isd:ied,k_size,JsdB:JedB))
1126	0	allocate(k0a_Lv(isd:ied,k_size,JsdB:JedB))
1127	0	allocate(k0a_Rv(isd:ied,k_size,JsdB:JedB))
1128	0	allocate(k0b_Lv(isd:ied,k_size,JsdB:JedB))
1129	0	allocate(k0b_Rv(isd:ied,k_size,JsdB:JedB))
1130		!$omp target enter data map(alloc: deep_wt_Lv, deep_wt_Rv, hP_Lv, hP_Rv, k0a_Lv, k0a_Rv, k0b_Lv, &
1131		!$omp k0b_Rv)
1132
1133		!$omp target enter data map(alloc: nPv)
1134
1135		!$OMP target teams loop collapse(2) &
1136		!$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
1137		!$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
1138		!$OMP wt_b,left_set,right_set,h_supply_frac_R, &
1139		!$OMP h_supply_frac_L,k) map(to: num_srt)
1140	0	do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
1141		! Set up the pairings for fluxes through the meridional faces.
1142
1143	0	do k=1,num_srt(i,j) ; h_demand_L(k) = 0.0 ; h_used_L(k) = 0.0 ; enddo
1144	0	do k=1,num_srt(i,j+1) ; h_demand_R(k) = 0.0 ; h_used_R(k) = 0.0 ; enddo
1145
1146		! First merge the left and right lists into a single, sorted list.
1147
1148		! Discard any layers that are lighter than the lightest in the other
1149		! column. They can only participate in mixing as the lighter part of a
1150		! pair of points.
1151	0	if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then
1152	0	kR = 1
1153	0	do kL=2,num_srt(i,j) ; if (rho_srt(i,kL,j) >= rho_srt(i,1,j+1)) exit ; enddo
1154	0	elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then
1155	0	kL = 1
1156	0	do kR=2,num_srt(i,j+1) ; if (rho_srt(i,kR,j+1) >= rho_srt(i,1,j)) exit ; enddo
1157		else
1158	0	kL = 1 ; kR = 1
1159		endif
1160	0	nP = 0
1161	0	do ! Loop to accumulate pairs of columns.
1162	0	if ((kL > num_srt(i,j)) .or. (kR > num_srt(i,j+1))) exit
1163
1164	0	if (rho_srt(i,kL,j) > rho_srt(i,kR,j+1)) then
1165		! The right point is lighter and defines the density for this trio.
1166	0	nP = nP+1 ; k = nP
1167	0	rho_pair = rho_srt(i,kR,j+1)
1168
1169	0	k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1)
1170	0	k0a_Lv(i,k,J) = k0_srt(i,kL-1,j) ; k0a_Rv(i,k,J) = k0b_Rv(i,k,J)
1171	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1172
1173	0	rho_a = rho_srt(i,kL-1,j) ; rho_b = rho_srt(i,kL,j)
1174	0	wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1175	0	wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1176	0	deep_wt_Lv(I,k,J) = wt_b ; deep_wt_Rv(I,k,J) = 1.0
1177
1178	0	h_demand_L(kL) = h_demand_L(kL) + 0.5h_srt(i,kR,j+1) wt_b
1179	0	h_demand_L(kL-1) = h_demand_L(kL-1) + 0.5h_srt(i,kR,j+1) (1.0-wt_b)
1180
1181	0	kR = kR+1 ; left_set(k) = .false. ; right_set(k) = .true.
1182	0	elseif (rho_srt(i,kL,j) < rho_srt(i,kR,j+1)) then
1183		! The left point is lighter and defines the density for this trio.
1184	0	nP = nP+1 ; k = nP
1185	0	rho_pair = rho_srt(i,kL,j)
1186	0	k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1)
1187	0	k0a_Lv(i,k,J) = k0b_Lv(i,k,J) ; k0a_Rv(i,k,J) = k0_srt(i,kR-1,j+1)
1188
1189	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1190
1191	0	rho_a = rho_srt(i,kR-1,j+1) ; rho_b = rho_srt(i,kR,j+1)
1192	0	wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1193	0	wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1194	0	deep_wt_Lv(I,k,J) = 1.0 ; deep_wt_Rv(I,k,J) = wt_b
1195
1196	0	h_demand_R(kR) = h_demand_R(kR) + 0.5h_srt(i,kL,j) wt_b
1197	0	h_demand_R(kR-1) = h_demand_R(kR-1) + 0.5h_srt(i,kL,j) (1.0-wt_b)
1198
1199	0	kL = kL+1 ; left_set(k) = .true. ; right_set(k) = .false.
1200	0	elseif ((k0_srt(i,kL,j) <= nkmb) .or. (k0_srt(i,kR,j+1) <= nkmb)) then
1201		! The densities are exactly equal and one layer is above the interior.
1202	0	nP = nP+1 ; k = nP
1203	0	k0b_Lv(i,k,J) = k0_srt(i,kL,j) ; k0b_Rv(i,k,J) = k0_srt(i,kR,j+1)
1204	0	k0a_Lv(i,k,J) = k0b_Lv(i,k,J) ; k0a_Rv(i,k,J) = k0b_Rv(i,k,J)
1205	0	kbs_Lp(k) = kL ; kbs_Rp(k) = kR
1206	0	deep_wt_Lv(I,k,J) = 1.0 ; deep_wt_Rv(I,k,J) = 1.0
1207
1208	0	h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i,kR,j+1)
1209	0	h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j)
1210
1211	0	kL = kL+1 ; kR = kR+1 ; left_set(k) = .true. ; right_set(k) = .true.
1212		else ! The densities are exactly equal and in the interior.
1213		! Mixing in this case has already occurred, so accumulate the thickness
1214		! demanded for that mixing and skip onward.
1215	0	h_demand_L(kL) = h_demand_L(kL) + 0.5*h_srt(i,kR,j+1)
1216	0	h_demand_R(kR) = h_demand_R(kR) + 0.5*h_srt(i,kL,j)
1217
1218	0	kL = kL+1 ; kR = kR+1
1219		endif
1220		enddo ! Loop to accumulate pairs of columns.
1221	0	nPv(i,J) = nP ! This is the number of active pairings.
1222
1223		! Determine what fraction of the thickness "demand" can be supplied.
1224	0	do k=1,num_srt(i,j+1)
1225	0	h_supply_frac_R(k) = 1.0
1226	0	if (h_demand_R(k) > 0.5*h_srt(i,k,j+1)) &
1227	0	h_supply_frac_R(k) = 0.5*h_srt(i,k,j+1) / h_demand_R(k)
1228		enddo
1229	0	do k=1,num_srt(i,j)
1230	0	h_supply_frac_L(k) = 1.0
1231	0	if (h_demand_L(k) > 0.5*h_srt(i,k,j)) &
1232	0	h_supply_frac_L(k) = 0.5*h_srt(i,k,j) / h_demand_L(k)
1233		enddo
1234
1235		! Distribute the "exported" thicknesses proportionately.
1236	0	do k=1,nPv(i,J)
1237	0	kL = kbs_Lp(k) ; kR = kbs_Rp(k)
1238	0	hP_Lv(I,k,J) = 0.0 ; hP_Rv(I,k,J) = 0.0
1239	0	if (left_set(k)) then ! Add the contributing thicknesses on the right.
1240	0	if (deep_wt_Rv(I,k,J) < 1.0) then
1241	0	hP_Rv(I,k,J) = 0.5h_srt(i,kL,j) min(h_supply_frac_R(kR), h_supply_frac_R(kR-1))
1242	0	wt_b = deep_wt_Rv(I,k,J)
1243	0	h_used_R(kR-1) = h_used_R(kR-1) + (1.0 - wt_b) * hP_Rv(I,k,J)
1244	0	h_used_R(kR) = h_used_R(kR) + wt_b * hP_Rv(I,k,J)
1245		else
1246	0	hP_Rv(I,k,J) = 0.5h_srt(i,kL,j) h_supply_frac_R(kR)
1247	0	h_used_R(kR) = h_used_R(kR) + hP_Rv(I,k,J)
1248		endif
1249		endif
1250	0	if (right_set(k)) then ! Add the contributing thicknesses on the left.
1251	0	if (deep_wt_Lv(I,k,J) < 1.0) then
1252	0	hP_Lv(I,k,J) = 0.5h_srt(i,kR,j+1) min(h_supply_frac_L(kL), h_supply_frac_L(kL-1))
1253	0	wt_b = deep_wt_Lv(I,k,J)
1254	0	h_used_L(kL-1) = h_used_L(kL-1) + (1.0 - wt_b) * hP_Lv(I,k,J)
1255	0	h_used_L(kL) = h_used_L(kL) + wt_b * hP_Lv(I,k,J)
1256		else
1257	0	hP_Lv(I,k,J) = 0.5h_srt(i,kR,j+1) h_supply_frac_L(kL)
1258	0	h_used_L(kL) = h_used_L(kL) + hP_Lv(I,k,J)
1259		endif
1260		endif
1261		enddo
1262
1263		! The left-over thickness (at least half the layer thickness) is now
1264		! added to the thicknesses of the importing columns.
1265	0	do k=1,nPv(i,J)
1266	0	if (left_set(k)) hP_Lv(I,k,J) = hP_Lv(I,k,J) + &
1267	0	(h_srt(i,kbs_Lp(k),j) - h_used_L(kbs_Lp(k)))
1268	0	if (right_set(k)) hP_Rv(I,k,J) = hP_Rv(I,k,J) + &
1269	0	(h_srt(i,kbs_Rp(k),j+1) - h_used_R(kbs_Rp(k)))
1270		enddo
1271
1272
1273		endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1274		!$omp target exit data map(release: h_srt, k0_srt)
1275		! The tracer-specific calculations start here.
1276
1277		!$omp target enter data map(alloc: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R, tr_flux_N, &
1278		!$omp tr_flux_S, tr_flux_E, tr_flux_W, tr_flux_conv)
1279
1280	0	do itt=1,max_itt
1281
1282	0	if (itt > 1) then ! The halos have already been filled if itt==1.
1283	0	call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1284		endif
1285	0	do m=1,ntr
1286		! Zero out tracer tendencies.
1287	0	if (CS%answer_date <= 20240330) then
1288	0	do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1289	0	tr_flux_conv(i,j,k) = 0.0
1290		enddo
1291		else
1292	0	do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1293	0	tr_flux_N(i,j,k) = 0.0 ; tr_flux_S(i,j,k) = 0.0
1294	0	tr_flux_E(i,j,k) = 0.0 ; tr_flux_W(i,j,k) = 0.0
1295		enddo
1296		endif
1297	0	do concurrent (k=1:2*nz, J=JsdB:JedB, i=isd:ied)
1298	0	tr_flux_3d(i,j,k) = 0.0
1299	0	tr_adj_vert_R(i,j,k) = 0.0 ; tr_adj_vert_L(i,j,k) = 0.0
1300		enddo
1301
1302		! collapse(2) is reproducible for CS%answer_date > 20240330 but not <= 20240330.
1303		! Do concurrent around j-loop doesn't seem to do the right thing.
1304		!$omp target teams loop collapse(2) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, &
1305		!$omp Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, &
1306		!$omp h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
1307	0	do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
1308		! Determine the fluxes through the zonal faces.
1309
1310		! Find the acceptable range of tracer concentration around this face.
1311	0	if (nPu(I,j) >= 1) then
1312	0	Tr_min_face = min(Tr(m)%t(i,j,1), Tr(m)%t(i+1,j,1))
1313	0	Tr_max_face = max(Tr(m)%t(i,j,1), Tr(m)%t(i+1,j,1))
1314	0	do k=2,nkmb
1315	0	Tr_min_face = min(Tr_min_face, Tr(m)%t(i,j,k), Tr(m)%t(i+1,j,k))
1316	0	Tr_max_face = max(Tr_max_face, Tr(m)%t(i,j,k), Tr(m)%t(i+1,j,k))
1317		enddo
1318
1319		! Include the next two layers denser than the densest buffer layer.
1320	0	kLa = nkmb+1 ; if (max_kRho(i,j) < nz+1) kLa = max_kRho(i,j)
1321	0	kLb = kLa ; if (max_kRho(i,j) < nz) kLb = max_kRho(i,j)+1
1322	0	kRa = nkmb+1 ; if (max_kRho(i+1,j) < nz+1) kRa = max_kRho(i+1,j)
1323	0	kRb = kRa ; if (max_kRho(i+1,j) < nz) kRb = max_kRho(i+1,j)+1
1324	0	Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La
1325	0	if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa)
1326	0	if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then
1327	0	if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
1328		else
1329	0	if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb)
1330		endif
1331	0	if (h(i+1,j,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i+1,j,kRa)
1332	0	if (h(i+1,j,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i+1,j,kRb)
1333	0	Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1334	0	Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1335
1336		! Include all points in diffusive pairings at this face.
1337	0	do k=1,nPu(I,j)
1338	0	Tr_Lb = Tr(m)%t(i,j,k0b_Lu(I,k,j))
1339	0	Tr_Rb = Tr(m)%t(i+1,j,k0b_Ru(I,k,j))
1340	0	Tr_La = Tr_Lb ; Tr_Ra = Tr_Rb
1341	0	if (deep_wt_Lu(I,k,j) < 1.0) Tr_La = Tr(m)%t(i,j,k0a_Lu(I,k,j))
1342	0	if (deep_wt_Ru(I,k,j) < 1.0) Tr_Ra = Tr(m)%t(i+1,j,k0a_Ru(I,k,j))
1343	0	Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1344	0	Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1345		enddo
1346		endif
1347
1348	0	do k=1,nPu(I,j)
1349	0	kLb = k0b_Lu(I,k,j) ; Tr_Lb = Tr(m)%t(i,j,kLb) ; Tr_av_L = Tr_Lb
1350	0	if (deep_wt_Lu(I,k,j) < 1.0) then
1351	0	kLa = k0a_Lu(I,k,j) ; Tr_La = Tr(m)%t(i,j,kLa)
1352	0	wt_b = deep_wt_Lu(I,k,j)
1353	0	Tr_av_L = wt_bTr_Lb + (1.0-wt_b)Tr_La
1354		endif
1355
1356	0	kRb = k0b_Ru(I,k,j) ; Tr_Rb = Tr(m)%t(i+1,j,kRb) ; Tr_av_R = Tr_Rb
1357	0	if (deep_wt_Ru(I,k,j) < 1.0) then
1358	0	kRa = k0a_Ru(I,k,j) ; Tr_Ra = Tr(m)%t(i+1,j,kRa)
1359	0	wt_b = deep_wt_Ru(I,k,j)
1360	0	Tr_av_R = wt_bTr_Rb + (1.0-wt_b)Tr_Ra
1361		endif
1362
1363	0	h_L = hP_Lu(I,k,j) ; h_R = hP_Ru(I,k,j)
1364	0	if (CS%answer_date <= 20240330) then
1365		Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * &
1366	0	((2.0 * h_L * h_R) / (h_L + h_R))
1367		else
1368		Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * &
1369	0	khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R)
1370		endif
1371
1372	0	if (deep_wt_Lu(I,k,j) >= 1.0) then
1373	0	if (CS%answer_date <= 20240330) then
1374	0	tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux
1375		else
1376	0	tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + Tr_flux
1377		endif
1378		else
1379	0	Tr_adj_vert = 0.0
1380	0	wt_b = deep_wt_Lu(I,k,j) ; wt_a = 1.0 - wt_b
1381	0	vol = hP_Lu(I,k,j) * G%areaT(i,j)
1382
1383		! Ensure that the tracer flux does not drive the tracer values
1384		! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1385		! does that the concentration in both contributing pieces exceed
1386		! this range equally. With down-gradient fluxes and the initial tracer
1387		! concentrations determining the valid range, the latter condition
1388		! only enters for large values of the effective diffusive CFL number.
1389	0	if (Tr_flux > 0.0) then
1390	0	if (Tr_La < Tr_Lb) then ; if (vol*(Tr_La-Tr_min_face) < Tr_flux) &
1391		Tr_adj_vert = -wt_a * min(Tr_flux - vol * (Tr_La-Tr_min_face), &
1392	0	(volwt_b) (Tr_Lb - Tr_La))
1393	0	else ; if (vol*(Tr_Lb-Tr_min_face) < Tr_flux) &
1394		Tr_adj_vert = wt_b * min(Tr_flux - vol * (Tr_Lb-Tr_min_face), &
1395	0	(volwt_a) (Tr_La - Tr_Lb))
1396		endif
1397	0	elseif (Tr_flux < 0.0) then
1398	0	if (Tr_La > Tr_Lb) then ; if (vol * (Tr_max_face-Tr_La) < -Tr_flux) &
1399		Tr_adj_vert = wt_a * min(-Tr_flux - vol * (Tr_max_face-Tr_La), &
1400	0	(volwt_b) (Tr_La - Tr_Lb))
1401	0	else ; if (vol*(Tr_max_face-Tr_Lb) < -Tr_flux) &
1402		Tr_adj_vert = -wt_b * min(-Tr_flux - vol * (Tr_max_face-Tr_Lb), &
1403	0	(volwt_a)(Tr_Lb - Tr_La))
1404		endif
1405		endif
1406
1407	0	if (CS%answer_date <= 20240330) then
1408	0	tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert)
1409	0	tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert)
1410		else
1411	0	tr_flux_E(i,j,kLa) = tr_flux_E(i,j,kLa) + (wt_a*Tr_flux + Tr_adj_vert)
1412	0	tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + (wt_b*Tr_flux - Tr_adj_vert)
1413		endif
1414		endif
1415
1416	0	if (deep_wt_Ru(I,k,j) >= 1.0) then
1417	0	if (CS%answer_date <= 20240330) then
1418	0	tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux
1419		else
1420	0	tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + Tr_flux
1421		endif
1422		else
1423	0	Tr_adj_vert = 0.0
1424	0	wt_b = deep_wt_Ru(I,k,j) ; wt_a = 1.0 - wt_b
1425	0	vol = hP_Ru(I,k,j) * G%areaT(i+1,j)
1426
1427		! Ensure that the tracer flux does not drive the tracer values
1428		! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1429		! does that the concentration in both contributing pieces exceed
1430		! this range equally. With down-gradient fluxes and the initial tracer
1431		! concentrations determining the valid range, the latter condition
1432		! only enters for large values of the effective diffusive CFL number.
1433	0	if (Tr_flux < 0.0) then
1434	0	if (Tr_Ra < Tr_Rb) then ; if (vol * (Tr_Ra-Tr_min_face) < -Tr_flux) &
1435		Tr_adj_vert = -wt_a * min(-Tr_flux - vol * (Tr_Ra-Tr_min_face), &
1436	0	(volwt_b) (Tr_Rb - Tr_Ra))
1437	0	else ; if (vol*(Tr_Rb-Tr_min_face) < (-Tr_flux)) &
1438		Tr_adj_vert = wt_b * min(-Tr_flux - vol * (Tr_Rb-Tr_min_face), &
1439	0	(volwt_a) (Tr_Ra - Tr_Rb))
1440		endif
1441	0	elseif (Tr_flux > 0.0) then
1442	0	if (Tr_Ra > Tr_Rb) then ; if (vol * (Tr_max_face-Tr_Ra) < Tr_flux) &
1443		Tr_adj_vert = wt_a * min(Tr_flux - vol * (Tr_max_face-Tr_Ra), &
1444	0	(volwt_b) (Tr_Ra - Tr_Rb))
1445	0	else ; if (vol*(Tr_max_face-Tr_Rb) < Tr_flux) &
1446		Tr_adj_vert = -wt_b * min(Tr_flux - vol * (Tr_max_face-Tr_Rb), &
1447	0	(volwt_a)(Tr_Rb - Tr_Ra))
1448		endif
1449		endif
1450
1451	0	if (CS%answer_date <= 20240330) then
1452	0	tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert)
1453	0	tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert)
1454		else
1455	0	tr_flux_W(i+1,j,kRa) = tr_flux_W(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert)
1456	0	tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert)
1457		endif
1458		endif
1459	0	if (associated(Tr(m)%df2d_x)) &
1460	0	Tr(m)%df2d_x(I,j) = Tr(m)%df2d_x(I,j) + Tr_flux * Idt
1461		enddo ! Loop over pairings at faces.
1462		endif ; enddo ; enddo ! i- & j- loops over zonal faces.
1463
1464		! this gives wrong result when using do concurrent on NVHPC 25.9
1465		!$omp target teams loop collapse(2) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, &
1466		!$omp Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, &
1467		!$omp h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
1468	0	do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
1469		! Determine the fluxes through the meridional faces.
1470
1471		! Find the acceptable range of tracer concentration around this face.
1472	0	if (nPv(i,J) >= 1) then
1473	0	Tr_min_face = min(Tr(m)%t(i,j,1), Tr(m)%t(i,j+1,1))
1474	0	Tr_max_face = max(Tr(m)%t(i,j,1), Tr(m)%t(i,j+1,1))
1475	0	do k=2,nkmb
1476	0	Tr_min_face = min(Tr_min_face, Tr(m)%t(i,j,k), Tr(m)%t(i,j+1,k))
1477	0	Tr_max_face = max(Tr_max_face, Tr(m)%t(i,j,k), Tr(m)%t(i,j+1,k))
1478		enddo
1479
1480		! Include the next two layers denser than the densest buffer layer.
1481	0	kLa = nkmb+1 ; if (max_kRho(i,j) < nz+1) kLa = max_kRho(i,j)
1482	0	kLb = kLa ; if (max_kRho(i,j) < nz) kLb = max_kRho(i,j)+1
1483	0	kRa = nkmb+1 ; if (max_kRho(i,j+1) < nz+1) kRa = max_kRho(i,j+1)
1484	0	kRb = kRa ; if (max_kRho(i,j+1) < nz) kRb = max_kRho(i,j+1)+1
1485	0	Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La
1486	0	if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa)
1487	0	if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then
1488	0	if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
1489		else
1490	0	if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb)
1491		endif
1492	0	if (h(i,j+1,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i,j+1,kRa)
1493	0	if (h(i,j+1,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i,j+1,kRb)
1494	0	Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1495	0	Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1496
1497		! Include all points in diffusive pairings at this face.
1498	0	do k=1,nPv(i,J)
1499	0	Tr_Lb = Tr(m)%t(i,j,k0b_Lv(i,k,J)) ; Tr_Rb = Tr(m)%t(i,j+1,k0b_Rv(i,k,J))
1500	0	Tr_La = Tr_Lb ; Tr_Ra = Tr_Rb
1501	0	if (deep_wt_Lv(I,k,J) < 1.0) Tr_La = Tr(m)%t(i,j,k0a_Lv(i,k,J))
1502	0	if (deep_wt_Rv(I,k,J) < 1.0) Tr_Ra = Tr(m)%t(i,j+1,k0a_Rv(i,k,J))
1503	0	Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1504	0	Tr_max_face = max(Tr_max_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
1505		enddo
1506		endif
1507
1508	0	do k=1,nPv(i,J)
1509	0	kLb = k0b_Lv(i,k,J) ; Tr_Lb = Tr(m)%t(i,j,kLb) ; Tr_av_L = Tr_Lb
1510	0	if (deep_wt_Lv(I,k,J) < 1.0) then
1511	0	kLa = k0a_Lv(i,k,J) ; Tr_La = Tr(m)%t(i,j,kLa)
1512	0	wt_b = deep_wt_Lv(I,k,J)
1513	0	Tr_av_L = wt_b * Tr_Lb + (1.0-wt_b) * Tr_La
1514		endif
1515
1516	0	kRb = k0b_Rv(i,k,J) ; Tr_Rb = Tr(m)%t(i,j+1,kRb) ; Tr_av_R = Tr_Rb
1517	0	if (deep_wt_Rv(I,k,J) < 1.0) then
1518	0	kRa = k0a_Rv(i,k,J) ; Tr_Ra = Tr(m)%t(i,j+1,kRa)
1519	0	wt_b = deep_wt_Rv(I,k,J)
1520	0	Tr_av_R = wt_b * Tr_Rb + (1.0-wt_b) * Tr_Ra
1521		endif
1522
1523	0	h_L = hP_Lv(I,k,J) ; h_R = hP_Rv(I,k,J)
1524		Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * &
1525	0	khdt_epi_y(i,J) * (Tr_av_L - Tr_av_R)
1526	0	Tr_flux_3d(i,J,k) = Tr_flux
1527
1528	0	if (deep_wt_Lv(I,k,J) < 1.0) then
1529	0	Tr_adj_vert = 0.0
1530	0	wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b
1531	0	vol = hP_Lv(I,k,J) * G%areaT(i,j)
1532
1533		! Ensure that the tracer flux does not drive the tracer values
1534		! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1535	0	if (Tr_flux > 0.0) then
1536	0	if (Tr_La < Tr_Lb) then ; if (vol * (Tr_La-Tr_min_face) < Tr_flux) &
1537		Tr_adj_vert = -wt_a * min(Tr_flux - vol * (Tr_La-Tr_min_face), &
1538	0	(volwt_b) (Tr_Lb - Tr_La))
1539	0	else ; if (vol*(Tr_Lb-Tr_min_face) < Tr_flux) &
1540		Tr_adj_vert = wt_b * min(Tr_flux - vol * (Tr_Lb-Tr_min_face), &
1541	0	(volwt_a) (Tr_La - Tr_Lb))
1542		endif
1543	0	elseif (Tr_flux < 0.0) then
1544	0	if (Tr_La > Tr_Lb) then ; if (vol * (Tr_max_face-Tr_La) < -Tr_flux) &
1545		Tr_adj_vert = wt_a * min(-Tr_flux - vol * (Tr_max_face-Tr_La), &
1546	0	(volwt_b) (Tr_La - Tr_Lb))
1547	0	else ; if (vol*(Tr_max_face-Tr_Lb) < -Tr_flux) &
1548		Tr_adj_vert = -wt_b * min(-Tr_flux - vol * (Tr_max_face-Tr_Lb), &
1549	0	(volwt_a)(Tr_Lb - Tr_La))
1550		endif
1551		endif
1552	0	Tr_adj_vert_L(i,J,k) = Tr_adj_vert
1553		endif
1554
1555	0	if (deep_wt_Rv(I,k,J) < 1.0) then
1556	0	Tr_adj_vert = 0.0
1557	0	wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b
1558	0	vol = hP_Rv(I,k,J) * G%areaT(i,j+1)
1559
1560		! Ensure that the tracer flux does not drive the tracer values
1561		! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1562	0	if (Tr_flux < 0.0) then
1563	0	if (Tr_Ra < Tr_Rb) then ; if (vol * (Tr_Ra-Tr_min_face) < -Tr_flux) &
1564		Tr_adj_vert = -wt_a * min(-Tr_flux - vol * (Tr_Ra-Tr_min_face), &
1565	0	(volwt_b) (Tr_Rb - Tr_Ra))
1566	0	else ; if (vol*(Tr_Rb-Tr_min_face) < (-Tr_flux)) &
1567		Tr_adj_vert = wt_b * min(-Tr_flux - vol * (Tr_Rb-Tr_min_face), &
1568	0	(volwt_a) (Tr_Ra - Tr_Rb))
1569		endif
1570	0	elseif (Tr_flux > 0.0) then
1571	0	if (Tr_Ra > Tr_Rb) then ; if (vol * (Tr_max_face-Tr_Ra) < Tr_flux) &
1572		Tr_adj_vert = wt_a * min(Tr_flux - vol * (Tr_max_face-Tr_Ra), &
1573	0	(volwt_b) (Tr_Ra - Tr_Rb))
1574	0	else ; if (vol*(Tr_max_face-Tr_Rb) < Tr_flux) &
1575		Tr_adj_vert = -wt_b * min(Tr_flux - vol * (Tr_max_face-Tr_Rb), &
1576	0	(volwt_a)(Tr_Rb - Tr_Ra))
1577		endif
1578		endif
1579	0	Tr_adj_vert_R(i,J,k) = Tr_adj_vert
1580		endif
1581	0	if (associated(Tr(m)%df2d_y)) &
1582	0	Tr(m)%df2d_y(i,J) = Tr(m)%df2d_y(i,J) + Tr_flux * Idt
1583		enddo ! Loop over pairings at faces.
1584		endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1585
1586		! The non-stride-1 loop order here is to facilitate openMP threading. However, it might be
1587		! suboptimal when openMP threading is not used, at which point it might be better to fuse
1588		! this loop with those that precede it and thereby eliminate the need for three 3-d arrays.
1589	0	if (CS%answer_date <= 20240330) then
1590		! KRa/b aren't guaranteed to be unique and update of tr_flux_conv(:, j/j+1, :) means this
1591		! loop must be serial in both j and k. i can be paralellised.
1592	0	do concurrent (i=is:ie) DO_LOCALITY(local(j, k, kLb, kRb, kLa, wt_b, wt_a, kRa))
1593	0	do J=js-1,je ; if (G%mask2dCv(i,J) > 0.0) then
1594	0	do k=1,nPv(i,J)
1595	0	kLb = k0b_Lv(i,k,J) ; kRb = k0b_Rv(i,k,J)
1596	0	if (deep_wt_Lv(I,k,J) >= 1.0) then
1597	0	tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k)
1598		else
1599	0	kLa = k0a_Lv(i,k,J)
1600	0	wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b
1601	0	tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k))
1602	0	tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k))
1603		endif
1604	0	if (deep_wt_Rv(I,k,J) >= 1.0) then
1605	0	tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k)
1606		else
1607	0	kRa = k0a_Rv(i,k,J)
1608	0	wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b
1609		tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + &
1610	0	(wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k))
1611		tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + &
1612	0	(wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k))
1613		endif
1614		enddo
1615		endif ; enddo
1616		enddo
1617		else
1618		! Update of tr_flux_N/S can be done independently so both i and j can be parallelised.
1619	0	do concurrent (J=js-1:je, i=is:ie, G%mask2dCv(i,J) > 0.0) DO_LOCALITY(local(k, kLb, kRb, kLa, wt_b, wt_a, kRa))
1620	0	do k=1,nPv(i,J)
1621	0	kLb = k0b_Lv(i,k,J) ; kRb = k0b_Rv(i,k,J)
1622	0	if (deep_wt_Lv(I,k,J) >= 1.0) then
1623	0	tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + Tr_flux_3d(i,J,k)
1624		else
1625	0	kLa = k0a_Lv(i,k,J)
1626	0	wt_b = deep_wt_Lv(I,k,J) ; wt_a = 1.0 - wt_b
1627	0	tr_flux_N(i,j,kLa) = tr_flux_N(i,j,kLa) + (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k))
1628	0	tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k))
1629		endif
1630	0	if (deep_wt_Rv(I,k,J) >= 1.0) then
1631	0	tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + Tr_flux_3d(i,J,k)
1632		else
1633	0	kRa = k0a_Rv(i,k,J)
1634	0	wt_b = deep_wt_Rv(I,k,J) ; wt_a = 1.0 - wt_b
1635	0	tr_flux_S(i,j+1,kRa) = tr_flux_S(i,j+1,kRa) + (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k))
1636	0	tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k))
1637		endif
1638		enddo
1639		enddo
1640		endif
1641
1642	0	if (CS%answer_date >= 20240331) then
1643	0	do concurrent (k=1:PEmax_kRho, j=js:je, i=is:ie)
1644		tr_flux_conv(i,j,k) = ((tr_flux_W(i,j,k) - tr_flux_E(i,j,k)) + &
1645	0	(tr_flux_S(i,j,k) - tr_flux_N(i,j,k)))
1646		enddo
1647		endif
1648
1649	0	do concurrent (k=1:PEmax_kRho, j=js:je, i=is:ie, (G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0))
1650	0	Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / (h(i,j,k)*G%areaT(i,j))
1651		enddo
1652
1653		! Do user controlled underflow of the tracer concentrations.
1654	0	if (Tr(m)%conc_underflow > 0.0) then
1655		!$omp target update from(Tr(m)%t)
1656		!$OMP parallel do default(shared)
1657	0	do k=1,nz ; do j=js,je ; do i=is,ie
1658	0	if (abs(Tr(m)%t(i,j,k)) < Tr(m)%conc_underflow) Tr(m)%t(i,j,k) = 0.0
1659		enddo ; enddo ; enddo
1660		!$omp target update to(Tr(m)%t)
1661		endif
1662
1663		enddo ! Loop over tracers
1664		enddo ! Loop over iterations
1665
1666		!$omp target exit data map(release: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R, tr_flux_N, &
1667		!$omp tr_flux_S, tr_flux_E, tr_flux_W, tr_flux_conv, nPv, nPu, max_kRho, rho_srt, num_srt, &
1668		!$omp max_srt)
1669
1670		!$omp target exit data map(release: deep_wt_Lu, deep_wt_Ru, hP_Lu, hP_Ru, k0a_Lu, k0a_Ru, &
1671		!$omp k0b_Lu, k0b_Ru)
1672
1673		!$omp target exit data map(release: deep_wt_Lv, deep_wt_Rv, hP_Lv, hP_Rv, k0a_Lv, k0a_Rv, &
1674		!$Omp k0b_Lv, k0b_Rv)
1675
1676	0	end subroutine tracer_epipycnal_ML_diff
1677
1678
1679		!> Initialize lateral tracer diffusion module
1680	1	subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic_CSp, CS)
1681		type(time_type), target, intent(in) :: Time !< current model time
1682		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1683		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1684		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1685		type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control
1686		type(EOS_type), target, intent(in) :: EOS !< Equation of state CS
1687		type(diabatic_CS), pointer, intent(in) :: diabatic_CSp !< Equation of state CS
1688		type(param_file_type), intent(in) :: param_file !< parameter file
1689		type(tracer_hor_diff_CS), pointer :: CS !< horz diffusion control structure
1690
1691		! This include declares and sets the variable "version".
1692		# include "version_variable.h"
1693		character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name.
1694		integer :: default_answer_date
1695
1696	1	if (associated(CS)) then
1697	0	call MOM_error(WARNING, "tracer_hor_diff_init called with associated control structure.")
1698	0	return
1699		endif
1700	1	allocate(CS)
1701
1702	1	CS%diag => diag
1703	1	CS%show_call_tree = callTree_showQuery()
1704
1705		! Read all relevant parameters and write them to the model log.
1706	1	call log_version(param_file, mdl, version, "")
1707		call get_param(param_file, mdl, "KHTR", CS%KhTr, &
1708		"The background along-isopycnal tracer diffusivity.", &
1709	1	units="m2 s-1", default=0.0, scale=US%m_to_L*2US%T_to_s)
1710		! call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", CS%KhTh_use_ebt_struct, &
1711		! "If true, uses the equivalent barotropic structure "//&
1712		! "as the vertical structure of the tracer diffusivity.",&
1713		! default=.false.)
1714		call get_param(param_file, mdl, "KHTR_SLOPE_CFF", CS%KhTr_Slope_Cff, &
1715		"The scaling coefficient for along-isopycnal tracer "//&
1716		"diffusivity using a shear-based (Visbeck-like) "//&
1717		"parameterization. A non-zero value enables this param.", &
1718	1	units="nondim", default=0.0)
1719		call get_param(param_file, mdl, "KHTR_MIN", CS%KhTr_Min, &
1720		"The minimum along-isopycnal tracer diffusivity.", &
1721	1	units="m2 s-1", default=0.0, scale=US%m_to_L*2US%T_to_s)
1722	1	if (CS%KhTr_Min > 0.0) then
1723		call get_param(param_file, mdl, "FULL_DEPTH_KHTR_MIN", CS%full_depth_khtr_min, &
1724		"If true, KHTR_MIN is enforced throughout the whole water column. "//&
1725		"Otherwise, KHTR_MIN is only enforced at the surface. This parameter "//&
1726		"is only available when KHTR_USE_EBT_STRUCT=True and KHTR_MIN>0.", &
1727	0	default=.false.)
1728		endif
1729		call get_param(param_file, mdl, "KHTR_MAX", CS%KhTr_Max, &
1730		"The maximum along-isopycnal tracer diffusivity.", &
1731	1	units="m2 s-1", default=0.0, scale=US%m_to_L*2US%T_to_s)
1732		call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", CS%KhTr_passivity_coeff, &
1733		"The coefficient that scales deformation radius over "//&
1734		"grid-spacing in passivity, where passivity is the ratio "//&
1735		"between along isopycnal mixing of tracers to thickness mixing. "//&
1736		"A non-zero value enables this parameterization.", &
1737	1	units="nondim", default=0.0)
1738		call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", CS%KhTr_passivity_min, &
1739		"The minimum passivity which is the ratio between "//&
1740		"along isopycnal mixing of tracers to thickness mixing.", &
1741	1	units="nondim", default=0.5)
1742		call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", CS%Diffuse_ML_interior, &
1743		"If true, enable epipycnal mixing between the surface "//&
1744	1	"boundary layer and the interior.", default=.false.)
1745		call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", CS%check_diffusive_CFL, &
1746		"If true, use enough iterations the diffusion to ensure "//&
1747		"that the diffusive equivalent of the CFL limit is not "//&
1748		"violated. If false, always use the greater of 1 or "//&
1749	1	"MAX_TR_DIFFUSION_CFL iteration.", default=.false.)
1750		call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", CS%max_diff_CFL, &
1751		"If positive, locally limit the along-isopycnal tracer "//&
1752		"diffusivity to keep the diffusive CFL locally at or "//&
1753		"below this value. The number of diffusive iterations "//&
1754		"is often this value or the next greater integer.", &
1755	1	units="nondim", default=-1.0)
1756		call get_param(param_File, mdl, "RECALC_NEUTRAL_SURF", CS%recalc_neutral_surf, &
1757		"If true, then recalculate the neutral surfaces if the \n"//&
1758		"diffusive CFL is exceeded. If false, assume that the \n"//&
1759	1	"positions of the surfaces do not change \n", default=.false.)
1760		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1761		"This sets the default value for the various _ANSWER_DATE parameters.", &
1762	1	default=99991231, do_not_log=.true.)
1763		call get_param(param_file, mdl, "HOR_DIFF_ANSWER_DATE", CS%answer_date, &
1764		"The vintage of the order of arithmetic to use for the tracer diffusion. "//&
1765		"Values of 20240330 or below recover the answers from the original form of the "//&
1766		"along-isopycnal mixed layer to interior mixing code, while higher values use "//&
1767		"mathematically equivalent expressions that recover rotational symmetry "//&
1768		"when DIFFUSE_ML_TO_INTERIOR is true.", &
1769	1	default=default_answer_date, do_not_log=.not.CS%Diffuse_ML_interior)
1770		call get_param(param_file, mdl, "HOR_DIFF_LIMIT_BUG", CS%limit_bug, &
1771		"If true and the answer date is 20240330 or below, use a rotational symmetry "//&
1772		"breaking bug when limiting the tracer properties in tracer_epipycnal_ML_diff.", &
1773	1	default=.false., do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331)))
1774	1	CS%ML_KhTR_scale = 1.0
1775	1	if (CS%Diffuse_ML_interior) then
1776		call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, &
1777		"With Diffuse_ML_interior, the ratio of the truly "//&
1778		"horizontal diffusivity in the mixed layer to the "//&
1779		"epipycnal diffusivity. The valid range is 0 to 1.", &
1780	0	units="nondim", default=1.0)
1781		endif
1782
1783		CS%use_neutral_diffusion = neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, &
1784	1	diabatic_CSp, CS%neutral_diffusion_CSp )
1785	1	if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// &
1786	0	"USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1787		CS%use_hor_bnd_diffusion = hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diabatic_CSp, &
1788	1	CS%hor_bnd_diffusion_CSp)
1789	1	if (CS%use_hor_bnd_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// &
1790	0	"USE_HORIZONTAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1791
1792	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
1793
1794	1	id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=CLOCK_MODULE)
1795	1	id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=CLOCK_MODULE)
1796	1	id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=CLOCK_ROUTINE)
1797	1	id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=CLOCK_ROUTINE)
1798
1799	1	CS%id_KhTr_u = -1
1800	1	CS%id_KhTr_v = -1
1801	1	CS%id_KhTr_h = -1
1802	1	CS%id_CFL = -1
1803
1804		CS%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCui, Time, &
1805	1	'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1806		CS%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCvi, Time, &
1807	1	'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T)
1808		CS%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesTi, Time, &
1809		'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1', conversion=US%L_to_m*2US%s_to_T, &
1810		cmor_field_name='diftrelo', &
1811		cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', &
1812	1	cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')
1813
1814		CS%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, Time, &
1815	1	'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'm2', conversion=US%L_to_m**2)
1816		CS%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, Time, &
1817	1	'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2', conversion=US%L_to_m**2)
1818	1	if (CS%check_diffusive_CFL) then
1819		CS%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, Time,&
1820	1	'Grid CFL number for lateral/neutral tracer diffusion', 'nondim')
1821		endif
1822
1823
1824		end subroutine tracer_hor_diff_init
1825
1826	1	subroutine tracer_hor_diff_end(CS)
1827		type(tracer_hor_diff_CS), pointer :: CS !< module control structure
1828
1829	1	call neutral_diffusion_end(CS%neutral_diffusion_CSp)
1830	1	call hor_bnd_diffusion_end(CS%hor_bnd_diffusion_CSp)
1831	1	if (associated(CS)) deallocate(CS)
1832
1833	1	end subroutine tracer_hor_diff_end
1834
1835
1836		!> \namespace mom_tracer_hor_diff
1837		!!
1838		!! \section section_intro Introduction to the module
1839		!!
1840		!! This module contains subroutines that handle horizontal
1841		!! diffusion (i.e., isoneutral or along layer) of tracers.
1842		!!
1843		!! Each of the tracers are subject to Fickian along-coordinate
1844		!! diffusion if Khtr is defined and positive. The tracer diffusion
1845		!! can use a suitable number of iterations to guarantee stability
1846		!! with an arbitrarily large time step.
1847
1848	0	end module MOM_tracer_hor_diff