MOM6
MOM_tracer_hor_diff.F90
Go to the documentation of this file.
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
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 ocntrol 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 velcities.
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
145  ! to time-integrated fluxes [H L2 ~> m3 or kg].
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
151  ! to time-integrated fluxes [H L2 ~> m3 or kg].
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  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
425  do j=js-1,je ; do i=is,ie
426  coef_y(i,j) = i_numitts * khdt_y(i,j)
427  enddo ; enddo
428  do j=js,je
429  do i=is-1,ie
430  coef_x(i,j) = i_numitts * khdt_x(i,j)
431  enddo
432  enddo
433 
434  do itt=1,num_itts
435  if (cs%show_call_tree) call calltree_waypoint("Calling neutral diffusion (tracer_hordiff)",itt)
436  if (itt>1) then ! Update halos for subsequent iterations
437  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
438  if (cs%recalc_neutral_surf) then
439  call neutral_diffusion_calc_coeffs(g, gv, us, h, tv%T, tv%S, cs%neutral_diffusion_CSp)
440  endif
441  endif
442  call neutral_diffusion(g, gv, h, coef_x, coef_y, i_numitts*dt, reg, us, cs%neutral_diffusion_CSp)
443  enddo ! itt
444 
445  else ! following if not using neutral diffusion, but instead along-surface diffusion
446 
447  if (cs%show_call_tree) call calltree_waypoint("Calculating horizontal diffusion (tracer_hordiff)")
448  do itt=1,num_itts
449  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
450  !$OMP parallel do default(shared) private(scale,Coef_y,Coef_x,Ihdxdy,dTr)
451  do k=1,nz
452  scale = i_numitts
453  if (cs%Diffuse_ML_interior) then
454  if (k<=gv%nkml) then
455  if (cs%ML_KhTr_scale <= 0.0) cycle
456  scale = i_numitts * cs%ML_KhTr_scale
457  endif
458  if ((k>gv%nkml) .and. (k<=gv%nk_rho_varies)) cycle
459  endif
460 
461  do j=js-1,je ; do i=is,ie
462  coef_y(i,j) = ((scale * khdt_y(i,j))*2.0*(h(i,j,k)*h(i,j+1,k))) / &
463  (h(i,j,k)+h(i,j+1,k)+h_neglect)
464  enddo ; enddo
465 
466  do j=js,je
467  do i=is-1,ie
468  coef_x(i,j) = ((scale * khdt_x(i,j))*2.0*(h(i,j,k)*h(i+1,j,k))) / &
469  (h(i,j,k)+h(i+1,j,k)+h_neglect)
470  enddo
471 
472  do i=is,ie
473  ihdxdy(i,j) = g%IareaT(i,j) / (h(i,j,k)+h_neglect)
474  enddo
475  enddo
476 
477  do m=1,ntr
478  do j=js,je ; do i=is,ie
479  dtr(i,j) = ihdxdy(i,j) * &
480  ((coef_x(i-1,j) * (reg%Tr(m)%t(i-1,j,k) - reg%Tr(m)%t(i,j,k)) - &
481  coef_x(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k))) + &
482  (coef_y(i,j-1) * (reg%Tr(m)%t(i,j-1,k) - reg%Tr(m)%t(i,j,k)) - &
483  coef_y(i,j) * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k))))
484  enddo ; enddo
485  if (associated(reg%Tr(m)%df_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
486  reg%Tr(m)%df_x(i,j,k) = reg%Tr(m)%df_x(i,j,k) + coef_x(i,j) &
487  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i+1,j,k)) * idt
488  enddo ; enddo ; endif
489  if (associated(reg%Tr(m)%df_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
490  reg%Tr(m)%df_y(i,j,k) = reg%Tr(m)%df_y(i,j,k) + coef_y(i,j) &
491  * (reg%Tr(m)%t(i,j,k) - reg%Tr(m)%t(i,j+1,k)) * idt
492  enddo ; enddo ; endif
493  if (associated(reg%Tr(m)%df2d_x)) then ; do j=js,je ; do i=g%IscB,g%IecB
494  reg%Tr(m)%df2d_x(i,j) = reg%Tr(m)%df2d_x(i,j) + 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)%df2d_y)) then ; do j=g%JscB,g%JecB ; do i=is,ie
498  reg%Tr(m)%df2d_y(i,j) = reg%Tr(m)%df2d_y(i,j) + 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  do j=js,je ; do i=is,ie
502  reg%Tr(m)%t(i,j,k) = reg%Tr(m)%t(i,j,k) + dtr(i,j)
503  enddo ; enddo
504  enddo
505 
506  enddo ! End of k loop.
507 
508  enddo ! End of "while" loop.
509 
510  endif ! endif for CS%use_neutral_diffusion
511  call cpu_clock_end(id_clock_diffuse)
512 
513 
514  if (cs%Diffuse_ML_interior) then
515  if (cs%show_call_tree) call calltree_waypoint("Calling epipycnal_ML_diff (tracer_hordiff)")
516  if (cs%debug) call mom_tracer_chksum("Before epipycnal diff ", reg%Tr, ntr, g)
517 
518  call cpu_clock_begin(id_clock_epimix)
519  call tracer_epipycnal_ml_diff(h, dt, reg%Tr, ntr, khdt_x, khdt_y, g, gv, us, &
520  cs, tv, num_itts)
521  call cpu_clock_end(id_clock_epimix)
522  endif
523 
524  if (cs%debug) call mom_tracer_chksum("After tracer diffusion ", reg%Tr, ntr, g)
525 
526  ! post diagnostics for 2d tracer diffusivity
527  if (cs%id_KhTr_u > 0) then
528  do j=js,je ; do i=is-1,ie
529  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
530  enddo ; enddo
531  call post_data(cs%id_KhTr_u, kh_u, cs%diag, mask=g%mask2dCu)
532  endif
533  if (cs%id_KhTr_v > 0) then
534  do j=js-1,je ; do i=is,ie
535  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
536  enddo ; enddo
537  call post_data(cs%id_KhTr_v, kh_v, cs%diag, mask=g%mask2dCv)
538  endif
539  if (cs%id_KhTr_h > 0) then
540  kh_h(:,:) = 0.0
541  do j=js,je ; do i=is-1,ie
542  kh_u(i,j) = g%mask2dCu(i,j)*kh_u(i,j)
543  enddo ; enddo
544  do j=js-1,je ; do i=is,ie
545  kh_v(i,j) = g%mask2dCv(i,j)*kh_v(i,j)
546  enddo ; enddo
547  do j=js,je ; do i=is,ie
548  normalize = 1.0 / ((g%mask2dCu(i-1,j)+g%mask2dCu(i,j)) + &
549  (g%mask2dCv(i,j-1)+g%mask2dCv(i,j)) + gv%H_subroundoff)
550  kh_h(i,j) = normalize*g%mask2dT(i,j)*((kh_u(i-1,j)+kh_u(i,j)) + &
551  (kh_v(i,j-1)+kh_v(i,j)))
552  enddo ; enddo
553  call post_data(cs%id_KhTr_h, kh_h, cs%diag, mask=g%mask2dT)
554  endif
555 
556 
557  if (cs%debug) then
558  call uvchksum("After tracer diffusion khdt_[xy]", khdt_x, khdt_y, &
559  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2)
560  if (cs%use_neutral_diffusion) then
561  call uvchksum("After tracer diffusion Coef_[xy]", coef_x, coef_y, &
562  g%HI, haloshift=0, symmetric=.true., scale=us%L_to_m**2)
563  endif
564  endif
565 
566  if (cs%id_khdt_x > 0) call post_data(cs%id_khdt_x, khdt_x, cs%diag)
567  if (cs%id_khdt_y > 0) call post_data(cs%id_khdt_y, khdt_y, cs%diag)
568 
569  if (cs%show_call_tree) call calltree_leave("tracer_hordiff()")
570 
571 end subroutine tracer_hordiff
572 
573 !> This subroutine does epipycnal diffusion of all tracers between the mixed
574 !! and buffer layers and the interior, using the diffusivity in CS%KhTr.
575 !! Multiple iterations are used (if necessary) so that there is no limit on the
576 !! acceptable time increment.
577 subroutine tracer_epipycnal_ml_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
578  GV, US, CS, tv, num_itts)
579  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
580  type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
581  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
582  real, intent(in) :: dt !< time step [T ~> s]
583  type(tracer_type), intent(inout) :: Tr(:) !< tracer array
584  integer, intent(in) :: ntr !< number of tracers
585  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: khdt_epi_x !< Zonal epipycnal diffusivity times
586  !! a time step and the ratio of the open face width over
587  !! the distance between adjacent tracer points [L2 ~> m2]
588  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: khdt_epi_y !< Meridional epipycnal diffusivity times
589  !! a time step and the ratio of the open face width over
590  !! the distance between adjacent tracer points [L2 ~> m2]
591  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
592  type(tracer_hor_diff_cs), intent(inout) :: CS !< module control structure
593  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic structure
594  integer, intent(in) :: num_itts !< number of iterations (usually=1)
595 
596 
597  real, dimension(SZI_(G), SZJ_(G)) :: &
598  Rml_max ! The maximum coordinate density within the mixed layer [R ~> kg m-3].
599  real, dimension(SZI_(G), SZJ_(G), max(1,GV%nk_rho_varies)) :: &
600  rho_coord ! The coordinate density that is used to mix along [R ~> kg m-3].
601 
602  ! The naming mnemnonic is a=above,b=below,L=Left,R=Right,u=u-point,v=v-point.
603  ! These are 1-D arrays of pointers to 2-d arrays to minimize memory usage.
604  type(p2d), dimension(SZJ_(G)) :: &
605  deep_wt_Lu, deep_wt_Ru, & ! The relative weighting of the deeper of a pair [nondim].
606  hP_Lu, hP_Ru ! The total thickness on each side for each pair [H ~> m or kg m-2].
607 
608  type(p2d), dimension(SZJB_(G)) :: &
609  deep_wt_Lv, deep_wt_Rv, & ! The relative weighting of the deeper of a pair [nondim].
610  hP_Lv, hP_Rv ! The total thickness on each side for each pair [H ~> m or kg m-2].
611 
612  type(p2di), dimension(SZJ_(G)) :: &
613  k0b_Lu, k0a_Lu, & ! The original k-indices of the layers that participate
614  k0b_Ru, k0a_Ru ! in each pair of mixing at u-faces.
615  type(p2di), dimension(SZJB_(G)) :: &
616  k0b_Lv, k0a_Lv, & ! The original k-indices of the layers that participate
617  k0b_Rv, k0a_Rv ! in each pair of mixing at v-faces.
618 
619  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
620  tr_flux_conv ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
621  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tr_flux_3d, Tr_adj_vert_L, Tr_adj_vert_R
622 
623  real, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
624  rho_srt, & ! The density of each layer of the sorted columns [R ~> kg m-3].
625  h_srt ! The thickness of each layer of the sorted columns [H ~> m or kg m-2].
626  integer, dimension(SZI_(G), SZK_(G), SZJ_(G)) :: &
627  k0_srt ! The original k-index that each layer of the sorted column
628  ! corresponds to.
629 
630  real, dimension(SZK_(G)) :: &
631  h_demand_L, & ! The thickness in the left (_L) or right (_R) column that
632  h_demand_R, & ! is demanded to match the thickness in the counterpart [H ~> m or kg m-2].
633  h_used_L, & ! The summed thickness from the left or right columns that
634  h_used_R, & ! have actually been used [H ~> m or kg m-2].
635  h_supply_frac_L, & ! The fraction of the demanded thickness that can
636  h_supply_frac_R ! actually be supplied from a layer.
637  integer, dimension(SZK_(G)) :: &
638  kbs_Lp, & ! The sorted indicies of the Left and Right columns for
639  kbs_Rp ! each pairing.
640 
641  integer, dimension(SZI_(G), SZJ_(G)) :: &
642  num_srt, & ! The number of layers that are sorted in each column.
643  k_end_srt, & ! The maximum index in each column that might need to be
644  ! sorted, based on neighboring values of max_kRho
645  max_krho ! The index of the layer whose target density is just denser
646  ! than the densest part of the mixed layer.
647  integer, dimension(SZJ_(G)) :: &
648  max_srt ! The maximum value of num_srt in a k-row.
649  integer, dimension(SZIB_(G), SZJ_(G)) :: &
650  nPu ! The number of epipycnal pairings at each u-point.
651  integer, dimension(SZI_(G), SZJB_(G)) :: &
652  nPv ! The number of epipycnal pairings at each v-point.
653  real :: h_exclude ! A thickness that layers must attain to be considered
654  ! for inclusion in mixing [H ~> m or kg m-2].
655  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
656  real :: I_maxitt ! The inverse of the maximum number of iterations.
657  real :: rho_pair, rho_a, rho_b ! Temporary densities [R ~> kg m-3].
658  real :: Tr_min_face ! The minimum and maximum tracer concentrations
659  real :: Tr_max_face ! associated with a pairing [Conc]
660  real :: Tr_La, Tr_Lb ! The 4 tracer concentrations that might be
661  real :: Tr_Ra, Tr_Rb ! associated with a pairing [Conc]
662  real :: Tr_av_L ! The average tracer concentrations on the left and right
663  real :: Tr_av_R ! sides of a pairing [Conc].
664  real :: Tr_flux ! The tracer flux from left to right in a pair [conc H L2 ~> conc m3 or conc kg].
665  real :: Tr_adj_vert ! A downward vertical adjustment to Tr_flux between the
666  ! two cells that make up one side of the pairing [conc H L2 ~> conc m3 or conc kg].
667  real :: h_L, h_R ! Thicknesses to the left and right [H ~> m or kg m-2].
668  real :: wt_a, wt_b ! Fractional weights of layers above and below [nondim].
669  real :: vol ! A cell volume or mass [H L2 ~> m3 or kg].
670  logical, dimension(SZK_(G)) :: &
671  left_set, & ! If true, the left or right point determines the density of
672  right_set ! of the trio. If densities are exactly equal, both are true.
673  real :: tmp
674  real :: p_ref_cv(SZI_(G))
675 
676  integer :: k_max, k_min, k_test, itmp
677  integer :: i, j, k, k2, m, is, ie, js, je, nz, nkmb
678  integer :: isd, ied, jsd, jed, IsdB, IedB, k_size
679  integer :: kL, kR, kLa, kLb, kRa, kRb, nP, itt, ns, max_itt
680  integer :: PEmax_kRho
681  integer :: isv, iev, jsv, jev ! The valid range of the indices.
682 
683  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
684  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
685  isdb = g%IsdB ; iedb = g%IedB
686  idt = 1.0 / dt
687  nkmb = gv%nk_rho_varies
688 
689  if (num_itts <= 1) then
690  max_itt = 1 ; i_maxitt = 1.0
691  else
692  max_itt = num_itts ; i_maxitt = 1.0 / (real(max_itt))
693  endif
694 
695  do i=is-2,ie+2 ; p_ref_cv(i) = tv%P_Ref ; enddo
696 
697  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
698  ! Determine which layers the mixed- and buffer-layers map into...
699  !$OMP parallel do default(shared)
700  do k=1,nkmb ; do j=js-2,je+2
701  call calculate_density(tv%T(:,j,k),tv%S(:,j,k), p_ref_cv, &
702  rho_coord(:,j,k), is-2, ie-is+5, tv%eqn_of_state, scale=us%kg_m3_to_R)
703  enddo ; enddo
704 
705  do j=js-2,je+2 ; do i=is-2,ie+2
706  rml_max(i,j) = rho_coord(i,j,1)
707  num_srt(i,j) = 0 ; max_krho(i,j) = 0
708  enddo ; enddo
709  do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2
710  if (rml_max(i,j) < rho_coord(i,j,k)) rml_max(i,j) = rho_coord(i,j,k)
711  enddo ; enddo ; enddo
712  ! Use bracketing and bisection to find the k-level that the densest of the
713  ! mixed and buffer layer corresponds to, such that:
714  ! GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
715 !$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,G,GV,Rml_max,max_kRho) &
716 !$OMP private(k_min,k_max,k_test)
717  do j=js-2,je+2 ; do i=is-2,ie+2 ; if (g%mask2dT(i,j) > 0.5) then
718  if (rml_max(i,j) > gv%Rlay(nz)) then ; max_krho(i,j) = nz+1
719  elseif (rml_max(i,j) <= gv%Rlay(nkmb+1)) then ; max_krho(i,j) = nkmb+1
720  else
721  k_min = nkmb+2 ; k_max = nz
722  do
723  k_test = (k_min + k_max) / 2
724  if (rml_max(i,j) <= gv%Rlay(k_test-1)) then ; k_max = k_test-1
725  elseif (gv%Rlay(k_test) < rml_max(i,j)) then ; k_min = k_test+1
726  else ; max_krho(i,j) = k_test ; exit ; endif
727 
728  if (k_min == k_max) then ; max_krho(i,j) = k_max ; exit ; endif
729  enddo
730  endif
731  endif ; enddo ; enddo
732 
733  pemax_krho = 0
734  do j=js-1,je+1 ; do i=is-1,ie+1
735  k_end_srt(i,j) = max(max_krho(i,j), max_krho(i-1,j), max_krho(i+1,j), &
736  max_krho(i,j-1), max_krho(i,j+1))
737  if (pemax_krho < k_end_srt(i,j)) pemax_krho = k_end_srt(i,j)
738  enddo ; enddo
739  if (pemax_krho > nz) pemax_krho = nz ! PEmax_kRho could have been nz+1.
740 
741  h_exclude = 10.0*(gv%Angstrom_H + gv%H_subroundoff)
742 !$OMP parallel default(none) shared(is,ie,js,je,nkmb,G,GV,h,h_exclude,num_srt,k0_srt, &
743 !$OMP rho_srt,h_srt,PEmax_kRho,k_end_srt,rho_coord,max_srt) &
744 !$OMP private(ns,tmp,itmp)
745 !$OMP do
746  do j=js-1,je+1
747  do k=1,nkmb ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
748  if (h(i,j,k) > h_exclude) then
749  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
750  k0_srt(i,ns,j) = k
751  rho_srt(i,ns,j) = rho_coord(i,j,k)
752  h_srt(i,ns,j) = h(i,j,k)
753  endif
754  endif ; enddo ; enddo
755  do k=nkmb+1,pemax_krho ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.5) then
756  if ((k<=k_end_srt(i,j)) .and. (h(i,j,k) > h_exclude)) then
757  num_srt(i,j) = num_srt(i,j) + 1 ; ns = num_srt(i,j)
758  k0_srt(i,ns,j) = k
759  rho_srt(i,ns,j) = gv%Rlay(k)
760  h_srt(i,ns,j) = h(i,j,k)
761  endif
762  endif ; enddo ; enddo
763  enddo
764  ! Sort each column by increasing density. This should already be close,
765  ! and the size of the arrays are small, so straight insertion is used.
766 !$OMP do
767  do j=js-1,je+1; do i=is-1,ie+1
768  do k=2,num_srt(i,j) ; if (rho_srt(i,k,j) < rho_srt(i,k-1,j)) then
769  ! The last segment needs to be shuffled earlier in the list.
770  do k2 = k,2,-1 ; if (rho_srt(i,k2,j) >= rho_srt(i,k2-1,j)) exit
771  itmp = k0_srt(i,k2-1,j) ; k0_srt(i,k2-1,j) = k0_srt(i,k2,j) ; k0_srt(i,k2,j) = itmp
772  tmp = rho_srt(i,k2-1,j) ; rho_srt(i,k2-1,j) = rho_srt(i,k2,j) ; rho_srt(i,k2,j) = tmp
773  tmp = h_srt(i,k2-1,j) ; h_srt(i,k2-1,j) = h_srt(i,k2,j) ; h_srt(i,k2,j) = tmp
774  enddo
775  endif ; enddo
776  enddo ; enddo
777 !$OMP do
778  do j=js-1,je+1
779  max_srt(j) = 0
780  do i=is-1,ie+1 ; max_srt(j) = max(max_srt(j), num_srt(i,j)) ; enddo
781  enddo
782 !$OMP end parallel
783 
784  do j=js,je
785  k_size = max(2*max_srt(j),1)
786  allocate(deep_wt_lu(j)%p(isdb:iedb,k_size))
787  allocate(deep_wt_ru(j)%p(isdb:iedb,k_size))
788  allocate(hp_lu(j)%p(isdb:iedb,k_size))
789  allocate(hp_ru(j)%p(isdb:iedb,k_size))
790  allocate(k0a_lu(j)%p(isdb:iedb,k_size))
791  allocate(k0a_ru(j)%p(isdb:iedb,k_size))
792  allocate(k0b_lu(j)%p(isdb:iedb,k_size))
793  allocate(k0b_ru(j)%p(isdb:iedb,k_size))
794  enddo
795 
796 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lu,k0_srt, &
797 !$OMP k0b_Ru,k0a_Lu,k0a_Ru,deep_wt_Lu,deep_wt_Ru, &
798 !$OMP h_srt,nkmb,nPu,hP_Lu,hP_Ru) &
799 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
800 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
801 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
802 !$OMP h_supply_frac_L)
803  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
804  ! Set up the pairings for fluxes through the zonal faces.
805 
806  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
807  do k=1,num_srt(i+1,j) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
808 
809  ! First merge the left and right lists into a single, sorted list.
810 
811  ! Discard any layers that are lighter than the lightest in the other
812  ! column. They can only participate in mixing as the lighter part of a
813  ! pair of points.
814  if (rho_srt(i,1,j) < rho_srt(i+1,1,j)) then
815  kr = 1
816  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i+1,1,j)) exit ; enddo
817  elseif (rho_srt(i+1,1,j) < rho_srt(i,1,j)) then
818  kl = 1
819  do kr=2,num_srt(i+1,j) ; if (rho_srt(i+1,kr,j) >= rho_srt(i,1,j)) exit ; enddo
820  else
821  kl = 1 ; kr = 1
822  endif
823  np = 0
824  do ! Loop to accumulate pairs of columns.
825  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i+1,j))) exit
826 
827  if (rho_srt(i,kl,j) > rho_srt(i+1,kr,j)) then
828  ! The right point is lighter and defines the density for this trio.
829  np = np+1 ; k = np
830  rho_pair = rho_srt(i+1,kr,j)
831 
832  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
833  k0a_lu(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
834  kbs_lp(k) = kl ; kbs_rp(k) = kr
835 
836  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
837  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
838  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
839  deep_wt_lu(j)%p(i,k) = wt_b ; deep_wt_ru(j)%p(i,k) = 1.0
840 
841  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j) * wt_b
842  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i+1,kr,j) * (1.0-wt_b)
843 
844  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
845  elseif (rho_srt(i,kl,j) < rho_srt(i+1,kr,j)) then
846  ! The left point is lighter and defines the density for this trio.
847  np = np+1 ; k = np
848  rho_pair = rho_srt(i,kl,j)
849  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
850  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)
851 
852  kbs_lp(k) = kl ; kbs_rp(k) = kr
853 
854  rho_a = rho_srt(i+1,kr-1,j) ; rho_b = rho_srt(i+1,kr,j)
855  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
856  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
857  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = wt_b
858 
859  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
860  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
861 
862  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
863  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i+1,kr,j) <= nkmb)) then
864  ! The densities are exactly equal and one layer is above the interior.
865  np = np+1 ; k = np
866  k0b_lu(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_ru(j)%p(i,k) = k0_srt(i+1,kr,j)
867  k0a_lu(j)%p(i,k) = k0b_lu(j)%p(i,k) ; k0a_ru(j)%p(i,k) = k0b_ru(j)%p(i,k)
868  kbs_lp(k) = kl ; kbs_rp(k) = kr
869  deep_wt_lu(j)%p(i,k) = 1.0 ; deep_wt_ru(j)%p(i,k) = 1.0
870 
871  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
872  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
873 
874  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
875  else ! The densities are exactly equal and in the interior.
876  ! Mixing in this case has already occurred, so accumulate the thickness
877  ! demanded for that mixing and skip onward.
878  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i+1,kr,j)
879  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
880 
881  kl = kl+1 ; kr = kr+1
882  endif
883  enddo ! Loop to accumulate pairs of columns.
884  npu(i,j) = np ! This is the number of active pairings.
885 
886  ! Determine what fraction of the thickness "demand" can be supplied.
887  do k=1,num_srt(i+1,j)
888  h_supply_frac_r(k) = 1.0
889  if (h_demand_r(k) > 0.5*h_srt(i+1,k,j)) &
890  h_supply_frac_r(k) = 0.5*h_srt(i+1,k,j) / h_demand_r(k)
891  enddo
892  do k=1,num_srt(i,j)
893  h_supply_frac_l(k) = 1.0
894  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
895  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
896  enddo
897 
898  ! Distribute the "exported" thicknesses proportionately.
899  do k=1,npu(i,j)
900  kl = kbs_lp(k) ; kr = kbs_rp(k)
901  hp_lu(j)%p(i,k) = 0.0 ; hp_ru(j)%p(i,k) = 0.0
902  if (left_set(k)) then ! Add the contributing thicknesses on the right.
903  if (deep_wt_ru(j)%p(i,k) < 1.0) then
904  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))
905  wt_b = deep_wt_ru(j)%p(i,k)
906  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b)*hp_ru(j)%p(i,k)
907  h_used_r(kr) = h_used_r(kr) + wt_b*hp_ru(j)%p(i,k)
908  else
909  hp_ru(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
910  h_used_r(kr) = h_used_r(kr) + hp_ru(j)%p(i,k)
911  endif
912  endif
913  if (right_set(k)) then ! Add the contributing thicknesses on the left.
914  if (deep_wt_lu(j)%p(i,k) < 1.0) then
915  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))
916  wt_b = deep_wt_lu(j)%p(i,k)
917  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b)*hp_lu(j)%p(i,k)
918  h_used_l(kl) = h_used_l(kl) + wt_b*hp_lu(j)%p(i,k)
919  else
920  hp_lu(j)%p(i,k) = 0.5*h_srt(i+1,kr,j) * h_supply_frac_l(kl)
921  h_used_l(kl) = h_used_l(kl) + hp_lu(j)%p(i,k)
922  endif
923  endif
924  enddo
925 
926  ! The left-over thickness (at least half the layer thickness) is now
927  ! added to the thicknesses of the importing columns.
928  do k=1,npu(i,j)
929  if (left_set(k)) hp_lu(j)%p(i,k) = hp_lu(j)%p(i,k) + &
930  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
931  if (right_set(k)) hp_ru(j)%p(i,k) = hp_ru(j)%p(i,k) + &
932  (h_srt(i+1,kbs_rp(k),j) - h_used_r(kbs_rp(k)))
933  enddo
934 
935  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
936 
937  do j=js-1,je
938  k_size = max(max_srt(j)+max_srt(j+1),1)
939  allocate(deep_wt_lv(j)%p(isd:ied,k_size))
940  allocate(deep_wt_rv(j)%p(isd:ied,k_size))
941  allocate(hp_lv(j)%p(isd:ied,k_size))
942  allocate(hp_rv(j)%p(isd:ied,k_size))
943  allocate(k0a_lv(j)%p(isd:ied,k_size))
944  allocate(k0a_rv(j)%p(isd:ied,k_size))
945  allocate(k0b_lv(j)%p(isd:ied,k_size))
946  allocate(k0b_rv(j)%p(isd:ied,k_size))
947  enddo
948 
949 !$OMP parallel do default(none) shared(is,ie,js,je,G,num_srt,rho_srt,k0b_Lv,k0b_Rv, &
950 !$OMP k0_srt,k0a_Lv,k0a_Rv,deep_wt_Lv,deep_wt_Rv, &
951 !$OMP h_srt,nkmb,nPv,hP_Lv,hP_Rv) &
952 !$OMP private(h_demand_L,h_used_L,h_demand_R,h_used_R, &
953 !$OMP kR,kL,nP,rho_pair,kbs_Lp,kbs_Rp,rho_a,rho_b, &
954 !$OMP wt_b,left_set,right_set,h_supply_frac_R, &
955 !$OMP h_supply_frac_L)
956  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
957  ! Set up the pairings for fluxes through the meridional faces.
958 
959  do k=1,num_srt(i,j) ; h_demand_l(k) = 0.0 ; h_used_l(k) = 0.0 ; enddo
960  do k=1,num_srt(i,j+1) ; h_demand_r(k) = 0.0 ; h_used_r(k) = 0.0 ; enddo
961 
962  ! First merge the left and right lists into a single, sorted list.
963 
964  ! Discard any layers that are lighter than the lightest in the other
965  ! column. They can only participate in mixing as the lighter part of a
966  ! pair of points.
967  if (rho_srt(i,1,j) < rho_srt(i,1,j+1)) then
968  kr = 1
969  do kl=2,num_srt(i,j) ; if (rho_srt(i,kl,j) >= rho_srt(i,1,j+1)) exit ; enddo
970  elseif (rho_srt(i,1,j+1) < rho_srt(i,1,j)) then
971  kl = 1
972  do kr=2,num_srt(i,j+1) ; if (rho_srt(i,kr,j+1) >= rho_srt(i,1,j)) exit ; enddo
973  else
974  kl = 1 ; kr = 1
975  endif
976  np = 0
977  do ! Loop to accumulate pairs of columns.
978  if ((kl > num_srt(i,j)) .or. (kr > num_srt(i,j+1))) exit
979 
980  if (rho_srt(i,kl,j) > rho_srt(i,kr,j+1)) then
981  ! The right point is lighter and defines the density for this trio.
982  np = np+1 ; k = np
983  rho_pair = rho_srt(i,kr,j+1)
984 
985  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
986  k0a_lv(j)%p(i,k) = k0_srt(i,kl-1,j) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
987  kbs_lp(k) = kl ; kbs_rp(k) = kr
988 
989  rho_a = rho_srt(i,kl-1,j) ; rho_b = rho_srt(i,kl,j)
990  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
991  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
992  deep_wt_lv(j)%p(i,k) = wt_b ; deep_wt_rv(j)%p(i,k) = 1.0
993 
994  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1) * wt_b
995  h_demand_l(kl-1) = h_demand_l(kl-1) + 0.5*h_srt(i,kr,j+1) * (1.0-wt_b)
996 
997  kr = kr+1 ; left_set(k) = .false. ; right_set(k) = .true.
998  elseif (rho_srt(i,kl,j) < rho_srt(i,kr,j+1)) then
999  ! The left point is lighter and defines the density for this trio.
1000  np = np+1 ; k = np
1001  rho_pair = rho_srt(i,kl,j)
1002  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1003  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)
1004 
1005  kbs_lp(k) = kl ; kbs_rp(k) = kr
1006 
1007  rho_a = rho_srt(i,kr-1,j+1) ; rho_b = rho_srt(i,kr,j+1)
1008  wt_b = 1.0 ; if (abs(rho_a - rho_b) > abs(rho_pair - rho_a)) &
1009  wt_b = (rho_pair - rho_a) / (rho_b - rho_a)
1010  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = wt_b
1011 
1012  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j) * wt_b
1013  h_demand_r(kr-1) = h_demand_r(kr-1) + 0.5*h_srt(i,kl,j) * (1.0-wt_b)
1014 
1015  kl = kl+1 ; left_set(k) = .true. ; right_set(k) = .false.
1016  elseif ((k0_srt(i,kl,j) <= nkmb) .or. (k0_srt(i,kr,j+1) <= nkmb)) then
1017  ! The densities are exactly equal and one layer is above the interior.
1018  np = np+1 ; k = np
1019  k0b_lv(j)%p(i,k) = k0_srt(i,kl,j) ; k0b_rv(j)%p(i,k) = k0_srt(i,kr,j+1)
1020  k0a_lv(j)%p(i,k) = k0b_lv(j)%p(i,k) ; k0a_rv(j)%p(i,k) = k0b_rv(j)%p(i,k)
1021  kbs_lp(k) = kl ; kbs_rp(k) = kr
1022  deep_wt_lv(j)%p(i,k) = 1.0 ; deep_wt_rv(j)%p(i,k) = 1.0
1023 
1024  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1025  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1026 
1027  kl = kl+1 ; kr = kr+1 ; left_set(k) = .true. ; right_set(k) = .true.
1028  else ! The densities are exactly equal and in the interior.
1029  ! Mixing in this case has already occurred, so accumulate the thickness
1030  ! demanded for that mixing and skip onward.
1031  h_demand_l(kl) = h_demand_l(kl) + 0.5*h_srt(i,kr,j+1)
1032  h_demand_r(kr) = h_demand_r(kr) + 0.5*h_srt(i,kl,j)
1033 
1034  kl = kl+1 ; kr = kr+1
1035  endif
1036  enddo ! Loop to accumulate pairs of columns.
1037  npv(i,j) = np ! This is the number of active pairings.
1038 
1039  ! Determine what fraction of the thickness "demand" can be supplied.
1040  do k=1,num_srt(i,j+1)
1041  h_supply_frac_r(k) = 1.0
1042  if (h_demand_r(k) > 0.5*h_srt(i,k,j+1)) &
1043  h_supply_frac_r(k) = 0.5*h_srt(i,k,j+1) / h_demand_r(k)
1044  enddo
1045  do k=1,num_srt(i,j)
1046  h_supply_frac_l(k) = 1.0
1047  if (h_demand_l(k) > 0.5*h_srt(i,k,j)) &
1048  h_supply_frac_l(k) = 0.5*h_srt(i,k,j) / h_demand_l(k)
1049  enddo
1050 
1051  ! Distribute the "exported" thicknesses proportionately.
1052  do k=1,npv(i,j)
1053  kl = kbs_lp(k) ; kr = kbs_rp(k)
1054  hp_lv(j)%p(i,k) = 0.0 ; hp_rv(j)%p(i,k) = 0.0
1055  if (left_set(k)) then ! Add the contributing thicknesses on the right.
1056  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1057  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))
1058  wt_b = deep_wt_rv(j)%p(i,k)
1059  h_used_r(kr-1) = h_used_r(kr-1) + (1.0 - wt_b) * hp_rv(j)%p(i,k)
1060  h_used_r(kr) = h_used_r(kr) + wt_b * hp_rv(j)%p(i,k)
1061  else
1062  hp_rv(j)%p(i,k) = 0.5*h_srt(i,kl,j) * h_supply_frac_r(kr)
1063  h_used_r(kr) = h_used_r(kr) + hp_rv(j)%p(i,k)
1064  endif
1065  endif
1066  if (right_set(k)) then ! Add the contributing thicknesses on the left.
1067  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1068  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))
1069  wt_b = deep_wt_lv(j)%p(i,k)
1070  h_used_l(kl-1) = h_used_l(kl-1) + (1.0 - wt_b) * hp_lv(j)%p(i,k)
1071  h_used_l(kl) = h_used_l(kl) + wt_b * hp_lv(j)%p(i,k)
1072  else
1073  hp_lv(j)%p(i,k) = 0.5*h_srt(i,kr,j+1) * h_supply_frac_l(kl)
1074  h_used_l(kl) = h_used_l(kl) + hp_lv(j)%p(i,k)
1075  endif
1076  endif
1077  enddo
1078 
1079  ! The left-over thickness (at least half the layer thickness) is now
1080  ! added to the thicknesses of the importing columns.
1081  do k=1,npv(i,j)
1082  if (left_set(k)) hp_lv(j)%p(i,k) = hp_lv(j)%p(i,k) + &
1083  (h_srt(i,kbs_lp(k),j) - h_used_l(kbs_lp(k)))
1084  if (right_set(k)) hp_rv(j)%p(i,k) = hp_rv(j)%p(i,k) + &
1085  (h_srt(i,kbs_rp(k),j+1) - h_used_r(kbs_rp(k)))
1086  enddo
1087 
1088 
1089  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1090 
1091 ! The tracer-specific calculations start here.
1092 
1093  ! Zero out tracer tendencies.
1094  do k=1,pemax_krho ; do j=js-1,je+1 ; do i=is-1,ie+1
1095  tr_flux_conv(i,j,k) = 0.0
1096  enddo ; enddo ; enddo
1097 
1098  do itt=1,max_itt
1099 
1100  if (itt > 1) then ! The halos have already been filled if itt==1.
1101  call do_group_pass(cs%pass_t, g%Domain, clock=id_clock_pass)
1102  endif
1103 
1104  do m=1,ntr
1105 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, &
1106 !$OMP k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru, &
1107 !$OMP hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) &
1108 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
1109 !$OMP Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
1110 !$OMP Tr_flux,Tr_adj_vert,wt_a,vol)
1111  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.5) then
1112  ! Determine the fluxes through the zonal faces.
1113 
1114  ! Find the acceptable range of tracer concentration around this face.
1115  if (npu(i,j) >= 1) then
1116  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1117  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i+1,j,1))
1118  do k=2,nkmb
1119  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1120  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i+1,j,k))
1121  enddo
1122 
1123  ! Include the next two layers denser than the densest buffer layer.
1124  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1125  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1126  kra = nkmb+1 ; if (max_krho(i+1,j) < nz+1) kra = max_krho(i+1,j)
1127  krb = kra ; if (max_krho(i+1,j) < nz) krb = max_krho(i+1,j)+1
1128  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1129  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1130  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1131  if (h(i+1,j,kra) > h_exclude) tr_ra = tr(m)%t(i+1,j,kra)
1132  if (h(i+1,j,krb) > h_exclude) tr_rb = tr(m)%t(i+1,j,krb)
1133  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1134  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1135 
1136  ! Include all points in diffusive pairings at this face.
1137  do k=1,npu(i,j)
1138  tr_lb = tr(m)%t(i,j,k0b_lu(j)%p(i,k))
1139  tr_rb = tr(m)%t(i+1,j,k0b_ru(j)%p(i,k))
1140  tr_la = tr_lb ; tr_ra = tr_rb
1141  if (deep_wt_lu(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lu(j)%p(i,k))
1142  if (deep_wt_ru(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i+1,j,k0a_ru(j)%p(i,k))
1143  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1144  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1145  enddo
1146  endif
1147 
1148  do k=1,npu(i,j)
1149  klb = k0b_lu(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1150  if (deep_wt_lu(j)%p(i,k) < 1.0) then
1151  kla = k0a_lu(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1152  wt_b = deep_wt_lu(j)%p(i,k)
1153  tr_av_l = wt_b*tr_lb + (1.0-wt_b)*tr_la
1154  endif
1155 
1156  krb = k0b_ru(j)%p(i,k) ; tr_rb = tr(m)%t(i+1,j,krb) ; tr_av_r = tr_rb
1157  if (deep_wt_ru(j)%p(i,k) < 1.0) then
1158  kra = k0a_ru(j)%p(i,k) ; tr_ra = tr(m)%t(i+1,j,kra)
1159  wt_b = deep_wt_ru(j)%p(i,k)
1160  tr_av_r = wt_b*tr_rb + (1.0-wt_b)*tr_ra
1161  endif
1162 
1163  h_l = hp_lu(j)%p(i,k) ; h_r = hp_ru(j)%p(i,k)
1164  tr_flux = i_maxitt * khdt_epi_x(i,j) * (tr_av_l - tr_av_r) * &
1165  ((2.0 * h_l * h_r) / (h_l + h_r))
1166 
1167 
1168  if (deep_wt_lu(j)%p(i,k) >= 1.0) then
1169  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux
1170  else
1171  tr_adj_vert = 0.0
1172  wt_b = deep_wt_lu(j)%p(i,k) ; wt_a = 1.0 - wt_b
1173  vol = hp_lu(j)%p(i,k) * g%areaT(i,j)
1174 
1175  ! Ensure that the tracer flux does not drive the tracer values
1176  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1177  ! does that the concentration in both contributing peices exceed
1178  ! this range equally. With downgradient fluxes and the initial tracer
1179  ! concentrations determining the valid range, the latter condition
1180  ! only enters for large values of the effective diffusive CFL number.
1181  if (tr_flux > 0.0) then
1182  if (tr_la < tr_lb) then ; if (vol*(tr_la-tr_min_face) < tr_flux) &
1183  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1184  (vol*wt_b) * (tr_lb - tr_la))
1185  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1186  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1187  (vol*wt_a) * (tr_la - tr_lb))
1188  endif
1189  elseif (tr_flux < 0.0) then
1190  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1191  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1192  (vol*wt_b) * (tr_la - tr_lb))
1193  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1194  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1195  (vol*wt_a)*(tr_lb - tr_la))
1196  endif
1197  endif
1198 
1199  tr_flux_conv(i,j,kla) = tr_flux_conv(i,j,kla) - (wt_a*tr_flux + tr_adj_vert)
1200  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - (wt_b*tr_flux - tr_adj_vert)
1201  endif
1202 
1203  if (deep_wt_ru(j)%p(i,k) >= 1.0) then
1204  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + tr_flux
1205  else
1206  tr_adj_vert = 0.0
1207  wt_b = deep_wt_ru(j)%p(i,k) ; wt_a = 1.0 - wt_b
1208  vol = hp_ru(j)%p(i,k) * g%areaT(i+1,j)
1209 
1210  ! Ensure that the tracer flux does not drive the tracer values
1211  ! outside of the range Tr_min_face <= Tr <= Tr_max_face, or if it
1212  ! does that the concentration in both contributing peices exceed
1213  ! this range equally. With downgradient fluxes and the initial tracer
1214  ! concentrations determining the valid range, the latter condition
1215  ! only enters for large values of the effective diffusive CFL number.
1216  if (tr_flux < 0.0) then
1217  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1218  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1219  (vol*wt_b) * (tr_rb - tr_ra))
1220  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1221  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1222  (vol*wt_a) * (tr_ra - tr_rb))
1223  endif
1224  elseif (tr_flux > 0.0) then
1225  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1226  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1227  (vol*wt_b) * (tr_ra - tr_rb))
1228  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1229  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1230  (vol*wt_a)*(tr_rb - tr_ra))
1231  endif
1232  endif
1233 
1234  tr_flux_conv(i+1,j,kra) = tr_flux_conv(i+1,j,kra) + &
1235  (wt_a*tr_flux - tr_adj_vert)
1236  tr_flux_conv(i+1,j,krb) = tr_flux_conv(i+1,j,krb) + &
1237  (wt_b*tr_flux + tr_adj_vert)
1238  endif
1239  if (associated(tr(m)%df2d_x)) &
1240  tr(m)%df2d_x(i,j) = tr(m)%df2d_x(i,j) + tr_flux * idt
1241  enddo ! Loop over pairings at faces.
1242  endif ; enddo ; enddo ! i- & j- loops over zonal faces.
1243 
1244 !$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, &
1245 !$OMP k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv, &
1246 !$OMP hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d, &
1247 !$OMP Tr_adj_vert_L,Tr_adj_vert_R,Idt) &
1248 !$OMP private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb, &
1249 !$OMP Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R, &
1250 !$OMP h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
1251  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
1252  ! Determine the fluxes through the meridional faces.
1253 
1254  ! Find the acceptable range of tracer concentration around this face.
1255  if (npv(i,j) >= 1) then
1256  tr_min_face = min(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1257  tr_max_face = max(tr(m)%t(i,j,1), tr(m)%t(i,j+1,1))
1258  do k=2,nkmb
1259  tr_min_face = min(tr_min_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1260  tr_max_face = max(tr_max_face, tr(m)%t(i,j,k), tr(m)%t(i,j+1,k))
1261  enddo
1262 
1263  ! Include the next two layers denser than the densest buffer layer.
1264  kla = nkmb+1 ; if (max_krho(i,j) < nz+1) kla = max_krho(i,j)
1265  klb = kla ; if (max_krho(i,j) < nz) klb = max_krho(i,j)+1
1266  kra = nkmb+1 ; if (max_krho(i,j+1) < nz+1) kra = max_krho(i,j+1)
1267  krb = kra ; if (max_krho(i,j+1) < nz) krb = max_krho(i,j+1)+1
1268  tr_la = tr_min_face ; tr_lb = tr_la ; tr_ra = tr_la ; tr_rb = tr_la
1269  if (h(i,j,kla) > h_exclude) tr_la = tr(m)%t(i,j,kla)
1270  if (h(i,j,klb) > h_exclude) tr_la = tr(m)%t(i,j,klb)
1271  if (h(i,j+1,kra) > h_exclude) tr_ra = tr(m)%t(i,j+1,kra)
1272  if (h(i,j+1,krb) > h_exclude) tr_rb = tr(m)%t(i,j+1,krb)
1273  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1274  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1275 
1276  ! Include all points in diffusive pairings at this face.
1277  do k=1,npv(i,j)
1278  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))
1279  tr_la = tr_lb ; tr_ra = tr_rb
1280  if (deep_wt_lv(j)%p(i,k) < 1.0) tr_la = tr(m)%t(i,j,k0a_lv(j)%p(i,k))
1281  if (deep_wt_rv(j)%p(i,k) < 1.0) tr_ra = tr(m)%t(i,j+1,k0a_rv(j)%p(i,k))
1282  tr_min_face = min(tr_min_face, tr_la, tr_lb, tr_ra, tr_rb)
1283  tr_max_face = max(tr_max_face, tr_la, tr_lb, tr_ra, tr_rb)
1284  enddo
1285  endif
1286 
1287  do k=1,npv(i,j)
1288  klb = k0b_lv(j)%p(i,k) ; tr_lb = tr(m)%t(i,j,klb) ; tr_av_l = tr_lb
1289  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1290  kla = k0a_lv(j)%p(i,k) ; tr_la = tr(m)%t(i,j,kla)
1291  wt_b = deep_wt_lv(j)%p(i,k)
1292  tr_av_l = wt_b * tr_lb + (1.0-wt_b) * tr_la
1293  endif
1294 
1295  krb = k0b_rv(j)%p(i,k) ; tr_rb = tr(m)%t(i,j+1,krb) ; tr_av_r = tr_rb
1296  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1297  kra = k0a_rv(j)%p(i,k) ; tr_ra = tr(m)%t(i,j+1,kra)
1298  wt_b = deep_wt_rv(j)%p(i,k)
1299  tr_av_r = wt_b * tr_rb + (1.0-wt_b) * tr_ra
1300  endif
1301 
1302  h_l = hp_lv(j)%p(i,k) ; h_r = hp_rv(j)%p(i,k)
1303  tr_flux = i_maxitt * ((2.0 * h_l * h_r) / (h_l + h_r)) * &
1304  khdt_epi_y(i,j) * (tr_av_l - tr_av_r)
1305  tr_flux_3d(i,j,k) = tr_flux
1306 
1307  if (deep_wt_lv(j)%p(i,k) < 1.0) then
1308  tr_adj_vert = 0.0
1309  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1310  vol = hp_lv(j)%p(i,k) * g%areaT(i,j)
1311 
1312  ! Ensure that the tracer flux does not drive the tracer values
1313  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1314  if (tr_flux > 0.0) then
1315  if (tr_la < tr_lb) then ; if (vol * (tr_la-tr_min_face) < tr_flux) &
1316  tr_adj_vert = -wt_a * min(tr_flux - vol * (tr_la-tr_min_face), &
1317  (vol*wt_b) * (tr_lb - tr_la))
1318  else ; if (vol*(tr_lb-tr_min_face) < tr_flux) &
1319  tr_adj_vert = wt_b * min(tr_flux - vol * (tr_lb-tr_min_face), &
1320  (vol*wt_a) * (tr_la - tr_lb))
1321  endif
1322  elseif (tr_flux < 0.0) then
1323  if (tr_la > tr_lb) then ; if (vol * (tr_max_face-tr_la) < -tr_flux) &
1324  tr_adj_vert = wt_a * min(-tr_flux - vol * (tr_max_face-tr_la), &
1325  (vol*wt_b) * (tr_la - tr_lb))
1326  else ; if (vol*(tr_max_face-tr_lb) < -tr_flux) &
1327  tr_adj_vert = -wt_b * min(-tr_flux - vol * (tr_max_face-tr_lb), &
1328  (vol*wt_a)*(tr_lb - tr_la))
1329  endif
1330  endif
1331  tr_adj_vert_l(i,j,k) = tr_adj_vert
1332  endif
1333 
1334  if (deep_wt_rv(j)%p(i,k) < 1.0) then
1335  tr_adj_vert = 0.0
1336  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1337  vol = hp_rv(j)%p(i,k) * g%areaT(i,j+1)
1338 
1339  ! Ensure that the tracer flux does not drive the tracer values
1340  ! outside of the range Tr_min_face <= Tr <= Tr_max_face.
1341  if (tr_flux < 0.0) then
1342  if (tr_ra < tr_rb) then ; if (vol * (tr_ra-tr_min_face) < -tr_flux) &
1343  tr_adj_vert = -wt_a * min(-tr_flux - vol * (tr_ra-tr_min_face), &
1344  (vol*wt_b) * (tr_rb - tr_ra))
1345  else ; if (vol*(tr_rb-tr_min_face) < (-tr_flux)) &
1346  tr_adj_vert = wt_b * min(-tr_flux - vol * (tr_rb-tr_min_face), &
1347  (vol*wt_a) * (tr_ra - tr_rb))
1348  endif
1349  elseif (tr_flux > 0.0) then
1350  if (tr_ra > tr_rb) then ; if (vol * (tr_max_face-tr_ra) < tr_flux) &
1351  tr_adj_vert = wt_a * min(tr_flux - vol * (tr_max_face-tr_ra), &
1352  (vol*wt_b) * (tr_ra - tr_rb))
1353  else ; if (vol*(tr_max_face-tr_rb) < tr_flux) &
1354  tr_adj_vert = -wt_b * min(tr_flux - vol * (tr_max_face-tr_rb), &
1355  (vol*wt_a)*(tr_rb - tr_ra))
1356  endif
1357  endif
1358  tr_adj_vert_r(i,j,k) = tr_adj_vert
1359  endif
1360  if (associated(tr(m)%df2d_y)) &
1361  tr(m)%df2d_y(i,j) = tr(m)%df2d_y(i,j) + tr_flux * idt
1362  enddo ! Loop over pairings at faces.
1363  endif ; enddo ; enddo ! i- & j- loops over meridional faces.
1364 !$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv, &
1365 !$OMP tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
1366 !$OMP deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
1367 !$OMP private(kLa,kLb,kRa,kRb,wt_b,wt_a)
1368  do i=is,ie ; do j=js-1,je ; if (g%mask2dCv(i,j) > 0.5) then
1369  do k=1,npv(i,j)
1370  klb = k0b_lv(j)%p(i,k); krb = k0b_rv(j)%p(i,k)
1371  if (deep_wt_lv(j)%p(i,k) >= 1.0) then
1372  tr_flux_conv(i,j,klb) = tr_flux_conv(i,j,klb) - tr_flux_3d(i,j,k)
1373  else
1374  kla = k0a_lv(j)%p(i,k)
1375  wt_b = deep_wt_lv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1376  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))
1377  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))
1378  endif
1379  if (deep_wt_rv(j)%p(i,k) >= 1.0) then
1380  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + tr_flux_3d(i,j,k)
1381  else
1382  kra = k0a_rv(j)%p(i,k)
1383  wt_b = deep_wt_rv(j)%p(i,k) ; wt_a = 1.0 - wt_b
1384  tr_flux_conv(i,j+1,kra) = tr_flux_conv(i,j+1,kra) + &
1385  (wt_a*tr_flux_3d(i,j,k) - tr_adj_vert_r(i,j,k))
1386  tr_flux_conv(i,j+1,krb) = tr_flux_conv(i,j+1,krb) + &
1387  (wt_b*tr_flux_3d(i,j,k) + tr_adj_vert_r(i,j,k))
1388  endif
1389  enddo
1390  endif ; enddo ; enddo
1391 !$OMP parallel do default(none) shared(PEmax_kRho,is,ie,js,je,G,h,Tr,tr_flux_conv,m)
1392  do k=1,pemax_krho ; do j=js,je ; do i=is,ie
1393  if ((g%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 0.0)) then
1394  tr(m)%t(i,j,k) = tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
1395  (h(i,j,k)*g%areaT(i,j))
1396  tr_flux_conv(i,j,k) = 0.0
1397  endif
1398  enddo ; enddo ; enddo
1399 
1400  enddo ! Loop over tracers
1401  enddo ! Loop over iterations
1402 
1403  do j=js,je
1404  deallocate(deep_wt_lu(j)%p) ; deallocate(deep_wt_ru(j)%p)
1405  deallocate(hp_lu(j)%p) ; deallocate(hp_ru(j)%p)
1406  deallocate(k0a_lu(j)%p) ; deallocate(k0a_ru(j)%p)
1407  deallocate(k0b_lu(j)%p) ; deallocate(k0b_ru(j)%p)
1408  enddo
1409 
1410  do j=js-1,je
1411  deallocate(deep_wt_lv(j)%p) ; deallocate(deep_wt_rv(j)%p)
1412  deallocate(hp_lv(j)%p) ; deallocate(hp_rv(j)%p)
1413  deallocate(k0a_lv(j)%p) ; deallocate(k0a_rv(j)%p)
1414  deallocate(k0b_lv(j)%p) ; deallocate(k0b_rv(j)%p)
1415  enddo
1416 
1417 end subroutine tracer_epipycnal_ml_diff
1418 
1419 
1420 !> Initialize lateral tracer diffusion module
1421 subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS)
1422  type(time_type), target, intent(in) :: time !< current model time
1423  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1424  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1425  type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control
1426  type(eos_type), target, intent(in) :: eos !< Equation of state CS
1427  type(diabatic_cs), pointer, intent(in) :: diabatic_csp !< Equation of state CS
1428  type(param_file_type), intent(in) :: param_file !< parameter file
1429  type(tracer_hor_diff_cs), pointer :: cs !< horz diffusion control structure
1430 
1431 ! This include declares and sets the variable "version".
1432 #include "version_variable.h"
1433  character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name.
1434  character(len=256) :: mesg ! Message for error messages.
1435 
1436  if (associated(cs)) then
1437  call mom_error(warning, "tracer_hor_diff_init called with associated control structure.")
1438  return
1439  endif
1440  allocate(cs)
1441 
1442  cs%diag => diag
1443  cs%show_call_tree = calltree_showquery()
1444 
1445  ! Read all relevant parameters and write them to the model log.
1446  call log_version(param_file, mdl, version, "")
1447  call get_param(param_file, mdl, "KHTR", cs%KhTr, &
1448  "The background along-isopycnal tracer diffusivity.", &
1449  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1450  call get_param(param_file, mdl, "KHTR_SLOPE_CFF", cs%KhTr_Slope_Cff, &
1451  "The scaling coefficient for along-isopycnal tracer "//&
1452  "diffusivity using a shear-based (Visbeck-like) "//&
1453  "parameterization. A non-zero value enables this param.", &
1454  units="nondim", default=0.0)
1455  call get_param(param_file, mdl, "KHTR_MIN", cs%KhTr_Min, &
1456  "The minimum along-isopycnal tracer diffusivity.", &
1457  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1458  call get_param(param_file, mdl, "KHTR_MAX", cs%KhTr_Max, &
1459  "The maximum along-isopycnal tracer diffusivity.", &
1460  units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s)
1461  call get_param(param_file, mdl, "KHTR_PASSIVITY_COEFF", cs%KhTr_passivity_coeff, &
1462  "The coefficient that scales deformation radius over "//&
1463  "grid-spacing in passivity, where passivity is the ratio "//&
1464  "between along isopycnal mixing of tracers to thickness mixing. "//&
1465  "A non-zero value enables this parameterization.", &
1466  units="nondim", default=0.0)
1467  call get_param(param_file, mdl, "KHTR_PASSIVITY_MIN", cs%KhTr_passivity_min, &
1468  "The minimum passivity which is the ratio between "//&
1469  "along isopycnal mixing of tracers to thickness mixing.", &
1470  units="nondim", default=0.5)
1471  call get_param(param_file, mdl, "DIFFUSE_ML_TO_INTERIOR", cs%Diffuse_ML_interior, &
1472  "If true, enable epipycnal mixing between the surface "//&
1473  "boundary layer and the interior.", default=.false.)
1474  call get_param(param_file, mdl, "CHECK_DIFFUSIVE_CFL", cs%check_diffusive_CFL, &
1475  "If true, use enough iterations the diffusion to ensure "//&
1476  "that the diffusive equivalent of the CFL limit is not "//&
1477  "violated. If false, always use the greater of 1 or "//&
1478  "MAX_TR_DIFFUSION_CFL iteration.", default=.false.)
1479  call get_param(param_file, mdl, "MAX_TR_DIFFUSION_CFL", cs%max_diff_CFL, &
1480  "If positive, locally limit the along-isopycnal tracer "//&
1481  "diffusivity to keep the diffusive CFL locally at or "//&
1482  "below this value. The number of diffusive iterations "//&
1483  "is often this value or the next greater integer.", &
1484  units="nondim", default=-1.0)
1485  call get_param(param_file, mdl, "RECALC_NEUTRAL_SURF", cs%recalc_neutral_surf, &
1486  "If true, then recalculate the neutral surfaces if the \n"//&
1487  "diffusive CFL is exceeded. If false, assume that the \n"//&
1488  "positions of the surfaces do not change \n", default = .false.)
1489  cs%ML_KhTR_scale = 1.0
1490  if (cs%Diffuse_ML_interior) then
1491  call get_param(param_file, mdl, "ML_KHTR_SCALE", cs%ML_KhTR_scale, &
1492  "With Diffuse_ML_interior, the ratio of the truly "//&
1493  "horizontal diffusivity in the mixed layer to the "//&
1494  "epipycnal diffusivity. The valid range is 0 to 1.", &
1495  units="nondim", default=1.0)
1496  endif
1497 
1498  cs%use_neutral_diffusion = neutral_diffusion_init(time, g, param_file, diag, eos, diabatic_csp, &
1499  cs%neutral_diffusion_CSp )
1500  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1501  "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1502  cs%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(time, g, param_file, diag, diabatic_csp, &
1503  cs%lateral_boundary_diffusion_CSp)
1504  if (cs%use_neutral_diffusion .and. cs%Diffuse_ML_interior) call mom_error(fatal, "MOM_tracer_hor_diff: "// &
1505  "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")
1506 
1507  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1508 
1509  id_clock_diffuse = cpu_clock_id('(Ocean diffuse tracer)', grain=clock_module)
1510  id_clock_epimix = cpu_clock_id('(Ocean epipycnal diffuse tracer)',grain=clock_module)
1511  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1512  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1513 
1514  cs%id_KhTr_u = -1
1515  cs%id_KhTr_v = -1
1516  cs%id_KhTr_h = -1
1517  cs%id_CFL = -1
1518 
1519  cs%id_KhTr_u = register_diag_field('ocean_model', 'KHTR_u', diag%axesCu1, time, &
1520  'Epipycnal tracer diffusivity at zonal faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1521  cs%id_KhTr_v = register_diag_field('ocean_model', 'KHTR_v', diag%axesCv1, time, &
1522  'Epipycnal tracer diffusivity at meridional faces of tracer cell', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1523  cs%id_KhTr_h = register_diag_field('ocean_model', 'KHTR_h', diag%axesT1, time, &
1524  'Epipycnal tracer diffusivity at tracer cell center', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
1525  cmor_field_name='diftrelo', &
1526  cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', &
1527  cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity')
1528 
1529  cs%id_khdt_x = register_diag_field('ocean_model', 'KHDT_x', diag%axesCu1, time, &
1530  'Epipycnal tracer diffusivity operator at zonal faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1531  cs%id_khdt_y = register_diag_field('ocean_model', 'KHDT_y', diag%axesCv1, time, &
1532  'Epipycnal tracer diffusivity operator at meridional faces of tracer cell', 'm2', conversion=us%L_to_m**2)
1533  if (cs%check_diffusive_CFL) then
1534  cs%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, time,&
1535  'Grid CFL number for lateral/neutral tracer diffusion', 'nondim')
1536  endif
1537 
1538 
1539 end subroutine tracer_hor_diff_init
1540 
1541 subroutine tracer_hor_diff_end(CS)
1542  type(tracer_hor_diff_cs), pointer :: cs !< module control structure
1543 
1544  call neutral_diffusion_end(cs%neutral_diffusion_CSp)
1545  if (associated(cs)) deallocate(cs)
1546 
1547 end subroutine tracer_hor_diff_end
1548 
1549 
1550 !> \namespace mom_tracer_hor_diff
1551 !!
1552 !! \section section_intro Introduction to the module
1553 !!
1554 !! This module contains subroutines that handle horizontal
1555 !! diffusion (i.e., isoneutral or along layer) of tracers.
1556 !!
1557 !! Each of the tracers are subject to Fickian along-coordinate
1558 !! diffusion if Khtr is defined and positive. The tracer diffusion
1559 !! can use a suitable number of iterations to guarantee stability
1560 !! with an arbitrarily large time step.
1561 
1562 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:38
mom_lateral_boundary_diffusion::lateral_boundary_diffusion
subroutine, public lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries....
Definition: MOM_lateral_boundary_diffusion.F90:123
mom_tracer_hor_diff::tracer_hordiff
subroutine, public tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y)
Compute along-coordinate diffusion of all tracers using the diffusivity in CSKhTr,...
Definition: MOM_tracer_hor_diff.F90:107
mom_tracer_hor_diff::tracer_epipycnal_ml_diff
subroutine tracer_epipycnal_ml_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, GV, US, CS, tv, num_itts)
This subroutine does epipycnal diffusion of all tracers between the mixed and buffer layers and the i...
Definition: MOM_tracer_hor_diff.F90:579
mom_neutral_diffusion::neutral_diffusion_init
logical function, public neutral_diffusion_init(Time, G, param_file, diag, EOS, diabatic_CSp, CS)
Read parameters and allocate control structure for neutral_diffusion module.
Definition: MOM_neutral_diffusion.F90:109
mom_tracer_hor_diff::id_clock_diffuse
integer id_clock_diffuse
CPU time clocks.
Definition: MOM_tracer_hor_diff.F90:97
mom_tracer_hor_diff::id_clock_epimix
integer id_clock_epimix
CPU time clocks.
Definition: MOM_tracer_hor_diff.F90:97
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 ocntrol 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:78
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_neutral_diffusion::neutral_diffusion_end
subroutine, public neutral_diffusion_end(CS)
Deallocates neutral_diffusion control structure.
Definition: MOM_neutral_diffusion.F90:2807
mom_tracer_hor_diff::id_clock_sync
integer id_clock_sync
CPU time clocks.
Definition: MOM_tracer_hor_diff.F90:97
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
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_meke_types
Definition: MOM_MEKE_types.F90:1
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:54
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_lateral_boundary_diffusion::lateral_boundary_diffusion_init
logical function, public lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS)
Initialization routine that reads runtime parameters and sets up pointers to other control structures...
Definition: MOM_lateral_boundary_diffusion.F90:62
mom_domains::do_group_pass
subroutine, public do_group_pass(group, MOM_dom, clock)
do_group_pass carries out a group halo update.
Definition: MOM_domains.F90:1090
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:26
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_neutral_diffusion::neutral_diffusion
subroutine, public neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)
Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.
Definition: MOM_neutral_diffusion.F90:493
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_neutral_diffusion::neutral_diffusion_calc_coeffs
subroutine, public neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS)
Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate spac...
Definition: MOM_neutral_diffusion.F90:264
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_error_handler::calltree_leave
subroutine, public calltree_leave(mesg)
Writes a message about leaving a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:151
mom_tracer_hor_diff::tracer_hor_diff_init
subroutine, public tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS)
Initialize lateral tracer diffusion module.
Definition: MOM_tracer_hor_diff.F90:1422
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_error_handler::mom_set_verbosity
subroutine, public mom_set_verbosity(verb)
This subroutine sets the level of verbosity filtering MOM error messages.
Definition: MOM_error_handler.F90:97
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_tracer_hor_diff::id_clock_pass
integer id_clock_pass
CPU time clocks.
Definition: MOM_tracer_hor_diff.F90:97
mom_tracer_hor_diff::tracer_hor_diff_end
subroutine, public tracer_hor_diff_end(CS)
Definition: MOM_tracer_hor_diff.F90:1542
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_error_handler::calltree_showquery
logical function, public calltree_showquery()
Returns True, if the verbosity>=6 indicating to show the call tree.
Definition: MOM_error_handler.F90:123
mom_tracer_registry::mom_tracer_chksum
subroutine, public mom_tracer_chksum(mesg, Tr, ntr, G)
This subroutine writes out chksums for tracers.
Definition: MOM_tracer_registry.F90:804
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_error_handler::calltree_waypoint
subroutine, public calltree_waypoint(mesg, n)
Writes a message about reaching a milestone if call tree reporting is active.
Definition: MOM_error_handler.F90:161
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:79
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_error_handler::calltree_enter
subroutine, public calltree_enter(mesg, n)
Writes a message about entering a subroutine if call tree reporting is active.
Definition: MOM_error_handler.F90:130
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:239
mom_neutral_diffusion::neutral_diffusion_cs
The control structure for the MOM_neutral_diffusion module.
Definition: MOM_neutral_diffusion.F90:40
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60