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