MOM6
Data Types
mom_tracer_advect Module Reference

Detailed Description

This module 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...
 
integer id_clock_advect
 CPU time clocks. More...
 
integer id_clock_pass
 CPU time clocks. More...
 
integer id_clock_sync
 CPU time clocks. More...
 
subroutine, public advect_tracer (h_end, uhtr, vhtr, OBC, dt, G, GV, US, 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, US, 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, US, usePPM, useHuynh)
 This subroutine does 1-d flux-form advection using a monotonic piecewise linear scheme. More...
 
subroutine, public tracer_advect_init (Time, G, US, param_file, diag, CS)
 Initialize lateral tracer advection module. More...
 
subroutine, public tracer_advect_end (CS)
 Close the tracer advection module. More...
 

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(unit_scale_type), intent(in)  US,
type(tracer_advect_cs), pointer  CS,
type(tracer_registry_type), pointer  Reg,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in), optional  h_prev_opt,
integer, intent(in), optional  max_iter_in,
logical, intent(in), 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), intent(out), 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 [H ~> m or kg m-2]
[in]uhtraccumulated volume/mass flux through zonal face [H L2 ~> m3 or kg]
[in]vhtraccumulated volume/mass flux through merid face [H L2 ~> m3 or kg]
obcspecifies whether, where, and what OBCs are used
[in]dttime increment [T ~> s]
[in]usA dimensional unit scaling type
cscontrol structure for module
regpointer to tracer registry
[in]h_prev_optlayer thickness before advection [H ~> m or kg m-2]
[in]max_iter_inThe maximum number of iterations
[in]x_first_inIf present, indicate whether to update first in the x- or y-direction.
[out]uhr_outaccumulated volume/mass flux through zonal face
[out]vhr_outaccumulated volume/mass flux through merid face
[out]h_outlayer thickness before advection [H ~> m or kg m-2]

Definition at line 52 of file MOM_tracer_advect.F90.

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

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(), mom_offline_main::offline_advection_layer(), and mom_offline_main::offline_redistribute_residual().

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,
type(unit_scale_type), intent(in)  US,
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
[in,out]trThe array of registered tracers to work on
[in,out]hprevcell volume at the end of previous tracer change [H L2 ~> m3 or kg]
[in,out]uhraccumulated volume/mass flux through the zonal face [H L2 ~> m3 or kg]
[in,out]uh_neglectA tiny zonal mass flux that can be neglected [H L2 ~> m3 or kg]
obcspecifies whether, where, and what OBCs are used
[in,out]domore_uIf true, there is more advection to be done in this u-row
[in]idtThe inverse of dt [T-1 ~> s-1]
[in]ntrThe number of tracers
[in]isThe starting tracer i-index to work on
[in]ieThe ending tracer i-index to work on
[in]jsThe starting tracer j-index to work on
[in]jeThe ending tracer j-index to work on
[in]kThe k-level to work on
[in]usA dimensional unit scaling type
[in]useppmIf true, use PPM instead of PLM
[in]usehuynhIf true, use the Huynh scheme for PPM interface values

Definition at line 330 of file MOM_tracer_advect.F90.

330  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
331  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
332  type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on
333  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hprev !< cell volume at the end of previous
334  !! tracer change [H L2 ~> m3 or kg]
335  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhr !< accumulated volume/mass flux through
336  !! the zonal face [H L2 ~> m3 or kg]
337  real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: uh_neglect !< A tiny zonal mass flux that can
338  !! be neglected [H L2 ~> m3 or kg]
339  type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used
340  logical, dimension(SZJ_(G),SZK_(G)), intent(inout) :: domore_u !< If true, there is more advection to be
341  !! done in this u-row
342  real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1]
343  integer, intent(in) :: ntr !< The number of tracers
344  integer, intent(in) :: is !< The starting tracer i-index to work on
345  integer, intent(in) :: ie !< The ending tracer i-index to work on
346  integer, intent(in) :: js !< The starting tracer j-index to work on
347  integer, intent(in) :: je !< The ending tracer j-index to work on
348  integer, intent(in) :: k !< The k-level to work on
349  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
350  logical, intent(in) :: usePPM !< If true, use PPM instead of PLM
351  logical, intent(in) :: useHuynh !< If true, use the Huynh scheme
352  !! for PPM interface values
353 
354  real, dimension(SZI_(G),ntr) :: &
355  slope_x ! The concentration slope per grid point [conc].
356  real, dimension(SZIB_(G),ntr) :: &
357  flux_x ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc].
358  real, dimension(SZI_(G),ntr) :: &
359  T_tmp ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].
360 
361  real :: maxslope ! The maximum concentration slope per grid point
362  ! consistent with monotonicity [conc].
363  real :: hup, hlos ! hup is the upwind volume, hlos is the
364  ! part of that volume that might be lost
365  ! due to advection out the other side of
366  ! the grid box, both in [H L2 ~> m3 or kg].
367  real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the
368  ! current iteration [H L2 ~> m3 or kg].
369  real, dimension(SZIB_(G)) :: &
370  hlst, & ! Work variable [H L2 ~> m3 or kg].
371  Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1].
372  CFL ! A nondimensional work variable [nondim].
373  real :: min_h ! The minimum thickness that can be realized during
374  ! any of the passes [H ~> m or kg m-2].
375  real :: h_neglect ! A thickness that is so small it is usually lost
376  ! in roundoff and can be neglected [H ~> m or kg m-2].
377  logical :: do_i(SZIB_(G)) ! If true, work on given points.
378  logical :: do_any_i
379  integer :: i, j, m, n, i_up, stencil
380  real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
381  real :: fac1,u_L_in,u_L_out ! terms used for time-stepping OBC reservoirs
382  type(OBC_segment_type), pointer :: segment=>null()
383  logical :: usePLMslope
384 
385  useplmslope = .not. (useppm .and. usehuynh)
386  ! stencil for calculating slope values
387  stencil = 1
388  if (useppm .and. .not. usehuynh) stencil = 2
389 
390  min_h = 0.1*gv%Angstrom_H
391  h_neglect = gv%H_subroundoff
392 
393 ! do I=is-1,ie ; ts2(I) = 0.0 ; enddo
394  do i=is-1,ie ; cfl(i) = 0.0 ; enddo
395 
396  do j=js,je ; if (domore_u(j,k)) then
397  domore_u(j,k) = .false.
398 
399  ! Calculate the i-direction profiles (slopes) of each tracer that
400  ! is being advected.
401  if (useplmslope) then
402  do m=1,ntr ; do i=is-stencil,ie+stencil
403  !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < &
404  ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))) then
405  ! maxslope = 4.0*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k))
406  !else
407  ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))
408  !endif
409  !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
410  ! slope_x(i,m) = 0.0
411  !elseif (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))<ABS(maxslope)) then
412  ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * &
413  ! 0.5*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))
414  !else
415  ! slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * 0.5*maxslope
416  !endif
417  tp = tr(m)%t(i+1,j,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i-1,j,k)
418  dmx = max( tp, tc, tm ) - tc
419  dmn= tc - min( tp, tc, tm )
420  slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
421  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
422  enddo ; enddo
423  endif ! usePLMslope
424 
425  ! make a copy of the tracers in case values need to be overridden for OBCs
426  do m = 1,ntr
427  do i=g%isd,g%ied
428  t_tmp(i,m) = tr(m)%t(i,j,k)
429  enddo
430  enddo
431  ! loop through open boundaries and recalculate flux terms
432  if (associated(obc)) then ; if (obc%OBC_pe) then
433  do n=1,obc%number_of_segments
434  segment=>obc%segment(n)
435  if (.not. associated(segment%tr_Reg)) cycle
436  if (segment%is_E_or_W) then
437  if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
438  i = segment%HI%IsdB
439  do m = 1,ntr ! replace tracers with OBC values
440  if (associated(segment%tr_Reg%Tr(m)%tres)) then
441  if (segment%direction == obc_direction_w) then
442  t_tmp(i,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
443  else
444  t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
445  endif
446  else
447  if (segment%direction == obc_direction_w) then
448  t_tmp(i,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
449  else
450  t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
451  endif
452  endif
453  enddo
454  do m = 1,ntr ! Apply update tracer values for slope calculation
455  do i=segment%HI%IsdB-1,segment%HI%IsdB+1
456  tp = t_tmp(i+1,m) ; tc = t_tmp(i,m) ; tm = t_tmp(i-1,m)
457  dmx = max( tp, tc, tm ) - tc
458  dmn= tc - min( tp, tc, tm )
459  slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
460  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
461  enddo
462  enddo
463 
464  endif
465  endif
466  enddo
467  endif; endif
468 
469 
470  ! Calculate the i-direction fluxes of each tracer, using as much
471  ! the minimum of the remaining mass flux (uhr) and the half the mass
472  ! in the cell plus whatever part of its half of the mass flux that
473  ! the flux through the other side does not require.
474  do i=is-1,ie
475  if (uhr(i,j,k) == 0.0) then
476  uhh(i) = 0.0
477  cfl(i) = 0.0
478  elseif (uhr(i,j,k) < 0.0) then
479  hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
480  hlos = max(0.0,uhr(i+1,j,k))
481  if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
482  ((0.5*hup + uhr(i,j,k)) < 0.0)) then
483  uhh(i) = min(-0.5*hup,-hup+hlos,0.0)
484  domore_u(j,k) = .true.
485  else
486  uhh(i) = uhr(i,j,k)
487  endif
488  !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)))
489  cfl(i) = - uhh(i) / (hprev(i+1,j,k) + h_neglect*g%areaT(i+1,j)) ! CFL is positive
490  else
491  hup = hprev(i,j,k) - g%areaT(i,j)*min_h
492  hlos = max(0.0,-uhr(i-1,j,k))
493  if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
494  ((0.5*hup - uhr(i,j,k)) < 0.0)) then
495  uhh(i) = max(0.5*hup,hup-hlos,0.0)
496  domore_u(j,k) = .true.
497  else
498  uhh(i) = uhr(i,j,k)
499  endif
500  !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
501  cfl(i) = uhh(i) / (hprev(i,j,k) + h_neglect*g%areaT(i,j)) ! CFL is positive
502  endif
503  enddo
504 
505 
506  if (useppm) then
507  do m=1,ntr ; do i=is-1,ie
508  ! centre cell depending on upstream direction
509  if (uhh(i) >= 0.0) then
510  i_up = i
511  else
512  i_up = i+1
513  endif
514 
515  ! Implementation of PPM-H3
516  tp = t_tmp(i_up+1,m) ; tc = t_tmp(i_up,m) ; tm = t_tmp(i_up-1,m)
517 
518  if (usehuynh) then
519  al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
520  al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
521  ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
522  ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
523  else
524  al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
525  ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
526  endif
527 
528  da = ar - al ; ma = 0.5*( ar + al )
529  if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.) then
530  al = tc ; ar = tc ! PCM for local extremum and bounadry cells
531  elseif ( da*(tc-ma) > (da*da)/6. ) then
532  al = 3.*tc - 2.*ar
533  elseif ( da*(tc-ma) < - (da*da)/6. ) then
534  ar = 3.*tc - 2.*al
535  endif
536 
537  a6 = 6.*tc - 3. * (ar + al) ! Curvature
538 
539  if (uhh(i) >= 0.0) then
540  flux_x(i,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
541  ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
542  else
543  flux_x(i,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
544  ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
545  endif
546  enddo ; enddo
547  else ! PLM
548  do m=1,ntr ; do i=is-1,ie
549  if (uhh(i) >= 0.0) then
550  ! Indirect implementation of PLM
551  !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m)
552  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
553  !flux_x(I,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) )
554  ! Alternative implementation of PLM
555  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
556  !flux_x(I,m) = uhh(I)*( aR - 0.5 * slope_x(i,m) * CFL(I) )
557  ! Alternative implementation of PLM
558  tc = t_tmp(i,m)
559  flux_x(i,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
560  ! Original implementation of PLM
561  !flux_x(I,m) = uhh(I)*(Tr(m)%t(i,j,k) + slope_x(i,m)*ts2(I))
562  else
563  ! Indirect implementation of PLM
564  !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
565  !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m)
566  !flux_x(I,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) )
567  ! Alternative implementation of PLM
568  !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
569  !flux_x(I,m) = uhh(I)*( aL + 0.5 * slope_x(i+1,m) * CFL(I) )
570  ! Alternative implementation of PLM
571  tc = t_tmp(i+1,m)
572  flux_x(i,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
573  ! Original implementation of PLM
574  !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I))
575  endif
576  !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j)))
577  enddo ; enddo
578  endif ! usePPM
579 
580  if (associated(obc)) then ; if (obc%OBC_pe) then
581  if (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally) then
582  do n=1,obc%number_of_segments
583  segment=>obc%segment(n)
584  if (.not. associated(segment%tr_Reg)) cycle
585  if (segment%is_E_or_W) then
586  if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
587  i = segment%HI%IsdB
588  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
589  ! Now changing to simply fixed inflows.
590  if ((uhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_w) .or. &
591  (uhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_e)) then
592  uhh(i) = uhr(i,j,k)
593  ! should the reservoir evolve for this case Kate ?? - Nope
594  do m=1,ntr
595  if (associated(segment%tr_Reg%Tr(m)%tres)) then
596  flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
597  else ; flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
598  enddo
599  endif
600  endif
601  endif
602  enddo
603  endif
604 
605  if (obc%open_u_BCs_exist_globally) then
606  do n=1,obc%number_of_segments
607  segment=>obc%segment(n)
608  i = segment%HI%IsdB
609  if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
610  if (segment%specified) cycle
611  if (.not. associated(segment%tr_Reg)) cycle
612 
613  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
614  if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
615  (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5)) then
616  uhh(i) = uhr(i,j,k)
617  do m=1,ntr
618  if (associated(segment%tr_Reg%Tr(m)%tres)) then
619  flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
620  else; flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif
621  enddo
622  endif
623  endif
624  enddo
625  endif
626  endif ; endif
627 
628  ! Calculate new tracer concentration in each cell after accounting
629  ! for the i-direction fluxes.
630  do i=is-1,ie
631  uhr(i,j,k) = uhr(i,j,k) - uhh(i)
632  if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
633  enddo
634  do i=is,ie
635  if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0)) then
636  do_i(i) = .true.
637  hlst(i) = hprev(i,j,k)
638  hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
639  if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false.
640  elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
641  hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
642  ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
643  else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
644  else
645  do_i(i) = .false.
646  endif
647  enddo
648 
649  ! update tracer concentration from i-flux and save some diagnostics
650  do m=1,ntr
651 
652  ! update tracer
653  do i=is,ie
654  if (do_i(i)) then
655  if (ihnew(i) > 0.0) then
656  tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
657  (flux_x(i,m) - flux_x(i-1,m))) * ihnew(i)
658  endif
659  endif
660  enddo
661 
662  ! diagnostics
663  if (associated(tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i)) then
664  tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,m)*idt
665  endif ; enddo ; endif
666  if (associated(tr(m)%ad2d_x)) then ; do i=is,ie ; if (do_i(i)) then
667  tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,m)*idt
668  endif ; enddo ; endif
669 
670  ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic).
671  ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
672  if (associated(tr(m)%advection_xy)) then
673  do i=is,ie ; if (do_i(i)) then
674  tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_x(i,m) - flux_x(i-1,m)) * &
675  idt * g%IareaT(i,j)
676  endif ; enddo
677  endif
678 
679  enddo
680 
681  endif ; enddo ! End of j-loop.
682 

Referenced by advect_tracer().

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,
type(unit_scale_type), intent(in)  US,
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
[in,out]trThe array of registered tracers to work on
[in,out]hprevcell volume at the end of previous tracer change [H L2 ~> m3 or kg]
[in,out]vhraccumulated volume/mass flux through the meridional face [H L2 ~> m3 or kg]
[in,out]vh_neglectA tiny meridional mass flux that can be neglected [H L2 ~> m3 or kg]
obcspecifies whether, where, and what OBCs are used
[in,out]domore_vIf true, there is more advection to be done in this v-row
[in]idtThe inverse of dt [T-1 ~> s-1]
[in]ntrThe number of tracers
[in]isThe starting tracer i-index to work on
[in]ieThe ending tracer i-index to work on
[in]jsThe starting tracer j-index to work on
[in]jeThe ending tracer j-index to work on
[in]kThe k-level to work on
[in]usA dimensional unit scaling type
[in]useppmIf true, use PPM instead of PLM
[in]usehuynhIf true, use the Huynh scheme for PPM interface values

Definition at line 689 of file MOM_tracer_advect.F90.

689  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
690  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
691  type(tracer_type), dimension(ntr), intent(inout) :: Tr !< The array of registered tracers to work on
692  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hprev !< cell volume at the end of previous
693  !! tracer change [H L2 ~> m3 or kg]
694  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhr !< accumulated volume/mass flux through
695  !! the meridional face [H L2 ~> m3 or kg]
696  real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vh_neglect !< A tiny meridional mass flux that can
697  !! be neglected [H L2 ~> m3 or kg]
698  type(ocean_OBC_type), pointer :: OBC !< specifies whether, where, and what OBCs are used
699  logical, dimension(SZJB_(G),SZK_(G)), intent(inout) :: domore_v !< If true, there is more advection to be
700  !! done in this v-row
701  real, intent(in) :: Idt !< The inverse of dt [T-1 ~> s-1]
702  integer, intent(in) :: ntr !< The number of tracers
703  integer, intent(in) :: is !< The starting tracer i-index to work on
704  integer, intent(in) :: ie !< The ending tracer i-index to work on
705  integer, intent(in) :: js !< The starting tracer j-index to work on
706  integer, intent(in) :: je !< The ending tracer j-index to work on
707  integer, intent(in) :: k !< The k-level to work on
708  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
709  logical, intent(in) :: usePPM !< If true, use PPM instead of PLM
710  logical, intent(in) :: useHuynh !< If true, use the Huynh scheme
711  !! for PPM interface values
712 
713  real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
714  slope_y ! The concentration slope per grid point [conc].
715  real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
716  flux_y ! The tracer flux across a boundary [H m2 conc ~> m3 conc or kg conc].
717  real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
718  T_tmp ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].
719  real :: maxslope ! The maximum concentration slope per grid point
720  ! consistent with monotonicity [conc].
721  real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
722  ! current iteration [H L2 ~> m3 or kg].
723  real :: hup, hlos ! hup is the upwind volume, hlos is the
724  ! part of that volume that might be lost
725  ! due to advection out the other side of
726  ! the grid box, both in [H L2 ~> m3 or kg].
727  real, dimension(SZIB_(G)) :: &
728  hlst, & ! Work variable [H L2 ~> m3 or kg].
729  Ihnew, & ! Work variable [H-1 L-2 ~> m-3 or kg-1].
730  CFL ! A nondimensional work variable.
731  real :: min_h ! The minimum thickness that can be realized during
732  ! any of the passes [H ~> m or kg m-2].
733  real :: h_neglect ! A thickness that is so small it is usually lost
734  ! in roundoff and can be neglected [H ~> m or kg m-2].
735  logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
736  logical :: do_i(SZIB_(G)) ! If true, work on given points.
737  logical :: do_any_i
738  integer :: i, j, j2, m, n, j_up, stencil
739  real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
740  real :: fac1,v_L_in,v_L_out ! terms used for time-stepping OBC reservoirs
741  type(OBC_segment_type), pointer :: segment=>null()
742  logical :: usePLMslope
743 
744  useplmslope = .not. (useppm .and. usehuynh)
745  ! stencil for calculating slope values
746  stencil = 1
747  if (useppm .and. .not. usehuynh) stencil = 2
748 
749  min_h = 0.1*gv%Angstrom_H
750  h_neglect = gv%H_subroundoff
751  !do i=is,ie ; ts2(i) = 0.0 ; enddo
752 
753  ! We conditionally perform work on tracer points: calculating the PLM slope,
754  ! and updating tracer concentration within a cell
755  ! this depends on whether there is a flux which would affect this tracer point,
756  ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires
757  ! slope calculations at the two tracer points on either side (as indicated by
758  ! the stencil variable), so we account for this with the do_j_tr flag array
759  !
760  ! Note: this does lead to unnecessary work in updating tracer concentrations,
761  ! since that doesn't need a wider stencil with the PPM advection scheme, but
762  ! this would require an additional loop, etc.
763  do_j_tr(:) = .false.
764  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
765 
766  ! Calculate the j-direction profiles (slopes) of each tracer that
767  ! is being advected.
768  if (useplmslope) then
769  do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr ; do i=is,ie
770  !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < &
771  ! ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then
772  ! maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))
773  !else
774  ! maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))
775  !endif
776  !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
777  ! slope_y(i,m,j) = 0.0
778  !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then
779  ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * &
780  ! 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))
781  !else
782  ! slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope
783  !endif
784  tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
785  dmx = max( tp, tc, tm ) - tc
786  dmn= tc - min( tp, tc, tm )
787  slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
788  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
789  enddo ; enddo ; endif ; enddo ! End of i-, m-, & j- loops.
790  endif ! usePLMslope
791 
792 
793  ! make a copy of the tracers in case values need to be overridden for OBCs
794 
795  do j=g%jsd,g%jed; do m=1,ntr; do i=g%isd,g%ied
796  t_tmp(i,m,j) = tr(m)%t(i,j,k)
797  enddo ; enddo ; enddo
798 
799  ! loop through open boundaries and recalculate flux terms
800  if (associated(obc)) then ; if (obc%OBC_pe) then
801  do n=1,obc%number_of_segments
802  segment=>obc%segment(n)
803  if (.not. associated(segment%tr_Reg)) cycle
804  do i=is,ie
805  if (segment%is_N_or_S) then
806  if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
807  j = segment%HI%JsdB
808  do m = 1,ntr ! replace tracers with OBC values
809  if (associated(segment%tr_Reg%Tr(m)%tres)) then
810  if (segment%direction == obc_direction_s) then
811  t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
812  else
813  t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
814  endif
815  else
816  if (segment%direction == obc_direction_s) then
817  t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
818  else
819  t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
820  endif
821  endif
822  enddo
823  do m = 1,ntr ! Apply update tracer values for slope calculation
824  do j=segment%HI%JsdB-1,segment%HI%JsdB+1
825  tp = t_tmp(i,m,j+1) ; tc = t_tmp(i,m,j) ; tm = t_tmp(i,m,j-1)
826  dmx = max( tp, tc, tm ) - tc
827  dmn= tc - min( tp, tc, tm )
828  slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
829  sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
830  enddo
831  enddo
832  endif
833  endif ! is_N_S
834  enddo ! i-loop
835  enddo ! segment loop
836  endif; endif
837 
838  ! Calculate the j-direction fluxes of each tracer, using as much
839  ! the minimum of the remaining mass flux (vhr) and the half the mass
840  ! in the cell plus whatever part of its half of the mass flux that
841  ! the flux through the other side does not require.
842  do j=js-1,je ; if (domore_v(j,k)) then
843  domore_v(j,k) = .false.
844 
845  do i=is,ie
846  if (vhr(i,j,k) == 0.0) then
847  vhh(i,j) = 0.0
848  cfl(i) = 0.0
849  elseif (vhr(i,j,k) < 0.0) then
850  hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
851  hlos = max(0.0,vhr(i,j+1,k))
852  if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
853  ((0.5*hup + vhr(i,j,k)) < 0.0)) then
854  vhh(i,j) = min(-0.5*hup,-hup+hlos,0.0)
855  domore_v(j,k) = .true.
856  else
857  vhh(i,j) = vhr(i,j,k)
858  endif
859  !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1))
860  cfl(i) = - vhh(i,j) / (hprev(i,j+1,k) + h_neglect*g%areaT(i,j+1)) ! CFL is positive
861  else
862  hup = hprev(i,j,k) - g%areaT(i,j)*min_h
863  hlos = max(0.0,-vhr(i,j-1,k))
864  if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
865  ((0.5*hup - vhr(i,j,k)) < 0.0)) then
866  vhh(i,j) = max(0.5*hup,hup-hlos,0.0)
867  domore_v(j,k) = .true.
868  else
869  vhh(i,j) = vhr(i,j,k)
870  endif
871  !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
872  cfl(i) = vhh(i,j) / (hprev(i,j,k) + h_neglect*g%areaT(i,j)) ! CFL is positive
873  endif
874  enddo
875 
876  if (useppm) then
877  do m=1,ntr ; do i=is,ie
878  ! centre cell depending on upstream direction
879  if (vhh(i,j) >= 0.0) then
880  j_up = j
881  else
882  j_up = j + 1
883  endif
884 
885  ! Implementation of PPM-H3
886  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)
887 
888  if (usehuynh) then
889  al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
890  al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
891  ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
892  ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
893  else
894  al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
895  ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
896  endif
897 
898  da = ar - al ; ma = 0.5*( ar + al )
899  if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.) then
900  al = tc ; ar = tc ! PCM for local extremum and bounadry cells
901  elseif ( da*(tc-ma) > (da*da)/6. ) then
902  al = 3.*tc - 2.*ar
903  elseif ( da*(tc-ma) < - (da*da)/6. ) then
904  ar = 3.*tc - 2.*al
905  endif
906 
907  a6 = 6.*tc - 3. * (ar + al) ! Curvature
908 
909  if (vhh(i,j) >= 0.0) then
910  flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
911  ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
912  else
913  flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
914  ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
915  endif
916  enddo ; enddo
917  else ! PLM
918  do m=1,ntr ; do i=is,ie
919  if (vhh(i,j) >= 0.0) then
920  ! Indirect implementation of PLM
921  !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j)
922  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
923  !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) )
924  ! Alternative implementation of PLM
925  !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
926  !flux_y(i,m,J) = vhh(i,J)*(aR - 0.5 * slope_y(i,m,j)*CFL(i))
927  ! Alternative implementation of PLM
928  tc = tr(m)%t(i,j,k)
929  flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
930  ! Original implementation of PLM
931  !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j,k) + slope_y(i,m,j)*ts2(i))
932  else
933  ! Indirect implementation of PLM
934  !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
935  !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1)
936  !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) )
937  ! Alternative implementation of PLM
938  !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
939  !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * slope_y(i,m,j+1)*CFL(i) )
940  ! Alternative implementation of PLM
941  tc = tr(m)%t(i,j+1,k)
942  flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
943  ! Original implementation of PLM
944  !flux_y(i,m,J) = vhh(i,J)*(Tr(m)%t(i,j+1,k) - slope_y(i,m,j+1)*ts2(i))
945  endif
946  enddo ; enddo
947  endif ! usePPM
948 
949  if (associated(obc)) then ; if (obc%OBC_pe) then
950  if (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) then
951  do n=1,obc%number_of_segments
952  segment=>obc%segment(n)
953  if (.not. segment%specified) cycle
954  if (.not. associated(segment%tr_Reg)) cycle
955  if (obc%segment(n)%is_N_or_S) then
956  if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB) then
957  do i=segment%HI%isd,segment%HI%ied
958  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
959  ! Now changing to simply fixed inflows.
960  if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
961  (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n)) then
962  vhh(i,j) = vhr(i,j,k)
963  do m=1,ntr
964  if (associated(segment%tr_Reg%Tr(m)%t)) then
965  flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
966  else ; flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
967  enddo
968  endif
969  enddo
970  endif
971  endif
972  enddo
973  endif
974 
975  if (obc%open_v_BCs_exist_globally) then
976  do n=1,obc%number_of_segments
977  segment=>obc%segment(n)
978  if (segment%specified) cycle
979  if (.not. associated(segment%tr_Reg)) cycle
980  if (segment%is_N_or_S .and. (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)) then
981  do i=segment%HI%isd,segment%HI%ied
982  ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
983  if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
984  (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5)) then
985  vhh(i,j) = vhr(i,j,k)
986  do m=1,ntr
987  if (associated(segment%tr_Reg%Tr(m)%t)) then
988  flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
989  else ; flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
990  enddo
991  endif
992  enddo
993  endif
994  enddo
995  endif
996  endif; endif
997 
998  else ! not domore_v.
999  do i=is,ie ; vhh(i,j) = 0.0 ; enddo
1000  do m=1,ntr ; do i=is,ie ; flux_y(i,m,j) = 0.0 ; enddo ; enddo
1001  endif ; enddo ! End of j-loop
1002 
1003  do j=js-1,je ; do i=is,ie
1004  vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
1005  if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
1006  enddo ; enddo
1007 
1008  ! Calculate new tracer concentration in each cell after accounting
1009  ! for the j-direction fluxes.
1010  do j=js,je ; if (do_j_tr(j)) then
1011  do i=is,ie
1012  if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0)) then
1013  do_i(i) = .true.
1014  hlst(i) = hprev(i,j,k)
1015  hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
1016  if (hprev(i,j,k) <= 0.0) then ; do_i(i) = .false.
1017  elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
1018  hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
1019  ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
1020  else ; ihnew(i) = 1.0 / hprev(i,j,k) ; endif
1021  else ; do_i(i) = .false. ; endif
1022  enddo
1023 
1024  ! update tracer and save some diagnostics
1025  do m=1,ntr
1026  do i=is,ie ; if (do_i(i)) then
1027  tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
1028  (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
1029  endif ; enddo
1030 
1031  ! diagnostics
1032  if (associated(tr(m)%ad_y)) then ; do i=is,ie ; if (do_i(i)) then
1033  tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
1034  endif ; enddo ; endif
1035  if (associated(tr(m)%ad2d_y)) then ; do i=is,ie ; if (do_i(i)) then
1036  tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
1037  endif ; enddo ; endif
1038 
1039  ! diagnose convergence of flux_y and add to convergence of flux_x.
1040  ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
1041  if (associated(tr(m)%advection_xy)) then
1042  do i=is,ie ; if (do_i(i)) then
1043  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 * &
1044  g%IareaT(i,j)
1045  endif ; enddo
1046  endif
1047 
1048  enddo
1049  endif ; enddo ! End of j-loop.
1050 

Referenced by advect_tracer().

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.

Parameters
csmodule control structure

Definition at line 1110 of file MOM_tracer_advect.F90.

1110  type(tracer_advect_CS), pointer :: CS !< module control structure
1111 
1112  if (associated(cs)) deallocate(cs)
1113 

Referenced by mom::mom_end().

Here is the caller graph for this function:

◆ 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(unit_scale_type), intent(in)  US,
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]usA dimensional unit scaling type
[in]param_fileopen file to parse for model parameters
[in,out]diagregulates diagnostic output
csmodule control structure

Definition at line 1055 of file MOM_tracer_advect.F90.

1055  type(time_type), target, intent(in) :: Time !< current model time
1056  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1057  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1058  type(param_file_type), intent(in) :: param_file !< open file to parse for model parameters
1059  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
1060  type(tracer_advect_CS), pointer :: CS !< module control structure
1061 
1062  integer, save :: init_calls = 0
1063 
1064  ! This include declares and sets the variable "version".
1065 # include "version_variable.h"
1066  character(len=40) :: mdl = "MOM_tracer_advect" ! This module's name.
1067  character(len=256) :: mesg ! Message for error messages.
1068 
1069  if (associated(cs)) then
1070  call mom_error(warning, "tracer_advect_init called with associated control structure.")
1071  return
1072  endif
1073  allocate(cs)
1074 
1075  cs%diag => diag
1076 
1077  ! Read all relevant parameters and write them to the model log.
1078  call log_version(param_file, mdl, version, "")
1079  call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
1080  desc="The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T)
1081  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1082  call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, &
1083  desc="The horizontal transport scheme for tracers:\n"//&
1084  " PLM - Piecewise Linear Method\n"//&
1085  " PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// &
1086  " PPM - Piecewise Parabolic Method (Colella-Woodward)" &
1087  , default='PLM')
1088  select case (trim(mesg))
1089  case ("PLM")
1090  cs%usePPM = .false.
1091  case ("PPM:H3")
1092  cs%usePPM = .true.
1093  cs%useHuynh = .true.
1094  case ("PPM")
1095  cs%usePPM = .true.
1096  cs%useHuynh = .false.
1097  case default
1098  call mom_error(fatal, "MOM_tracer_advect, tracer_advect_init: "//&
1099  "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
1100  end select
1101 
1102  id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=clock_module)
1103  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
1104  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)
1105 

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

Here is the call graph for this function:

Variable Documentation

◆ id_clock_advect

integer mom_tracer_advect::id_clock_advect

CPU time clocks.

Definition at line 41 of file MOM_tracer_advect.F90.

41 integer :: id_clock_advect

Referenced by advect_tracer(), and tracer_advect_init().

◆ id_clock_pass

integer mom_tracer_advect::id_clock_pass
private

CPU time clocks.

Definition at line 42 of file MOM_tracer_advect.F90.

42 integer :: id_clock_pass

Referenced by advect_tracer(), and tracer_advect_init().

◆ id_clock_sync

integer mom_tracer_advect::id_clock_sync
private

CPU time clocks.

Definition at line 43 of file MOM_tracer_advect.F90.

43 integer :: id_clock_sync

Referenced by advect_tracer(), and tracer_advect_init().