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