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