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