MOM6
Data Types | Functions/Subroutines | Variables
mom_tracer_advect Module Reference

Detailed Description

This program contains the subroutines that advect tracers along coordinate surfaces.

This program contains the subroutines that advect tracers horizontally (i.e. along layers).

section_mom_advect_intro

  • advect_tracer advects tracer concentrations using a combination of the modified flux advection scheme from Easter (Mon. Wea. Rev., 1993) with tracer distributions given by the monotonic modified van Leer scheme proposed by Lin et al. (Mon. Wea. Rev., 1994). This scheme conserves the total amount of tracer while avoiding spurious maxima and minima of the tracer concentration. If a higher order accuracy scheme is needed, suggest monotonic piecewise parabolic method, as described in Carpenter et al. (MWR, 1990).
  • advect_tracer has 4 arguments, described below. This subroutine determines the volume of a layer in a grid cell at the previous instance when the tracer concentration was changed, so it is essential that the volume fluxes should be correct. It is also important that the tracer advection occurs before each calculation of the diabatic forcing.

Data Types

type  tracer_advect_cs
 Control structure for this module. More...
 

Functions/Subroutines

subroutine, public advect_tracer (h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg, h_prev_opt, max_iter_in, x_first_in, uhr_out, vhr_out, h_out)
 This routine time steps the tracer concentration using a monotonic, conservative, weakly diffusive scheme. More...
 
subroutine advect_x (Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, is, ie, js, je, k, G, GV, usePPM, useHuynh)
 This subroutine does 1-d flux-form advection in the zonal direction using a monotonic piecewise linear scheme. More...
 
subroutine advect_y (Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, is, ie, js, je, k, G, GV, usePPM, useHuynh)
 This subroutine does 1-d flux-form advection using a monotonic piecewise linear scheme. More...
 
subroutine, public tracer_advect_init (Time, G, param_file, diag, CS)
 Initialize lateral tracer advection module. More...
 
subroutine, public tracer_advect_end (CS)
 Close the tracer advection module. More...
 

Variables

integer id_clock_advect
 
integer id_clock_pass
 
integer id_clock_sync
 

Function/Subroutine Documentation

◆ advect_tracer()

subroutine, public mom_tracer_advect::advect_tracer ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h_end,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  uhtr,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  vhtr,
type(ocean_obc_type), pointer  OBC,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(tracer_advect_cs), pointer  CS,
type(tracer_registry_type), pointer  Reg,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), optional  h_prev_opt,
integer, optional  max_iter_in,
logical, optional  x_first_in,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(out), optional  uhr_out,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(out), optional  vhr_out,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), optional  h_out 
)

This routine time steps the tracer concentration using a monotonic, conservative, weakly diffusive scheme.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]h_endlayer thickness after advection (m or kg m-2)
[in]uhtraccumulated volume/mass flux through zonal face (m3 or kg)
[in]vhtraccumulated volume/mass flux through merid face (m3 or kg)
obcspecifies whether, where, and what OBCs are used
[in]dttime increment (seconds)
cscontrol structure for module
regpointer to tracer registry
h_prev_optlayer thickness before advection (m or kg m-2)
[out]uhr_outaccumulated volume/mass flux through zonal face (m3 or kg)
[out]vhr_outaccumulated volume/mass flux through merid face (m3 or kg)
h_outlayer thickness before advection (m or kg m-2)

Definition at line 49 of file MOM_tracer_advect.F90.

References advect_x(), advect_y(), mom_domains::do_group_pass(), id_clock_advect, id_clock_pass, id_clock_sync, and mom_error_handler::mom_error().

Referenced by mom_offline_main::offline_advection_ale(), and mom_offline_main::offline_redistribute_residual().

49  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
50  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
51  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_end !< layer thickness after advection (m or kg m-2)
52  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhtr !< accumulated volume/mass flux through zonal face (m3 or kg)
53  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhtr !< accumulated volume/mass flux through merid face (m3 or kg)
54  type(ocean_obc_type), pointer :: obc !< specifies whether, where, and what OBCs are used
55  real, intent(in) :: dt !< time increment (seconds)
56  type(tracer_advect_cs), pointer :: cs !< control structure for module
57  type(tracer_registry_type), pointer :: reg !< pointer to tracer registry
58  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional :: h_prev_opt !< layer thickness before advection (m or kg m-2)
59  integer, optional :: max_iter_in
60  logical, optional :: x_first_in
61  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: uhr_out !< accumulated volume/mass flux through zonal face (m3 or kg)
62  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), optional, intent(out) :: vhr_out !< accumulated volume/mass flux through merid face (m3 or kg)
63  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional :: h_out !< layer thickness before advection (m or kg m-2)
64 
65  type(tracer_type) :: tr(max_fields_) ! The array of registered tracers
66  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
67  hprev ! cell volume at the end of previous tracer change (m3)
68  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
69  uhr ! The remaining zonal thickness flux (m3)
70  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
71  vhr ! The remaining meridional thickness fluxes (m3)
72  real :: uh_neglect(szib_(g),szj_(g)) ! uh_neglect and vh_neglect are the
73  real :: vh_neglect(szi_(g),szjb_(g)) ! magnitude of remaining transports that
74  ! can be simply discarded (m3 or kg).
75 
76  real :: landvolfill ! An arbitrary? nonzero cell volume, m3.
77  real :: idt ! 1/dt in s-1.
78  logical :: domore_u(szj_(g),szk_(g)) ! domore__ indicate whether there is more
79  logical :: domore_v(szjb_(g),szk_(g)) ! advection to be done in the corresponding
80  ! row or column.
81  logical :: x_first ! If true, advect in the x-direction first.
82  integer :: max_iter ! maximum number of iterations in each layer
83  integer :: domore_k(szk_(g))
84  integer :: stencil ! stencil of the advection scheme
85  integer :: nsten_halo ! number of stencils that fit in the halos
86  integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any
87  integer :: isv, iev, jsv, jev ! The valid range of the indices.
88  integer :: isdb, iedb, jsdb, jedb
89 
90  domore_u(:,:) = .false.
91  domore_v(:,:) = .false.
92  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
93  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
94  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
95  landvolfill = 1.0e-20 ! This is arbitrary, but must be positive.
96  stencil = 2 ! The scheme's stencil; 2 for PLM and PPM:H3
97 
98  if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_advect: "// &
99  "tracer_advect_init must be called before advect_tracer.")
100  if (.not. associated(reg)) call mom_error(fatal, "MOM_tracer_advect: "// &
101  "register_tracer must be called before advect_tracer.")
102  if (reg%ntr==0) return
103  call cpu_clock_begin(id_clock_advect)
104  x_first = (mod(g%first_direction,2) == 0)
105 
106  ! increase stencil size for Colella & Woodward PPM
107  if (cs%usePPM .and. .not. cs%useHuynh) stencil = 3
108 
109  ntr = reg%ntr
110  do m=1,ntr ; tr(m) = reg%Tr(m) ; enddo
111  idt = 1.0/dt
112 
113  max_iter = 2*int(ceiling(dt/cs%dt)) + 1
114 
115  if(present(max_iter_in)) max_iter = max_iter_in
116  if(present(x_first_in)) x_first = x_first_in
117  call cpu_clock_begin(id_clock_pass)
118  call create_group_pass(cs%pass_uhr_vhr_t_hprev, uhr, vhr, g%Domain)
119  call create_group_pass(cs%pass_uhr_vhr_t_hprev, hprev, g%Domain)
120  do m=1,ntr
121  call create_group_pass(cs%pass_uhr_vhr_t_hprev, tr(m)%t, g%Domain)
122  enddo
123  call cpu_clock_end(id_clock_pass)
124 
125 !$OMP parallel default(none) shared(nz,jsd,jed,IsdB,IedB,uhr,jsdB,jedB,Isd,Ied,vhr, &
126 !$OMP hprev,domore_k,js,je,is,ie,uhtr,vhtr,G,GV,h_end,&
127 !$OMP uh_neglect,vh_neglect,ntr,Tr,h_prev_opt)
128 
129 ! This initializes the halos of uhr and vhr because pass_vector might do
130 ! calculations on them, even though they are never used.
131 !$OMP do
132 
133  do k = 1, nz
134  do j = jsd, jed; do i = isdb, iedb; uhr(i,j,k) = 0.0; enddo ; enddo
135  do j = jsdb, jedb; do i = isd, ied; vhr(i,j,k) = 0.0; enddo ; enddo
136  do j = jsd, jed; do i = isd, ied; hprev(i,j,k) = 0.0; enddo ; enddo
137  domore_k(k)=1
138 ! Put the remaining (total) thickness fluxes into uhr and vhr.
139  do j=js,je ; do i=is-1,ie ; uhr(i,j,k) = uhtr(i,j,k) ; enddo ; enddo
140  do j=js-1,je ; do i=is,ie ; vhr(i,j,k) = vhtr(i,j,k) ; enddo ; enddo
141  if (.not. present(h_prev_opt)) then
142  ! This loop reconstructs the thickness field the last time that the
143  ! tracers were updated, probably just after the diabatic forcing. A useful
144  ! diagnostic could be to compare this reconstruction with that older value.
145  do i=is,ie ; do j=js,je
146  hprev(i,j,k) = max(0.0, g%areaT(i,j)*h_end(i,j,k) + &
147  ((uhr(i,j,k) - uhr(i-1,j,k)) + (vhr(i,j,k) - vhr(i,j-1,k))))
148  ! In the case that the layer is now dramatically thinner than it was previously,
149  ! add a bit of mass to avoid truncation errors. This will lead to
150  ! non-conservation of tracers
151  hprev(i,j,k) = hprev(i,j,k) + &
152  max(0.0, 1.0e-13*hprev(i,j,k) - g%areaT(i,j)*h_end(i,j,k))
153  enddo ; enddo
154  else
155  do i=is,ie ; do j=js,je
156  hprev(i,j,k) = h_prev_opt(i,j,k);
157  enddo ; enddo
158  endif
159  enddo
160 
161 
162 !$OMP do
163  do j=jsd,jed ; do i=isd,ied-1
164  uh_neglect(i,j) = gv%H_subroundoff*min(g%areaT(i,j),g%areaT(i+1,j))
165  enddo ; enddo
166 !$OMP do
167  do j=jsd,jed-1 ; do i=isd,ied
168  vh_neglect(i,j) = gv%H_subroundoff*min(g%areaT(i,j),g%areaT(i,j+1))
169  enddo ; enddo
170 
171 !$OMP do
172  ! initialize diagnostic fluxes and tendencies
173  do m=1,ntr
174  if (associated(tr(m)%ad_x)) then
175  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
176  tr(m)%ad_x(i,j,k) = 0.0
177  enddo ; enddo ; enddo
178  endif
179  if (associated(tr(m)%ad_y)) then
180  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
181  tr(m)%ad_y(i,j,k) = 0.0
182  enddo ; enddo ; enddo
183  endif
184  if (associated(tr(m)%advection_xy)) then
185  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
186  tr(m)%advection_xy(i,j,k) = 0.0
187  enddo ; enddo ; enddo
188  endif
189  if (associated(tr(m)%ad2d_x)) then
190  do j=jsd,jed ; do i=isd,ied ; tr(m)%ad2d_x(i,j) = 0.0 ; enddo ; enddo
191  endif
192  if (associated(tr(m)%ad2d_y)) then
193  do j=jsd,jed ; do i=isd,ied ; tr(m)%ad2d_y(i,j) = 0.0 ; enddo ; enddo
194  endif
195  enddo
196 !$OMP end parallel
197 
198  isv = is ; iev = ie ; jsv = js ; jev = je
199 
200  do itt=1,max_iter
201 
202  if (isv > is-stencil) then
203  call cpu_clock_begin(id_clock_pass)
204  call do_group_pass(cs%pass_uhr_vhr_t_hprev, g%Domain)
205  call cpu_clock_end(id_clock_pass)
206 
207  nsten_halo = min(is-isd,ied-ie,js-jsd,jed-je)/stencil
208  isv = is-nsten_halo*stencil ; jsv = js-nsten_halo*stencil
209  iev = ie+nsten_halo*stencil ; jev = je+nsten_halo*stencil
210  ! Reevaluate domore_u & domore_v unless the valid range is the same size as
211  ! before. Also, do this if there is Strang splitting.
212  if ((nsten_halo > 1) .or. (itt==1)) then
213 !$OMP parallel do default(none) shared(nz,domore_k,jsv,jev,domore_u,isv,iev,stencil, &
214 !$OMP uhr,domore_v,vhr)
215  do k=1,nz ; if (domore_k(k) > 0) then
216  do j=jsv,jev ; if (.not.domore_u(j,k)) then
217  do i=isv+stencil-1,iev-stencil; if (uhr(i,j,k) /= 0.0) then
218  domore_u(j,k) = .true. ; exit
219  endif ; enddo ! i-loop
220  endif ; enddo
221  do j=jsv+stencil-1,jev-stencil ; if (.not.domore_v(j,k)) then
222  do i=isv+stencil,iev-stencil; if (vhr(i,j,k) /= 0.0) then
223  domore_v(j,k) = .true. ; exit
224  endif ; enddo ! i-loop
225  endif ; enddo
226 
227  ! At this point, domore_k is global. Change it so that it indicates
228  ! whether any work is needed on a layer on this processor.
229  domore_k(k) = 0
230  do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
231  do j=jsv+stencil-1,jev-stencil ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
232 
233  endif ; enddo ! k-loop
234  endif
235  endif
236 
237  ! Set the range of valid points after this iteration.
238  isv = isv + stencil ; iev = iev - stencil
239  jsv = jsv + stencil ; jev = jev - stencil
240 
241 !$OMP parallel do default(none) shared(nz,domore_k,x_first,Tr,hprev,uhr,uh_neglect, &
242 !$OMP OBC,domore_u,ntr,Idt,isv,iev,jsv,jev,stencil, &
243 !$OMP G,GV,CS,vhr,vh_neglect,domore_v)
244 
245  ! To ensure positive definiteness of the thickness at each iteration, the
246  ! mass fluxes out of each layer are checked each step, and limited to keep
247  ! the thicknesses positive. This means that several iterations may be required
248  ! for all the transport to happen. The sum over domore_k keeps the processors
249  ! synchronized. This may not be very efficient, but it should be reliable.
250  do k=1,nz ; if (domore_k(k) > 0) then
251 
252  if (x_first) then
253 
254  ! First, advect zonally.
255  call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
256  isv, iev, jsv-stencil, jev+stencil, k, g, gv, cs%usePPM, cs%useHuynh)
257 
258  ! Next, advect meridionally.
259  call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
260  isv, iev, jsv, jev, k, g, gv, cs%usePPM, cs%useHuynh)
261 
262  domore_k(k) = 0
263  do j=jsv-stencil,jev+stencil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
264  do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
265 
266  else
267 
268  ! First, advect meridionally.
269  call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
270  isv-stencil, iev+stencil, jsv, jev, k, g, gv, cs%usePPM, cs%useHuynh)
271 
272  ! Next, advect zonally.
273  call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
274  isv, iev, jsv, jev, k, g, gv, cs%usePPM, cs%useHuynh)
275 
276  domore_k(k) = 0
277  do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
278  do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
279 
280  endif
281 
282 
283  endif ; enddo ! End of k-loop
284 
285  ! If the advection just isn't finishing after max_iter, move on.
286  if (itt >= max_iter) then
287  exit
288  endif
289 
290  ! Exit if there are no layers that need more iterations.
291  if (isv > is-stencil) then
292  do_any = 0
293  call cpu_clock_begin(id_clock_sync)
294  call sum_across_pes(domore_k(:), nz)
295  call cpu_clock_end(id_clock_sync)
296  do k=1,nz ; do_any = do_any + domore_k(k) ; enddo
297  if (do_any == 0) then
298  exit
299  endif
300 
301  endif
302 
303  enddo ! Iterations loop
304 
305  if(present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:)
306  if(present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:)
307  if(present(h_out)) h_out(:,:,:) = hprev(:,:,:)
308 
309  call cpu_clock_end(id_clock_advect)
310 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ advect_x()

subroutine mom_tracer_advect::advect_x ( type(tracer_type), dimension(ntr), intent(inout)  Tr,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  hprev,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)  uhr,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(inout)  uh_neglect,
type(ocean_obc_type), pointer  OBC,
logical, dimension( g %jsd: g %jed, g %ke), intent(inout)  domore_u,
integer, intent(in)  ntr,
real, intent(in)  Idt,
integer, intent(in)  is,
integer, intent(in)  ie,
integer, intent(in)  js,
integer, intent(in)  je,
integer, intent(in)  k,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
logical, intent(in)  usePPM,
logical, intent(in)  useHuynh 
)
private

This subroutine does 1-d flux-form advection in the zonal direction using a monotonic piecewise linear scheme.

Parameters
[in,out]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure

Definition at line 318 of file MOM_tracer_advect.F90.

Referenced by advect_tracer().

318  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
319  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
320  type(tracer_type), dimension(ntr), intent(inout) :: tr
321  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hprev
322  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhr
323  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: uh_neglect
324  type(ocean_obc_type), pointer :: obc
325  logical, dimension(SZJ_(G),SZK_(G)), intent(inout) :: domore_u
326  real, intent(in) :: idt
327  integer, intent(in) :: ntr, is, ie, js, je,k
328  logical, intent(in) :: useppm, usehuynh
329 
330  real, dimension(SZI_(G),ntr) :: &
331  slope_x ! The concentration slope per grid point in units of
332  ! concentration (nondim.).
333  real, dimension(SZIB_(G),ntr) :: &
334  flux_x ! The tracer flux across a boundary in m3*conc or kg*conc.
335  real :: maxslope ! The maximum concentration slope per grid point
336  ! consistent with monotonicity, in conc. (nondim.).
337  real :: hup, hlos ! hup is the upwind volume, hlos is the
338  ! part of that volume that might be lost
339  ! due to advection out the other side of
340  ! the grid box, both in m3 or kg.
341  real :: uhh(szib_(g)) ! The zonal flux that occurs during the
342  ! current iteration, in m3 or kg.
343  real, dimension(SZIB_(G)) :: &
344  hlst, ihnew, & ! Work variables with units of m3 or kg and m-3 or kg-1.
345  cfl ! A nondimensional work variable.
346  real :: min_h ! The minimum thickness that can be realized during
347  ! any of the passes, in m or kg m-2.
348  real :: h_neglect ! A thickness that is so small it is usually lost
349  ! in roundoff and can be neglected, in m.
350  logical :: do_i(szib_(g)) ! If true, work on given points.
351  logical :: do_any_i
352  integer :: i, j, m, n, i_up, stencil
353  real :: ar, al, dmx, dmn, tp, tc, tm, da, ma, a6
354  logical :: useplmslope
355 
356  useplmslope = .not. (useppm .and. usehuynh)
357  ! stencil for calculating slope values
358  stencil = 1
359  if (useppm .and. .not. usehuynh) stencil = 2
360 
361  min_h = 0.1*gv%Angstrom
362  h_neglect = gv%H_subroundoff
363 
364 ! do I=is-1,ie ; ts2(I) = 0.0 ; enddo
365  do i=is-1,ie ; cfl(i) = 0.0 ; enddo
366 
367  do j=js,je ; if (domore_u(j,k)) then
368  domore_u(j,k) = .false.
369 
370  ! Calculate the i-direction profiles (slopes) of each tracer that
371  ! is being advected.
372  if (useplmslope) then
373  do m=1,ntr ; do i=is-stencil,ie+stencil
374  !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < &
375  ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))) then
376  ! maxslope = 4.0*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k))
377  !else
378  ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))
379  !endif
380  !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
381  ! slope_x(i,m) = 0.0
382  !elseif (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))<ABS(maxslope)) then
383  ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * &
384  ! 0.5*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))
385  !else
386  ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * 0.5*maxslope
387  !endif
388  tp = tr(m)%t(i+1,j,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i-1,j,k)
389  dmx = max( tp, tc, tm ) - tc
390  dmn= tc - min( tp, tc, tm )
391  slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
392  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
393  enddo ; enddo
394  endif ! usePLMslope
395 
396  ! Calculate the i-direction fluxes of each tracer, using as much
397  ! the minimum of the remaining mass flux (uhr) and the half the mass
398  ! in the cell plus whatever part of its half of the mass flux that
399  ! the flux through the other side does not require.
400  do i=is-1,ie
401  if (uhr(i,j,k) == 0.0) then
402  uhh(i) = 0.0
403  cfl(i) = 0.0
404  elseif (uhr(i,j,k) < 0.0) then
405  hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
406  hlos = max(0.0,uhr(i+1,j,k))
407  if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
408  ((0.5*hup + uhr(i,j,k)) < 0.0)) then
409  uhh(i) = min(-0.5*hup,-hup+hlos,0.0)
410  domore_u(j,k) = .true.
411  else
412  uhh(i) = uhr(i,j,k)
413  endif
414  !ts2(I) = 0.5*(1.0 + uhh(I)/(hprev(i+1,j,k)+h_neglect))
415  cfl(i) = - uhh(i)/(hprev(i+1,j,k)+h_neglect) ! CFL is positive
416  else
417  hup = hprev(i,j,k) - g%areaT(i,j)*min_h
418  hlos = max(0.0,-uhr(i-1,j,k))
419  if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
420  ((0.5*hup - uhr(i,j,k)) < 0.0)) then
421  uhh(i) = max(0.5*hup,hup-hlos,0.0)
422  domore_u(j,k) = .true.
423  else
424  uhh(i) = uhr(i,j,k)
425  endif
426  !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
427  cfl(i) = uhh(i)/(hprev(i,j,k)+h_neglect) ! CFL is positive
428  endif
429  enddo
430 
431 
432  if (useppm) then
433  do m=1,ntr ; do i=is-1,ie
434  ! centre cell depending on upstream direction
435  if (uhh(i) >= 0.0) then
436  i_up = i
437  else
438  i_up = i+1
439  endif
440 
441  ! Implementation of PPM-H3
442  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)
443 
444  if (usehuynh) then
445  al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
446  al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
447  ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
448  ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
449  else
450  al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
451  ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
452  endif
453 
454  da = ar - al ; ma = 0.5*( ar + al )
455  if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.) then
456  al = tc ; ar = tc ! PCM for local extremum and bounadry cells
457  elseif ( da*(tc-ma) > (da*da)/6. ) then
458  al = 3.*tc - 2.*ar
459  elseif ( da*(tc-ma) < - (da*da)/6. ) then
460  ar = 3.*tc - 2.*al
461  endif
462 
463  a6 = 6.*tc - 3. * (ar + al) ! Curvature
464 
465  if (uhh(i) >= 0.0) then
466  flux_x(i,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
467  ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
468  else
469  flux_x(i,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
470  ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
471  endif
472  enddo ; enddo
473  else ! PLM
474  do m=1,ntr ; do i=is-1,ie
475  if (uhh(i) >= 0.0) then
476  ! Indirect implementation of PLM
477  !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m)
478  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
479  !flux_x(I,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) )
480  ! Alternative implementation of PLM
481  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
482  !flux_x(I,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) )
483  ! Alternative implementation of PLM
484  tc = tr(m)%t(i,j,k)
485  flux_x(i,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
486  ! Original implementation of PLM
487  !flux_x(I,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I))
488  else
489  ! Indirect implementation of PLM
490  !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
491  !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m)
492  !flux_x(I,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) )
493  ! Alternative implementation of PLM
494  !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
495  !flux_x(I,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) )
496  ! Alternative implementation of PLM
497  tc = tr(m)%t(i+1,j,k)
498  flux_x(i,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
499  ! Original implementation of PLM
500  !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I))
501  endif
502  !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
503  enddo ; enddo
504  endif ! usePPM
505 
506  if (associated(obc)) then ; if (obc%OBC_pe) then ; if (obc%specified_u_BCs_exist_globally) then
507  do n=1,obc%number_of_segments
508  if (obc%segment(n)%is_E_or_W) then
509  if (j >= obc%segment(n)%HI%jsd .and. j<= obc%segment(n)%HI%jed) then
510  i = obc%segment(n)%HI%IsdB
511  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
512  if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
513  (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5)) then
514  uhh(i) = uhr(i,j,k)
515  do m=1,ntr
516  if (associated(tr(m)%OBC_in_u)) then
517  flux_x(i,m) = uhh(i)*tr(m)%OBC_in_u(i,j,k)
518  else ; flux_x(i,m) = uhh(i)*tr(m)%OBC_inflow_conc ; endif
519  enddo
520  endif
521  endif
522  endif
523  enddo
524  endif ; endif ; endif
525 
526  ! Calculate new tracer concentration in each cell after accounting
527  ! for the i-direction fluxes.
528  do i=is-1,ie
529  uhr(i,j,k) = uhr(i,j,k) - uhh(i)
530  if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
531  enddo
532  do i=is,ie
533  if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0)) then
534  do_i(i) = .true.
535  hlst(i) = hprev(i,j,k)
536  hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
537  if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false.
538  elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
539  hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
540  ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
541  else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
542  else
543  do_i(i) = .false.
544  endif
545  enddo
546 
547  ! update tracer concentration from i-flux and save some diagnostics
548  do m=1,ntr
549 
550  ! update tracer
551  do i=is,ie ; if ((do_i(i)) .and. (ihnew(i) > 0.0)) then
552  tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
553  (flux_x(i,m) - flux_x(i-1,m))) * ihnew(i)
554  endif ; enddo
555 
556  ! diagnostics
557  if (associated(tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i)) then
558  tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,m)*idt
559  endif ; enddo ; endif
560  if (associated(tr(m)%ad2d_x)) then ; do i=is,ie ; if (do_i(i)) then
561  tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,m)*idt
562  endif ; enddo ; endif
563 
564  ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic).
565  ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
566  if (associated(tr(m)%advection_xy)) then
567  do i=is,ie ; if (do_i(i)) then
568  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)
569  endif ; enddo
570  endif
571 
572  enddo
573 
574  endif ; enddo ! End of j-loop.
575 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the caller graph for this function:

◆ advect_y()

subroutine mom_tracer_advect::advect_y ( type(tracer_type), dimension(ntr), intent(inout)  Tr,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  hprev,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)  vhr,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(inout)  vh_neglect,
type(ocean_obc_type), pointer  OBC,
logical, dimension( g %jsdb: g %jedb, g %ke), intent(inout)  domore_v,
integer, intent(in)  ntr,
real, intent(in)  Idt,
integer, intent(in)  is,
integer, intent(in)  ie,
integer, intent(in)  js,
integer, intent(in)  je,
integer, intent(in)  k,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
logical, intent(in)  usePPM,
logical, intent(in)  useHuynh 
)
private

This subroutine does 1-d flux-form advection using a monotonic piecewise linear scheme.

Parameters
[in,out]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure

Definition at line 582 of file MOM_tracer_advect.F90.

Referenced by advect_tracer().

582  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
583  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
584  type(tracer_type), dimension(ntr), intent(inout) :: tr
585  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hprev
586  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhr
587  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect
588  type(ocean_obc_type), pointer :: obc
589  logical, dimension(SZJB_(G),SZK_(G)), intent(inout) :: domore_v
590  real, intent(in) :: idt
591  integer, intent(in) :: ntr, is, ie, js, je,k
592  logical, intent(in) :: useppm, usehuynh
593 
594  real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
595  slope_y ! The concentration slope per grid point in units of
596  ! concentration (nondim.).
597  real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
598  flux_y ! The tracer flux across a boundary in m3 * conc or kg*conc.
599  real :: maxslope ! The maximum concentration slope per grid point
600  ! consistent with monotonicity, in conc. (nondim.).
601  real :: vhh(szi_(g),szjb_(g)) ! The meridional flux that occurs during the
602  ! current iteration, in m3 or kg.
603  real :: hup, hlos ! hup is the upwind volume, hlos is the
604  ! part of that volume that might be lost
605  ! due to advection out the other side of
606  ! the grid box, both in m3 or kg.
607  real, dimension(SZIB_(G)) :: &
608  hlst, ihnew, & ! Work variables with units of m3 or kg and m-3 or kg-1.
609  cfl ! A nondimensional work variable.
610  real :: min_h ! The minimum thickness that can be realized during
611  ! any of the passes, in m or kg m-2.
612  real :: h_neglect ! A thickness that is so small it is usually lost
613  ! in roundoff and can be neglected, in m.
614  logical :: do_j_tr(szj_(g)) ! If true, calculate the tracer profiles.
615  logical :: do_i(szib_(g)) ! If true, work on given points.
616  logical :: do_any_i
617  integer :: i, j, j2, m, n, j_up, stencil
618  real :: ar, al, dmx, dmn, tp, tc, tm, da, ma, a6
619  logical :: useplmslope
620 
621  useplmslope = .not. (useppm .and. usehuynh)
622  ! stencil for calculating slope values
623  stencil = 1
624  if (useppm .and. .not. usehuynh) stencil = 2
625 
626  min_h = 0.1*gv%Angstrom
627  h_neglect = gv%H_subroundoff
628 
629  !do i=is,ie ; ts2(i) = 0.0 ; enddo
630 
631  ! We conditionally perform work on tracer points: calculating the PLM slope,
632  ! and updating tracer concentration within a cell
633  ! this depends on whether there is a flux which would affect this tracer point,
634  ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires
635  ! slope calculations at the two tracer points on either side (as indicated by
636  ! the stencil variable), so we account for this with the do_j_tr flag array
637  !
638  ! Note: this does lead to unnecessary work in updating tracer concentrations,
639  ! since that doesn't need a wider stencil with the PPM advection scheme, but
640  ! this would require an additional loop, etc.
641  do_j_tr(:) = .false.
642  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
643 
644  ! Calculate the j-direction profiles (slopes) of each tracer that
645  ! is being advected.
646  if (useplmslope) then
647  do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr ; do i=is,ie
648  !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < &
649  ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then
650  ! maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))
651  !else
652  ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))
653  !endif
654  !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
655  ! slope_y(i,m,j) = 0.0
656  !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then
657  ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * &
658  ! 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))
659  !else
660  ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope
661  !endif
662  tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
663  dmx = max( tp, tc, tm ) - tc
664  dmn= tc - min( tp, tc, tm )
665  slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
666  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
667  enddo ; enddo ; endif ; enddo ! End of i-, m-, & j- loops.
668  endif ! usePLMslope
669 
670  ! Calculate the j-direction fluxes of each tracer, using as much
671  ! the minimum of the remaining mass flux (vhr) and the half the mass
672  ! in the cell plus whatever part of its half of the mass flux that
673  ! the flux through the other side does not require.
674  do j=js-1,je ; if (domore_v(j,k)) then
675  domore_v(j,k) = .false.
676 
677  do i=is,ie
678  if (vhr(i,j,k) == 0.0) then
679  vhh(i,j) = 0.0
680  cfl(i) = 0.0
681  elseif (vhr(i,j,k) < 0.0) then
682  hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
683  hlos = max(0.0,vhr(i,j+1,k))
684  if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
685  ((0.5*hup + vhr(i,j,k)) < 0.0)) then
686  vhh(i,j) = min(-0.5*hup,-hup+hlos,0.0)
687  domore_v(j,k) = .true.
688  else
689  vhh(i,j) = vhr(i,j,k)
690  endif
691  !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k)+h_neglect))
692  cfl(i) = - vhh(i,j) / (hprev(i,j+1,k)+h_neglect) ! CFL is positive
693  else
694  hup = hprev(i,j,k) - g%areaT(i,j)*min_h
695  hlos = max(0.0,-vhr(i,j-1,k))
696  if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
697  ((0.5*hup - vhr(i,j,k)) < 0.0)) then
698  vhh(i,j) = max(0.5*hup,hup-hlos,0.0)
699  domore_v(j,k) = .true.
700  else
701  vhh(i,j) = vhr(i,j,k)
702  endif
703  !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k)+h_neglect))
704  cfl(i) = vhh(i,j) / (hprev(i,j,k)+h_neglect) ! CFL is positive
705  endif
706  enddo
707 
708  if (useppm) then
709  do m=1,ntr ; do i=is,ie
710  ! centre cell depending on upstream direction
711  if (vhh(i,j) >= 0.0) then
712  j_up = j
713  else
714  j_up = j + 1
715  endif
716 
717  ! Implementation of PPM-H3
718  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)
719 
720  if (usehuynh) then
721  al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
722  al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
723  ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
724  ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
725  else
726  al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
727  ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
728  endif
729 
730  da = ar - al ; ma = 0.5*( ar + al )
731  if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.) then
732  al = tc ; ar = tc ! PCM for local extremum and bounadry cells
733  elseif ( da*(tc-ma) > (da*da)/6. ) then
734  al = 3.*tc - 2.*ar
735  elseif ( da*(tc-ma) < - (da*da)/6. ) then
736  ar = 3.*tc - 2.*al
737  endif
738 
739  a6 = 6.*tc - 3. * (ar + al) ! Curvature
740 
741  if (vhh(i,j) >= 0.0) then
742  flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
743  ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
744  else
745  flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
746  ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
747  endif
748  enddo ; enddo
749  else ! PLM
750  do m=1,ntr ; do i=is,ie
751  if (vhh(i,j) >= 0.0) then
752  ! Indirect implementation of PLM
753  !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j)
754  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
755  !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) )
756  ! Alternative implementation of PLM
757  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
758  !flux_y(i,m,J) = vhh(i,J)*(aR - 0.5 * slope_y(i,m,j)*CFL(i))
759  ! Alternative implementation of PLM
760  tc = tr(m)%t(i,j,k)
761  flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
762  ! Original implementation of PLM
763  !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j,k) + slope_y(i,m,j)*ts2(i))
764  else
765  ! Indirect implementation of PLM
766  !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
767  !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1)
768  !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) )
769  ! Alternative implementation of PLM
770  !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
771  !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * slope_y(i,m,j+1)*CFL(i) )
772  ! Alternative implementation of PLM
773  tc = tr(m)%t(i,j+1,k)
774  flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
775  ! Original implementation of PLM
776  !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j+1,k) - slope_y(i,m,j+1)*ts2(i))
777  endif
778  enddo ; enddo
779  endif ! usePPM
780 
781  if (associated(obc)) then ; if (obc%OBC_pe) then ; if (obc%specified_v_BCs_exist_globally) then
782  do n=1,obc%number_of_segments
783  if (obc%segment(n)%is_N_or_S) then
784  if (j >= obc%segment(n)%HI%JsdB .and. j<= obc%segment(n)%HI%JedB) then
785  i = obc%segment(n)%HI%isd
786  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
787  if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
788  (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5)) then
789  vhh(i,j) = vhr(i,j,k)
790  do m=1,ntr
791  if (associated(tr(m)%OBC_in_v)) then
792  flux_y(i,m,j) = vhh(i,j)*tr(m)%OBC_in_v(i,j,k)
793  else ; flux_y(i,m,j) = vhh(i,j)*tr(m)%OBC_inflow_conc ; endif
794  enddo
795  endif
796  endif
797  endif
798  enddo
799  endif ; endif ; endif
800  else ! not domore_v.
801  do i=is,ie ; vhh(i,j) = 0.0 ; enddo
802  do m=1,ntr ; do i=is,ie ; flux_y(i,m,j) = 0.0 ; enddo ; enddo
803  endif ; enddo ! End of j-loop
804 
805  do j=js-1,je ; do i=is,ie
806  vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
807  if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
808  enddo ; enddo
809 
810  ! Calculate new tracer concentration in each cell after accounting
811  ! for the j-direction fluxes.
812  do j=js,je ; if (do_j_tr(j)) then
813  do i=is,ie
814  if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0)) then
815  do_i(i) = .true.
816  hlst(i) = hprev(i,j,k)
817  hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
818  if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false.
819  elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
820  hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
821  ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
822  else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
823  else ; do_i(i) = .false. ; endif
824  enddo
825 
826  ! update tracer and save some diagnostics
827  do m=1,ntr
828  do i=is,ie ; if (do_i(i)) then
829  tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
830  (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
831  endif ; enddo
832 
833  ! diagnostics
834  if (associated(tr(m)%ad_y)) then ; do i=is,ie ; if (do_i(i)) then
835  tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
836  endif ; enddo ; endif
837  if (associated(tr(m)%ad2d_y)) then ; do i=is,ie ; if (do_i(i)) then
838  tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
839  endif ; enddo ; endif
840 
841  ! diagnose convergence of flux_y and add to convergence of flux_x.
842  ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
843  if (associated(tr(m)%advection_xy)) then
844  do i=is,ie ; if (do_i(i)) then
845  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)
846  endif ; enddo
847  endif
848 
849 
850  enddo
851  endif ; enddo ! End of j-loop.
852 
853 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the caller graph for this function:

◆ tracer_advect_end()

subroutine, public mom_tracer_advect::tracer_advect_end ( type(tracer_advect_cs), pointer  CS)

Close the tracer advection module.

Definition at line 912 of file MOM_tracer_advect.F90.

912  type(tracer_advect_cs), pointer :: cs
913 
914  if (associated(cs)) deallocate(cs)
915 

◆ tracer_advect_init()

subroutine, public mom_tracer_advect::tracer_advect_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(tracer_advect_cs), pointer  CS 
)

Initialize lateral tracer advection module.

Parameters
[in]timecurrent model time
[in]gocean grid structure
[in]param_fileopen file to parse for model parameters
[in,out]diagregulates diagnostic output
csmodule control structure

Definition at line 858 of file MOM_tracer_advect.F90.

References id_clock_advect, id_clock_pass, id_clock_sync, and mom_error_handler::mom_error().

858  type(time_type), target, intent(in) :: time !< current model time
859  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
860  type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
861  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
862  type(tracer_advect_cs), pointer :: cs !< module control structure
863 
864  integer, save :: init_calls = 0
865 
866 ! This include declares and sets the variable "version".
867 #include "version_variable.h"
868  character(len=40) :: mdl = "MOM_tracer_advect" ! This module's name.
869  character(len=256) :: mesg ! Message for error messages.
870 
871  if (associated(cs)) then
872  call mom_error(warning, "tracer_advect_init called with associated control structure.")
873  return
874  endif
875  allocate(cs)
876 
877  cs%diag => diag
878 
879  ! Read all relevant parameters and write them to the model log.
880  call log_version(param_file, mdl, version, "")
881  call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
882  desc="The (baroclinic) dynamics time step.", units="s")
883  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
884  call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, &
885  desc="The horizontal transport scheme for tracers:\n"//&
886  " PLM - Piecewise Linear Method\n"//&
887  " PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// &
888  " PPM - Piecewise Parabolic Method (Colella-Woodward)" &
889  , default='PLM')
890  select case (trim(mesg))
891  case ("PLM")
892  cs%usePPM = .false.
893  case ("PPM:H3")
894  cs%usePPM = .true.
895  cs%useHuynh = .true.
896  case ("PPM")
897  cs%usePPM = .true.
898  cs%useHuynh = .false.
899  case default
900  call mom_error(fatal, "MOM_tracer_advect, tracer_advect_init: "//&
901  "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
902  end select
903 
904  id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=clock_module)
905  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
906  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
907 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:

Variable Documentation

◆ id_clock_advect

integer mom_tracer_advect::id_clock_advect

Definition at line 39 of file MOM_tracer_advect.F90.

Referenced by advect_tracer(), and tracer_advect_init().

39 integer :: id_clock_advect

◆ id_clock_pass

integer mom_tracer_advect::id_clock_pass
private

Definition at line 40 of file MOM_tracer_advect.F90.

Referenced by advect_tracer(), and tracer_advect_init().

40 integer :: id_clock_pass

◆ id_clock_sync

integer mom_tracer_advect::id_clock_sync
private

Definition at line 41 of file MOM_tracer_advect.F90.

Referenced by advect_tracer(), and tracer_advect_init().

41 integer :: id_clock_sync