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 (segment%direction == obc_direction_w) then
482 t_tmp(i,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
483 else
484 t_tmp(i+1,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
485 endif
486 enddo
487 do m = 1,ntr ! Apply update tracer values for slope calculation
488 do i=segment%HI%IsdB-1,segment%HI%IsdB+1
489 tp = t_tmp(i+1,m) ; tc = t_tmp(i,m) ; tm = t_tmp(i-1,m)
490 dmx = max( tp, tc, tm ) - tc
491 dmn= tc - min( tp, tc, tm )
492 slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
493 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
494 enddo
495 enddo
496
497 endif
498 endif
499 enddo
500 endif ; endif
501
502
503 ! Calculate the i-direction fluxes of each tracer, using as much
504 ! the minimum of the remaining mass flux (uhr) and the half the mass
505 ! in the cell plus whatever part of its half of the mass flux that
506 ! the flux through the other side does not require.
507 do i=is-1,ie
508 if ((uhr(i,j,k) == 0.0) .or. &
509 ((uhr(i,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. &
510 ((uhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
511 uhh(i) = 0.0
512 cfl(i) = 0.0
513 elseif (uhr(i,j,k) < 0.0) then
514 hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
515 hlos = max(0.0, uhr(i+1,j,k))
516 if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
517 ((0.5*hup + uhr(i,j,k)) < 0.0)) then
518 uhh(i) = min(-0.5*hup, -hup+hlos, 0.0)
519 domore_u(j,k) = .true.
520 else
521 uhh(i) = uhr(i,j,k)
522 endif
523 cfl(i) = - uhh(i) / (hprev(i+1,j,k)) ! CFL is positive
524 else
525 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
526 hlos = max(0.0, -uhr(i-1,j,k))
527 if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
528 ((0.5*hup - uhr(i,j,k)) < 0.0)) then
529 uhh(i) = max(0.5*hup, hup-hlos, 0.0)
530 domore_u(j,k) = .true.
531 else
532 uhh(i) = uhr(i,j,k)
533 endif
534 cfl(i) = uhh(i) / (hprev(i,j,k)) ! CFL is positive
535 endif
536 enddo
537
538 do m=1,ntr
539
540 if ((advect_schemes(m) == advect_ppm) .or. (advect_schemes(m) == advect_ppmh3)) then
541 do i=is-1,ie
542 ! centre cell depending on upstream direction
543 if (uhh(i) >= 0.0) then
544 i_up = i
545 else
546 i_up = i+1
547 endif
548
549 ! Implementation of PPM-H3
550 tp = t_tmp(i_up+1,m) ; tc = t_tmp(i_up,m) ; tm = t_tmp(i_up-1,m)
551
552 if (advect_schemes(m) == advect_ppmh3) then
553 al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
554 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
555 ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
556 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
557 else
558 al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
559 ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
560 endif
561
562 da = ar - al ; ma = 0.5*( ar + al )
563 if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.) then
564 al = tc ; ar = tc ! PCM for local extrema and boundary cells
565 elseif ( da*(tc-ma) > (da*da)/6. ) then
566 al = (3.*tc) - 2.*ar
567 elseif ( da*(tc-ma) < - (da*da)/6. ) then
568 ar = (3.*tc) - 2.*al
569 endif
570
571 a6 = 6.*tc - 3. * (ar + al) ! Curvature
572
573 if (uhh(i) >= 0.0) then
574 flux_x(i,j,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
575 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
576 else
577 flux_x(i,j,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
578 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
579 endif
580 enddo
581 else ! PLM
582 do i=is-1,ie
583 if (uhh(i) >= 0.0) then
584 ! Indirect implementation of PLM
585 !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m)
586 !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
587 !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) )
588 ! Alternative implementation of PLM
589 tc = t_tmp(i,m)
590 flux_x(i,j,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
591 else
592 ! Indirect implementation of PLM
593 !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
594 !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m)
595 !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) )
596 ! Alternative implementation of PLM
597 tc = t_tmp(i+1,m)
598 flux_x(i,j,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
599 endif
600 enddo
601 endif ! usePPM
602 enddo
603
604 if (associated(obc)) then ; if (obc%OBC_pe) then
605 if (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally) then
606 do n=1,obc%number_of_segments
607 segment=>obc%segment(n)
608 if (.not. associated(segment%tr_Reg)) cycle
609 if (segment%is_E_or_W) then
610 if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
611 i = segment%HI%IsdB
612 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
613 ! Now changing to simply fixed inflows.
614 if ((uhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_w) .or. &
615 (uhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_e)) then
616 uhh(i) = uhr(i,j,k)
617 ! should the reservoir evolve for this case Kate ?? - Nope
618 do m=1,segment%tr_Reg%ntseg
619 ntr_id = segment%tr_reg%Tr(m)%ntr_index
620 flux_x(i,j,ntr_id) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
621 enddo
622 endif
623 endif
624 endif
625 enddo
626 endif
627
628 if (obc%open_u_BCs_exist_globally) then
629 do n=1,obc%number_of_segments
630 segment=>obc%segment(n)
631 i = segment%HI%IsdB
632 if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
633 if (segment%specified) cycle
634 if (.not. associated(segment%tr_Reg)) cycle
635
636 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
637 if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
638 (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5)) then
639 uhh(i) = uhr(i,j,k)
640 do m=1,segment%tr_Reg%ntseg
641 ntr_id = segment%tr_reg%Tr(m)%ntr_index
642 flux_x(i,j,ntr_id) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
643 enddo
644 endif
645 endif
646 enddo
647 endif
648 endif ; endif
649
650 ! Calculate new tracer concentration in each cell after accounting
651 ! for the i-direction fluxes.
652 do i=is-1,ie
653 uhr(i,j,k) = uhr(i,j,k) - uhh(i)
654 if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
655 enddo
656 do i=is,ie
657 if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0)) then
658 do_i(i,j) = .true.
659 hlst(i) = hprev(i,j,k)
660 hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
661 if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
662 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
663 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
664 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
665 else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
666 else
667 do_i(i,j) = .false.
668 endif
669 enddo
670
671 ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
672 if (associated(obc)) then
673 if ((.not.obc%exterior_OBC_bug) .and. (obc%OBC_pe) .and. &
674 (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally)) then
675 ! OBC_DIRECTION_E / OBC_DIRECTION_W on the west / east edge
676 do i=is,ie ; if ((obc%segnum_u(i-1,j) > 0) .or. (obc%segnum_u(i,j) < 0)) &
677 do_i(i,j) = .false.
678 enddo
679 endif
680 endif
681
682 ! update tracer concentration from i-flux and save some diagnostics
683 do m=1,ntr
684
685 ! update tracer
686 do i=is,ie
687 if (do_i(i,j)) then
688 if (ihnew(i) > 0.0) then
689 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
690 (flux_x(i,j,m) - flux_x(i-1,j,m))) * ihnew(i)
691 endif
692 endif
693 enddo
694
695 ! diagnostics
696 if (associated(tr(m)%ad_x)) then ; do i=is-1,ie
697 tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,j,m)*idt
698 enddo ; endif
699
700 ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic).
701 ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
702 if (associated(tr(m)%advection_xy)) then
703 do i=is,ie ; if (do_i(i,j)) then
704 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - &
705 (flux_x(i,j,m) - flux_x(i-1,j,m)) * &
706 idt * g%IareaT(i,j)
707 endif ; enddo
708 endif
709
710 enddo
711
712 endif ; enddo ! End of j-loop.
713
714 ! Do user controlled underflow of the tracer concentrations.
715 do m=1,ntr ; if (tr(m)%conc_underflow > 0.0) then
716 do j=js,je ; do i=is,ie
717 if (abs(tr(m)%t(i,j,k)) < tr(m)%conc_underflow) tr(m)%t(i,j,k) = 0.0
718 enddo ; enddo
719 endif ; enddo
720
721 ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active.
722
723 !$OMP ordered
724 do m=1,ntr ; if (associated(tr(m)%ad2d_x)) then
725 do j=js,je ; if (domore_u_initial(j,k)) then
726 do i=is-1,ie
727 tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,j,m)*idt
728 enddo
729 endif ; enddo
730 endif ; enddo ! End of m-loop.
731 !$OMP end ordered
732
733end subroutine advect_x
734
735!> This subroutine does 1-d flux-form advection using a monotonic piecewise
736!! linear scheme.
737subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
738 is, ie, js, je, k, G, GV, US, advect_schemes)
739 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
740 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
741 integer, intent(in) :: ntr !< The number of tracers
742 type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on
743 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: hprev !< cell volume at the end of previous
744 !! tracer change [H L2 ~> m3 or kg]
745 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhr !< accumulated volume/mass flux through
746 !! the meridional face [H L2 ~> m3 or kg]
747 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect !< A tiny meridional mass flux that can
748 !! be neglected [H L2 ~> m3 or kg]
749 type(ocean_obc_type), pointer :: OBC !< specifies whether, where, and what OBCs are used
750 logical, dimension(SZJB_(G),SZK_(GV)), intent(inout) :: domore_v !< If true, there is more advection to be
751 !! done in this v-row
752 real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1]
753 integer, intent(in) :: is !< The starting tracer i-index to work on
754 integer, intent(in) :: ie !< The ending tracer i-index to work on
755 integer, intent(in) :: js !< The starting tracer j-index to work on
756 integer, intent(in) :: je !< The ending tracer j-index to work on
757 integer, intent(in) :: k !< The k-level to work on
758 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
759 integer, dimension(ntr), intent(in) :: advect_schemes !< list of advection schemes to use
760
761 real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
762 slope_y ! The concentration slope per grid point [conc].
763 real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
764 flux_y ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc].
765 real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
766 T_tmp ! The copy of the tracer concentration at constant i,k [conc].
767 real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
768 ! current iteration [H L2 ~> m3 or kg].
769 real :: hup, hlos ! hup is the upwind volume, hlos is the
770 ! part of that volume that might be lost
771 ! due to advection out the other side of
772 ! the grid box, both in [H L2 ~> m3 or kg].
773 real, dimension(SZIB_(G)) :: &
774 hlst, & ! Work variable [H L2 ~> m3 or kg].
775 Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1].
776 CFL ! The absolute value of the advective upwind-cell CFL number [nondim].
777 real :: min_h ! The minimum thickness that can be realized during
778 ! any of the passes [H ~> m or kg m-2].
779 real :: tiny_h ! The smallest numerically invertible thickness [H ~> m or kg m-2].
780 real :: h_neglect ! A thickness that is so small it is usually lost
781 ! in roundoff and can be neglected [H ~> m or kg m-2].
782 real :: aR, aL ! Reconstructed tracer concentrations at the right and left edges [conc]
783 real :: dMx ! Difference between the maximum of the surrounding cell concentrations and
784 ! the value in the cell whose reconstruction is being found [conc]
785 real :: dMn ! Difference between the tracer average in the cell whose reconstruction
786 ! is being found and the minimum of the surrounding values [conc]
787 real :: Tp, Tc, Tm ! Tracer concentrations around the upstream cell [conc]
788 real :: dA ! Difference between the reconstruction tracer edge values [conc]
789 real :: mA ! Average of the reconstruction tracer edge values [conc]
790 real :: a6 ! Curvature of the reconstruction tracer values [conc]
791 logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
792 logical :: do_i(SZI_(G), SZJ_(G)) ! If true, work on given points.
793 logical :: usePLMslope
794 integer :: i, j, j2, m, n, j_up, stencil, ntr_id
795 type(obc_segment_type), pointer :: segment=>null()
796 logical :: domore_v_initial(SZJB_(G)) ! Initial state of domore_v
797
798 useplmslope = .false.
799 ! stencil for calculating slope values
800 stencil = 1
801 do m = 1,ntr
802 if ((advect_schemes(m) == advect_plm) .or. (advect_schemes(m) == advect_ppm)) &
803 useplmslope = .true.
804 if (advect_schemes(m) == advect_ppm) stencil = 2
805 enddo
806
807 min_h = 0.1*gv%Angstrom_H
808 tiny_h = tiny(min_h)
809 h_neglect = gv%H_subroundoff
810
811 ! We conditionally perform work on tracer points: calculating the PLM slope,
812 ! and updating tracer concentration within a cell
813 ! this depends on whether there is a flux which would affect this tracer point,
814 ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires
815 ! slope calculations at the two tracer points on either side (as indicated by
816 ! the stencil variable), so we account for this with the do_j_tr flag array
817 !
818 ! Note: this does lead to unnecessary work in updating tracer concentrations,
819 ! since that doesn't need a wider stencil with the PPM advection scheme, but
820 ! this would require an additional loop, etc.
821 do_j_tr(:) = .false.
822 do j=js-1,je
823 if (domore_v(j,k)) then ; do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ; enddo ; endif
824 enddo
825 domore_v_initial(:) = domore_v(:,k)
826
827 ! Calculate the j-direction profiles (slopes) of each tracer that
828 ! is being advected.
829 if (useplmslope) then
830 do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr ; do i=is,ie
831 !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < &
832 ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then
833 ! maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))
834 !else
835 ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))
836 !endif
837 !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
838 ! slope_y(i,m,j) = 0.0
839 !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then
840 ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * &
841 ! 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))
842 !else
843 ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope
844 !endif
845 tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
846 dmx = max( tp, tc, tm ) - tc
847 dmn = tc - min( tp, tc, tm )
848 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
849 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
850 enddo ; enddo ; endif ; enddo ! End of i-, m-, & j- loops.
851 endif ! usePLMslope
852
853
854 ! make a copy of the tracers in case values need to be overridden for OBCs
855
856 do j=g%jsd,g%jed ; do m=1,ntr ; do i=g%isd,g%ied
857 t_tmp(i,m,j) = tr(m)%t(i,j,k)
858 enddo ; enddo ; enddo
859
860 ! loop through open boundaries and recalculate flux terms
861 if (associated(obc)) then ; if (obc%OBC_pe) then
862 do n=1,obc%number_of_segments
863 segment=>obc%segment(n)
864 if (.not. associated(segment%tr_Reg)) cycle
865 do i=is,ie
866 if (segment%is_N_or_S) then
867 if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
868 j = segment%HI%JsdB
869 do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values
870 ntr_id = segment%tr_reg%Tr(m)%ntr_index
871 if (segment%direction == obc_direction_s) then
872 t_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
873 else
874 t_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
875 endif
876 enddo
877 do m = 1,ntr ! Apply update tracer values for slope calculation
878 do j=segment%HI%JsdB-1,segment%HI%JsdB+1
879 tp = t_tmp(i,m,j+1) ; tc = t_tmp(i,m,j) ; tm = t_tmp(i,m,j-1)
880 dmx = max( tp, tc, tm ) - tc
881 dmn= tc - min( tp, tc, tm )
882 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
883 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
884 enddo
885 enddo
886 endif
887 endif ! is_N_S
888 enddo ! i-loop
889 enddo ! segment loop
890 endif ; endif
891
892 ! Calculate the j-direction fluxes of each tracer, using as much
893 ! the minimum of the remaining mass flux (vhr) and the half the mass
894 ! in the cell plus whatever part of its half of the mass flux that
895 ! the flux through the other side does not require.
896 do j=js-1,je ; if (domore_v(j,k)) then
897 domore_v(j,k) = .false.
898
899 do i=is,ie
900 if ((vhr(i,j,k) == 0.0) .or. &
901 ((vhr(i,j,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. &
902 ((vhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
903 vhh(i,j) = 0.0
904 cfl(i) = 0.0
905 elseif (vhr(i,j,k) < 0.0) then
906 hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
907 hlos = max(0.0, vhr(i,j+1,k))
908 if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
909 ((0.5*hup + vhr(i,j,k)) < 0.0)) then
910 vhh(i,j) = min(-0.5*hup, -hup+hlos, 0.0)
911 domore_v(j,k) = .true.
912 else
913 vhh(i,j) = vhr(i,j,k)
914 endif
915 cfl(i) = - vhh(i,j) / hprev(i,j+1,k) ! CFL is positive
916 else
917 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
918 hlos = max(0.0, -vhr(i,j-1,k))
919 if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
920 ((0.5*hup - vhr(i,j,k)) < 0.0)) then
921 vhh(i,j) = max(0.5*hup, hup-hlos, 0.0)
922 domore_v(j,k) = .true.
923 else
924 vhh(i,j) = vhr(i,j,k)
925 endif
926 cfl(i) = vhh(i,j) / hprev(i,j,k) ! CFL is positive
927 endif
928 enddo
929
930 do m=1,ntr
931
932 if ((advect_schemes(m) == advect_ppm) .or. (advect_schemes(m) == advect_ppmh3)) then
933 do i=is,ie
934 ! centre cell depending on upstream direction
935 if (vhh(i,j) >= 0.0) then
936 j_up = j
937 else
938 j_up = j + 1
939 endif
940
941 ! Implementation of PPM-H3
942 tp = t_tmp(i,m,j_up+1) ; tc = t_tmp(i,m,j_up) ; tm = t_tmp(i,m,j_up-1)
943
944 if (advect_schemes(m) == advect_ppmh3) then
945 al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
946 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
947 ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
948 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
949 else
950 al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
951 ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
952 endif
953
954 da = ar - al ; ma = 0.5*( ar + al )
955 if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.) then
956 al = tc ; ar = tc ! PCM for local extrema and boundary cells
957 elseif ( da*(tc-ma) > (da*da)/6. ) then
958 al = (3.*tc) - 2.*ar
959 elseif ( da*(tc-ma) < - (da*da)/6. ) then
960 ar = (3.*tc) - 2.*al
961 endif
962
963 a6 = 6.*tc - 3. * (ar + al) ! Curvature
964
965 if (vhh(i,j) >= 0.0) then
966 flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
967 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
968 else
969 flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
970 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
971 endif
972 enddo
973 else ! PLM
974 do i=is,ie
975 if (vhh(i,j) >= 0.0) then
976 ! Indirect implementation of PLM
977 !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j)
978 !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
979 !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) )
980 ! Alternative implementation of PLM
981 tc = t_tmp(i,m,j)
982 flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
983 else
984 ! Indirect implementation of PLM
985 !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
986 !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1)
987 !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) )
988 ! Alternative implementation of PLM
989 tc = t_tmp(i,m,j+1)
990 flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
991 endif
992 enddo
993 endif ! usePPM
994 enddo
995
996 if (associated(obc)) then ; if (obc%OBC_pe) then
997 if (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) then
998 do n=1,obc%number_of_segments
999 segment=>obc%segment(n)
1000 if (.not. segment%specified) cycle
1001 if (.not. associated(segment%tr_Reg)) cycle
1002 if (obc%segment(n)%is_N_or_S) then
1003 if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB) then
1004 do i=segment%HI%isd,segment%HI%ied
1005 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
1006 ! Now changing to simply fixed inflows.
1007 if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
1008 (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n)) then
1009 vhh(i,j) = vhr(i,j,k)
1010 do m=1,segment%tr_Reg%ntseg
1011 ntr_id = segment%tr_reg%Tr(m)%ntr_index
1012 flux_y(i,ntr_id,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
1013 enddo
1014 endif
1015 enddo
1016 endif
1017 endif
1018 enddo
1019 endif
1020
1021 if (obc%open_v_BCs_exist_globally) then
1022 do n=1,obc%number_of_segments
1023 segment=>obc%segment(n)
1024 if (segment%specified) cycle
1025 if (.not. associated(segment%tr_Reg)) cycle
1026 if (segment%is_N_or_S .and. (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)) then
1027 do i=segment%HI%isd,segment%HI%ied
1028 ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
1029 if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
1030 (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5)) then
1031 vhh(i,j) = vhr(i,j,k)
1032 do m=1,segment%tr_Reg%ntseg
1033 ntr_id = segment%tr_reg%Tr(m)%ntr_index
1034 flux_y(i,ntr_id,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
1035 enddo
1036 endif
1037 enddo
1038 endif
1039 enddo
1040 endif
1041 endif ; endif
1042
1043 else ! not domore_v.
1044 do i=is,ie ; vhh(i,j) = 0.0 ; enddo
1045 do m=1,ntr ; do i=is,ie ; flux_y(i,m,j) = 0.0 ; enddo ; enddo
1046 endif ; enddo ! End of j-loop
1047
1048 do j=js-1,je ; do i=is,ie
1049 vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
1050 if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
1051 enddo ; enddo
1052
1053 ! Calculate new tracer concentration in each cell after accounting
1054 ! for the j-direction fluxes.
1055 do j=js,je ; if (do_j_tr(j)) then
1056 do i=is,ie
1057 if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0)) then
1058 do_i(i,j) = .true.
1059 hlst(i) = hprev(i,j,k)
1060 hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
1061 if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
1062 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
1063 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
1064 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
1065 else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
1066 else ; do_i(i,j) = .false. ; endif
1067 enddo
1068
1069 ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
1070 if (associated(obc)) then
1071 if ((.not.obc%exterior_OBC_bug) .and. (obc%OBC_pe) .and. &
1072 (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) then
1073 ! OBC_DIRECTION_N / OBC_DIRECTION_S on the south / north edge
1074 do i=is,ie ; if ((obc%segnum_v(i,j-1) > 0) .or. (obc%segnum_v(i,j) < 0)) &
1075 do_i(i,j) = .false.
1076 enddo
1077 endif
1078 endif
1079
1080 ! update tracer and save some diagnostics
1081 do m=1,ntr
1082 do i=is,ie ; if (do_i(i,j)) then
1083 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
1084 (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
1085 endif ; enddo
1086
1087 ! diagnose convergence of flux_y and add to convergence of flux_x.
1088 ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
1089 if (associated(tr(m)%advection_xy)) then
1090 do i=is,ie ; if (do_i(i,j)) then
1091 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - &
1092 (flux_y(i,m,j) - flux_y(i,m,j-1))* idt * &
1093 g%IareaT(i,j)
1094 endif ; enddo
1095 endif
1096
1097 enddo
1098 endif ; enddo ! End of j-loop.
1099
1100 ! Do user controlled underflow of the tracer concentrations.
1101 do m=1,ntr ; if (tr(m)%conc_underflow > 0.0) then
1102 do j=js,je ; do i=is,ie
1103 if (abs(tr(m)%t(i,j,k)) < tr(m)%conc_underflow) tr(m)%t(i,j,k) = 0.0
1104 enddo ; enddo
1105 endif ; enddo
1106
1107 ! compute ad_y and ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active.
1108 !$OMP ordered
1109 do m=1,ntr ; if (associated(tr(m)%ad_y)) then
1110 do j=js-1,je ; if (domore_v_initial(j)) then
1111 do i=is,ie
1112 tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
1113 enddo
1114 endif ; enddo
1115 endif ; enddo ! End of m-loop.
1116
1117 do m=1,ntr ; if (associated(tr(m)%ad2d_y)) then
1118 do j=js-1,je ; if (domore_v_initial(j)) then
1119 do i=is,ie
1120 tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
1121 enddo
1122 endif ; enddo
1123 endif ; enddo ! End of m-loop.
1124 !$OMP end ordered
1125
1126end subroutine advect_y
1127
1128!> Initialize lateral tracer advection module
1129subroutine tracer_advect_init(Time, G, US, param_file, diag, CS)
1130 type(time_type), target, intent(in) :: time !< current model time
1131 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
1132 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1133 type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
1134 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
1135 type(tracer_advect_cs), pointer :: cs !< module control structure
1136
1137 ! This include declares and sets the variable "version".
1138# include "version_variable.h"
1139 character(len=40) :: mdl = "MOM_tracer_advect" ! This module's name.
1140 character(len=256) :: mesg ! Message for error messages.
1141
1142 if (associated(cs)) then
1143 call mom_error(warning, "tracer_advect_init called with associated control structure.")
1144 return
1145 endif
1146 allocate(cs)
1147
1148 cs%diag => diag
1149
1150 ! Read all relevant parameters and write them to the model log.
1151 call log_version(param_file, mdl, version, "")
1152 call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
1153 desc="The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T)
1154 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1155 call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, &
1156 desc="The horizontal transport scheme for tracers:\n"//&
1157 trim(traceradvectionschemedoc), default='PLM')
1158
1159 ! Get the integer value of the tracer scheme
1160 call set_tracer_advect_scheme(cs%default_advect_scheme, mesg)
1161
1162 if (cs%default_advect_scheme == advect_ppmh3) then
1163 call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", &
1164 cs%useHuynhStencilBug, &
1165 desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " &
1166 // "This is incorrect and will produce regressions in certain " &
1167 // "configurations, but may be required to reproduce results in " &
1168 // "legacy simulations.", &
1169 default=.false.)
1170 endif
1171
1172 id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=clock_module)
1173 id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1174 id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1175
1176end subroutine tracer_advect_init
1177
1178!> Close the tracer advection module
1179subroutine tracer_advect_end(CS)
1180 type(tracer_advect_cs), pointer :: cs !< module control structure
1181
1182 if (associated(cs)) deallocate(cs)
1183
1184end subroutine tracer_advect_end
1185
1186
1187!> \namespace mom_tracer_advect
1188!!
1189!! This program contains the subroutines that advect tracers
1190!! horizontally (i.e. along layers).
1191!!
1192!! \section section_mom_advect_intro
1193!!
1194!! * advect_tracer advects tracer concentrations using a combination
1195!! of the modified flux advection scheme from Easter (Mon. Wea. Rev.,
1196!! 1993) with tracer distributions given by the monotonic piecewise
1197!! parabolic method, as described in Carpenter et al. (MWR, 1990).
1198!! This scheme conserves the total amount of tracer while avoiding
1199!! spurious maxima and minima of the tracer concentration.
1200!!
1201!! * advect_tracer subroutine determines the volume of a layer in
1202!! a grid cell at the previous instance when the tracer concentration
1203!! was changed, so it is essential that the volume fluxes should be
1204!! correct. It is also important that the tracer advection occurs
1205!! before each calculation of the diabatic forcing.
1206!!
1207!! The advection scheme of some tracers can be set to be different
1208!! to that used by active tracers.
1209
1210
1211end module mom_tracer_advect