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