MOM_tracer_hor_diff.F90

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