MOM6
MOM_tracer_hor_diff.F90
1 !> Main routine for lateral (along surface or neutral) diffusion of tracers
2 module mom_tracer_hor_diff
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
7 use mom_cpu_clock, only : clock_module, clock_routine
8 use mom_diag_mediator, only : post_data, diag_ctrl
9 use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
10 use mom_domains, only : sum_across_pes, max_across_pes
11 use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
12 use mom_domains, only : pass_vector
13 use mom_debugging, only : hchksum, uvchksum
14 use mom_diabatic_driver, only : diabatic_cs
15 use mom_eos, only : calculate_density, eos_type, eos_domain
16 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
17 use mom_error_handler, only : mom_set_verbosity, calltree_showquery
18 use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint
19 use mom_file_parser, only : get_param, log_version, param_file_type
20 use mom_grid, only : ocean_grid_type
21 use mom_lateral_mixing_coeffs, only : varmix_cs
22 use mom_meke_types, only : meke_type
23 use mom_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end
24 use mom_neutral_diffusion, only : neutral_diffusion_cs
25 use mom_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion
26 use mom_lateral_boundary_diffusion, only : lateral_boundary_diffusion_cs, lateral_boundary_diffusion_init
27 use mom_lateral_boundary_diffusion, only : lateral_boundary_diffusion
28 use mom_tracer_registry, only : tracer_registry_type, tracer_type, mom_tracer_chksum
29 use mom_unit_scaling, only : unit_scale_type
30 use mom_variables, only : thermo_var_ptrs
31 use mom_verticalgrid, only : verticalgrid_type
32 
33 implicit none ; private
34 
35 #include <MOM_memory.h>
36 
37 public tracer_hordiff, tracer_hor_diff_init, tracer_hor_diff_end
38 
39 !> The control structure for along-layer and epineutral tracer diffusion
40 type, public :: tracer_hor_diff_cs ; private
41  real :: khtr !< The along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
42  real :: khtr_slope_cff !< The non-dimensional coefficient in KhTr formula [nondim]
43  real :: khtr_min !< Minimum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
44  real :: khtr_max !< Maximum along-isopycnal tracer diffusivity [L2 T-1 ~> m2 s-1].
45  real :: khtr_passivity_coeff !< Passivity coefficient that scales Rd/dx (default = 0)
46  !! where passivity is the ratio between along-isopycnal
47  !! tracer mixing and thickness mixing [nondim]
48  real :: khtr_passivity_min !< Passivity minimum (default = 1/2) [nondim]
49  real :: ml_khtr_scale !< With Diffuse_ML_interior, the ratio of the
50  !! truly horizontal diffusivity in the mixed
51  !! layer to the epipycnal diffusivity [nondim].
52  real :: max_diff_cfl !< If positive, locally limit the along-isopycnal
53  !! tracer diffusivity to keep the diffusive CFL
54  !! locally at or below this value [nondim].
55  logical :: diffuse_ml_interior !< If true, diffuse along isopycnals between
56  !! the mixed layer and the interior.
57  logical :: check_diffusive_cfl !< If true, automatically iterate the diffusion
58  !! to ensure that the diffusive equivalent of
59  !! the CFL limit is not violated.
60  logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within
61  !! tracer_hor_diff.
62  logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within
63  !! tracer_hor_diff.
64  logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been
65  !! exceeded
66  type(neutral_diffusion_cs), pointer :: neutral_diffusion_csp => null() !< Control structure for neutral diffusion.
67  type(lateral_boundary_diffusion_cs), pointer :: lateral_boundary_diffusion_csp => null() !< Control structure for
68  !! lateral boundary mixing.
69  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
70  !! regulate the timing of diagnostic output.
71  logical :: debug !< If true, write verbose checksums for debugging purposes.
72  logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level.
73  logical :: first_call = .true. !< This is true until after the first call
74  !>@{ Diagnostic IDs
75  integer :: id_khtr_u = -1
76  integer :: id_khtr_v = -1
77  integer :: id_khtr_h = -1
78  integer :: id_cfl = -1
79  integer :: id_khdt_x = -1
80  integer :: id_khdt_y = -1
81  !>@}
82 
83  type(group_pass_type) :: pass_t !< For group halo pass, used in both
84  !! tracer_hordiff and tracer_epipycnal_ML_diff
85 end type tracer_hor_diff_cs
86 
87 !> A type that can be used to create arrays of pointers to 2D arrays
88 type p2d
89  real, dimension(:,:), pointer :: p => null() !< A pointer to a 2D array of reals
90 end type p2d
91 !> A type that can be used to create arrays of pointers to 2D integer arrays
92 type p2di
93  integer, dimension(:,:), pointer :: p => null() !< A pointer to a 2D array of integers
94 end type p2di
95 
96 !>@{ CPU time clocks
97 integer :: id_clock_diffuse, id_clock_epimix, id_clock_pass, id_clock_sync
98 !>@}
99 
100 contains
101 
102 !> Compute along-coordinate diffusion of all tracers
103 !! using the diffusivity in CS%KhTr, or using space-dependent diffusivity.
104 !! Multiple iterations are used (if necessary) so that there is no limit
105 !! on the acceptable time increment.
106 subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y)
107  type(ocean_grid_type), intent(inout) :: g !< Grid type
108  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
109  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
110  real, intent(in) :: dt !< time step [T ~> s]
111  type(meke_type), pointer :: meke !< MEKE type
112  type(varmix_cs), pointer :: varmix !< Variable mixing type
113  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
114  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
115  type(tracer_hor_diff_cs), pointer :: cs !< module control structure
116  type(tracer_registry_type), pointer :: reg !< registered tracers
117  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
118  !! thermodynamic fields, including potential temp and
119  !! salinity or mixed layer density. Absent fields have
120  !! NULL ptrs, and these may (probably will) point to
121  !! some of the same arrays as Tr does. tv is required
122  !! for epipycnal mixing between mixed layer and the interior.
123  ! Optional inputs for offline tracer transport
124  logical, optional, intent(in) :: do_online_flag !< If present and true, do online
125  !! tracer transport with stored velocities.
126  real, dimension(SZIB_(G),SZJ_(G)), &
127  optional, intent(in) :: read_khdt_x !< If present, these are the zonal
128  !! diffusivities from previous run.
129  real, dimension(SZI_(G),SZJB_(G)), &
130  optional, intent(in) :: read_khdt_y !< If present, these are the meridional
131  !! diffusivities from previous run.
132 
133 
134  real, dimension(SZI_(G),SZJ_(G)) :: &
135  ihdxdy, & ! The inverse of the volume or mass of fluid in a layer in a
136  ! grid cell [H-1 L-2 ~> m-3 or kg-1].
137  kh_h, & ! The tracer diffusivity averaged to tracer points [L2 T-1 ~> m2 s-1].
138  cfl, & ! A diffusive CFL number for each cell [nondim].
139  dtr ! The change in a tracer's concentration, in units of concentration [Conc].
140 
141  real, dimension(SZIB_(G),SZJ_(G)) :: &
142  khdt_x, & ! The value of Khtr*dt times the open face width divided by
143  ! the distance between adjacent tracer points [L2 ~> m2].
144  coef_x, & ! The coefficients relating zonal tracer differences to time-integrated
145  ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others.
146  kh_u ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
147  real, dimension(SZI_(G),SZJB_(G)) :: &
148  khdt_y, & ! The value of Khtr*dt times the open face width divided by
149  ! the distance between adjacent tracer points [L2 ~> m2].
150  coef_y, & ! The coefficients relating meridional tracer differences to time-integrated
151  ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others.
152  kh_v ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1].
153 
154  real :: khdt_max ! The local limiting value of khdt_x or khdt_y [L2 ~> m2].
155  real :: max_cfl ! The global maximum of the diffusive CFL number.
156  logical :: use_varmix, resoln_scaled, do_online, use_eady
157  integer :: s_idx, t_idx ! Indices for temperature and salinity if needed
158  integer :: i, j, k, m, is, ie, js, je, nz, ntr, itt, num_itts
159  real :: i_numitts ! The inverse of the number of iterations, num_itts.
160  real :: scale ! The fraction of khdt_x or khdt_y that is applied in this
161  ! layer for this iteration [nondim].
162  real :: idt ! The inverse of the time step [T-1 ~> s-1].
163  real :: h_neglect ! A thickness that is so small it is usually lost
164  ! in roundoff and can be neglected [H ~> m or kg m-2].
165  real :: kh_loc ! The local value of Kh [L2 T-1 ~> m2 s-1].
166  real :: res_fn ! The local value of the resolution function [nondim].
167  real :: rd_dx ! The local value of deformation radius over grid-spacing [nondim].
168  real :: normalize ! normalization used for diagnostic Kh_h; diffusivity averaged to h-points.
169 
170  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
171 
172  do_online = .true.
173  if (present(do_online_flag)) do_online = do_online_flag
174 
175  if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
176  "register_tracer must be called before tracer_hordiff.")
177  if (loc(reg)==0) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
178  "register_tracer must be called before tracer_hordiff.")
179  if ((reg%ntr==0) .or. ((cs%KhTr <= 0.0) .and. .not.associated(varmix)) ) return
180 
181  if (cs%show_call_tree) call calltree_enter("tracer_hordiff(), MOM_tracer_hor_diff.F90")
182 
183  call cpu_clock_begin(id_clock_diffuse)
184 
185  ntr = reg%ntr
186  idt = 1.0 / dt
187  h_neglect = gv%H_subroundoff
188 
189  if (cs%Diffuse_ML_interior .and. cs%first_call) then ; if (is_root_pe()) then
190  do m=1,ntr ; if (associated(reg%Tr(m)%df_x) .or. associated(reg%Tr(m)%df_y)) &
191  call mom_error(warning, "tracer_hordiff: Tracer "//trim(reg%Tr(m)%name)// &
192  " has associated 3-d diffusive flux diagnostics. These are not "//&
193  "valid when DIFFUSE_ML_TO_INTERIOR is defined. Use 2-d tracer "//&
194  "diffusion diagnostics instead to get accurate total fluxes.")
195  enddo
196  endif ; endif
197  cs%first_call = .false.
198 
199  if (cs%debug) call mom_tracer_chksum("Before tracer diffusion ", reg%Tr, ntr, g)
200 
201  use_varmix = .false. ; resoln_scaled = .false. ; use_eady = .false.
202  if (Associated(varmix)) then
203  use_varmix = varmix%use_variable_mixing
204  resoln_scaled = varmix%Resoln_scaled_KhTr
205  use_eady = cs%KhTr_Slope_Cff > 0.
206  endif
207 
208  call cpu_clock_begin(id_clock_pass)
209  do m=1,ntr
210  call create_group_pass(cs%pass_t, reg%Tr(m)%t(:,:,:), g%Domain)
211  enddo
212  call cpu_clock_end(id_clock_pass)
213 
214  if (cs%show_call_tree) call calltree_waypoint("Calculating diffusivity (tracer_hordiff)")
215 
216  if (do_online) then
217  if (use_varmix) then
218  !$OMP parallel do default(shared) private(Kh_loc,Rd_dx)
219  do j=js,je ; do i=is-1,ie
220  kh_loc = cs%KhTr
221  if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2u(i,j)*varmix%SN_u(i,j)
222  if (associated(meke%Kh)) &
223  kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
224  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
225  if (resoln_scaled) &
226  kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
227  kh_u(i,j) = max(kh_loc, cs%KhTr_min)
228  if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
229  rd_dx=0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i+1,j) ) ! Rd/dx at u-points
230  kh_loc = kh_u(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
231  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
232  kh_u(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
233  endif
234  enddo ; enddo
235  !$OMP parallel do default(shared) private(Kh_loc,Rd_dx)
236  do j=js-1,je ; do i=is,ie
237  kh_loc = cs%KhTr
238  if (use_eady) kh_loc = kh_loc + cs%KhTr_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
239  if (associated(meke%Kh)) &
240  kh_loc = kh_loc + meke%KhTr_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
241  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max)
242  if (resoln_scaled) &
243  kh_loc = kh_loc * 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
244  kh_v(i,j) = max(kh_loc, cs%KhTr_min)
245  if (cs%KhTr_passivity_coeff>0.) then ! Apply passivity
246  rd_dx = 0.5*( varmix%Rd_dx_h(i,j)+varmix%Rd_dx_h(i,j+1) ) ! Rd/dx at v-points
247  kh_loc = kh_v(i,j)*max( cs%KhTr_passivity_min, cs%KhTr_passivity_coeff*rd_dx )
248  if (cs%KhTr_max > 0.) kh_loc = min(kh_loc, cs%KhTr_max) ! Re-apply max
249  kh_v(i,j) = max(kh_loc, cs%KhTr_min) ! Re-apply min
250  endif
251  enddo ; enddo
252 
253  !$OMP parallel do default(shared)
254  do j=js,je ; do i=is-1,ie
255  khdt_x(i,j) = dt*(kh_u(i,j)*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
256  enddo ; enddo
257  !$OMP parallel do default(shared)
258  do j=js-1,je ; do i=is,ie
259  khdt_y(i,j) = dt*(kh_v(i,j)*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
260  enddo ; enddo
261  elseif (resoln_scaled) then
262  !$OMP parallel do default(shared) private(Res_fn)
263  do j=js,je ; do i=is-1,ie
264  res_fn = 0.5 * (varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i+1,j))
265  kh_u(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
266  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j))) * res_fn
267  enddo ; enddo
268  !$OMP parallel do default(shared) private(Res_fn)
269  do j=js-1,je ; do i=is,ie
270  res_fn = 0.5*(varmix%Res_fn_h(i,j) + varmix%Res_fn_h(i,j+1))
271  kh_v(i,j) = max(cs%KhTr * res_fn, cs%KhTr_min)
272  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j))) * res_fn
273  enddo ; enddo
274  else ! Use a simple constant diffusivity.
275  if (cs%id_KhTr_u > 0) then
276  !$OMP parallel do default(shared)
277  do j=js,je ; do i=is-1,ie
278  kh_u(i,j) = cs%KhTr
279  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
280  enddo ; enddo
281  else
282  !$OMP parallel do default(shared)
283  do j=js,je ; do i=is-1,ie
284  khdt_x(i,j) = dt*(cs%KhTr*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
285  enddo ; enddo
286  endif
287  if (cs%id_KhTr_v > 0) then
288  !$OMP parallel do default(shared)
289  do j=js-1,je ; do i=is,ie
290  kh_v(i,j) = cs%KhTr
291  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
292  enddo ; enddo
293  else
294  !$OMP parallel do default(shared)
295  do j=js-1,je ; do i=is,ie
296  khdt_y(i,j) = dt*(cs%KhTr*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
297  enddo ; enddo
298  endif
299  endif ! VarMix
300 
301  if (cs%max_diff_CFL > 0.0) then
302  if ((cs%id_KhTr_u > 0) .or. (cs%id_KhTr_h > 0)) then
303  !$OMP parallel do default(shared) private(khdt_max)
304  do j=js,je ; do i=is-1,ie
305  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
306  if (khdt_x(i,j) > khdt_max) then
307  khdt_x(i,j) = khdt_max
308  if (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)) > 0.0) &
309  kh_u(i,j) = khdt_x(i,j) / (dt*(g%dy_Cu(i,j)*g%IdxCu(i,j)))
310  endif
311  enddo ; enddo
312  else
313  !$OMP parallel do default(shared) private(khdt_max)
314  do j=js,je ; do i=is-1,ie
315  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i+1,j))
316  khdt_x(i,j) = min(khdt_x(i,j), khdt_max)
317  enddo ; enddo
318  endif
319  if ((cs%id_KhTr_v > 0) .or. (cs%id_KhTr_h > 0)) then
320  !$OMP parallel do default(shared) private(khdt_max)
321  do j=js-1,je ; do i=is,ie
322  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
323  if (khdt_y(i,j) > khdt_max) then
324  khdt_y(i,j) = khdt_max
325  if (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)) > 0.0) &
326  kh_v(i,j) = khdt_y(i,j) / (dt*(g%dx_Cv(i,j)*g%IdyCv(i,j)))
327  endif
328  enddo ; enddo
329  else
330  !$OMP parallel do default(shared) private(khdt_max)
331  do j=js-1,je ; do i=is,ie
332  khdt_max = 0.125*cs%max_diff_CFL * min(g%areaT(i,j), g%areaT(i,j+1))
333  khdt_y(i,j) = min(khdt_y(i,j), khdt_max)
334  enddo ; enddo
335  endif
336  endif
337 
338  else ! .not. do_online
339  !$OMP parallel do default(shared)
340  do j=js,je ; do i=is-1,ie
341  khdt_x(i,j) = us%m_to_L**2*read_khdt_x(i,j)
342  enddo ; enddo
343  !$OMP parallel do default(shared)
344  do j=js-1,je ; do i=is,ie
345  khdt_y(i,j) = us%m_to_L**2*read_khdt_y(i,j)
346  enddo ; enddo
347  call pass_vector(khdt_x, khdt_y, g%Domain)
348  endif ! do_online
349 
350  if (cs%check_diffusive_CFL) then
351  if (cs%show_call_tree) call calltree_waypoint("Checking diffusive CFL (tracer_hordiff)")
352  max_cfl = 0.0
353  do j=js,je ; do i=is,ie
354  cfl(i,j) = 2.0*((khdt_x(i-1,j) + khdt_x(i,j)) + &
355  (khdt_y(i,j-1) + khdt_y(i,j))) * g%IareaT(i,j)
356  if (max_cfl < cfl(i,j)) max_cfl = cfl(i,j)
357  enddo ; enddo
358  call cpu_clock_begin(id_clock_sync)
359  call max_across_pes(max_cfl)
360  call cpu_clock_end(id_clock_sync)
361  num_itts = max(1, ceiling(max_cfl - 4.0*epsilon(max_cfl)))
362  i_numitts = 1.0 / (real(num_itts))
363  if (cs%id_CFL > 0) call post_data(cs%id_CFL, cfl, cs%diag, mask=g%mask2dT)
364  elseif (cs%max_diff_CFL > 0.0) then
365  num_itts = max(1, ceiling(cs%max_diff_CFL - 4.0*epsilon(cs%max_diff_CFL)))
366  i_numitts = 1.0 / (real(num_itts))
367  else
368  num_itts = 1 ; i_numitts = 1.0
369  endif
370 
371  do m=1,ntr
372  if (associated(reg%Tr(m)%df_x)) then
373  do k=1,nz ; do j=js,je ; do i=is-1,ie
374  reg%Tr(m)%df_x(i,j,k) = 0.0
375  enddo ; enddo ; enddo
376  endif
377  if (associated(reg%Tr(m)%df_y)) then
378  do k=1,nz ; do j=js-1,je ; do i=is,ie
379  reg%Tr(m)%df_y(i,j,k) = 0.0
380  enddo ; enddo ; enddo
381  endif
382  if (associated(reg%Tr(m)%df2d_x)) then
383  do j=js,je ; do i=is-1,ie ; reg%Tr(m)%df2d_x(i,j) = 0.0 ; enddo ; enddo
384  endif
385  if (associated(reg%Tr(m)%df2d_y)) then
386  do j=js-1,je ; do i=is,ie ; reg%Tr(m)%df2d_y(i,j) = 0.0 ; enddo ; enddo
387  endif
388  enddo
389 
390  if (cs%use_lateral_boundary_diffusion) then
391 
392  if (cs%show_call_tree) call calltree_waypoint("Calling lateral boundary mixing (tracer_hordiff)")
393 
394  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
395 
396  do j=js-1,je ; do i=is,ie
397  coef_y(i,j) = i_numitts * khdt_y(i,j)
398  enddo ; enddo
399  do j=js,je
400  do i=is-1,ie
401  coef_x(i,j) = i_numitts * khdt_x(i,j)
402  enddo
403  enddo
404 
405  do itt=1,num_itts
406  if (cs%show_call_tree) call calltree_waypoint("Calling lateral boundary diffusion (tracer_hordiff)",itt)
407  if (itt>1) then ! Update halos for subsequent iterations
408  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
409  endif
410  call lateral_boundary_diffusion(g, gv, us, h, coef_x, coef_y, i_numitts*dt, reg, &
411  cs%lateral_boundary_diffusion_CSp)
412  enddo ! itt
413  endif
414 
415  if (cs%use_neutral_diffusion) then
416 
417  if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion coeffs (tracer_hordiff)")
418 
419  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
420  ! We are assuming that neutral surfaces do not evolve (much) as a result of multiple
421  ! lateral diffusion iterations. Otherwise the call to neutral_diffusion_calc_coeffs()
422  ! would be inside the itt-loop. -AJA
423 
424  if (associated(tv%p_surf)) then
425  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp, p_surf=tv%p_surf)
426  else
427  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
428  endif
429  do j=js-1,je ; do i=is,ie
430  coef_y(i,j) = i_numitts * khdt_y(i,j)
431  enddo ; enddo
432  do j=js,je
433  do i=is-1,ie
434  coef_x(i,j) = i_numitts * khdt_x(i,j)
435  enddo
436  enddo
437 
438  do itt=1,num_itts
439  if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt)
440  if (itt>1) then ! Update halos for subsequent iterations
441  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
442  if (cs%recalc_neutral_surf) then
443  if (associated(tv%p_surf)) then
444  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp, p_surf=tv%p_surf)
445  else
446  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
447  endif
448  endif
449  endif
450  call neutral_diffusion(g, gv, h, coef_x, coef_y, i_numitts*dt, reg, us, cs%neutral_diffusion_CSp)
451  enddo ! itt
452 
453  else ! following if not using neutral diffusion, but instead along-surface diffusion
454 
455  if (cs%show_call_tree) call calltree_waypoint("Calculating horizontal diffusion (tracer_hordiff)")
456  do itt=1,num_itts
457  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
458  !$OMP parallel do default(shared) private(scale,Coef_y,Coef_x,Ihdxdy,dTr)
459  do k=1,nz
460  scale = i_numitts
461  if (cs%Diffuse_ML_interior) then
462  if (k<=gv%nkml) then
463  if (cs%ML_KhTr_scale <= 0.0) cycle
464  scale = i_numitts * cs%ML_KhTr_scale
465  endif
466  if ((k>gv%nkml) .and. (k<=gv%nk_rho_varies)) cycle
467  endif
468 
469  do j=js-1,je ; do i=is,ie
470  coef_y(i,j) = ((scale * khdt_y(i,j))*2.0*(h(i,j,k)*h(i,j+1,k))) / &
471  (h(i,j,k)+h(i,j+1,k)+h_neglect)
472  enddo ; enddo
473 
474  do j=js,je
475  do i=is-1,ie
476  coef_x(i,j) = ((scale * khdt_x(i,j))*2.0*(h(i,j,k)*h(i+1,j,k))) / &
477  (h(i,j,k)+h(i+1,j,k)+h_neglect)
478  enddo
479 
480  do i=is,ie
481  ihdxdy(i,j) = g%IareaT(i,j) / (h(i,j,k)+h_neglect)
482  enddo
483  enddo
484 
485  do m=1,ntr
486  do j=js,je ; do i=is,ie
487  dtr(i,j) = ihdxdy(i,j) * &
488  ((coef_x(i-1,j) * (reg%Tr(m)%t(i-1,j,k) - reg%Tr(m)%t(i,j,k)) - &
489  coef_x(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))) + &
490  (coef_y(i,j-1) * (reg%Tr(m)%t(i,j-1,k) - reg%Tr(m)%t(i,j,k)) - &
491  coef_y(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))))
492  enddo ; enddo
493  if (associated(reg%Tr(m)%df_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
494  reg%Tr(m)%df_x(i,j,k) = reg%Tr(m)%df_x(i,j,k) + coef_x(i,j) &
495  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k)) * idt
496  enddo ; enddo ; endif
497  if (associated(reg%Tr(m)%df_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
498  reg%Tr(m)%df_y(i,j,k) = reg%Tr(m)%df_y(i,j,k) + coef_y(i,j) &
499  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k)) * idt
500  enddo ; enddo ; endif
501  if (associated(reg%Tr(m)%df2d_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
502  reg%Tr(m)%df2d_x(i,j) = reg%Tr(m)%df2d_x(i,j) + coef_x(i,j) &
503  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k)) * idt
504  enddo ; enddo ; endif
505  if (associated(reg%Tr(m)%df2d_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
506  reg%Tr(m)%df2d_y(i,j) = reg%Tr(m)%df2d_y(i,j) + coef_y(i,j) &
507  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k)) * idt
508  enddo ; enddo ; endif
509  do j=js,je ; do i=is,ie
510  reg%Tr(m)%t(i,j,k) = reg%Tr(m)%t(i,j,k) + dtr(i,j)
511  enddo ; enddo
512  enddo
513 
514  enddo ! End of k loop.
515 
516  enddo ! End of "while" loop.
517 
518  endif ! endif for CS%use_neutral_diffusion
519  call cpu_clock_end(id_clock_diffuse)
520 
521 
522  if (cs%Diffuse_ML_interior) then
523  if (cs%show_call_tree) call calltree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)")
524  if (cs%debug) call mom_tracer_chksum("Before epipycnal diff ", reg%Tr, ntr, g)
525 
526  call cpu_clock_begin(id_clock_epimix)
527  call tracer_epipycnal_ml_diff(h, dt, reg%Tr, ntr, khdt_x, khdt_y, g, gv, us, &
528  cs, tv, num_itts)
529  call cpu_clock_end(id_clock_epimix)
530  endif
531 
532  if (cs%debug) call mom_tracer_chksum("After tracer diffusion ", reg%Tr, ntr, g)
533 
534  ! post diagnostics for 2d tracer diffusivity
535  if (cs%id_KhTr_u > 0) then
536  do j=js,je ; do i=is-1,ie
537  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
538  enddo ; enddo
539  call post_data(cs%id_KhTr_u, kh_u, cs%diag, mask=g%mask2dCu)
540  endif
541  if (cs%id_KhTr_v > 0) then
542  do j=js-1,je ; do i=is,ie
543  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
544  enddo ; enddo
545  call post_data(cs%id_KhTr_v, kh_v, cs%diag, mask=g%mask2dCv)
546  endif
547  if (cs%id_KhTr_h > 0) then
548  kh_h(:,:) = 0.0
549  do j=js,je ; do i=is-1,ie
550  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
551  enddo ; enddo
552  do j=js-1,je ; do i=is,ie
553  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
554  enddo ; enddo
555  do j=js,je ; do i=is,ie
556  normalize = 1.0 / ((g%mask2dCu(i-1,j)+g%mask2dCu(i,j)) + &
557  (g%mask2dCv(i,j-1)+g%mask2dCv(i,j)) + gv%H_subroundoff)
558  kh_h(i,j) = normalize*g%mask2dT(i,j)*((kh_u(i-1,j)+kh_u(i,j)) + &
559  (kh_v(i,j-1)+kh_v(i,j)))
560  enddo ; enddo
561  call post_data(cs%id_KhTr_h, kh_h, cs%diag, mask=g%mask2dT)
562  endif
563 
564 
565  if (cs%debug) then
566  call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, &
567  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2, &
568  scalar_pair=.true.)
569  if (cs%use_neutral_diffusion) then
570  call uvchksum("After tracer diffusion Coef_[xy]", coef_x, coef_y, &
571  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2, &
572  scalar_pair=.true.)
573  endif
574  endif
575 
576  if (cs%id_khdt_x > 0) call post_data(cs%id_khdt_x, khdt_x, cs%diag)
577  if (cs%id_khdt_y > 0) call post_data(cs%id_khdt_y, khdt_y, cs%diag)
578 
579  if (cs%show_call_tree) call calltree_leave("tracer_hordiff()")
580 
581 end subroutine tracer_hordiff
582 
583 !> This subroutine does epipycnal diffusion of all tracers between the mixed
584 !! and buffer layers and the interior, using the diffusivity in CS%KhTr.
585 !! Multiple iterations are used (if necessary) so that there is no limit on the
586 !! acceptable time increment.
587 subroutine tracer_epipycnal_ml_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
588  GV, US, CS, tv, num_itts)
589  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
590  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
591  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
592  real, intent(in) :: dt !< time step [T ~> s]
593  type(tracer_type), intent(inout) :: Tr(:) !< tracer array
594  integer, intent(in) :: ntr !< number of tracers
595  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: khdt_epi_x !< Zonal epipycnal diffusivity times
596  !! a time step and the ratio of the open face width over
597  !! the distance between adjacent tracer points [L2 ~> m2]
598  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times
599  !! a time step and the ratio of the open face width over
600  !! the distance between adjacent tracer points [L2 ~> m2]
601  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
602  type(tracer_hor_diff_cs), intent(inout) :: CS !< module control structure
603  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure
604  integer, intent(in) :: num_itts !< number of iterations (usually=1)
605 
606 
607  real, dimension(SZI_(G), SZJ_(G)) :: &
608  Rml_max ! The maximum coordinate density within the mixed layer [R ~> kg m-3].
609  real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: &
610  rho_coord ! The coordinate density that is used to mix along [R ~> kg m-3].
611 
612  ! The naming mnemonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point.
613  ! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage.
614  type(p2d), dimension(SZJ_(G)) :: &
615  deep_wt_Lu, deep_wt_Ru, & ! The relative weighting of the deeper of a pair [nondim].
616  hP_Lu, hP_Ru ! The total thickness on each side for each pair [H ~> m or kg m-2].
617 
618  type(p2d), dimension(SZJB_(G)) :: &
619  deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair [nondim].
620  hP_Lv, hP_Rv ! The total thickness on each side for each pair [H ~> m or kg m-2].
621 
622  type(p2di), dimension(SZJ_(G)) :: &
623  k0b_Lu, k0a_Lu, & ! The original k-indices of the layers that participate
624  k0b_Ru, k0a_Ru ! in each pair of mixing at u-faces.
625  type(p2di), dimension(SZJB_(G)) :: &
626  k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate
627  k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces.
628 
629  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
630  tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
631  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R
632 
633  real, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
634  rho_srt, & ! The density of each layer of the sorted columns [R ~> kg m-3].
635  h_srt ! The thickness of each layer of the sorted columns [H ~> m or kg m-2].
636  integer, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
637  k0_srt ! The original k-index that each layer of the sorted column
638  ! corresponds to.
639 
640  real, dimension(SZK_(G)) :: &
641  h_demand_L, & ! The thickness in the left (_L) or right (_R) column that
642  h_demand_R, & ! is demanded to match the thickness in the counterpart [H ~> m or kg m-2].
643  h_used_L, & ! The summed thickness from the left or right columns that
644  h_used_R, & ! have actually been used [H ~> m or kg m-2].
645  h_supply_frac_L, & ! The fraction of the demanded thickness that can
646  h_supply_frac_R ! actually be supplied from a layer.
647  integer, dimension(SZI_(G), SZJ_(G)) :: &
648  num_srt, & ! The number of layers that are sorted in each column.
649  k_end_srt, & ! The maximum index in each column that might need to be
650  ! sorted, based on neighboring values of max_kRho
651  max_krho ! The index of the layer whose target density is just denser
652  ! than the densest part of the mixed layer.
653  integer, dimension(SZJ_(G)) :: &
654  max_srt ! The maximum value of num_srt in a k-row.
655  integer, dimension(SZIB_(G), SZJ_(G)) :: &
656  nPu ! The number of epipycnal pairings at each u-point.
657  integer, dimension(SZI_(G), SZJB_(G)) :: &
658  nPv ! The number of epipycnal pairings at each v-point.
659  real :: h_exclude ! A thickness that layers must attain to be considered
660  ! for inclusion in mixing [H ~> m or kg m-2].
661  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
662  real :: I_maxitt ! The inverse of the maximum number of iterations.
663  real :: rho_pair, rho_a, rho_b ! Temporary densities [R ~> kg m-3].
664  real :: Tr_min_face ! The minimum and maximum tracer concentrations
665  real :: Tr_max_face ! associated with a pairing [Conc]
666  real :: Tr_La, Tr_Lb ! The 4 tracer concentrations that might be
667  real :: Tr_Ra, Tr_Rb ! associated with a pairing [Conc]
668  real :: Tr_av_L ! The average tracer concentrations on the left and right
669  real :: Tr_av_R ! sides of a pairing [Conc].
670  real :: Tr_flux ! The tracer flux from left to right in a pair [conc H L2 ~> conc m3 or conc kg].
671  real :: Tr_adj_vert ! A downward vertical adjustment to Tr_flux between the
672  ! two cells that make up one side of the pairing [conc H L2 ~> conc m3 or conc kg].
673  real :: h_L, h_R ! Thicknesses to the left and right [H ~> m or kg m-2].
674  real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim].
675  real :: vol ! A cell volume or mass [H L2 ~> m3 or kg].
676 
677  ! The total number of pairings is usually much less than twice the number of layers, but
678  ! the memory in these 1-d columns of pairings can be allocated generously for safety.
679  integer, dimension(SZK_(G)*2) :: &
680  kbs_Lp, & ! The sorted indices of the Left and Right columns for
681  kbs_Rp ! each pairing.
682  logical, dimension(SZK_(G)*2) :: &
683  left_set, & ! If true, the left or right point determines the density of
684  right_set ! of the trio. If densities are exactly equal, both are true.
685 
686  real :: tmp
687  real :: p_ref_cv(SZI_(G)) ! The reference pressure for the coordinate density [R L2 T-2 ~> Pa]
688 
689  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
690  integer :: k_max, k_min, k_test, itmp
691  integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb
692  integer :: isd, ied, jsd, jed, IsdB, IedB, k_size
693  integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt
694  integer :: PEmax_kRho
695  integer :: isv, iev, jsv, jev ! The valid range of the indices.
696 
697  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
698  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
699  isdb = g%IsdB ; iedb = g%IedB
700  idt = 1.0 / dt
701  nkmb = gv%nk_rho_varies
702 
703  if (num_itts <= 1) then
704  max_itt = 1 ; i_maxitt = 1.0
705  else
706  max_itt = num_itts ; i_maxitt = 1.0 / (real(max_itt))
707  endif
708 
709  do i=is-2,ie+2 ; p_ref_cv(i) = tv%P_Ref ; enddo
710  eosdom(:) = eos_domain(g%HI,halo=2)
711 
712  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
713  ! Determine which layers the mixed- and buffer-layers map into...
714  !$OMP parallel do default(shared)
715  do k=1,nkmb ; do j=js-2,je+2
716  call calculate_density(tv%T(:,j,k),tv%S(:,j,k), p_ref_cv, rho_coord(:,j,k), &
717  tv%eqn_of_state, eosdom)
718  enddo ; enddo
719 
720  do j=js-2,je+2 ; do i=is-2,ie+2
721  rml_max(i,j) = rho_coord(i,j,1)
722  num_srt(i,j) = 0 ; max_krho(i,j) = 0
723  enddo ; enddo
724  do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2
725  if (rml_max(i,j) < rho_coord(i,j,k)) rml_max(i,j) = rho_coord(i,j,k)
726  enddo ; enddo ; enddo
727  ! Use bracketing and bisection to find the k-level that the densest of the
728  ! mixed and buffer layer corresponds to, such that:
729  ! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
730  !$OMP parallel do default(shared) private(k_min,k_max,k_test)
731  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (g%mask2dT(i,j) > 0.5) then
732  if ((rml_max(i,j) > gv%Rlay(nz)) .or. (nkmb+1 > nz)) then ; max_krho(i,j) = nz+1
733  elseif ((rml_max(i,j) <= gv%Rlay(nkmb+1)) .or. (nkmb+2 > nz)) then ; max_krho(i,j) = nkmb+1
734  else
735  k_min = nkmb+2 ; k_max = nz
736  do
737  k_test = (k_min + k_max) / 2
738  if (rml_max(i,j) <= gv%Rlay(k_test-1)) then ; k_max = k_test-1
739  elseif (gv%Rlay(k_test) < rml_max(i,j)) then ; k_min = k_test+1
740  else ; max_krho(i,j) = k_test ; exit ; endif
741 
742  if (k_min == k_max) then ; max_krho(i,j) = k_max ; exit ; endif
743  enddo
744  endif
745  endif ; enddo ; enddo
746 
747  pemax_krho = 0
748  do j=js-1,je+1 ; do i=is-1,ie+1
749  k_end_srt(i,j) = max(max_krho(i,j), max_krho(i-1,j), max_krho(i+1,j), &
750  max_krho(i,j-1), max_krho(i,j+1))
751  if (pemax_krho < k_end_srt(i,j)) pemax_krho = k_end_srt(i,j)
752  enddo ; enddo
753  if (pemax_krho > nz) pemax_krho = nz ! PEmax_kRho could have been nz+1.
754 
755  h_exclude = 10.0*(gv%Angstrom_H + gv%H_subroundoff)
756  !$OMP parallel default(shared) private(ns,tmp,itmp)
757  !$OMP do
758  do j=js-1,je+1
759  do k=1,nkmb ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
760  if (h(i,j,k) > h_exclude) then
761  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
762  k0_srt(i,ns,j) = k
763  rho_srt(i,ns,j) = rho_coord(i,j,k)
764  h_srt(i,ns,j) = h(i,j,k)
765  endif
766  endif ; enddo ; enddo
767  do k=nkmb+1,pemax_krho ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
768  if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
769  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
770  k0_srt(i,ns,j) = k
771  rho_srt(i,ns,j) = gv%Rlay(k)
772  h_srt(i,ns,j) = h(i,j,k)
773  endif
774  endif ; enddo ; enddo
775  enddo
776  ! Sort each column by increasing density. This should already be close,
777  ! and the size of the arrays are small, so straight insertion is used.
778  !$OMP do
779  do j=js-1,je+1; do i=is-1,ie+1
780  do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then
781  ! The last segment needs to be shuffled earlier in the list.
782  do k2 = k,2,-1 ; if (rho_srt(i,k2,j) >= rho_srt(i,k2-1,j)) exit
783  itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp
784  tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp
785  tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp
786  enddo
787  endif ; enddo
788  enddo ; enddo
789  !$OMP do
790  do j=js-1,je+1
791  max_srt(j) = 0
792  do i=is-1,ie+1 ; max_srt(j) = max(max_srt(j), num_srt(i,j)) ; enddo
793  enddo
794  !$OMP end parallel
795 
796  do j=js,je
797  k_size = max(2*max_srt(j),1)
798  allocate(deep_wt_lu(j)%p(isdb:iedb,k_size))
799  allocate(deep_wt_ru(j)%p(isdb:iedb,k_size))
800  allocate(hp_lu(j)%p(isdb:iedb,k_size))
801  allocate(hp_ru(j)%p(isdb:iedb,k_size))
802  allocate(k0a_lu(j)%p(isdb:iedb,k_size))
803  allocate(k0a_ru(j)%p(isdb:iedb,k_size))
804  allocate(k0b_lu(j)%p(isdb:iedb,k_size))
805  allocate(k0b_ru(j)%p(isdb:iedb,k_size))
806  enddo
807 
808 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lu,k0_srt, &
809 !$OMP k0b_Ru,k0a_Lu,k0a_Ru,deep_wt_Lu,deep_wt_Ru, &
810 !$OMP h_srt,nkmb,nPu,hP_Lu,hP_Ru) &
811 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
812 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
813 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
814 !$OMP h_supply_frac_L)
815  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
816  ! Set up the pairings for fluxes through the zonal faces.
817 
818  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
819  do k=1,num_srt(i+1,j) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
820 
821  ! First merge the left and right lists into a single, sorted list.
822 
823  ! Discard any layers that are lighter than the lightest in the other
824  ! column. They can only participate in mixing as the lighter part of a
825  ! pair of points.
826  if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then
827  kr = 1
828  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i+1,1,j)) exit ; enddo
829  elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then
830  kl = 1
831  do kr=2,num_srt(i+1,j) ; if (rho_srt(i+1,kr,j) >= rho_srt(i,1,j)) exit ; enddo
832  else
833  kl = 1 ; kr = 1
834  endif
835  np = 0
836  do ! Loop to accumulate pairs of columns.
837  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i+1,j))) exit
838 
839  if (rho_srt(i,kl,j) > rho_srt(i+1,kr,j)) then
840  ! The right point is lighter and defines the density for this trio.
841  np = np+1 ; k = np
842  rho_pair = rho_srt(i+1,kr,j)
843 
844  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
845  k0a_lu(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
846  kbs_lp(k) = kl ; kbs_rp(k) = kr
847 
848  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
849  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
850  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
851  deep_wt_lu(j)%p(i,k) = wt_b ; deep_wt_ru(j)%p(i,k) = 1.0
852 
853  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j) * wt_b
854  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i+1,kr,j) * (1.0-wt_b)
855 
856  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
857  elseif (rho_srt(i,kl,j) < rho_srt(i+1,kr,j)) then
858  ! The left point is lighter and defines the density for this trio.
859  np = np+1 ; k = np
860  rho_pair = rho_srt(i,kl,j)
861  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
862  k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0_srt(i+1,kr-1,j)
863 
864  kbs_lp(k) = kl ; kbs_rp(k) = kr
865 
866  rho_a = rho_srt(i+1,kr-1,j) ; rho_b = rho_srt(i+1,kr,j)
867  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
868  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
869  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = wt_b
870 
871  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
872  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
873 
874  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
875  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i+1,kr,j) <= nkmb)) then
876  ! The densities are exactly equal and one layer is above the interior.
877  np = np+1 ; k = np
878  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
879  k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
880  kbs_lp(k) = kl ; kbs_rp(k) = kr
881  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = 1.0
882 
883  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
884  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
885 
886  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
887  else ! The densities are exactly equal and in the interior.
888  ! Mixing in this case has already occurred, so accumulate the thickness
889  ! demanded for that mixing and skip onward.
890  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
891  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
892 
893  kl = kl+1 ; kr = kr+1
894  endif
895  enddo ! Loop to accumulate pairs of columns.
896  npu(i,j) = np ! This is the number of active pairings.
897 
898  ! Determine what fraction of the thickness "demand" can be supplied.
899  do k=1,num_srt(i+1,j)
900  h_supply_frac_r(k) = 1.0
901  if (h_demand_r(k) > 0.5*h_srt(i+1,k,j)) &
902  h_supply_frac_r(k) = 0.5*h_srt(i+1,k,j) / h_demand_r(k)
903  enddo
904  do k=1,num_srt(i,j)
905  h_supply_frac_l(k) = 1.0
906  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
907  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
908  enddo
909 
910  ! Distribute the "exported" thicknesses proportionately.
911  do k=1,npu(i,j)
912  kl = kbs_lp(k) ; kr = kbs_rp(k)
913  hp_lu(j)%p(i,k) = 0.0 ; hp_ru(j)%p(i,k) = 0.0
914  if (left_set(k)) then ! Add the contributing thicknesses on the right.
915  if (deep_wt_ru(j)%p(i,k) < 1.0) then
916  hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * min(h_supply_frac_r(kr), h_supply_frac_r(kr-1))
917  wt_b = deep_wt_ru(j)%p(i,k)
918  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b)*hp_ru(j)%p(i,k)
919  h_used_r(kr) = h_used_r(kr) + wt_b*hp_ru(j)%p(i,k)
920  else
921  hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
922  h_used_r(kr) = h_used_r(kr) + hp_ru(j)%p(i,k)
923  endif
924  endif
925  if (right_set(k)) then ! Add the contributing thicknesses on the left.
926  if (deep_wt_lu(j)%p(i,k) < 1.0) then
927  hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * min(h_supply_frac_l(kl), h_supply_frac_l(kl-1))
928  wt_b = deep_wt_lu(j)%p(i,k)
929  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b)*hp_lu(j)%p(i,k)
930  h_used_l(kl) = h_used_l(kl) + wt_b*hp_lu(j)%p(i,k)
931  else
932  hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * h_supply_frac_l(kl)
933  h_used_l(kl) = h_used_l(kl) + hp_lu(j)%p(i,k)
934  endif
935  endif
936  enddo
937 
938  ! The left-over thickness (at least half the layer thickness) is now
939  ! added to the thicknesses of the importing columns.
940  do k=1,npu(i,j)
941  if (left_set(k)) hp_lu(j)%p(i,k) = hp_lu(j)%p(i,k) + &
942  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
943  if (right_set(k)) hp_ru(j)%p(i,k) = hp_ru(j)%p(i,k) + &
944  (h_srt(i+1,kbs_rp(k),j) - h_used_r(kbs_rp(k)))
945  enddo
946 
947  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
948 
949  do j=js-1,je
950  k_size = max(max_srt(j)+max_srt(j+1),1)
951  allocate(deep_wt_lv(j)%p(isd:ied,k_size))
952  allocate(deep_wt_rv(j)%p(isd:ied,k_size))
953  allocate(hp_lv(j)%p(isd:ied,k_size))
954  allocate(hp_rv(j)%p(isd:ied,k_size))
955  allocate(k0a_lv(j)%p(isd:ied,k_size))
956  allocate(k0a_rv(j)%p(isd:ied,k_size))
957  allocate(k0b_lv(j)%p(isd:ied,k_size))
958  allocate(k0b_rv(j)%p(isd:ied,k_size))
959  enddo
960 
961 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lv,k0b_Rv, &
962 !$OMP k0_srt,k0a_Lv,k0a_Rv,deep_wt_Lv,deep_wt_Rv, &
963 !$OMP h_srt,nkmb,nPv,hP_Lv,hP_Rv) &
964 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
965 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
966 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
967 !$OMP h_supply_frac_L)
968  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
969  ! Set up the pairings for fluxes through the meridional faces.
970 
971  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
972  do k=1,num_srt(i,j+1) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
973 
974  ! First merge the left and right lists into a single, sorted list.
975 
976  ! Discard any layers that are lighter than the lightest in the other
977  ! column. They can only participate in mixing as the lighter part of a
978  ! pair of points.
979  if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then
980  kr = 1
981  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i,1,j+1)) exit ; enddo
982  elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then
983  kl = 1
984  do kr=2,num_srt(i,j+1) ; if (rho_srt(i,kr,j+1) >= rho_srt(i,1,j)) exit ; enddo
985  else
986  kl = 1 ; kr = 1
987  endif
988  np = 0
989  do ! Loop to accumulate pairs of columns.
990  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i,j+1))) exit
991 
992  if (rho_srt(i,kl,j) > rho_srt(i,kr,j+1)) then
993  ! The right point is lighter and defines the density for this trio.
994  np = np+1 ; k = np
995  rho_pair = rho_srt(i,kr,j+1)
996 
997  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
998  k0a_lv(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
999  kbs_lp(k) = kl ; kbs_rp(k) = kr
1000 
1001  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
1002  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1003  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1004  deep_wt_lv(j)%p(i,k) = wt_b ; deep_wt_rv(j)%p(i,k) = 1.0
1005 
1006  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1) * wt_b
1007  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i,kr,j+1) * (1.0-wt_b)
1008 
1009  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
1010  elseif (rho_srt(i,kl,j) < rho_srt(i,kr,j+1)) then
1011  ! The left point is lighter and defines the density for this trio.
1012  np = np+1 ; k = np
1013  rho_pair = rho_srt(i,kl,j)
1014  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1015  k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k) ; k0a_rv(j)%p(i,k) = k0_srt(i,kr-1,j+1)
1016 
1017  kbs_lp(k) = kl ; kbs_rp(k) = kr
1018 
1019  rho_a = rho_srt(i,kr-1,j+1) ; rho_b = rho_srt(i,kr,j+1)
1020  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1021  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1022  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = wt_b
1023 
1024  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
1025  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
1026 
1027  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
1028  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i,kr,j+1) <= nkmb)) then
1029  ! The densities are exactly equal and one layer is above the interior.
1030  np = np+1 ; k = np
1031  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1032  k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
1033  kbs_lp(k) = kl ; kbs_rp(k) = kr
1034  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = 1.0
1035 
1036  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1037  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1038 
1039  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
1040  else ! The densities are exactly equal and in the interior.
1041  ! Mixing in this case has already occurred, so accumulate the thickness
1042  ! demanded for that mixing and skip onward.
1043  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1044  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1045 
1046  kl = kl+1 ; kr = kr+1
1047  endif
1048  enddo ! Loop to accumulate pairs of columns.
1049  npv(i,j) = np ! This is the number of active pairings.
1050 
1051  ! Determine what fraction of the thickness "demand" can be supplied.
1052  do k=1,num_srt(i,j+1)
1053  h_supply_frac_r(k) = 1.0
1054  if (h_demand_r(k) > 0.5*h_srt(i,k,j+1)) &
1055  h_supply_frac_r(k) = 0.5*h_srt(i,k,j+1) / h_demand_r(k)
1056  enddo
1057  do k=1,num_srt(i,j)
1058  h_supply_frac_l(k) = 1.0
1059  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
1060  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
1061  enddo
1062 
1063  ! Distribute the "exported" thicknesses proportionately.
1064  do k=1,npv(i,j)
1065  kl = kbs_lp(k) ; kr = kbs_rp(k)
1066  hp_lv(j)%p(i,k) = 0.0 ; hp_rv(j)%p(i,k) = 0.0
1067  if (left_set(k)) then ! Add the contributing thicknesses on the right.
1068  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1069  hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * min(h_supply_frac_r(kr), h_supply_frac_r(kr-1))
1070  wt_b = deep_wt_rv(j)%p(i,k)
1071  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b) * hp_rv(j)%p(i,k)
1072  h_used_r(kr) = h_used_r(kr) + wt_b * hp_rv(j)%p(i,k)
1073  else
1074  hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
1075  h_used_r(kr) = h_used_r(kr) + hp_rv(j)%p(i,k)
1076  endif
1077  endif
1078  if (right_set(k)) then ! Add the contributing thicknesses on the left.
1079  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1080  hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * min(h_supply_frac_l(kl), h_supply_frac_l(kl-1))
1081  wt_b = deep_wt_lv(j)%p(i,k)
1082  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b) * hp_lv(j)%p(i,k)
1083  h_used_l(kl) = h_used_l(kl) + wt_b * hp_lv(j)%p(i,k)
1084  else
1085  hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * h_supply_frac_l(kl)
1086  h_used_l(kl) = h_used_l(kl) + hp_lv(j)%p(i,k)
1087  endif
1088  endif
1089  enddo
1090 
1091  ! The left-over thickness (at least half the layer thickness) is now
1092  ! added to the thicknesses of the importing columns.
1093  do k=1,npv(i,j)
1094  if (left_set(k)) hp_lv(j)%p(i,k) = hp_lv(j)%p(i,k) + &
1095  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
1096  if (right_set(k)) hp_rv(j)%p(i,k) = hp_rv(j)%p(i,k) + &
1097  (h_srt(i,kbs_rp(k),j+1) - h_used_r(kbs_rp(k)))
1098  enddo
1099 
1100 
1101  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1102 
1103 ! The tracer-specific calculations start here.
1104 
1105  ! Zero out tracer tendencies.
1106  do k=1,pemax_krho ; do j=js-1,je+1 ; do i=is-1,ie+1
1107  tr_flux_conv(i,j,k) = 0.0
1108  enddo ; enddo ; enddo
1109 
1110  do itt=1,max_itt
1111 
1112  if (itt > 1) then ! The halos have already been filled if itt==1.
1113  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
1114  endif
1115 
1116  do m=1,ntr
1117 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, &
1118 !$OMP k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru, &
1119 !$OMP hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) &
1120 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
1121 !$OMP Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
1122 !$OMP Tr_flux,Tr_adj_vert,wt_a,vol)
1123  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
1124  ! Determine the fluxes through the zonal faces.
1125 
1126  ! Find the acceptable range of tracer concentration around this face.
1127  if (npu(i,j) >= 1) then
1128  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1129  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1130  do k=2,nkmb
1131  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1132  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1133  enddo
1134 
1135  ! Include the next two layers denser than the densest buffer layer.
1136  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1137  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1138  kra = nkmb+1 ; if (max_krho(i+1,j) < nz+1) kra = max_krho(i+1,j)
1139  krb = kra ; if (max_krho(i+1,j) < nz) krb = max_krho(i+1,j)+1
1140  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1141  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1142  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1143  if (h(i+1,j,kra) > h_exclude) tr_ra = tr(m)%t(i+1,j,kra)
1144  if (h(i+1,j,krb) > h_exclude) tr_rb = tr(m)%t(i+1,j,krb)
1145  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1146  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1147 
1148  ! Include all points in diffusive pairings at this face.
1149  do k=1,npu(i,j)
1150  tr_lb = tr(m)%t(i,j,k0b_lu(j)%p(i,k))
1151  tr_rb = tr(m)%t(i+1,j,k0b_ru(j)%p(i,k))
1152  tr_la = tr_lb ; tr_ra = tr_rb
1153  if (deep_wt_lu(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lu(j)%p(i,k))
1154  if (deep_wt_ru(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i+1,j,k0a_ru(j)%p(i,k))
1155  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1156  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1157  enddo
1158  endif
1159 
1160  do k=1,npu(i,j)
1161  klb = k0b_lu(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1162  if (deep_wt_lu(j)%p(i,k) < 1.0) then
1163  kla = k0a_lu(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1164  wt_b = deep_wt_lu(j)%p(i,k)
1165  tr_av_l = wt_b*tr_lb + (1.0-wt_b)*tr_la
1166  endif
1167 
1168  krb = k0b_ru(j)%p(i,k) ; tr_rb = tr(m)%t(i+1,j,krb) ; tr_av_r = tr_rb
1169  if (deep_wt_ru(j)%p(i,k) < 1.0) then
1170  kra = k0a_ru(j)%p(i,k) ; tr_ra = tr(m)%t(i+1,j,kra)
1171  wt_b = deep_wt_ru(j)%p(i,k)
1172  tr_av_r = wt_b*tr_rb + (1.0-wt_b)*tr_ra
1173  endif
1174 
1175  h_l = hp_lu(j)%p(i,k) ; h_r = hp_ru(j)%p(i,k)
1176  tr_flux = i_maxitt * khdt_epi_x(i,j) * (tr_av_l - tr_av_r) * &
1177  ((2.0 * h_l * h_r) / (h_l + h_r))
1178 
1179 
1180  if (deep_wt_lu(j)%p(i,k) >= 1.0) then
1181  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux
1182  else
1183  tr_adj_vert = 0.0
1184  wt_b = deep_wt_lu(j)%p(i,k) ; wt_a = 1.0 - wt_b
1185  vol = hp_lu(j)%p(i,k) * g%areaT(i,j)
1186 
1187  ! Ensure that the tracer flux does not drive the tracer values
1188  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1189  ! does that the concentration in both contributing pieces exceed
1190  ! this range equally. With down-gradient fluxes and the initial tracer
1191  ! concentrations determining the valid range, the latter condition
1192  ! only enters for large values of the effective diffusive CFL number.
1193  if (tr_flux > 0.0) then
1194  if (tr_la < tr_lb) then ; if (vol*(tr_la-tr_min_face) < tr_flux) &
1195  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1196  (vol*wt_b) * (tr_lb - tr_la))
1197  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1198  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1199  (vol*wt_a) * (tr_la - tr_lb))
1200  endif
1201  elseif (tr_flux < 0.0) then
1202  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1203  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1204  (vol*wt_b) * (tr_la - tr_lb))
1205  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1206  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1207  (vol*wt_a)*(tr_lb - tr_la))
1208  endif
1209  endif
1210 
1211  tr_flux_conv(i,j,kla) = tr_flux_conv(i,j,kla) - (wt_a*tr_flux + tr_adj_vert)
1212  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - (wt_b*tr_flux - tr_adj_vert)
1213  endif
1214 
1215  if (deep_wt_ru(j)%p(i,k) >= 1.0) then
1216  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + tr_flux
1217  else
1218  tr_adj_vert = 0.0
1219  wt_b = deep_wt_ru(j)%p(i,k) ; wt_a = 1.0 - wt_b
1220  vol = hp_ru(j)%p(i,k) * g%areaT(i+1,j)
1221 
1222  ! Ensure that the tracer flux does not drive the tracer values
1223  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1224  ! does that the concentration in both contributing pieces exceed
1225  ! this range equally. With down-gradient fluxes and the initial tracer
1226  ! concentrations determining the valid range, the latter condition
1227  ! only enters for large values of the effective diffusive CFL number.
1228  if (tr_flux < 0.0) then
1229  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1230  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1231  (vol*wt_b) * (tr_rb - tr_ra))
1232  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1233  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1234  (vol*wt_a) * (tr_ra - tr_rb))
1235  endif
1236  elseif (tr_flux > 0.0) then
1237  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1238  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1239  (vol*wt_b) * (tr_ra - tr_rb))
1240  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1241  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1242  (vol*wt_a)*(tr_rb - tr_ra))
1243  endif
1244  endif
1245 
1246  tr_flux_conv(i+1,j,kra) = tr_flux_conv(i+1,j,kra) + &
1247  (wt_a*tr_flux - tr_adj_vert)
1248  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + &
1249  (wt_b*tr_flux + tr_adj_vert)
1250  endif
1251  if (associated(tr(m)%df2d_x)) &
1252  tr(m)%df2d_x(i,j) = tr(m)%df2d_x(i,j) + tr_flux * idt
1253  enddo ! Loop over pairings at faces.
1254  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
1255 
1256 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, &
1257 !$OMP k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv, &
1258 !$OMP hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d, &
1259 !$OMP Tr_adj_vert_L,Tr_adj_vert_R,Idt) &
1260 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, &
1261 !$OMP Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, &
1262 !$OMP h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
1263  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
1264  ! Determine the fluxes through the meridional faces.
1265 
1266  ! Find the acceptable range of tracer concentration around this face.
1267  if (npv(i,j) >= 1) then
1268  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1269  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1270  do k=2,nkmb
1271  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1272  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1273  enddo
1274 
1275  ! Include the next two layers denser than the densest buffer layer.
1276  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1277  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1278  kra = nkmb+1 ; if (max_krho(i,j+1) < nz+1) kra = max_krho(i,j+1)
1279  krb = kra ; if (max_krho(i,j+1) < nz) krb = max_krho(i,j+1)+1
1280  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1281  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1282  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1283  if (h(i,j+1,kra) > h_exclude) tr_ra = tr(m)%t(i,j+1,kra)
1284  if (h(i,j+1,krb) > h_exclude) tr_rb = tr(m)%t(i,j+1,krb)
1285  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1286  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1287 
1288  ! Include all points in diffusive pairings at this face.
1289  do k=1,npv(i,j)
1290  tr_lb = tr(m)%t(i,j,k0b_lv(j)%p(i,k)) ; tr_rb = tr(m)%t(i,j+1,k0b_rv(j)%p(i,k))
1291  tr_la = tr_lb ; tr_ra = tr_rb
1292  if (deep_wt_lv(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lv(j)%p(i,k))
1293  if (deep_wt_rv(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i,j+1,k0a_rv(j)%p(i,k))
1294  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1295  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1296  enddo
1297  endif
1298 
1299  do k=1,npv(i,j)
1300  klb = k0b_lv(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1301  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1302  kla = k0a_lv(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1303  wt_b = deep_wt_lv(j)%p(i,k)
1304  tr_av_l = wt_b * tr_lb + (1.0-wt_b) * tr_la
1305  endif
1306 
1307  krb = k0b_rv(j)%p(i,k) ; tr_rb = tr(m)%t(i,j+1,krb) ; tr_av_r = tr_rb
1308  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1309  kra = k0a_rv(j)%p(i,k) ; tr_ra = tr(m)%t(i,j+1,kra)
1310  wt_b = deep_wt_rv(j)%p(i,k)
1311  tr_av_r = wt_b * tr_rb + (1.0-wt_b) * tr_ra
1312  endif
1313 
1314  h_l = hp_lv(j)%p(i,k) ; h_r = hp_rv(j)%p(i,k)
1315  tr_flux = i_maxitt * ((2.0 * h_l * h_r) / (h_l + h_r)) * &
1316  khdt_epi_y(i,j) * (tr_av_l - tr_av_r)
1317  tr_flux_3d(i,j,k) = tr_flux
1318 
1319  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1320  tr_adj_vert = 0.0
1321  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1322  vol = hp_lv(j)%p(i,k) * g%areaT(i,j)
1323 
1324  ! Ensure that the tracer flux does not drive the tracer values
1325  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1326  if (tr_flux > 0.0) then
1327  if (tr_la < tr_lb) then ; if (vol * (tr_la-tr_min_face) < tr_flux) &
1328  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1329  (vol*wt_b) * (tr_lb - tr_la))
1330  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1331  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1332  (vol*wt_a) * (tr_la - tr_lb))
1333  endif
1334  elseif (tr_flux < 0.0) then
1335  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1336  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1337  (vol*wt_b) * (tr_la - tr_lb))
1338  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1339  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1340  (vol*wt_a)*(tr_lb - tr_la))
1341  endif
1342  endif
1343  tr_adj_vert_l(i,j,k) = tr_adj_vert
1344  endif
1345 
1346  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1347  tr_adj_vert = 0.0
1348  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1349  vol = hp_rv(j)%p(i,k) * g%areaT(i,j+1)
1350 
1351  ! Ensure that the tracer flux does not drive the tracer values
1352  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1353  if (tr_flux < 0.0) then
1354  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1355  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1356  (vol*wt_b) * (tr_rb - tr_ra))
1357  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1358  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1359  (vol*wt_a) * (tr_ra - tr_rb))
1360  endif
1361  elseif (tr_flux > 0.0) then
1362  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1363  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1364  (vol*wt_b) * (tr_ra - tr_rb))
1365  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1366  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1367  (vol*wt_a)*(tr_rb - tr_ra))
1368  endif
1369  endif
1370  tr_adj_vert_r(i,j,k) = tr_adj_vert
1371  endif
1372  if (associated(tr(m)%df2d_y)) &
1373  tr(m)%df2d_y(i,j) = tr(m)%df2d_y(i,j) + tr_flux * idt
1374  enddo ! Loop over pairings at faces.
1375  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1376 !$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv, &
1377 !$OMP tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
1378 !$OMP deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
1379 !$OMP private(kLa,kLb,kRa,kRb,wt_b,wt_a)
1380  do i=is,ie ; do j=js-1,je ; if (g%mask2dCv(i,j) > 0.5) then
1381  do k=1,npv(i,j)
1382  klb = k0b_lv(j)%p(i,k); krb = k0b_rv(j)%p(i,k)
1383  if (deep_wt_lv(j)%p(i,k) >= 1.0) then
1384  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux_3d(i,j,k)
1385  else
1386  kla = k0a_lv(j)%p(i,k)
1387  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1388  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))
1389  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))
1390  endif
1391  if (deep_wt_rv(j)%p(i,k) >= 1.0) then
1392  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + tr_flux_3d(i,j,k)
1393  else
1394  kra = k0a_rv(j)%p(i,k)
1395  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1396  tr_flux_conv(i,j+1,kra) = tr_flux_conv(i,j+1,kra) + &
1397  (wt_a*tr_flux_3d(i,j,k) - tr_adj_vert_r(i,j,k))
1398  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + &
1399  (wt_b*tr_flux_3d(i,j,k) + tr_adj_vert_r(i,j,k))
1400  endif
1401  enddo
1402  endif ; enddo ; enddo
1403 !$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m)
1404  do k=1,pemax_krho ; do j=js,je ; do i=is,ie
1405  if ((g%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 0.0)) then
1406  tr(m)%t(i,j,k) = tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
1407  (h(i,j,k)*g%areaT(i,j))
1408  tr_flux_conv(i,j,k) = 0.0
1409  endif
1410  enddo ; enddo ; enddo
1411 
1412  enddo ! Loop over tracers
1413  enddo ! Loop over iterations
1414 
1415  do j=js,je
1416  deallocate(deep_wt_lu(j)%p) ; deallocate(deep_wt_ru(j)%p)
1417  deallocate(hp_lu(j)%p) ; deallocate(hp_ru(j)%p)
1418  deallocate(k0a_lu(j)%p) ; deallocate(k0a_ru(j)%p)
1419  deallocate(k0b_lu(j)%p) ; deallocate(k0b_ru(j)%p)
1420  enddo
1421 
1422  do j=js-1,je
1423  deallocate(deep_wt_lv(j)%p) ; deallocate(deep_wt_rv(j)%p)
1424  deallocate(hp_lv(j)%p) ; deallocate(hp_rv(j)%p)
1425  deallocate(k0a_lv(j)%p) ; deallocate(k0a_rv(j)%p)
1426  deallocate(k0b_lv(j)%p) ; deallocate(k0b_rv(j)%p)
1427  enddo
1428 
1429 end subroutine tracer_epipycnal_ml_diff
1430 
1431 
1432 !> Initialize lateral tracer diffusion module
1433 subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS)
1434  type(time_type), target, intent(in) :: time !< current model time
1435  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1436  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1437  type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control
1438  type(eos_type), target, intent(in) :: eos !< Equation of state CS
1439  type(diabatic_cs), pointer, intent(in) :: diabatic_csp !< Equation of state CS
1440  type(param_file_type), intent(in) :: param_file !< parameter file
1441  type(tracer_hor_diff_cs), pointer :: cs !< horz diffusion control structure
1442 
1443 ! This include declares and sets the variable "version".
1444 #include "version_variable.h"
1445  character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name.
1446  character(len=256) :: mesg ! Message for error messages.
1447 
1448  if (associated(cs)) then
1449  call mom_error(warning, "tracer_hor_diff_init called with associated control structure.")
1450  return
1451  endif
1452  allocate(cs)
1453 
1454  cs%diag => diag
1455  cs%show_call_tree = calltree_showquery()
1456 
1457  ! Read all relevant parameters and write them to the model log.
1458  call log_version(param_file, mdl, version, "")
1459  call get_param(param_file, mdl, "KHTR", cs%KhTr, &
1460  "The background along-isopycnal tracer diffusivity.", &
1461  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1462  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", cs%KhTr_Slope_Cff, &
1463  "The scaling coefficient for along-isopycnal tracer "//&
1464  "diffusivity using a shear-based (Visbeck-like) "//&
1465  "parameterization. A non-zero value enables this param.", &
1466  units="nondim", default=0.0)
1467  call get_param(param_file, mdl, "KHTR_MIN", cs%KhTr_Min, &
1468  "The minimum along-isopycnal tracer diffusivity.", &
1469  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1470  call get_param(param_file, mdl, "KHTR_MAX", cs%KhTr_Max, &
1471  "The maximum along-isopycnal tracer diffusivity.", &
1472  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1473  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", cs%KhTr_passivity_coeff, &
1474  "The coefficient that scales deformation radius over "//&
1475  "grid-spacing in passivity, where passivity is the ratio "//&
1476  "between along isopycnal mixing of tracers to thickness mixing. "//&
1477  "A non-zero value enables this parameterization.", &
1478  units="nondim", default=0.0)
1479  call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", cs%KhTr_passivity_min, &
1480  "The minimum passivity which is the ratio between "//&
1481  "along isopycnal mixing of tracers to thickness mixing.", &
1482  units="nondim", default=0.5)
1483  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
1484  "If true, enable epipycnal mixing between the surface "//&
1485  "boundary layer and the interior.", default=.false.)
1486  call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", cs%check_diffusive_CFL, &
1487  "If true, use enough iterations the diffusion to ensure "//&
1488  "that the diffusive equivalent of the CFL limit is not "//&
1489  "violated. If false, always use the greater of 1 or "//&
1490  "MAX_TR_DIFFUSION_CFL iteration.", default=.false.)
1491  call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", cs%max_diff_CFL, &
1492  "If positive, locally limit the along-isopycnal tracer "//&
1493  "diffusivity to keep the diffusive CFL locally at or "//&
1494  "below this value. The number of diffusive iterations "//&
1495  "is often this value or the next greater integer.", &
1496  units="nondim", default=-1.0)
1497  call get_param(param_file, mdl, "RECALC_NEUTRAL_SURF", cs%recalc_neutral_surf, &
1498  "If true, then recalculate the neutral surfaces if the \n"//&
1499  "diffusive CFL is exceeded. If false, assume that the \n"//&
1500  "positions of the surfaces do not change \n", default = .false.)
1501  cs%ML_KhTR_scale = 1.0
1502  if (cs%Diffuse_ML_interior) then
1503  call get_param(param_file, mdl, "ML_KHTR_SCALE", cs%ML_KhTR_scale, &
1504  "With Diffuse_ML_interior, the ratio of the truly "//&
1505  "horizontal diffusivity in the mixed layer to the "//&
1506  "epipycnal diffusivity. The valid range is 0 to 1.", &
1507  units="nondim", default=1.0)
1508  endif
1509 
1510  cs%use_neutral_diffusion = neutral_diffusion_init(time, g, us, param_file, diag, eos, &
1511  diabatic_csp, cs%neutral_diffusion_CSp )
1512  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1513  "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1514  cs%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(time, g, param_file, diag, diabatic_csp, &
1515  cs%lateral_boundary_diffusion_CSp)
1516  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1517  "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1518 
1519  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1520 
1521  id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=clock_module)
1522  id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=clock_module)
1523  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1524  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1525 
1526  cs%id_KhTr_u = -1
1527  cs%id_KhTr_v = -1
1528  cs%id_KhTr_h = -1
1529  cs%id_CFL = -1
1530 
1531  cs%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCu1, time, &
1532  'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1533  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
1534  'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1535  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time, &
1536  'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
1537  cmor_field_name='diftrelo', &
1538  cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', &
1539  cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')
1540 
1541  cs%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, time, &
1542  'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1543  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
1544  'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1545  if (cs%check_diffusive_CFL) then
1546  cs%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, time,&
1547  'Grid CFL number for lateral/neutral tracer diffusion', 'nondim')
1548  endif
1549 
1550 
1551 end subroutine tracer_hor_diff_init
1552 
1553 subroutine tracer_hor_diff_end(CS)
1554  type(tracer_hor_diff_cs), pointer :: cs !< module control structure
1555 
1556  call neutral_diffusion_end(cs%neutral_diffusion_CSp)
1557  if (associated(cs)) deallocate(cs)
1558 
1559 end subroutine tracer_hor_diff_end
1560 
1561 
1562 !> \namespace mom_tracer_hor_diff
1563 !!
1564 !! \section section_intro Introduction to the module
1565 !!
1566 !! This module contains subroutines that handle horizontal
1567 !! diffusion (i.e., isoneutral or along layer) of tracers.
1568 !!
1569 !! Each of the tracers are subject to Fickian along-coordinate
1570 !! diffusion if Khtr is defined and positive. The tracer diffusion
1571 !! can use a suitable number of iterations to guarantee stability
1572 !! with an arbitrarily large time step.
1573 
1574 end module mom_tracer_hor_diff
mom_tracer_hor_diff::p2di
A type that can be used to create arrays of pointers to 2D integer arrays.
Definition: MOM_tracer_hor_diff.F90:92
mom_lateral_boundary_diffusion::lateral_boundary_diffusion_cs
Sets parameters for lateral boundary mixing module.
Definition: MOM_lateral_boundary_diffusion.F90:39
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_tracer_hor_diff
Main routine for lateral (along surface or neutral) diffusion of tracers.
Definition: MOM_tracer_hor_diff.F90:2
mom_tracer_registry::tracer_type
The tracer type.
Definition: MOM_tracer_registry.F90:38
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_tracer_hor_diff::tracer_hor_diff_cs
The control structure for along-layer and epineutral tracer diffusion.
Definition: MOM_tracer_hor_diff.F90:40
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:59
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:108
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_meke_types::meke_type
This type is used to exchange information related to the MEKE calculations.
Definition: MOM_MEKE_types.F90:8
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_lateral_mixing_coeffs::varmix_cs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:27
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_tracer_hor_diff::p2d
A type that can be used to create arrays of pointers to 2D arrays.
Definition: MOM_tracer_hor_diff.F90:88
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:138
mom_lateral_boundary_diffusion
Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by meso...
Definition: MOM_lateral_boundary_diffusion.F90:4
mom_diabatic_driver
This routine drives the diabatic/dianeutral physics for MOM.
Definition: MOM_diabatic_driver.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_lateral_mixing_coeffs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:2
mom_domains::create_group_pass
Set up a group of halo updates.
Definition: MOM_domains.F90:84
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_neutral_diffusion
A column-wise toolbox for implementing neutral diffusion.
Definition: MOM_neutral_diffusion.F90:2
mom_diabatic_driver::diabatic_cs
Control structure for this module.
Definition: MOM_diabatic_driver.F90:92
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240
mom_neutral_diffusion::neutral_diffusion_cs
The control structure for the MOM_neutral_diffusion module.
Definition: MOM_neutral_diffusion.F90:43
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:68