MOM_tracer_advect.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!> This module contains the subroutines that advect tracers along coordinate surfaces.
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9use mom_cpu_clock, only : clock_module, clock_routine
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, pass_var
14use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
15use mom_file_parser, only : get_param, log_version, param_file_type
16use mom_grid, only : ocean_grid_type
17use mom_open_boundary, only : ocean_obc_type, obc_none, obc_direction_e
18use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s
20use mom_tracer_registry, only : tracer_registry_type, tracer_type
23use mom_tracer_advect_schemes, only : advect_plm, advect_ppmh3, advect_ppm
24use mom_tracer_advect_schemes, only : set_tracer_advect_scheme, traceradvectionschemedoc
25implicit none ; private
26
27#include <MOM_memory.h>
28
29public advect_tracer
32
33!> Control structure for this module
34type, public :: tracer_advect_cs ; private
35 real :: dt !< The baroclinic dynamics time step [T ~> s].
36 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
37 !< timing of diagnostic output.
38 logical :: debug !< If true, write verbose checksums for debugging purposes.
39 logical :: usehuynhstencilbug = .false. !< If true, use the incorrect stencil width.
40 !! This is provided for compatibility with legacy simuations.
41 type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structure used for group passes
42 integer :: default_advect_scheme = -1 !< Determines which reconstruction to use
43end type tracer_advect_cs
44
45!>@{ CPU time clocks
46integer :: id_clock_advect
47integer :: id_clock_pass
48integer :: id_clock_sync
49!>@}
50
51contains
52
53!> This routine time steps the tracer concentration using a
54!! monotonic, conservative, weakly diffusive scheme.
55subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first_in, &
56 vol_prev, max_iter_in, update_vol_prev, uhr_out, vhr_out, &
57 flux_type)
58 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
59 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
60 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
61 intent(in) :: h_end !< Layer thickness after advection [H ~> m or kg m-2]
62 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
63 intent(in) :: uhtr !< Accumulated volume or mass flux through the
64 !! zonal faces [H L2 ~> m3 or kg]
65 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
66 intent(in) :: vhtr !< Accumulated volume or mass flux through the
67 !! meridional faces [H L2 ~> m3 or kg]
68 type(ocean_obc_type), pointer :: obc !< specifies whether, where, and what OBCs are used
69 real, intent(in) :: dt !< time increment [T ~> s]
70 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
71 type(tracer_advect_cs), pointer :: cs !< control structure for module
72 type(tracer_registry_type), pointer :: reg !< pointer to tracer registry
73 logical, optional, intent(in) :: x_first_in !< If present, indicate whether to update
74 !! first in the x- or y-direction.
75 ! The next four optional arguments are only used in offline tracer mode.
76 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
77 optional, intent(inout) :: vol_prev !< Cell volume before advection [H L2 ~> m3 or kg].
78 !! If update_vol_prev is true, the returned value is
79 !! the cell volume after the transport that was done
80 !! by this call, and if all the transport could be
81 !! accommodated it should be close to h_end*G%areaT.
82 integer, optional, intent(in) :: max_iter_in !< The maximum number of iterations
83 logical, optional, intent(in) :: update_vol_prev !< If present and true, update vol_prev to
84 !! return its value after the tracer have been updated.
85 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
86 optional, intent(out) :: uhr_out !< Remaining accumulated volume or mass fluxes
87 !! through the zonal faces [H L2 ~> m3 or kg]
88 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
89 optional, intent(out) :: vhr_out !< Remaining accumulated volume or mass fluxes
90 !! through the meridional faces [H L2 ~> m3 or kg]
91 ! The next optional argument is for diagnosing resolved vs parameterized tracer flux and control
92 ! which diagnostics are written. The tracers are only updated if flux_type = 0 (the default). Otherwise
93 ! the routines are dry run to collect diagnostics.
94 integer, optional, intent(in) :: flux_type !< Indicates whether uhtr, vhtr are the flux due to
95 !! the residual (= 0), resolved (= 1), or parameterized (= 2)
96 !! flow
97
98 ! local variables
99 integer :: flux_type_ctrl !< To allow setting a default value for flux_type
100 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
101 hprev ! cell volume at the end of previous tracer change [H L2 ~> m3 or kg]
102 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
103 uhr ! The remaining zonal thickness flux [H L2 ~> m3 or kg]
104 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
105 vhr ! The remaining meridional thickness fluxes [H L2 ~> m3 or kg]
106 real :: uh_neglect(szib_(g),szj_(g)) ! uh_neglect and vh_neglect are the
107 real :: vh_neglect(szi_(g),szjb_(g)) ! magnitude of remaining transports that
108 ! can be simply discarded [H L2 ~> m3 or kg].
109
110 real :: landvolfill ! An arbitrary? nonzero cell volume [H L2 ~> m3 or kg].
111 logical :: use_ppm_stencil ! If true, use the correct PPM stencil width.
112 real :: idt ! 1/dt [T-1 ~> s-1].
113 logical :: domore_u(szj_(g),szk_(gv)) ! domore_u and domore_v indicate whether there is more
114 logical :: domore_v(szjb_(g),szk_(gv)) ! advection to be done in the corresponding row or column.
115 logical :: x_first ! If true, advect in the x-direction first.
116 logical :: advect_this_tracer(reg%ntr) ! If true, advect the mth tracer. Diagnostics of advection due to the
117 ! resolved and parameterized flow are collected by re-running the advection
118 ! routines with different advecting fluxes without updating the tracer.
119 ! This can be expensive if there are lots of tracers and only a few you
120 ! want diagnostics about. We therefore only calculate advection on the
121 ! tracers for which there are active resolved/parameterized diagnostics.
122 integer :: max_iter ! maximum number of iterations in each layer
123 integer :: domore_k(szk_(gv))
124 integer :: stencil ! stencil of the advection scheme
125 integer :: nsten_halo ! number of stencils that fit in the halos
126 integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any
127 integer :: isv, iev, jsv, jev ! The valid range of the indices.
128 integer :: isdb, iedb, jsdb, jedb
129 integer :: stencil_local ! Stencil for the local adection scheme
130 integer :: local_advect_scheme(reg%ntr) ! contains the list of the advection for each tracer
131
132 domore_u(:,:) = .false.
133 domore_v(:,:) = .false.
134 advect_this_tracer(:) = .false.
135 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
136 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
137 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
138 landvolfill = 1.0e-20 ! This is arbitrary, but must be positive.
139 stencil = 2 ! The scheme's stencil; 2 for PLM
140
141 ntr = reg%ntr
142 idt = 1.0 / dt
143
144 if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_advect: "// &
145 "tracer_advect_init must be called before advect_tracer.")
146 if (.not. associated(reg)) call mom_error(fatal, "MOM_tracer_advect: "// &
147 "register_tracer must be called before advect_tracer.")
148 if (reg%ntr==0) return
149 call cpu_clock_begin(id_clock_advect)
150 x_first = (mod(g%first_direction,2) == 0)
151
152 ! Choose the maximum stencil from all the local advection scheme
153 do m = 1,ntr
154
155 local_advect_scheme(m) = reg%Tr(m)%advect_scheme
156 if (local_advect_scheme(m) < 0) local_advect_scheme(m) = cs%default_advect_scheme
157
158 if (local_advect_scheme(m) == advect_plm) then
159 stencil_local = 2
160 elseif (local_advect_scheme(m) == advect_ppm) then
161 stencil_local = 3
162 elseif (local_advect_scheme(m) == advect_ppmh3) then
163 if (cs%useHuynhStencilBug) then
164 stencil_local = 2
165 else
166 stencil_local = 3
167 endif
168 endif
169 stencil = max(stencil, stencil_local)
170 enddo
171
172 if (min(is-isd,ied-ie,js-jsd,jed-je) < stencil) then
173 call mom_error(fatal, "MOM_tracer_advect: "//&
174 "stencil is wider than the halo.")
175 endif
176
177 max_iter = 2*int(ceiling(dt/cs%dt)) + 1
178
179 if (present(max_iter_in)) max_iter = max_iter_in
180 if (present(x_first_in)) x_first = x_first_in
181
182 flux_type_ctrl = 0
183 if (present(flux_type)) flux_type_ctrl = flux_type ! default to residual flow
184
185 call cpu_clock_begin(id_clock_pass)
186 call create_group_pass(cs%pass_uhr_vhr_t_hprev, uhr, vhr, g%Domain)
187 call create_group_pass(cs%pass_uhr_vhr_t_hprev, hprev, g%Domain)
188 do m=1,ntr
189 call create_group_pass(cs%pass_uhr_vhr_t_hprev, reg%Tr(m)%t, g%Domain)
190 enddo
191 call cpu_clock_end(id_clock_pass)
192
193 !$OMP parallel default(shared)
194
195 ! This initializes the halos of uhr and vhr because pass_vector might do
196 ! calculations on them, even though they are never used.
197 !$OMP do
198 do k=1,nz
199 do j=jsd,jed ; do i=isdb,iedb ; uhr(i,j,k) = 0.0 ; enddo ; enddo
200 do j=jsdb,jedb ; do i=isd,ied ; vhr(i,j,k) = 0.0 ; enddo ; enddo
201 do j=jsd,jed ; do i=isd,ied ; hprev(i,j,k) = 0.0 ; enddo ; enddo
202 domore_k(k)=1
203 ! Put the remaining (total) thickness fluxes into uhr and vhr.
204 do j=js,je ; do i=is-1,ie ; uhr(i,j,k) = uhtr(i,j,k) ; enddo ; enddo
205 do j=js-1,je ; do i=is,ie ; vhr(i,j,k) = vhtr(i,j,k) ; enddo ; enddo
206 if (.not. present(vol_prev)) then
207 ! This loop reconstructs the thickness field the last time that the
208 ! tracers were updated, probably just after the diabatic forcing. A useful
209 ! diagnostic could be to compare this reconstruction with that older value.
210 do j=js,je ; do i=is,ie
211 hprev(i,j,k) = max(0.0, g%areaT(i,j)*h_end(i,j,k) + &
212 ((uhr(i,j,k) - uhr(i-1,j,k)) + (vhr(i,j,k) - vhr(i,j-1,k))))
213 ! In the case that the layer is now dramatically thinner than it was previously,
214 ! add a bit of mass to avoid truncation errors. This will lead to
215 ! non-conservation of tracers
216 hprev(i,j,k) = hprev(i,j,k) + &
217 max(0.0, 1.0e-13*hprev(i,j,k) - g%areaT(i,j)*h_end(i,j,k))
218 enddo ; enddo
219 else
220 do j=js,je ; do i=is,ie
221 hprev(i,j,k) = vol_prev(i,j,k)
222 enddo ; enddo
223 endif
224 enddo
225
226
227 !$OMP do
228 do j=jsd,jed ; do i=isd,ied-1
229 uh_neglect(i,j) = gv%H_subroundoff * min(g%areaT(i,j), g%areaT(i+1,j))
230 enddo ; enddo
231 !$OMP do
232 do j=jsd,jed-1 ; do i=isd,ied
233 vh_neglect(i,j) = gv%H_subroundoff * min(g%areaT(i,j), g%areaT(i,j+1))
234 enddo ; enddo
235
236 ! initialize diagnostic fluxes and tendencies and determine which tracers to advect
237 if (flux_type_ctrl == 0) then ! Flux is residual
238 !$OMP do
239 do m=1,ntr
240 advect_this_tracer(m) = .true. ! Advect all the tracers regardless of diagnostic output
241 if (associated(reg%Tr(m)%ad_x)) reg%Tr(m)%ad_x(:,:,:) = 0.0
242 if (associated(reg%Tr(m)%ad_y)) reg%Tr(m)%ad_y(:,:,:) = 0.0
243 if (associated(reg%Tr(m)%advection_xy)) reg%Tr(m)%advection_xy(:,:,:) = 0.0
244 if (associated(reg%Tr(m)%ad2d_x)) reg%Tr(m)%ad2d_x(:,:) = 0.0
245 if (associated(reg%Tr(m)%ad2d_y)) reg%Tr(m)%ad2d_y(:,:) = 0.0
246 enddo
247 elseif (flux_type_ctrl == 1) then ! Flux is resolved
248 do m=1,ntr
249 if (associated(reg%Tr(m)%ad_x_resolved)) then
250 reg%Tr(m)%ad_x_resolved(:,:,:) = 0.0
251 advect_this_tracer(m) = .true. ! advect this tracer
252 endif
253 if (associated(reg%Tr(m)%ad_y_resolved)) then
254 reg%Tr(m)%ad_y_resolved(:,:,:) = 0.0
255 advect_this_tracer(m) = .true. ! advect this tracer
256 endif
257 enddo
258 elseif (flux_type_ctrl == 2) then ! Flux is parameterized
259 do m=1,ntr
260 if (associated(reg%Tr(m)%ad_x_param)) then
261 reg%Tr(m)%ad_x_param(:,:,:) = 0.0
262 advect_this_tracer(m) = .true. ! advect this tracer
263 endif
264 if (associated(reg%Tr(m)%ad_y_param)) then
265 reg%Tr(m)%ad_y_param(:,:,:) = 0.0
266 advect_this_tracer(m) = .true. ! advect this tracer
267 endif
268 enddo
269 else
270 call mom_error(fatal, &
271 "Inconsistent flux type in advect_tracer. Must be of 0 (residual), 1 (resolved), or 2 (parameterized)")
272 endif ! flux_type_ctrl
273 !$OMP end parallel
274
275 isv = is ; iev = ie ; jsv = js ; jev = je
276 nsten_halo = min(is - isd, ied - ie, js - jsd, jed - je) / stencil
277
278 do itt=1,max_iter
279
280 if (isv > is-stencil) then
281 call do_group_pass(cs%pass_uhr_vhr_t_hprev, g%Domain, clock=id_clock_pass)
282
283 isv = is - nsten_halo * stencil ; jsv = js - nsten_halo * stencil
284 iev = ie + nsten_halo * stencil ; jev = je + nsten_halo * stencil
285 ! Reevaluate domore_u & domore_v unless the valid range is the same size as
286 ! before. Also, do this if there is Strang splitting.
287 if ((nsten_halo > 1) .or. (itt==1)) then
288 !$OMP parallel do default(shared)
289 do k=1,nz ; if (domore_k(k) > 0) then
290 do j=jsv,jev ; if (.not.domore_u(j,k)) then
291 do i=isv+stencil-1,iev-stencil ; if (uhr(i,j,k) /= 0.0) then
292 domore_u(j,k) = .true. ; exit
293 endif ; enddo ! i-loop
294 endif ; enddo
295 do j=jsv+stencil-1,jev-stencil ; if (.not.domore_v(j,k)) then
296 do i=isv+stencil,iev-stencil ; if (vhr(i,j,k) /= 0.0) then
297 domore_v(j,k) = .true. ; exit
298 endif ; enddo ! i-loop
299 endif ; enddo
300
301 ! At this point, domore_k is global. Change it so that it indicates
302 ! whether any work is needed on a layer on this processor.
303 domore_k(k) = 0
304 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
305 do j=jsv+stencil-1,jev-stencil ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
306
307 endif ; enddo ! k-loop
308 endif
309 endif
310
311 ! Set the range of valid points after this iteration.
312 isv = isv + stencil ; iev = iev - stencil
313 jsv = jsv + stencil ; jev = jev - stencil
314
315 ! To ensure positive definiteness of the thickness at each iteration, the
316 ! mass fluxes out of each layer are checked each step, and limited to keep
317 ! the thicknesses positive. This means that several iterations may be required
318 ! for all the transport to happen. The sum over domore_k keeps the processors
319 ! synchronized. This may not be very efficient, but it should be reliable.
320
321 !$OMP parallel default(shared)
322
323 if (x_first) then
324
325 !$OMP do ordered
326 do k=1,nz ; if (domore_k(k) > 0) then
327 ! First, advect zonally.
328 call advect_x(reg%Tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
329 isv, iev, jsv-stencil, jev+stencil, k, g, gv, us, &
330 flux_type_ctrl, advect_this_tracer, local_advect_scheme)
331 endif ; enddo
332
333 !$OMP do ordered
334 do k=1,nz ; if (domore_k(k) > 0) then
335 ! Next, advect meridionally.
336 call advect_y(reg%Tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
337 isv, iev, jsv, jev, k, g, gv, us, flux_type_ctrl, advect_this_tracer, &
338 local_advect_scheme)
339
340 ! Update domore_k(k) for the next iteration
341 domore_k(k) = 0
342 do j=jsv-stencil,jev+stencil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
343 do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
344
345 endif ; enddo
346
347 else
348
349 !$OMP do ordered
350 do k=1,nz ; if (domore_k(k) > 0) then
351 ! First, advect meridionally.
352 call advect_y(reg%Tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
353 isv-stencil, iev+stencil, jsv, jev, k, g, gv, us, &
354 flux_type_ctrl, advect_this_tracer, local_advect_scheme)
355 endif ; enddo
356
357 !$OMP do ordered
358 do k=1,nz ; if (domore_k(k) > 0) then
359 ! Next, advect zonally.
360 call advect_x(reg%Tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
361 isv, iev, jsv, jev, k, g, gv, us, flux_type_ctrl, advect_this_tracer, &
362 local_advect_scheme)
363
364 ! Update domore_k(k) for the next iteration
365 domore_k(k) = 0
366 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
367 do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
368 endif ; enddo
369
370 endif ! x_first
371
372 !$OMP end parallel
373
374 ! If the advection just isn't finishing after max_iter, move on.
375 if (itt >= max_iter) then
376 exit
377 endif
378
379 ! Exit if there are no layers that need more iterations.
380 if (isv > is-stencil) then
381 do_any = 0
382 call cpu_clock_begin(id_clock_sync)
383 call sum_across_pes(domore_k(:), nz)
384 call cpu_clock_end(id_clock_sync)
385 do k=1,nz ; do_any = do_any + domore_k(k) ; enddo
386 if (do_any == 0) then
387 exit
388 endif
389
390 endif
391
392 enddo ! Iterations loop
393
394 if (present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:)
395 if (present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:)
396 if (present(vol_prev) .and. present(update_vol_prev)) then
397 if (update_vol_prev) vol_prev(:,:,:) = hprev(:,:,:)
398 endif
399
400 call cpu_clock_end(id_clock_advect)
401
402end subroutine advect_tracer
403
404
405!> This subroutine does 1-d flux-form advection in the zonal direction using
406!! a monotonic piecewise linear scheme.
407subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
408 is, ie, js, je, k, G, GV, US, flux_type, advect_this_tracer, &
409 advect_schemes)
410 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
411 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
412 integer, intent(in) :: ntr !< The number of tracers
413 type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on
414 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous
415 !! tracer change [H L2 ~> m3 or kg]
416 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhr !< accumulated volume/mass flux through
417 !! the zonal face [H L2 ~> m3 or kg]
418 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: uh_neglect !< A tiny zonal mass flux that can
419 !! be neglected [H L2 ~> m3 or kg]
420 type(ocean_obc_type), pointer :: OBC !< specifies whether, where, and what OBCs are used
421 logical, dimension(SZJ_(G),SZK_(GV)), intent(inout) :: domore_u !< If true, there is more advection to be
422 !! done in this u-row
423 real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1]
424 integer, intent(in) :: is !< The starting tracer i-index to work on
425 integer, intent(in) :: ie !< The ending tracer i-index to work on
426 integer, intent(in) :: js !< The starting tracer j-index to work on
427 integer, intent(in) :: je !< The ending tracer j-index to work on
428 integer, intent(in) :: k !< The k-level to work on
429 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
430 integer, intent(in) :: flux_type !< Indicates whether uhtr, vhtr are the flux
431 !! due to the residual (= 0), resolved (= 1),
432 !! or parameterized (= 2) flow
433 logical, dimension(ntr), intent(in) :: advect_this_tracer !< If true, advect this tracer
434 integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use
435
436 real, dimension(SZI_(G),ntr) :: &
437 slope_x ! The concentration slope per grid point [conc].
438 real, dimension(SZIB_(G),SZJ_(G),ntr) :: &
439 flux_x ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc].
440 real, dimension(SZI_(G),ntr) :: &
441 T_tmp ! The copy of the tracer concentration at constant i,k [conc].
442
443 real :: hup, hlos ! hup is the upwind volume, hlos is the
444 ! part of that volume that might be lost
445 ! due to advection out the other side of
446 ! the grid box, both in [H L2 ~> m3 or kg].
447 real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the
448 ! current iteration [H L2 ~> m3 or kg].
449 real, dimension(SZIB_(G)) :: &
450 hlst, & ! Work variable [H L2 ~> m3 or kg].
451 Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1].
452 CFL ! The absolute value of the advective upwind-cell CFL number [nondim].
453 real :: min_h ! The minimum thickness that can be realized during
454 ! any of the passes [H ~> m or kg m-2].
455 real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2].
456 real :: h_neglect ! A thickness that is so small it is usually lost
457 ! in roundoff and can be neglected [H ~> m or kg m-2].
458 real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc]
459 real :: dMx ! Difference between the maximum of the surrounding cell concentrations and
460 ! the value in the cell whose reconstruction is being found [conc]
461 real :: dMn ! Difference between the tracer concentration in the cell whose reconstruction
462 ! is being found and the minimum of the surrounding values [conc]
463 real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc]
464 real :: dA ! Difference between the reconstruction tracer edge values [conc]
465 real :: mA ! Average of the reconstruction tracer edge values [conc]
466 real :: a6 ! Curvature of the reconstruction tracer values [conc]
467 logical :: do_i(SZI_(G),SZJ_(G)) ! If true, work on given points.
468 logical :: usePLMslope
469 integer :: i, j, m, n, i_up, stencil, ntr_id
470 type(obc_segment_type), pointer :: segment=>null()
471 logical, dimension(SZJ_(G),SZK_(GV)) :: domore_u_initial
472
473 ! keep a local copy of the initial values of domore_u, which is to be used when computing ad2d_x
474 ! diagnostic at the end of this subroutine.
475 domore_u_initial = domore_u
476
477 useplmslope = .false.
478 ! stencil for calculating slope values
479 stencil = 1
480 do m = 1,ntr
481 if ((advect_schemes(m) == advect_plm) .or. (advect_schemes(m) == advect_ppm)) &
482 useplmslope = .true.
483 if (advect_schemes(m) == advect_ppm) stencil = 2
484 enddo
485
486 min_h = 0.1*gv%Angstrom_H
487 tiny_h = tiny(min_h)
488 h_neglect = gv%H_subroundoff
489
490 do i=is-1,ie ; cfl(i) = 0.0 ; enddo
491
492 do j=js,je ; if (domore_u(j,k)) then
493 domore_u(j,k) = .false.
494
495 ! Calculate the i-direction profiles (slopes) of each tracer that is being advected.
496 if (useplmslope) then
497 do m=1,ntr
498 if (advect_this_tracer(m)) then
499 do i=is-stencil,ie+stencil
500 !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < &
501 ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))) then
502 ! maxslope = 4.0*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k))
503 !else
504 ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))
505 !endif
506 !if ((Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) * (Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k)) < 0.0) then
507 ! slope_x(i,m) = 0.0
508 !elseif (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))<ABS(maxslope)) then
509 ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * &
510 ! 0.5*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))
511 !else
512 ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * 0.5*maxslope
513 !endif
514 tp = tr(m)%t(i+1,j,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i-1,j,k)
515 dmx = max( tp, tc, tm ) - tc
516 dmn= tc - min( tp, tc, tm )
517 slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
518 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
519 enddo
520 endif ! advect_this_tracer
521 enddo
522 endif ! usePLMslope
523
524 ! make a copy of the tracers in case values need to be overridden for OBCs
525 do m = 1,ntr
526 if (advect_this_tracer(m)) then
527 do i=g%isd,g%ied
528 t_tmp(i,m) = tr(m)%t(i,j,k)
529 enddo
530 endif ! advect_this_tracer
531 enddo
532 ! loop through open boundaries and recalculate flux terms
533 if (associated(obc)) then ; if (obc%OBC_pe) then
534 do n=1,obc%number_of_segments
535 segment=>obc%segment(n)
536 if (.not. associated(segment%tr_Reg)) cycle
537 if (segment%is_E_or_W) then
538 if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
539 i = segment%HI%IsdB
540 do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values
541 ntr_id = segment%tr_reg%Tr(m)%ntr_index
542 if (advect_this_tracer(ntr_id)) then
543 if (segment%direction == obc_direction_w) then
544 t_tmp(i,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
545 else
546 t_tmp(i+1,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
547 endif
548 endif ! advect_this_tracer
549 enddo
550 do m = 1,ntr ! Apply update tracer values for slope calculation
551 if (advect_this_tracer(m)) then
552 do i=segment%HI%IsdB-1,segment%HI%IsdB+1
553 tp = t_tmp(i+1,m) ; tc = t_tmp(i,m) ; tm = t_tmp(i-1,m)
554 dmx = max( tp, tc, tm ) - tc
555 dmn= tc - min( tp, tc, tm )
556 slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
557 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
558 enddo
559 endif ! advect_this_tracer
560 enddo
561
562 endif
563 endif
564 enddo
565 endif ; endif
566
567
568 ! Calculate the i-direction fluxes of each tracer, using as much
569 ! the minimum of the remaining mass flux (uhr) and the half the mass
570 ! in the cell plus whatever part of its half of the mass flux that
571 ! the flux through the other side does not require.
572 do i=is-1,ie
573 if ((uhr(i,j,k) == 0.0) .or. &
574 ((uhr(i,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. &
575 ((uhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
576 uhh(i) = 0.0
577 cfl(i) = 0.0
578 elseif (uhr(i,j,k) < 0.0) then
579 hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
580 hlos = max(0.0, uhr(i+1,j,k))
581 if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
582 ((0.5*hup + uhr(i,j,k)) < 0.0)) then
583 uhh(i) = min(-0.5*hup, -hup+hlos, 0.0)
584 domore_u(j,k) = .true.
585 else
586 uhh(i) = uhr(i,j,k)
587 endif
588 cfl(i) = - uhh(i) / (hprev(i+1,j,k)) ! CFL is positive
589 else
590 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
591 hlos = max(0.0, -uhr(i-1,j,k))
592 if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
593 ((0.5*hup - uhr(i,j,k)) < 0.0)) then
594 uhh(i) = max(0.5*hup, hup-hlos, 0.0)
595 domore_u(j,k) = .true.
596 else
597 uhh(i) = uhr(i,j,k)
598 endif
599 cfl(i) = uhh(i) / (hprev(i,j,k)) ! CFL is positive
600 endif
601 enddo
602
603 do m=1,ntr
604 if (advect_this_tracer(m)) then
605 if ((advect_schemes(m) == advect_ppm) .or. (advect_schemes(m) == advect_ppmh3)) then
606 do i=is-1,ie
607 ! centre cell depending on upstream direction
608 if (uhh(i) >= 0.0) then
609 i_up = i
610 else
611 i_up = i+1
612 endif
613
614 ! Implementation of PPM-H3
615 tp = t_tmp(i_up+1,m) ; tc = t_tmp(i_up,m) ; tm = t_tmp(i_up-1,m)
616
617 if (advect_schemes(m) == advect_ppmh3) then
618 al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
619 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
620 ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
621 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
622 else
623 al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
624 ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
625 endif
626
627 da = ar - al ; ma = 0.5*( ar + al )
628 if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.) then
629 al = tc ; ar = tc ! PCM for local extrema and boundary cells
630 elseif ( da*(tc-ma) > (da*da)/6. ) then
631 al = (3.*tc) - 2.*ar
632 elseif ( da*(tc-ma) < - (da*da)/6. ) then
633 ar = (3.*tc) - 2.*al
634 endif
635
636 a6 = 6.*tc - 3. * (ar + al) ! Curvature
637
638 if (uhh(i) >= 0.0) then
639 flux_x(i,j,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
640 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
641 else
642 flux_x(i,j,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
643 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
644 endif
645 enddo
646 else ! PLM
647 do i=is-1,ie
648 if (uhh(i) >= 0.0) then
649 ! Indirect implementation of PLM
650 !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m)
651 !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
652 !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) )
653 ! Alternative implementation of PLM
654 tc = t_tmp(i,m)
655 flux_x(i,j,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
656 else
657 ! Indirect implementation of PLM
658 !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
659 !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m)
660 !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) )
661 ! Alternative implementation of PLM
662 tc = t_tmp(i+1,m)
663 flux_x(i,j,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
664 endif
665 enddo
666 endif ! usePPM
667 endif ! advect_this_tracer
668 enddo
669
670 if (associated(obc)) then ; if (obc%OBC_pe) then
671 if (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally) then
672 do n=1,obc%number_of_segments
673 segment=>obc%segment(n)
674 if (.not. associated(segment%tr_Reg)) cycle
675 if (segment%is_E_or_W) then
676 if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
677 i = segment%HI%IsdB
678 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
679 ! Now changing to simply fixed inflows.
680 if ((uhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_w) .or. &
681 (uhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_e)) then
682 uhh(i) = uhr(i,j,k)
683 ! should the reservoir evolve for this case Kate ?? - Nope
684 do m=1,segment%tr_Reg%ntseg
685 ntr_id = segment%tr_reg%Tr(m)%ntr_index
686 if (advect_this_tracer(ntr_id)) then
687 flux_x(i,j,ntr_id) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
688 endif ! advect_this_tracer
689 enddo
690 endif
691 endif
692 endif
693 enddo
694 endif
695
696 if (obc%open_u_BCs_exist_globally) then
697 do n=1,obc%number_of_segments
698 segment=>obc%segment(n)
699 i = segment%HI%IsdB
700 if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
701 if (segment%specified) cycle
702 if (.not. associated(segment%tr_Reg)) cycle
703
704 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
705 if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
706 (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5)) then
707 uhh(i) = uhr(i,j,k)
708 do m=1,segment%tr_Reg%ntseg
709 ntr_id = segment%tr_reg%Tr(m)%ntr_index
710 if (advect_this_tracer(ntr_id)) then
711 flux_x(i,j,ntr_id) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
712 endif ! advect_this_tracer
713 enddo
714 endif
715 endif
716 enddo
717 endif
718 endif ; endif
719
720 ! Calculate new tracer concentration in each cell after accounting
721 ! for the i-direction fluxes.
722 do i=is-1,ie
723 uhr(i,j,k) = uhr(i,j,k) - uhh(i)
724 if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
725 enddo
726 do i=is,ie
727 if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0)) then
728 do_i(i,j) = .true.
729 hlst(i) = hprev(i,j,k)
730 hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
731 if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
732 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
733 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
734 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
735 else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
736 else
737 do_i(i,j) = .false.
738 endif
739 enddo
740
741 ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
742 if (associated(obc)) then
743 if ((.not.obc%exterior_OBC_bug) .and. (obc%OBC_pe) .and. &
744 (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally)) then
745 ! OBC_DIRECTION_E / OBC_DIRECTION_W on the west / east edge
746 do i=is,ie ; if ((obc%segnum_u(i-1,j) > 0) .or. (obc%segnum_u(i,j) < 0)) &
747 do_i(i,j) = .false.
748 enddo
749 endif
750 endif
751
752 ! update tracer concentration from i-flux and save some diagnostics
753 do m=1,ntr
754 if (advect_this_tracer(m)) then
755
756 ! update tracer
757 if (flux_type == 0) then ! Only update tracer if using residual flux
758 do i=is,ie
759 if (do_i(i,j)) then
760 if (ihnew(i) > 0.0) then
761 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
762 (flux_x(i,j,m) - flux_x(i-1,j,m))) * ihnew(i)
763 endif
764 endif
765 enddo
766 endif ! flux_type == 0
767
768 ! diagnostics
769 if (flux_type == 0) then
770 if (associated(tr(m)%ad_x)) then ; do i=is-1,ie
771 tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,j,m)*idt
772 enddo ; endif
773
774 ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic).
775 ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
776 if (associated(tr(m)%advection_xy)) then
777 do i=is,ie ; if (do_i(i,j)) then
778 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_x(i,j,m) - flux_x(i-1,j,m)) * &
779 idt * g%IareaT(i,j)
780 endif ; enddo
781 endif
782 elseif (flux_type == 1) then
783 if (associated(tr(m)%ad_x_resolved)) then ; do i=is-1,ie
784 tr(m)%ad_x_resolved(i,j,k) = tr(m)%ad_x_resolved(i,j,k) + flux_x(i,j,m)*idt
785 enddo ; endif
786 elseif (flux_type == 2) then
787 if (associated(tr(m)%ad_x_param)) then ; do i=is-1,ie
788 tr(m)%ad_x_param(i,j,k) = tr(m)%ad_x_param(i,j,k) + flux_x(i,j,m)*idt
789 enddo ; endif
790 endif ! the case of flux_type not equal 0, 1, or 2 is caught in advect_tracer above.
791 endif ! advect_this_tracer
792 enddo
793
794 endif ; enddo ! End of j-loop.
795
796 ! Do user controlled underflow of the tracer concentrations.
797 if (flux_type == 0) then ! Only update tracer if using residual flux
798 do m=1,ntr ; if (tr(m)%conc_underflow > 0.0) then
799 do j=js,je ; do i=is,ie
800 if (abs(tr(m)%t(i,j,k)) < tr(m)%conc_underflow) tr(m)%t(i,j,k) = 0.0
801 enddo ; enddo
802 endif ; enddo
803 endif
804
805 ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active.
806
807 if (flux_type == 0) then ! Only update tracer if using residual flux
808 !$OMP ordered
809 do m=1,ntr ; if (associated(tr(m)%ad2d_x)) then
810 do j=js,je ; if (domore_u_initial(j,k)) then
811 do i=is-1,ie
812 tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,j,m)*idt
813 enddo
814 endif ; enddo
815 endif ; enddo ! End of m-loop.
816 !$OMP end ordered
817 endif
818
819end subroutine advect_x
820
821!> This subroutine does 1-d flux-form advection using a monotonic piecewise
822!! linear scheme.
823subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
824 is, ie, js, je, k, G, GV, US, flux_type, advect_this_tracer, &
825 advect_schemes)
826 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
827 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
828 integer, intent(in) :: ntr !< The number of tracers
829 type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on
830 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous
831 !! tracer change [H L2 ~> m3 or kg]
832 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhr !< accumulated volume/mass flux through
833 !! the meridional face [H L2 ~> m3 or kg]
834 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect !< A tiny meridional mass flux that can
835 !! be neglected [H L2 ~> m3 or kg]
836 type(ocean_obc_type), pointer :: OBC !< specifies whether, where, and what OBCs are used
837 logical, dimension(SZJB_(G),SZK_(GV)), intent(inout) :: domore_v !< If true, there is more advection to be
838 !! done in this v-row
839 real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1]
840 integer, intent(in) :: is !< The starting tracer i-index to work on
841 integer, intent(in) :: ie !< The ending tracer i-index to work on
842 integer, intent(in) :: js !< The starting tracer j-index to work on
843 integer, intent(in) :: je !< The ending tracer j-index to work on
844 integer, intent(in) :: k !< The k-level to work on
845 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
846 integer, intent(in) :: flux_type !< Indicates whether uhtr, vhtr are the flux
847 !! due to the residual (= 0), resolved (= 1),
848 !! or parameterized (= 2) flow
849 logical, dimension(ntr), intent(in) :: advect_this_tracer !< If true, advect this tracer
850 integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use
851
852 real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
853 slope_y ! The concentration slope per grid point [conc].
854 real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
855 flux_y ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc].
856 real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
857 T_tmp ! The copy of the tracer concentration at constant i,k [conc].
858 real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
859 ! current iteration [H L2 ~> m3 or kg].
860 real :: hup, hlos ! hup is the upwind volume, hlos is the
861 ! part of that volume that might be lost
862 ! due to advection out the other side of
863 ! the grid box, both in [H L2 ~> m3 or kg].
864 real, dimension(SZIB_(G)) :: &
865 hlst, & ! Work variable [H L2 ~> m3 or kg].
866 Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1].
867 CFL ! The absolute value of the advective upwind-cell CFL number [nondim].
868 real :: min_h ! The minimum thickness that can be realized during
869 ! any of the passes [H ~> m or kg m-2].
870 real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2].
871 real :: h_neglect ! A thickness that is so small it is usually lost
872 ! in roundoff and can be neglected [H ~> m or kg m-2].
873 real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc]
874 real :: dMx ! Difference between the maximum of the surrounding cell concentrations and
875 ! the value in the cell whose reconstruction is being found [conc]
876 real :: dMn ! Difference between the tracer average in the cell whose reconstruction
877 ! is being found and the minimum of the surrounding values [conc]
878 real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc]
879 real :: dA ! Difference between the reconstruction tracer edge values [conc]
880 real :: mA ! Average of the reconstruction tracer edge values [conc]
881 real :: a6 ! Curvature of the reconstruction tracer values [conc]
882 logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
883 logical :: do_i(SZI_(G), SZJ_(G)) ! If true, work on given points.
884 logical :: usePLMslope
885 integer :: i, j, j2, m, n, j_up, stencil, ntr_id
886 type(obc_segment_type), pointer :: segment=>null()
887 logical :: domore_v_initial(SZJB_(G)) ! Initial state of domore_v
888
889 useplmslope = .false.
890 ! stencil for calculating slope values
891 stencil = 1
892 do m = 1,ntr
893 if ((advect_schemes(m) == advect_plm) .or. (advect_schemes(m) == advect_ppm)) &
894 useplmslope = .true.
895 if (advect_schemes(m) == advect_ppm) stencil = 2
896 enddo
897
898 min_h = 0.1*gv%Angstrom_H
899 tiny_h = tiny(min_h)
900 h_neglect = gv%H_subroundoff
901
902 ! We conditionally perform work on tracer points: calculating the PLM slope,
903 ! and updating tracer concentration within a cell
904 ! this depends on whether there is a flux which would affect this tracer point,
905 ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires
906 ! slope calculations at the two tracer points on either side (as indicated by
907 ! the stencil variable), so we account for this with the do_j_tr flag array
908 !
909 ! Note: this does lead to unnecessary work in updating tracer concentrations,
910 ! since that doesn't need a wider stencil with the PPM advection scheme, but
911 ! this would require an additional loop, etc.
912 do_j_tr(:) = .false.
913 do j=js-1,je
914 if (domore_v(j,k)) then ; do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ; enddo ; endif
915 enddo
916 domore_v_initial(:) = domore_v(:,k)
917
918 ! Calculate the j-direction profiles (slopes) of each tracer that
919 ! is being advected.
920 if (useplmslope) then
921 do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr
922 if (advect_this_tracer(m)) then
923 do i=is,ie
924 !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < &
925 ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then
926 ! maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))
927 !else
928 ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))
929 !endif
930 !if ((Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k)) < 0.0) then
931 ! slope_y(i,m,j) = 0.0
932 !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then
933 ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * &
934 ! 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))
935 !else
936 ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope
937 !endif
938 tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
939 dmx = max( tp, tc, tm ) - tc
940 dmn = tc - min( tp, tc, tm )
941 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
942 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
943 enddo
944 endif ! advect_this_tracer
945 enddo ; endif ; enddo ! End of i-, m-, & j- loops.
946 endif ! usePLMslope
947
948
949 ! make a copy of the tracers in case values need to be overridden for OBCs
950
951 do j=g%jsd,g%jed ; do m=1,ntr
952 if (advect_this_tracer(m)) then
953 do i=g%isd,g%ied
954 t_tmp(i,m,j) = tr(m)%t(i,j,k)
955 enddo
956 endif ! advect_this_tracer
957 enddo ; enddo
958
959 ! loop through open boundaries and recalculate flux terms
960 if (associated(obc)) then ; if (obc%OBC_pe) then
961 do n=1,obc%number_of_segments
962 segment=>obc%segment(n)
963 if (.not. associated(segment%tr_Reg)) cycle
964 do i=is,ie
965 if (segment%is_N_or_S) then
966 if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
967 j = segment%HI%JsdB
968 do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values
969 ntr_id = segment%tr_reg%Tr(m)%ntr_index
970 if (advect_this_tracer(ntr_id)) then
971 if (segment%direction == obc_direction_s) then
972 t_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
973 else
974 t_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
975 endif
976 endif ! advect_this_tracer
977 enddo
978 do m = 1,ntr ! Apply update tracer values for slope calculation
979 if (advect_this_tracer(m)) then
980 do j=segment%HI%JsdB-1,segment%HI%JsdB+1
981 tp = t_tmp(i,m,j+1) ; tc = t_tmp(i,m,j) ; tm = t_tmp(i,m,j-1)
982 dmx = max( tp, tc, tm ) - tc
983 dmn= tc - min( tp, tc, tm )
984 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
985 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
986 enddo
987 endif ! advect_this_tracer
988 enddo
989 endif
990 endif ! is_N_S
991 enddo ! i-loop
992 enddo ! segment loop
993 endif ; endif
994
995 ! Calculate the j-direction fluxes of each tracer, using as much
996 ! the minimum of the remaining mass flux (vhr) and the half the mass
997 ! in the cell plus whatever part of its half of the mass flux that
998 ! the flux through the other side does not require.
999 do j=js-1,je ; if (domore_v(j,k)) then
1000 domore_v(j,k) = .false.
1001
1002 do i=is,ie
1003 if ((vhr(i,j,k) == 0.0) .or. &
1004 ((vhr(i,j,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. &
1005 ((vhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
1006 vhh(i,j) = 0.0
1007 cfl(i) = 0.0
1008 elseif (vhr(i,j,k) < 0.0) then
1009 hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
1010 hlos = max(0.0, vhr(i,j+1,k))
1011 if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
1012 ((0.5*hup + vhr(i,j,k)) < 0.0)) then
1013 vhh(i,j) = min(-0.5*hup, -hup+hlos, 0.0)
1014 domore_v(j,k) = .true.
1015 else
1016 vhh(i,j) = vhr(i,j,k)
1017 endif
1018 cfl(i) = - vhh(i,j) / hprev(i,j+1,k) ! CFL is positive
1019 else
1020 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
1021 hlos = max(0.0, -vhr(i,j-1,k))
1022 if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
1023 ((0.5*hup - vhr(i,j,k)) < 0.0)) then
1024 vhh(i,j) = max(0.5*hup, hup-hlos, 0.0)
1025 domore_v(j,k) = .true.
1026 else
1027 vhh(i,j) = vhr(i,j,k)
1028 endif
1029 cfl(i) = vhh(i,j) / hprev(i,j,k) ! CFL is positive
1030 endif
1031 enddo
1032
1033 do m=1,ntr
1034 if (advect_this_tracer(m)) then
1035 if ((advect_schemes(m) == advect_ppm) .or. (advect_schemes(m) == advect_ppmh3)) then
1036 do i=is,ie
1037 ! centre cell depending on upstream direction
1038 if (vhh(i,j) >= 0.0) then
1039 j_up = j
1040 else
1041 j_up = j + 1
1042 endif
1043
1044 ! Implementation of PPM-H3
1045 tp = t_tmp(i,m,j_up+1) ; tc = t_tmp(i,m,j_up) ; tm = t_tmp(i,m,j_up-1)
1046
1047 if (advect_schemes(m) == advect_ppmh3) then
1048 al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
1049 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
1050 ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
1051 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
1052 else
1053 al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
1054 ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
1055 endif
1056
1057 da = ar - al ; ma = 0.5*( ar + al )
1058 if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.) then
1059 al = tc ; ar = tc ! PCM for local extrema and boundary cells
1060 elseif ( da*(tc-ma) > (da*da)/6. ) then
1061 al = (3.*tc) - 2.*ar
1062 elseif ( da*(tc-ma) < - (da*da)/6. ) then
1063 ar = (3.*tc) - 2.*al
1064 endif
1065
1066 a6 = 6.*tc - 3. * (ar + al) ! Curvature
1067
1068 if (vhh(i,j) >= 0.0) then
1069 flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
1070 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
1071 else
1072 flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
1073 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
1074 endif
1075 enddo
1076 else ! PLM
1077 do i=is,ie
1078 if (vhh(i,j) >= 0.0) then
1079 ! Indirect implementation of PLM
1080 !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j)
1081 !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
1082 !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) )
1083 ! Alternative implementation of PLM
1084 tc = t_tmp(i,m,j)
1085 flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
1086 else
1087 ! Indirect implementation of PLM
1088 !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
1089 !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1)
1090 !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) )
1091 ! Alternative implementation of PLM
1092 tc = t_tmp(i,m,j+1)
1093 flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
1094 endif
1095 enddo
1096 endif ! usePPM
1097 endif ! advect_this_tracer
1098 enddo
1099
1100 if (associated(obc)) then ; if (obc%OBC_pe) then
1101 if (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) then
1102 do n=1,obc%number_of_segments
1103 segment=>obc%segment(n)
1104 if (.not. segment%specified) cycle
1105 if (.not. associated(segment%tr_Reg)) cycle
1106 if (obc%segment(n)%is_N_or_S) then
1107 if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB) then
1108 do i=segment%HI%isd,segment%HI%ied
1109 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
1110 ! Now changing to simply fixed inflows.
1111 if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
1112 (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n)) then
1113 vhh(i,j) = vhr(i,j,k)
1114 do m=1,segment%tr_Reg%ntseg
1115 ntr_id = segment%tr_reg%Tr(m)%ntr_index
1116 if (advect_this_tracer(ntr_id)) then
1117 flux_y(i,ntr_id,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
1118 endif ! advect_this_tracer
1119 enddo
1120 endif
1121 enddo
1122 endif
1123 endif
1124 enddo
1125 endif
1126
1127 if (obc%open_v_BCs_exist_globally) then
1128 do n=1,obc%number_of_segments
1129 segment=>obc%segment(n)
1130 if (segment%specified) cycle
1131 if (.not. associated(segment%tr_Reg)) cycle
1132 if (segment%is_N_or_S .and. (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)) then
1133 do i=segment%HI%isd,segment%HI%ied
1134 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
1135 if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
1136 (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5)) then
1137 vhh(i,j) = vhr(i,j,k)
1138 do m=1,segment%tr_Reg%ntseg
1139 ntr_id = segment%tr_reg%Tr(m)%ntr_index
1140 if (advect_this_tracer(ntr_id)) then
1141 flux_y(i,ntr_id,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
1142 endif ! advect_this_tracer
1143 enddo
1144 endif
1145 enddo
1146 endif
1147 enddo
1148 endif
1149 endif ; endif
1150
1151 else ! not domore_v.
1152 do i=is,ie ; vhh(i,j) = 0.0 ; enddo
1153 do m=1,ntr
1154 if (advect_this_tracer(m)) then
1155 do i=is,ie ; flux_y(i,m,j) = 0.0 ; enddo
1156 endif ! advect_this_tracer
1157 enddo
1158 endif ; enddo ! End of j-loop
1159
1160 do j=js-1,je ; do i=is,ie
1161 vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
1162 if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
1163 enddo ; enddo
1164
1165 ! Calculate new tracer concentration in each cell after accounting
1166 ! for the j-direction fluxes.
1167 do j=js,je ; if (do_j_tr(j)) then
1168 do i=is,ie
1169 if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0)) then
1170 do_i(i,j) = .true.
1171 hlst(i) = hprev(i,j,k)
1172 hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
1173 if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
1174 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
1175 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
1176 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
1177 else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
1178 else ; do_i(i,j) = .false. ; endif
1179 enddo
1180
1181 ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
1182 if (associated(obc)) then
1183 if ((.not.obc%exterior_OBC_bug) .and. (obc%OBC_pe) .and. &
1184 (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) then
1185 ! OBC_DIRECTION_N / OBC_DIRECTION_S on the south / north edge
1186 do i=is,ie ; if ((obc%segnum_v(i,j-1) > 0) .or. (obc%segnum_v(i,j) < 0)) &
1187 do_i(i,j) = .false.
1188 enddo
1189 endif
1190 endif
1191
1192 ! update tracer and save some diagnostics
1193 do m=1,ntr
1194 if (advect_this_tracer(m)) then
1195 if (flux_type == 0) then ! Only update tracer if using residual flux
1196 do i=is,ie ; if (do_i(i,j)) then
1197 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
1198 (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
1199 endif ; enddo
1200
1201 ! diagnose convergence of flux_y and add to convergence of flux_x.
1202 ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
1203 if (associated(tr(m)%advection_xy)) then
1204 do i=is,ie ; if (do_i(i,j)) then
1205 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_y(i,m,j) - flux_y(i,m,j-1))* idt * &
1206 g%IareaT(i,j)
1207 endif ; enddo
1208 endif
1209 endif ! flux_type == 0
1210 endif ! advect_this_tracer
1211 enddo
1212 endif ; enddo ! End of j-loop.
1213
1214 ! Do user controlled underflow of the tracer concentrations.
1215 if (flux_type == 0) then ! Only update tracer if using residual flux
1216 do m=1,ntr ; if (tr(m)%conc_underflow > 0.0) then
1217 do j=js,je ; do i=is,ie
1218 if (abs(tr(m)%t(i,j,k)) < tr(m)%conc_underflow) tr(m)%t(i,j,k) = 0.0
1219 enddo ; enddo
1220 endif ; enddo
1221 endif
1222
1223 ! compute ad_y and ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active.
1224 if (flux_type == 0) then
1225 !$OMP ordered
1226 do m=1,ntr ; if (associated(tr(m)%ad_y)) then
1227 do j=js-1,je ; if (domore_v_initial(j)) then
1228 do i=is,ie
1229 tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
1230 enddo
1231 endif ; enddo
1232 endif ; enddo ! End of m-loop.
1233
1234 do m=1,ntr ; if (associated(tr(m)%ad2d_y)) then
1235 do j=js-1,je ; if (domore_v_initial(j)) then
1236 do i=is,ie
1237 tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
1238 enddo
1239 endif ; enddo
1240 endif ; enddo ! End of m-loop.
1241 !$OMP end ordered
1242 elseif (flux_type == 1) then
1243 !$OMP ordered
1244 do m=1,ntr ; if (associated(tr(m)%ad_y_resolved)) then
1245 do j=js-1,je ; if (domore_v_initial(j)) then
1246 do i=is,ie
1247 tr(m)%ad_y_resolved(i,j,k) = tr(m)%ad_y_resolved(i,j,k) + flux_y(i,m,j)*idt
1248 enddo
1249 endif ; enddo
1250 endif ; enddo ! End of m-loop.
1251 !$OMP end ordered
1252 elseif (flux_type == 2) then
1253 !$OMP ordered
1254 do m=1,ntr ; if (associated(tr(m)%ad_y_param)) then
1255 do j=js-1,je ; if (domore_v_initial(j)) then
1256 do i=is,ie
1257 tr(m)%ad_y_param(i,j,k) = tr(m)%ad_y_param(i,j,k) + flux_y(i,m,j)*idt
1258 enddo
1259 endif ; enddo
1260 endif ; enddo ! End of m-loop.
1261 !$OMP end ordered
1262 endif ! the case of flux_type not equal 0, 1, or 2 is caught in advect_tracer above.
1263
1264end subroutine advect_y
1265
1266!> Initialize lateral tracer advection module
1267subroutine tracer_advect_init(Time, G, US, param_file, diag, CS)
1268 type(time_type), target, intent(in) :: time !< current model time
1269 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1270 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1271 type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
1272 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
1273 type(tracer_advect_cs), pointer :: cs !< module control structure
1274
1275 ! This include declares and sets the variable "version".
1276# include "version_variable.h"
1277 character(len=40) :: mdl = "MOM_tracer_advect" ! This module's name.
1278 character(len=256) :: mesg ! Message for error messages.
1279
1280 if (associated(cs)) then
1281 call mom_error(warning, "tracer_advect_init called with associated control structure.")
1282 return
1283 endif
1284 allocate(cs)
1285
1286 cs%diag => diag
1287
1288 ! Read all relevant parameters and write them to the model log.
1289 call log_version(param_file, mdl, version, "")
1290 call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
1291 desc="The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T)
1292 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1293 call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, &
1294 desc="The horizontal transport scheme for tracers:\n"//&
1295 trim(traceradvectionschemedoc), default='PLM')
1296
1297 ! Get the integer value of the tracer scheme
1298 call set_tracer_advect_scheme(cs%default_advect_scheme, mesg)
1299
1300 if (cs%default_advect_scheme == advect_ppmh3) then
1301 call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", &
1302 cs%useHuynhStencilBug, &
1303 desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " &
1304 // "This is incorrect and will produce regressions in certain " &
1305 // "configurations, but may be required to reproduce results in " &
1306 // "legacy simulations.", &
1307 default=.false.)
1308 endif
1309
1310 id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=clock_module)
1311 id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1312 id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1313
1314end subroutine tracer_advect_init
1315
1316!> Close the tracer advection module
1317subroutine tracer_advect_end(CS)
1318 type(tracer_advect_cs), pointer :: cs !< module control structure
1319
1320 if (associated(cs)) deallocate(cs)
1321
1322end subroutine tracer_advect_end
1323
1324
1325!> \namespace mom_tracer_advect
1326!!
1327!! This program contains the subroutines that advect tracers
1328!! horizontally (i.e. along layers).
1329!!
1330!! \section section_mom_advect_intro
1331!!
1332!! * advect_tracer advects tracer concentrations using a combination
1333!! of the modified flux advection scheme from Easter (Mon. Wea. Rev.,
1334!! 1993) with tracer distributions given by the monotonic piecewise
1335!! parabolic method, as described in Carpenter et al. (MWR, 1990).
1336!! This scheme conserves the total amount of tracer while avoiding
1337!! spurious maxima and minima of the tracer concentration.
1338!!
1339!! * advect_tracer subroutine determines the volume of a layer in
1340!! a grid cell at the previous instance when the tracer concentration
1341!! was changed, so it is essential that the volume fluxes should be
1342!! correct. It is also important that the tracer advection occurs
1343!! before each calculation of the diabatic forcing.
1344!!
1345!! The advection scheme of some tracers can be set to be different
1346!! to that used by active tracers.
1347
1348
1349end module mom_tracer_advect