ROMS
Loading...
Searching...
No Matches
step3d_t.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE step3d_t_mod
3#if !defined TS_FIXED && (defined NONLINEAR && defined SOLVE3D)
4!
5!git $Id$
6!=======================================================================
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md Hernan G. Arango !
10!========================================== Alexander F. Shchepetkin ===
11! !
12! This routine time-steps tracer equations. Notice that advective !
13! and diffusive terms are time-stepped differently. It applies the !
14! corrector time-step for horizontal/vertical advection, vertical !
15! diffusion, nudging if necessary, and lateral boundary conditions. !
16! !
17! A different horizontal/vertical advection scheme is allowed for !
18! each tracer. If the MPDATA or HSIMT monotonic scheme, it is applied !
19! to both horizontal and vertical advective fluxes. !
20! !
21! Notice that at input the tracer arrays have: !
22! !
23! t(:,:,:,nnew,:) m Tunits n+1 horizontal/vertical diffusion !
24! terms plus source/sink terms !
25! (biology, sediment), if any !
26! !
27! t(:,:,:,3 ,:) Tunits n+1/2 advective terms and vertical !
28! diffusion predictor step !
29! !
30!=======================================================================
31!
32 implicit none
33!
34 PRIVATE
35 PUBLIC :: step3d_t
36!
37 CONTAINS
38!
39!***********************************************************************
40 SUBROUTINE step3d_t (ng, tile)
41!***********************************************************************
42!
43 USE mod_param
44# ifdef DIAGNOSTICS_TS
45 USE mod_diags
46# endif
47 USE mod_grid
48 USE mod_mixing
49 USE mod_ocean
50# if defined SEDIMENT && defined SED_MORPH
51 USE mod_sedbed
52# endif
53 USE mod_stepping
54!
55! Imported variable declarations.
56!
57 integer, intent(in) :: ng, tile
58!
59! Local variable declarations.
60!
61 character (len=*), parameter :: myfile = &
62 & __FILE__
63!
64# include "tile.h"
65!
66# ifdef PROFILE
67 CALL wclock_on (ng, inlm, 35, __line__, myfile)
68# endif
69 CALL step3d_t_tile (ng, tile, &
70 & lbi, ubi, lbj, ubj, &
71 & imins, imaxs, jmins, jmaxs, &
72 & nrhs(ng), nstp(ng), nnew(ng), &
73# ifdef MASKING
74 & grid(ng) % rmask, &
75 & grid(ng) % umask, &
76 & grid(ng) % vmask, &
77# endif
78# ifdef WET_DRY
79 & grid(ng) % rmask_wet, &
80 & grid(ng) % umask_wet, &
81 & grid(ng) % vmask_wet, &
82# endif
83 & grid(ng) % omn, &
84 & grid(ng) % om_u, &
85 & grid(ng) % om_v, &
86 & grid(ng) % on_u, &
87 & grid(ng) % on_v, &
88 & grid(ng) % pm, &
89 & grid(ng) % pn, &
90 & grid(ng) % Hz, &
91 & grid(ng) % Huon, &
92 & grid(ng) % Hvom, &
93 & grid(ng) % z_r, &
94 & mixing(ng) % Akt, &
95 & ocean(ng) % W, &
96# ifdef OMEGA_IMPLICIT
97 & ocean(ng) % Wi, &
98# endif
99# ifdef WEC_VF
100 & ocean(ng) % W_stokes, &
101# endif
102# if defined SEDIMENT && defined SED_MORPH
103 & sedbed(ng) % bed_thick, &
104# endif
105# if defined FLOATS && defined FLOAT_VWALK
106 & mixing(ng) % dAktdz, &
107# endif
108# ifdef DIAGNOSTICS_TS
109 & diags(ng) % DiaTwrk, &
110# endif
111 & ocean(ng) % t)
112# ifdef PROFILE
113 CALL wclock_off (ng, inlm, 35, __line__, myfile)
114# endif
115!
116 RETURN
117 END SUBROUTINE step3d_t
118!
119!***********************************************************************
120 SUBROUTINE step3d_t_tile (ng, tile, &
121 & LBi, UBi, LBj, UBj, &
122 & IminS, ImaxS, JminS, JmaxS, &
123 & nrhs, nstp, nnew, &
124# ifdef MASKING
125 & rmask, umask, vmask, &
126# endif
127# ifdef WET_DRY
128 & rmask_wet, umask_wet, vmask_wet, &
129# endif
130 & omn, om_u, om_v, on_u, on_v, &
131 & pm, pn, &
132 & Hz, Huon, Hvom, &
133 & z_r, &
134 & Akt, &
135 & W, &
136# ifdef OMEGA_IMPLICIT
137 & Wi, &
138# endif
139# ifdef WEC_VF
140 & W_stokes, &
141# endif
142# if defined SEDIMENT && defined SED_MORPH
143 & bed_thick, &
144# endif
145# if defined FLOATS && defined FLOAT_VWALK
146 & dAktdz, &
147# endif
148# ifdef DIAGNOSTICS_TS
149 & DiaTwrk, &
150# endif
151 & t)
152!***********************************************************************
153!
154 USE mod_param
155 USE mod_clima
156 USE mod_ncparam
157# if defined NESTING && !defined ONE_WAY
158 USE mod_nesting
159# endif
160 USE mod_scalars
161 USE mod_sources
162!
163 USE exchange_3d_mod, ONLY : exchange_r3d_tile
164# ifdef DISTRIBUTE
165 USE mp_exchange_mod, ONLY : mp_exchange3d
166 USE mp_exchange_mod, ONLY : mp_exchange4d
167# endif
168 USE mpdata_adiff_mod
169# ifdef NESTING
170 USE nesting_mod, ONLY : bry_fluxes
171# endif
172 USE t3dbc_mod, ONLY : t3dbc_tile
173!
174! Imported variable declarations.
175!
176 integer, intent(in) :: ng, tile
177 integer, intent(in) :: LBi, UBi, LBj, UBj
178 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
179 integer, intent(in) :: nrhs, nstp, nnew
180!
181# ifdef ASSUMED_SHAPE
182# ifdef MASKING
183 real(r8), intent(in) :: rmask(LBi:,LBj:)
184 real(r8), intent(in) :: umask(LBi:,LBj:)
185 real(r8), intent(in) :: vmask(LBi:,LBj:)
186# endif
187# ifdef WET_DRY
188 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
189 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
190 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
191# endif
192 real(r8), intent(in) :: omn(LBi:,LBj:)
193 real(r8), intent(in) :: om_u(LBi:,LBj:)
194 real(r8), intent(in) :: om_v(LBi:,LBj:)
195 real(r8), intent(in) :: on_u(LBi:,LBj:)
196 real(r8), intent(in) :: on_v(LBi:,LBj:)
197 real(r8), intent(in) :: pm(LBi:,LBj:)
198 real(r8), intent(in) :: pn(LBi:,LBj:)
199 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
200 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
201 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
202 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
203# ifdef SUN
204 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
205# else
206 real(r8), intent(in) :: Akt(LBi:,LBj:,0:,:)
207# endif
208 real(r8), intent(in) :: W(LBi:,LBj:,0:)
209# ifdef OMEGA_IMPLICIT
210 real(r8), intent(in) :: Wi(LBi:,LBj:,0:)
211# endif
212# ifdef WEC_VF
213 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
214# endif
215# if defined SEDIMENT && defined SED_MORPH
216 real(r8), intent(in) :: bed_thick(LBi:,LBj:,:)
217# endif
218# ifdef DIAGNOSTICS_TS
219 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
220# endif
221# ifdef SUN
222 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
223# else
224 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
225# endif
226# if defined FLOATS && defined FLOAT_VWALK
227 real(r8), intent(out) :: dAktdz(LBi:,LBj:,:)
228# endif
229
230# else
231
232# ifdef MASKING
233 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
234 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
235 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
236# endif
237# ifdef WET_DRY
238 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
239 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
240 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
241# endif
242 real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
243 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
244 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
245 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
246 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
247 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
248 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
249 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
250 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
251 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
252 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
253 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
254 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
255# ifdef OMEGA_IMPLICIT
256 real(r8), intent(in) :: Wi(LBi:UBi,LBj:UBj,0:N(ng))
257# endif
258# ifdef WEC_VF
259 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
260# endif
261# if defined SEDIMENT && defined SED_MORPH
262 real(r8), intent(in) :: bed_thick(LBi:UBi,LBj:UBj,3)
263# endif
264# ifdef DIAGNOSTICS_TS
265 real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
266 & NDT)
267# endif
268 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
269
270# if defined FLOATS && defined FLOAT_VWALK
271 real(r8), intent(out) :: dAktdz(LBi:UBi,LBj:UBj,N(ng))
272# endif
273# endif
274!
275! Local variable declarations.
276!
277 logical :: LapplySrc, Lhsimt, Lmpdata
278!
279# ifdef NESTING
280 integer :: ILB, IUB, JLB, JUB
281 integer :: dg, cr, rg
282# endif
283 integer :: IminT, ImaxT, JminT, JmaxT
284 integer :: Isrc, Jsrc
285 integer :: i, ic, ii, is, itrc, j, jj, k, ltrc
286# if defined AGE_MEAN && defined T_PASSIVE
287 integer :: iage
288# endif
289# ifdef DIAGNOSTICS_TS
290 integer :: idiag
291# endif
292 real(r8) :: eps = 1.0e-16_r8
293 real(r8) :: eps1 = 1.0e-12_r8
294
295 real(r8) :: cff, cff1, cff2, cff3
296 real(r8) :: betaL, betaR, betaD, betaU
297 real(r8) :: rL, rR, rD, rU, rkaL, rkaR, rkaD, rkaU
298 real(r8) :: a1, b1, sw, sw_eta, sw_xi
299
300 real(r8), dimension(IminS:ImaxS) :: gradX, KaX, oKaX
301 real(r8), dimension(JminS:JmaxS) :: gradE, KaE, oKaE
302 real(r8), dimension(0:N(ng)) :: gradZ, KaZ, oKaZ
303
304 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
305 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
306 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
307 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
308# ifdef OMEGA_IMPLICIT
309 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmin
310 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmax
311# endif
312
313 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
314 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
315 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curv
316 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
317
318 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
319
320# ifdef DIAGNOSTICS_TS
321 real(r8), allocatable :: Dhadv(:,:,:)
322 real(r8), allocatable :: Dvadv(:,:,:,:)
323# endif
324 real(r8), allocatable :: Ta(:,:,:,:)
325 real(r8), allocatable :: Ua(:,:,:)
326 real(r8), allocatable :: Va(:,:,:)
327 real(r8), allocatable :: Wa(:,:,:)
328
329# include "set_bounds.h"
330
331# ifdef NESTING
332!
333! Notice that the trace flux boundary arrays are dimensioned with the
334! global dimensions of grid to facilitate processing.
335!
336 ilb=bounds(ng)%LBi(-1)
337 iub=bounds(ng)%UBi(-1)
338 jlb=bounds(ng)%LBj(-1)
339 jub=bounds(ng)%UBj(-1)
340# endif
341!
342!-----------------------------------------------------------------------
343! Time-step horizontal advection term.
344!-----------------------------------------------------------------------
345!
346 lhsimt =any(hadvection(:,ng)%HSIMT).and. &
347 & any(vadvection(:,ng)%HSIMT)
348 lmpdata=any(hadvection(:,ng)%MPDATA).and. &
349 & any(vadvection(:,ng)%MPDATA)
350!
351! Allocate local arrays for MPDATA.
352!
353 IF (lmpdata) THEN
354# ifdef DIAGNOSTICS_TS
355 IF (.not.allocated(dhadv)) THEN
356 allocate ( dhadv(imins:imaxs,jmins:jmaxs,3) )
357 dhadv=0.0_r8
358 END IF
359 IF (.not.allocated(dvadv)) THEN
360 allocate ( dvadv(imins:imaxs,jmins:jmaxs,n(ng),nt(ng)) )
361 dvadv=0.0_r8
362 END IF
363# endif
364 IF (.not.allocated(ta)) THEN
365 allocate ( ta(imins:imaxs,jmins:jmaxs,n(ng),nt(ng)) )
366 ta=0.0_r8
367 END IF
368 IF (.not.allocated(ua)) THEN
369 allocate ( ua(imins:imaxs,jmins:jmaxs,n(ng)) )
370 ua=0.0_r8
371 END IF
372 IF (.not.allocated(va)) THEN
373 allocate ( va(imins:imaxs,jmins:jmaxs,n(ng)) )
374 va=0.0_r8
375 END IF
376 IF (.not.allocated(wa)) THEN
377 allocate ( wa(imins:imaxs,jmins:jmaxs,0:n(ng)) )
378 wa=0.0_r8
379 END IF
380 END IF
381!
382! Compute reciprocal thickness, 1/Hz.
383!
384 IF (lmpdata.or.lhsimt) THEN
385 DO k=1,n(ng)
386 DO j=jstrm2,jendp2
387 DO i=istrm2,iendp2
388 ohz(i,j,k)=1.0_r8/hz(i,j,k)
389 END DO
390 END DO
391 END DO
392 ELSE
393 DO k=1,n(ng)
394 DO j=jstr,jend
395 DO i=istr,iend
396 ohz(i,j,k)=1.0_r8/hz(i,j,k)
397 END DO
398 END DO
399 END DO
400 END IF
401!
402! Horizontal tracer advection. It is possible to have a different
403! advection schme for each tracer.
404!
405 t_loop1 : DO itrc=1,nt(ng)
406!
407! The MPDATA and HSIMT algorithms requires a three-point footprint, so
408! exchange boundary data on t(:,:,:,nnew,:) so other processes computed
409! earlier (horizontal diffusion, biology, or sediment) are accounted.
410!
411 IF ((hadvection(itrc,ng)%MPDATA).or. &
412 & (hadvection(itrc,ng)%HSIMT)) THEN
413 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
414 CALL exchange_r3d_tile (ng, tile, &
415 & lbi, ubi, lbj, ubj, 1, n(ng), &
416 & t(:,:,:,nnew,itrc))
417 END IF
418
419# ifdef DISTRIBUTE
420 CALL mp_exchange3d (ng, tile, inlm, 1, &
421 & lbi, ubi, lbj, ubj, 1, n(ng), &
422 & nghostpoints, &
423 & ewperiodic(ng), nsperiodic(ng), &
424 & t(:,:,:,nnew,itrc))
425# endif
426 END IF
427!
428! Compute horizontal tracer advection fluxes.
429!
430 k_loop : DO k=1,n(ng)
431!
432 hadv_flux : IF (hadvection(itrc,ng)%CENTERED2) THEN
433!
434! Second-order, centered differences horizontal advective fluxes.
435!
436 DO j=jstr,jend
437 DO i=istr,iend+1
438 fx(i,j)=huon(i,j,k)* &
439 & 0.5_r8*(t(i-1,j,k,3,itrc)+ &
440 & t(i ,j,k,3,itrc))
441 END DO
442 END DO
443 DO j=jstr,jend+1
444 DO i=istr,iend
445 fe(i,j)=hvom(i,j,k)* &
446 & 0.5_r8*(t(i,j-1,k,3,itrc)+ &
447 & t(i,j ,k,3,itrc))
448 END DO
449 END DO
450!
451 ELSE IF (hadvection(itrc,ng)%MPDATA) THEN
452!
453! First-order, upstream differences horizontal advective fluxes.
454!
455 DO j=jstrvm2,jendp2i
456 DO i=istrum2,iendp3
457 cff1=max(huon(i,j,k),0.0_r8)
458 cff2=min(huon(i,j,k),0.0_r8)
459 fx(i,j)=cff1*t(i-1,j,k,3,itrc)+ &
460 & cff2*t(i ,j,k,3,itrc)
461 END DO
462 END DO
463 DO j=jstrvm2,jendp3
464 DO i=istrum2,iendp2i
465 cff1=max(hvom(i,j,k),0.0_r8)
466 cff2=min(hvom(i,j,k),0.0_r8)
467 fe(i,j)=cff1*t(i,j-1,k,3,itrc)+ &
468 & cff2*t(i,j ,k,3,itrc)
469 END DO
470 END DO
471!
472 ELSE IF (hadvection(itrc,ng)%HSIMT) THEN
473!
474! Third High-order Spatial Interpolation at the Middle Temporal level
475! (HSIMT; Wu and Zhu, 2010) with a Total Variation Diminishing (TVD)
476! limiter horizontal advection fluxes.
477!
478! Hui Wu and Jianrong Zhu, 2010: Advection scheme with 3rd high-order
479! spatial interpolation at the middle temporal level and its
480! application to saltwater intrusion in the Changjiang Estuary,
481! Ocean Modelling 33, 33-51, doi:10.1016/j.ocemod.2009.12.001
482!
483 DO j=jstr,jend
484 DO i=istru-1,iendp2
485 cff=0.125_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
486 & dt(ng)
487 cff1=cff*(ohz(i-1,j,k)+ohz(i,j,k))
488 gradx(i)=t(i,j,k,3,itrc)-t(i-1,j,k,3,itrc)
489 kax(i)=1.0_r8-abs(huon(i,j,k)*cff1)
490# ifdef MASKING
491 gradx(i)=gradx(i)*umask(i,j)
492 kax(i)=kax(i)*umask(i,j)
493# endif
494 END DO
495 IF (.not.ewperiodic(ng)) THEN
496 IF (domain(ng)%Western_Edge(tile)) THEN
497 IF (huon(istr,j,k).ge.0.0_r8) THEN
498 gradx(istr-1)=0.0_r8
499 kax(istr-1)=0.0_r8
500 END IF
501 END IF
502 IF (domain(ng)%Eastern_Edge(tile)) THEN
503 IF (huon(iend+1,j,k).lt.0.0_r8) THEN
504 gradx(iend+2)=0.0_r8
505 kax(iend+2)=0.0_r8
506 END IF
507 END IF
508 END IF
509 DO i=istr,iend+1
510 IF (kax(i).le.eps1) THEN
511 okax(i)=0.0_r8
512 ELSE
513 okax(i)=1.0_r8/max(kax(i),eps1)
514 END IF
515 IF (huon(i,j,k).ge.0.0_r8) THEN
516 IF (abs(gradx(i)).le.eps1) THEN
517 rl=0.0_r8
518 rkal=0.0_r8
519 ELSE
520 rl=gradx(i-1)/gradx(i)
521 rkal=kax(i-1)*okax(i)
522 END IF
523 a1= cc1*kax(i)+cc2-cc3*okax(i)
524 b1=-cc1*kax(i)+cc2+cc3*okax(i)
525 betal=a1+b1*rl
526 cff=0.5_r8*max(0.0_r8, &
527 & min(2.0_r8, 2.0_r8*rl*rkal, betal))* &
528 & gradx(i)*kax(i)
529# ifdef MASKING
530 ii=max(i-2,0)
531 cff=cff*rmask(ii,j)
532# endif
533 sw_xi=t(i-1,j,k,3,itrc)+cff
534 ELSE
535 IF (abs(gradx(i)).le.eps1) THEN
536 rr=0.0_r8
537 rkar=0.0_r8
538 ELSE
539 rr=gradx(i+1)/gradx(i)
540 rkar=kax(i+1)*okax(i)
541 END IF
542 a1= cc1*kax(i)+cc2-cc3*okax(i)
543 b1=-cc1*kax(i)+cc2+cc3*okax(i)
544 betar=a1+b1*rr
545 cff=0.5_r8*max(0.0_r8, &
546 & min(2.0_r8, 2.0_r8*rr*rkar, betar))* &
547 & gradx(i)*kax(i)
548# ifdef MASKING
549 ii=min(i+1,lm(ng)+1)
550 cff=cff*rmask(ii,j)
551# endif
552 sw_xi=t(i,j,k,3,itrc)-cff
553 END IF
554 fx(i,j)=sw_xi*huon(i,j,k)
555 END DO
556 END DO
557!
558 DO i=istr,iend
559 DO j=jstrv-1,jendp2
560 cff=0.125_r8*(pn(i,j)+pn(i,j-1))*(pm(i,j)+pm(i,j-1))* &
561 & dt(ng)
562 cff1=cff*(ohz(i,j,k)+ohz(i,j-1,k))
563 grade(j)=t(i,j,k,3,itrc)-t(i,j-1,k,3,itrc)
564 kae(j)=1.0_r8-abs(hvom(i,j,k)*cff1)
565# ifdef MASKING
566 grade(j)=grade(j)*vmask(i,j)
567 kae(j)=kae(j)*vmask(i,j)
568# endif
569 END DO
570 IF (.not.nsperiodic(ng)) THEN
571 IF (domain(ng)%Southern_Edge(tile)) THEN
572 IF (hvom(i,jstr,k).ge.0.0_r8) THEN
573 grade(jstr-1)=0.0_r8
574 kae(jstr-1)=0.0_r8
575 END IF
576 END IF
577 IF (domain(ng)%Northern_Edge(tile)) THEN
578 IF (hvom(i,jend+1,k).lt.0.0_r8) THEN
579 grade(jend+2)=0.0_r8
580 kae(jend+2)=0.0_r8
581 END IF
582 END IF
583 END IF
584 DO j=jstr,jend+1
585 IF (kae(j).le.eps1) THEN
586 okae(j)=0.0_r8
587 ELSE
588 okae(j)=1.0_r8/max(kae(j),eps1)
589 END IF
590 IF (hvom(i,j,k).ge.0.0_r8) THEN
591 IF (abs(grade(j)).le.eps1) THEN
592 rd=0.0_r8
593 rkad=0.0_r8
594 ELSE
595 rd=grade(j-1)/grade(j)
596 rkad=kae(j-1)*okae(j)
597 END IF
598 a1= cc1*kae(j)+cc2-cc3*okae(j)
599 b1=-cc1*kae(j)+cc2+cc3*okae(j)
600 betad=a1+b1*rd
601 cff=0.5_r8*max(0.0_r8, &
602 & min(2.0_r8, 2.0_r8*rd*rkad, betad))* &
603 & grade(j)*kae(j)
604# ifdef MASKING
605 jj=max(j-2,0)
606 cff=cff*rmask(i,jj)
607# endif
608 sw_eta=t(i,j-1,k,3,itrc)+cff
609 ELSE
610 IF (abs(grade(j)).le.eps1) THEN
611 ru=0.0_r8
612 rkau=0.0_r8
613 ELSE
614 ru=grade(j+1)/grade(j)
615 rkau=kae(j+1)*okae(j)
616 END IF
617 a1= cc1*kae(j)+cc2-cc3*okae(j)
618 b1=-cc1*kae(j)+cc2+cc3*okae(j)
619 betau=a1+b1*ru
620 cff=0.5*max(0.0_r8, &
621 & min(2.0_r8, 2.0_r8*ru*rkau, betau))* &
622 & grade(j)*kae(j)
623# ifdef MASKING
624 jj=min(j+1,mm(ng)+1)
625 cff=cff*rmask(i,jj)
626# endif
627 sw_eta=t(i,j,k,3,itrc)-cff
628 END IF
629 fe(i,j)=sw_eta*hvom(i,j,k)
630 END DO
631 END DO
632!
633 ELSE IF ((hadvection(itrc,ng)%AKIMA4).or. &
634 & (hadvection(itrc,ng)%CENTERED4).or. &
635 & (hadvection(itrc,ng)%SPLIT_U3).or. &
636 & (hadvection(itrc,ng)%UPSTREAM3)) THEN
637!
638! Fourth-order Akima, fourth-order centered differences, or third-order
639! upstream-biased horizontal advective fluxes.
640!
641 DO j=jstr,jend
642 DO i=istrm1,iendp2
643 fx(i,j)=t(i ,j,k,3,itrc)- &
644 & t(i-1,j,k,3,itrc)
645# ifdef MASKING
646 fx(i,j)=fx(i,j)*umask(i,j)
647# endif
648 END DO
649 END DO
650 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
651 IF (domain(ng)%Western_Edge(tile)) THEN
652 DO j=jstr,jend
653 fx(istr-1,j)=fx(istr,j)
654 END DO
655 END IF
656 END IF
657 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
658 IF (domain(ng)%Eastern_Edge(tile)) THEN
659 DO j=jstr,jend
660 fx(iend+2,j)=fx(iend+1,j)
661 END DO
662 END IF
663 END IF
664!
665 DO j=jstr,jend
666 DO i=istr-1,iend+1
667 IF (hadvection(itrc,ng)%UPSTREAM3) THEN
668 curv(i,j)=fx(i+1,j)-fx(i,j)
669 ELSE IF (hadvection(itrc,ng)%AKIMA4) THEN
670 cff=2.0_r8*fx(i+1,j)*fx(i,j)
671 IF (cff.gt.eps) THEN
672 grad(i,j)=cff/(fx(i+1,j)+fx(i,j))
673 ELSE
674 grad(i,j)=0.0_r8
675 END IF
676 ELSE IF ((hadvection(itrc,ng)%CENTERED4).or. &
677 & (hadvection(itrc,ng)%SPLIT_U3)) THEN
678 grad(i,j)=0.5_r8*(fx(i+1,j)+fx(i,j))
679 END IF
680 END DO
681 END DO
682!
683 cff1=1.0_r8/6.0_r8
684 cff2=1.0_r8/3.0_r8
685 DO j=jstr,jend
686 DO i=istr,iend+1
687 IF (hadvection(itrc,ng)%UPSTREAM3) THEN
688 fx(i,j)=huon(i,j,k)*0.5_r8* &
689 & (t(i-1,j,k,3,itrc)+ &
690 & t(i ,j,k,3,itrc))- &
691 & cff1*(curv(i-1,j)*max(huon(i,j,k),0.0_r8)+ &
692 & curv(i ,j)*min(huon(i,j,k),0.0_r8))
693 ELSE IF ((hadvection(itrc,ng)%AKIMA4).or. &
694 & (hadvection(itrc,ng)%CENTERED4).or. &
695 & (hadvection(itrc,ng)%SPLIT_U3)) THEN
696 fx(i,j)=huon(i,j,k)*0.5_r8* &
697 & (t(i-1,j,k,3,itrc)+ &
698 & t(i ,j,k,3,itrc)- &
699 & cff2*(grad(i ,j)- &
700 & grad(i-1,j)))
701 END IF
702 END DO
703 END DO
704!
705 DO j=jstrm1,jendp2
706 DO i=istr,iend
707 fe(i,j)=t(i,j ,k,3,itrc)- &
708 & t(i,j-1,k,3,itrc)
709# ifdef MASKING
710 fe(i,j)=fe(i,j)*vmask(i,j)
711# endif
712 END DO
713 END DO
714 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
715 IF (domain(ng)%Southern_Edge(tile)) THEN
716 DO i=istr,iend
717 fe(i,jstr-1)=fe(i,jstr)
718 END DO
719 END IF
720 END IF
721 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
722 IF (domain(ng)%Northern_Edge(tile)) THEN
723 DO i=istr,iend
724 fe(i,jend+2)=fe(i,jend+1)
725 END DO
726 END IF
727 END IF
728!
729 DO j=jstr-1,jend+1
730 DO i=istr,iend
731 IF (hadvection(itrc,ng)%UPSTREAM3) THEN
732 curv(i,j)=fe(i,j+1)-fe(i,j)
733 ELSE IF (hadvection(itrc,ng)%AKIMA4) THEN
734 cff=2.0_r8*fe(i,j+1)*fe(i,j)
735 IF (cff.gt.eps) THEN
736 grad(i,j)=cff/(fe(i,j+1)+fe(i,j))
737 ELSE
738 grad(i,j)=0.0_r8
739 END IF
740 ELSE IF ((hadvection(itrc,ng)%CENTERED4).or. &
741 & (hadvection(itrc,ng)%SPLIT_U3)) THEN
742 grad(i,j)=0.5_r8*(fe(i,j+1)+fe(i,j))
743 END IF
744 END DO
745 END DO
746!
747 cff1=1.0_r8/6.0_r8
748 cff2=1.0_r8/3.0_r8
749 DO j=jstr,jend+1
750 DO i=istr,iend
751 IF (hadvection(itrc,ng)%UPSTREAM3) THEN
752 fe(i,j)=hvom(i,j,k)*0.5_r8* &
753 & (t(i,j-1,k,3,itrc)+ &
754 & t(i,j ,k,3,itrc))- &
755 & cff1*(curv(i,j-1)*max(hvom(i,j,k),0.0_r8)+ &
756 & curv(i,j )*min(hvom(i,j,k),0.0_r8))
757 ELSE IF ((hadvection(itrc,ng)%AKIMA4).or. &
758 & (hadvection(itrc,ng)%CENTERED4).or. &
759 & (hadvection(itrc,ng)%SPLIT_U3)) THEN
760 fe(i,j)=hvom(i,j,k)*0.5_r8* &
761 & (t(i,j-1,k,3,itrc)+ &
762 & t(i,j ,k,3,itrc)- &
763 & cff2*(grad(i,j )- &
764 & grad(i,j-1)))
765 END IF
766 END DO
767 END DO
768 END IF hadv_flux
769!
770! Apply tracers point sources to the horizontal advection terms,
771! if any.
772!
773! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
774! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
775!
776 IF (luvsrc(ng)) THEN
777 DO is=1,nsrc(ng)
778 isrc=sources(ng)%Isrc(is)
779 jsrc=sources(ng)%Jsrc(is)
780 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
781 IF ((hadvection(itrc,ng)%MPDATA).or. &
782 & (hadvection(itrc,ng)%HSIMT)) THEN
783 lapplysrc=(istrum2.le.isrc).and. &
784 & (isrc.le.iendp3).and. &
785 & (jstrvm2.le.jsrc).and. &
786 & (jsrc.le.jendp2i)
787 ELSE
788 lapplysrc=(istr.le.isrc).and. &
789 & (isrc.le.iend+1).and. &
790 & (jstr.le.jsrc).and. &
791 & (jsrc.le.jend)
792 END IF
793 IF (lapplysrc) THEN
794 IF (ltracersrc(itrc,ng)) THEN
795 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
796 & sources(ng)%Tsrc(is,k,itrc)
797# ifdef MASKING
798 ELSE
799 IF ((rmask(isrc ,jsrc).eq.0.0_r8).and. &
800 & (rmask(isrc-1,jsrc).eq.1.0_r8)) THEN
801 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
802 & t(isrc-1,jsrc,k,3,itrc)
803 ELSE IF ((rmask(isrc ,jsrc).eq.1.0_r8).and. &
804 & (rmask(isrc-1,jsrc).eq.0.0_r8)) THEN
805 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
806 & t(isrc ,jsrc,k,3,itrc)
807 END IF
808# endif
809 END IF
810 END IF
811 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
812 IF ((hadvection(itrc,ng)%MPDATA).or. &
813 & (hadvection(itrc,ng)%HSIMT)) THEN
814 lapplysrc=(istrum2.le.isrc).and. &
815 & (isrc.le.iendp2i).and. &
816 & (jstrvm2.le.jsrc).and. &
817 & (jsrc.le.jendp3)
818 ELSE
819 lapplysrc=(istr.le.isrc).and. &
820 & (isrc.le.iend).and. &
821 & (jstr.le.jsrc).and. &
822 & (jsrc.le.jend+1)
823 END IF
824 IF (lapplysrc) THEN
825 IF (ltracersrc(itrc,ng)) THEN
826 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
827 & sources(ng)%Tsrc(is,k,itrc)
828# ifdef MASKING
829 ELSE
830 IF ((rmask(isrc,jsrc ).eq.0.0_r8).and. &
831 & (rmask(isrc,jsrc-1).eq.1.0_r8)) THEN
832 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
833 & t(isrc,jsrc-1,k,3,itrc)
834 ELSE IF ((rmask(isrc,jsrc ).eq.1.0_r8).and. &
835 & (rmask(isrc,jsrc-1).eq.0.0_r8)) THEN
836 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
837 & t(isrc,jsrc ,k,3,itrc)
838 END IF
839# endif
840 END IF
841 END IF
842 END IF
843 END DO
844 END IF
845
846# if defined NESTING && !defined ONE_WAY
847!
848! If refinement grids, extract tracer horizontal advection fluxes
849! (Hz*u*T/n, Hz*v*T/m) at the grid contact boundary (physical
850! domain perimeter) to be used in two-way nesting.
851!
852 IF (refinedgrid(ng)) THEN
853 DO cr=1,ncontact
854 dg=rcontact(cr)%donor_grid
855 rg=rcontact(cr)%receiver_grid
856 IF (ng.eq.rg) THEN
857 CALL bry_fluxes (dg, rg, cr, inlm, tile, &
858 & imins, imaxs, jmins, jmaxs, &
859 & ilb, iub, jlb, jub, &
860 & dt(ng), fx, fe, &
861 & bry_contact(iwest, cr)%Tflux(:,k,itrc), &
862 & bry_contact(ieast, cr)%Tflux(:,k,itrc), &
863 & bry_contact(isouth,cr)%Tflux(:,k,itrc), &
864 & bry_contact(inorth,cr)%Tflux(:,k,itrc))
865 END IF
866 END DO
867 END IF
868# endif
869!
870! If MPDATA, time-step horizontal advection for intermediate diffusive
871! tracer, Ta (m Tunits).
872!
873 hadv_stepping : IF (hadvection(itrc,ng)%MPDATA) THEN
874 DO j=jstrvm2,jendp2i
875 DO i=istrum2,iendp2i
876 cff=dt(ng)*pm(i,j)*pn(i,j)
877 cff1=cff*(fx(i+1,j)-fx(i,j))
878 cff2=cff*(fe(i,j+1)-fe(i,j))
879 cff3=cff1+cff2
880 ta(i,j,k,itrc)=t(i,j,k,nnew,itrc)-cff3
881# ifdef DIAGNOSTICS_TS
882 dhadv(i,j,itxadv)=-cff1
883 dhadv(i,j,ityadv)=-cff2
884 dhadv(i,j,ithadv)=-cff3
885# endif
886 END DO
887 END DO
888# ifdef DIAGNOSTICS_TS
889!
890 DO j=jstr,jend
891 DO i=istr,iend
892 diatwrk(i,j,k,itrc,itxadv)=dhadv(i,j,itxadv)
893 diatwrk(i,j,k,itrc,ityadv)=dhadv(i,j,ityadv)
894 diatwrk(i,j,k,itrc,ithadv)=dhadv(i,j,ithadv)
895 END DO
896 END DO
897# endif
898!
899! OTHERWISE, time-step horizontal advection term. Advective fluxes
900! have units of Tunits m3/s. The new tracer has units of m Tunits.
901!
902 ELSE
903 DO j=jstr,jend
904 DO i=istr,iend
905 cff=dt(ng)*pm(i,j)*pn(i,j)
906 cff1=cff*(fx(i+1,j)-fx(i,j))
907 cff2=cff*(fe(i,j+1)-fe(i,j))
908 cff3=cff1+cff2
909 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff3
910# ifdef DIAGNOSTICS_TS
911 diatwrk(i,j,k,itrc,itxadv)=-cff1
912 diatwrk(i,j,k,itrc,ityadv)=-cff2
913 diatwrk(i,j,k,itrc,ithadv)=-cff3
914# endif
915 END DO
916 END DO
917 END IF hadv_stepping
918 END DO k_loop
919 END DO t_loop1
920!
921!-----------------------------------------------------------------------
922! Time-step vertical advection term.
923!-----------------------------------------------------------------------
924!
925 t_loop2 : DO itrc=1,nt(ng)
926 IF (vadvection(itrc,ng)%MPDATA) THEN
927 jmint=jstrvm2
928 jmaxt=jendp2i
929 ELSE
930 jmint=jstr
931 jmaxt=jend
932 END IF
933!
934 j_loop1 : DO j=jmint,jmaxt ! start pipelined J-loop
935!
936 vadv_flux : IF (vadvection(itrc,ng)%SPLINES) THEN
937!
938! Build conservative parabolic splines for the vertical derivatives
939! "FC" of the tracer. Then, the interfacial "FC" values are
940! converted to vertical advective flux (Tunits m3/s).
941!
942 DO i=istr,iend
943# ifdef NEUMANN
944 fc(i,0)=1.5_r8*t(i,j,1,3,itrc)
945 cf(i,1)=0.5_r8
946# else
947 fc(i,0)=2.0_r8*t(i,j,1,3,itrc)
948 cf(i,1)=1.0_r8
949# endif
950 END DO
951 DO k=1,n(ng)-1
952 DO i=istr,iend
953 cff=1.0_r8/(2.0_r8*hz(i,j,k)+ &
954 & hz(i,j,k+1)*(2.0_r8-cf(i,k)))
955 cf(i,k+1)=cff*hz(i,j,k)
956 fc(i,k)=cff*(3.0_r8*(hz(i,j,k )*t(i,j,k+1,3,itrc)+ &
957 & hz(i,j,k+1)*t(i,j,k ,3,itrc))- &
958 & hz(i,j,k+1)*fc(i,k-1))
959 END DO
960 END DO
961 DO i=istr,iend
962# ifdef NEUMANN
963 fc(i,n(ng))=(3.0_r8*t(i,j,n(ng),3,itrc)-fc(i,n(ng)-1))/ &
964 & (2.0_r8-cf(i,n(ng)))
965# else
966 fc(i,n(ng))=(2.0_r8*t(i,j,n(ng),3,itrc)-fc(i,n(ng)-1))/ &
967 & (1.0_r8-cf(i,n(ng)))
968# endif
969 END DO
970 DO k=n(ng)-1,0,-1
971 DO i=istr,iend
972 fc(i,k)=fc(i,k)-cf(i,k+1)*fc(i,k+1)
973# ifdef WEC_VF
974 fc(i,k+1)=(w(i,j,k+1)+w_stokes(i,j,k+1))*fc(i,k+1)
975# else
976 fc(i,k+1)=w(i,j,k+1)*fc(i,k+1)
977# endif
978 END DO
979 END DO
980 DO i=istr,iend
981 fc(i,n(ng))=0.0_r8
982 fc(i,0)=0.0_r8
983 END DO
984!
985 ELSE IF (vadvection(itrc,ng)%AKIMA4) THEN
986!
987! Fourth-order, Akima vertical advective flux (Tunits m3/s).
988!
989 DO k=1,n(ng)-1
990 DO i=istr,iend
991 fc(i,k)=t(i,j,k+1,3,itrc)- &
992 & t(i,j,k ,3,itrc)
993 END DO
994 END DO
995 DO i=istr,iend
996 fc(i,0)=fc(i,1)
997 fc(i,n(ng))=fc(i,n(ng)-1)
998 END DO
999 DO k=1,n(ng)
1000 DO i=istr,iend
1001 cff=2.0_r8*fc(i,k)*fc(i,k-1)
1002 IF (cff.gt.eps) THEN
1003 cf(i,k)=cff/(fc(i,k)+fc(i,k-1))
1004 ELSE
1005 cf(i,k)=0.0_r8
1006 END IF
1007 END DO
1008 END DO
1009 cff1=1.0_r8/3.0_r8
1010 DO k=1,n(ng)-1
1011 DO i=istr,iend
1012# ifdef WEC_VF
1013 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1014# else
1015 fc(i,k)=w(i,j,k)* &
1016# endif
1017 & 0.5_r8*(t(i,j,k ,3,itrc)+ &
1018 & t(i,j,k+1,3,itrc)- &
1019 & cff1*(cf(i,k+1)-cf(i,k)))
1020 END DO
1021 END DO
1022 DO i=istr,iend
1023 fc(i,0)=0.0_r8
1024 fc(i,n(ng))=0.0_r8
1025 END DO
1026!
1027 ELSE IF (vadvection(itrc,ng)%CENTERED2) THEN
1028!
1029! Second-order, central differences vertical advective flux
1030! (Tunits m3/s).
1031!
1032 DO k=1,n(ng)-1
1033 DO i=istr,iend
1034# ifdef WEC_VF
1035 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1036# else
1037 fc(i,k)=w(i,j,k)* &
1038# endif
1039 & 0.5_r8*(t(i,j,k ,3,itrc)+ &
1040 & t(i,j,k+1,3,itrc))
1041 END DO
1042 END DO
1043 DO i=istr,iend
1044 fc(i,0)=0.0_r8
1045 fc(i,n(ng))=0.0_r8
1046 END DO
1047!
1048 ELSE IF (vadvection(itrc,ng)%MPDATA) THEN
1049!
1050! First_order, upstream differences vertical advective flux
1051! (Tunits m3/s).
1052!
1053 DO i=istrum2,iendp2i
1054 DO k=1,n(ng)-1
1055# ifdef WEC_VF
1056 cff1=max(w(i,j,k)+w_stokes(i,j,k),0.0_r8)
1057 cff2=min(w(i,j,k)+w_stokes(i,j,k),0.0_r8)
1058# else
1059 cff1=max(w(i,j,k),0.0_r8)
1060 cff2=min(w(i,j,k),0.0_r8)
1061# endif
1062 fc(i,k)=cff1*t(i,j,k ,3,itrc)+ &
1063 & cff2*t(i,j,k+1,3,itrc)
1064 END DO
1065 fc(i,0)=0.0_r8
1066 fc(i,n(ng))=0.0_r8
1067 END DO
1068!
1069 ELSE IF (vadvection(itrc,ng)%HSIMT) THEN
1070!
1071! Third High-order Spatial Interpolation at the Middle Temporal level
1072! (HSIMT; Wu and Zhu, 2010) with a Total Variation Diminishing (TVD)
1073! limiter vertical advection flux (Tunits m3/s).
1074!
1075 DO i=istr,iend
1076 kaz(0)=0.0_r8
1077 okaz(0)=0.0_r8
1078 gradz(0)=0.0_r8
1079 DO k=1,n(ng)-1
1080 cff=pm(i,j)*pn(i,j)*dt(ng)
1081# ifdef WEC_VF
1082 kaz(k)=1.0_r8-abs(cff*(w(i,j,k)+w_stokes(i,j,k))/ &
1083 & (z_r(i,j,k+1)-z_r(i,j,k)))
1084# else
1085 kaz(k)=1.0_r8-abs(cff*w(i,j,k)/ &
1086 & (z_r(i,j,k+1)-z_r(i,j,k)))
1087# endif
1088 okaz(k)=1.0_r8/kaz(k)
1089 gradz(k)=t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)
1090 END DO
1091 kaz(n(ng))=0.0_r8
1092 okaz(n(ng))=0.0_r8
1093 gradz(n(ng))=0.0_r8
1094!
1095 DO k=1,n(ng)-1
1096# ifdef WEC_VF
1097 cff1=w(i,j,k)+w_stokes(i,j,k)
1098# else
1099 cff1=w(i,j,k)
1100# endif
1101 IF ((k.eq.1).and.(cff1.ge.0.0_r8)) THEN
1102 fc(i,k)=cff1*t(i,j,k,3,itrc)
1103 ELSE IF ((k.eq.n(ng)-1).and.(cff1.lt.0.0_r8)) THEN
1104 fc(i,k)=cff1*t(i,j,k+1,3,itrc)
1105 ELSE
1106 IF (cff1.ge.0) THEN
1107 IF (abs(gradz(k)).le.eps1) THEN
1108 rd=0.0_r8
1109 rkad=0.0_r8
1110 ELSE
1111 rd=gradz(k-1)/gradz(k)
1112 rkad=kaz(k-1)*okaz(k)
1113 END IF
1114 a1= cc1*kaz(k)+cc2-cc3*okaz(k)
1115 b1=-cc1*kaz(k)+cc2+cc3*okaz(k)
1116 betad=a1+b1*rd
1117 cff=0.5_r8*max(0.0_r8, &
1118 & min(2.0_r8, 2.0_r8*rd*rkad, betad))* &
1119 & gradz(k)*kaz(k)
1120 sw=t(i,j,k,3,itrc)+cff
1121 ELSE
1122 IF (abs(gradz(k)).le.eps1) THEN
1123 ru=0.0_r8
1124 rkau=0.0_r8
1125 ELSE
1126 ru=gradz(k+1)/gradz(k)
1127 rkau=kaz(k+1)*okaz(k)
1128 END IF
1129 a1= cc1*kaz(k)+cc2-cc3*okaz(k)
1130 b1=-cc1*kaz(k)+cc2+cc3*okaz(k)
1131 betau=a1+b1*ru
1132 cff=0.5_r8*max(0.0_r8, &
1133 & min(2.0_r8, 2.0_r8*ru*rkau, betau))* &
1134 & gradz(k)*kaz(k)
1135 sw=t(i,j,k+1,3,itrc)-cff
1136 END IF
1137 fc(i,k)=cff1*sw
1138 END IF
1139 END DO
1140 fc(i,0)=0.0_r8
1141 fc(i,n(ng))=0.0_r8
1142 END DO
1143!
1144 ELSE IF ((vadvection(itrc,ng)%CENTERED4).or. &
1145 & (vadvection(itrc,ng)%SPLIT_U3)) THEN
1146!
1147! Fourth-order, central differences vertical advective flux
1148! (Tunits m3/s).
1149!
1150 cff1=0.5_r8
1151 cff2=7.0_r8/12.0_r8
1152 cff3=1.0_r8/12.0_r8
1153 DO k=2,n(ng)-2
1154 DO i=istr,iend
1155# ifdef WEC_VF
1156 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1157# else
1158 fc(i,k)=w(i,j,k)* &
1159# endif
1160 & (cff2*(t(i,j,k ,3,itrc)+ &
1161 & t(i,j,k+1,3,itrc))- &
1162 & cff3*(t(i,j,k-1,3,itrc)+ &
1163 & t(i,j,k+2,3,itrc)))
1164 END DO
1165 END DO
1166 DO i=istr,iend
1167 fc(i,0)=0.0_r8
1168# ifdef WEC_VF
1169 fc(i,1)=(w(i,j,1)+w_stokes(i,j,1))* &
1170# else
1171 fc(i,1)=w(i,j,1)* &
1172# endif
1173 & (cff1*t(i,j,1,3,itrc)+ &
1174 & cff2*t(i,j,2,3,itrc)- &
1175 & cff3*t(i,j,3,3,itrc))
1176# ifdef WEC_VF
1177 fc(i,n(ng)-1)=(w(i,j,n(ng)-1)+w_stokes(i,j,n(ng)-1))* &
1178# else
1179 fc(i,n(ng)-1)=w(i,j,n(ng)-1)* &
1180# endif
1181 & (cff1*t(i,j,n(ng) ,3,itrc)+ &
1182 & cff2*t(i,j,n(ng)-1,3,itrc)- &
1183 & cff3*t(i,j,n(ng)-2,3,itrc))
1184 fc(i,n(ng))=0.0_r8
1185 END DO
1186 END IF vadv_flux
1187!
1188! If MPDATA and cell-centered point source (LwSrc), augment the
1189! intermediate tracer, Ta, with the tracer divergence. For other
1190! advection schemes LwSrc is applied after the vertical advection
1191! is completed (J. Wilkin).
1192!
1193! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
1194!
1195 IF (lwsrc(ng)) THEN
1196 IF (vadvection(itrc,ng)%MPDATA) THEN
1197 DO is=1,nsrc(ng)
1198 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
1199 isrc=sources(ng)%Isrc(is)
1200 jsrc=sources(ng)%Jsrc(is)
1201 IF (((istr.le.isrc).and.(isrc.le.iend+1)).and. &
1202 & ((jstr.le.jsrc).and.(jsrc.le.jend+1)).and. &
1203 & (j.eq.jsrc)) THEN
1204 DO k=1,n(ng)
1205 cff=dt(ng)*pm(isrc,jsrc)*pn(isrc,jsrc)
1206 IF (ltracersrc(itrc,ng)) THEN
1207 cff3=sources(ng)%Tsrc(is,k,itrc)
1208 ELSE
1209 cff3=t(isrc,jsrc,k,3,itrc)
1210 END IF
1211 ta(isrc,jsrc,k,itrc)=ta(isrc,jsrc,k,itrc)+ &
1212 & cff*sources(ng)%Qsrc(is,k)* &
1213 & cff3
1214 END DO
1215 END IF
1216 END IF
1217 END DO
1218 END IF
1219 END IF
1220
1221# ifdef SED_MORPH
1222!
1223! If MPDATA, augment the intermediate tracer, Ta, with the tracer
1224! divergence. For other advection schemes the bed change is applied
1225! after the vertical advection is completed (similar to Lwsrc).
1226! Use 1/N to distribute the bed change in each water column layer.
1227!
1228 IF (vadvection(itrc,ng)%MPDATA) THEN
1229 cff1=1.0_r8/n(ng)
1230 DO k=1,n(ng)
1231 DO i=istrum2,iendp2i
1232 cff=cff1*(bed_thick(i,j,nstp)-bed_thick(i,j,3))
1233 cff3=t(i,j,k,3,itrc)
1234 ta(i,j,k,itrc)=ta(i,j,k,itrc)-cff*cff3
1235 END DO
1236 END DO
1237 END IF
1238# endif
1239!
1240! If MPDATA, time-step vertical advection for intermediate diffusive
1241! tracer, Ta (Tunits).
1242# ifdef DIAGNOSTICS_TS
1243! Convert units of tracer diagnostic terms to Tunits.
1244# endif
1245!
1246 vadv_stepping : IF (vadvection(itrc,ng)%MPDATA) THEN
1247 DO i=istrum2,iendp2i
1248 cf(i,0)=dt(ng)*pm(i,j)*pn(i,j)
1249 END DO
1250 DO k=1,n(ng)
1251 DO i=istrum2,iendp2i
1252 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1253 ta(i,j,k,itrc)=(ta(i,j,k,itrc)-cff1)*ohz(i,j,k)
1254# ifdef DIAGNOSTICS_TS
1255 dvadv(i,j,k,itrc)=-cff1
1256# endif
1257 END DO
1258 END DO
1259
1260# ifdef OMEGA_IMPLICIT
1261!
1262! If MPDATA, compute off-diagonal coefficients FC [dt*Wi*pm*pn] for the
1263! implicit vertical advection term located at horizontal RHO-points and
1264! vertical W-points. Also, set FC at the top and bottom levels.
1265!
1266! It needs to be done after the Wexplicit and before the anti-diffusion
1267! steps.
1268!
1269 cff=dt(ng)
1270 DO k=1,n(ng)-1
1271 DO i=istr,iend
1272 cff1=cff*pm(i,j)*pn(i,j)
1273 fcmax(i,k)=max(wi(i,j,k),0.0_r8)*cff1
1274 fcmin(i,k)=min(wi(i,j,k),0.0_r8)*cff1
1275 END DO
1276 END DO
1277 DO i=istr,iend
1278 fcmax(i,0)=0.0_r8
1279 fcmin(i,0)=0.0_r8
1280 fcmax(i,n(ng))=0.0_r8
1281 fcmin(i,n(ng))=0.0_r8
1282 END DO
1283!
1284! Compute diagonal matrix coefficients BC and load right-hand-side
1285! terms for the tracer equation into DC.
1286!
1287 DO k=1,n(ng)
1288 DO i=istr,iend
1289 bc(i,k)=hz(i,j,k)+fcmax(i,k)-fcmin(i,k-1)
1290 dc(i,k)=ta(i,j,k,itrc)*hz(i,j,k)
1291 END DO
1292 END DO
1293!
1294! Solve the tridiagonal system.
1295!
1296 DO i=istr,iend
1297 cff=1.0_r8/bc(i,1)
1298 cf(i,1)=cff*fcmin(i,1)
1299 dc(i,1)=cff*dc(i,1)
1300 END DO
1301 DO k=2,n(ng)-1
1302 DO i=istr,iend
1303 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
1304 cf(i,k)=cff*fcmin(i,k)
1305 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
1306 END DO
1307 END DO
1308!
1309! Compute new solution by back substitution.
1310!
1311 DO i=istr,iend
1312# ifdef DIAGNOSTICS_TS
1313 cff1=ta(i,j,n(ng),itrc)*ohz(i,j,n(ng))
1314# endif
1315 cff=1.0_r8/(bc(i,n(ng))+fcmax(i,n(ng)-1)*cf(i,n(ng)-1))
1316 dc(i,n(ng))=cff*(dc(i,n(ng))+ &
1317 & fcmax(i,n(ng)-1)*dc(i,n(ng)-1))
1318 ta(i,j,n(ng),itrc)=dc(i,n(ng))
1319# ifdef DIAGNOSTICS_TS
1320 diatwrk(i,j,n(ng),itrc,itvadv)= &
1321 & diatwrk(i,j,n(ng),itrc,itvadv)+ &
1322 & ta(i,j,n(ng),itrc)-cff1
1323# endif
1324 END DO
1325 DO k=n(ng)-1,1,-1
1326 DO i=istr,iend
1327# ifdef DIAGNOSTICS_TS
1328 cff1=ta(i,j,k,itrc)*ohz(i,j,k)
1329# endif
1330 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1331 ta(i,j,k,itrc)=dc(i,k)
1332# ifdef DIAGNOSTICS_TS
1333 diatwrk(i,j,k,itrc,itvadv)=diatwrk(i,j,k,itrc,itvadv)+&
1334 & ta(i,j,k,itrc)-cff1
1335# endif
1336 END DO
1337 END DO
1338# endif
1339!
1340! OTHERWISE, time-step vertical advection term (m Tunits, or Tunits if
1341! conservative, parabolic splines diffusion).
1342# ifdef DIAGNOSTICS_TS
1343! Convert units of tracer diagnostic terms to Tunits.
1344# endif
1345!
1346 ELSE
1347 DO i=istr,iend
1348 cf(i,0)=dt(ng)*pm(i,j)*pn(i,j)
1349 END DO
1350 DO k=1,n(ng)
1351 DO i=istr,iend
1352 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1353 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
1354# ifdef SPLINES_VDIFF
1355 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1356# endif
1357# ifdef DIAGNOSTICS_TS
1358 diatwrk(i,j,k,itrc,itvadv)=-cff1
1359 DO idiag=1,ndt
1360 diatwrk(i,j,k,itrc,idiag)=diatwrk(i,j,k,itrc,idiag)* &
1361 & ohz(i,j,k)
1362 END DO
1363# endif
1364 END DO
1365 END DO
1366 END IF vadv_stepping
1367 END DO j_loop1
1368 END DO t_loop2
1369!
1370!-----------------------------------------------------------------------
1371! Compute anti-diffusive velocities to corrected advected tracers
1372! using MPDATA recursive method. Notice that pipelined J-loop ended.
1373!-----------------------------------------------------------------------
1374!
1375 t_loop3 : DO itrc=1,nt(ng)
1376 mpdata : IF ((hadvection(itrc,ng)%MPDATA).and. &
1377 & (vadvection(itrc,ng)%MPDATA)) THEN
1378 CALL mpdata_adiff_tile (ng, tile, &
1379 & lbi, ubi, lbj, ubj, &
1380 & imins, imaxs, jmins, jmaxs, &
1381# ifdef MASKING
1382 & rmask, umask, vmask, &
1383# endif
1384# ifdef WET_DRY
1385 & rmask_wet, umask_wet, vmask_wet, &
1386# endif
1387 & pm, pn, omn, om_u, on_v, &
1388 & z_r, ohz, &
1389 & huon, hvom, w, &
1390# ifdef WEC_VF
1391 & w_stokes, &
1392# endif
1393# ifdef OMEGA_IMPLICIT
1394 & wi, &
1395# endif
1396 & t(:,:,:,3,itrc), &
1397 & ta(:,:,:,itrc), ua, va, wa)
1398!
1399! Compute anti-diffusive corrected advection fluxes.
1400!
1401 DO k=1,n(ng)
1402 DO j=jstr,jend
1403 DO i=istr,iend+1
1404 cff1=max(ua(i,j,k),0.0_r8)
1405 cff2=min(ua(i,j,k),0.0_r8)
1406 fx(i,j)=(cff1*ta(i-1,j,k,itrc)+ &
1407 & cff2*ta(i ,j,k,itrc))* &
1408 & 0.5_r8*(hz(i,j,k)+hz(i-1,j,k))*on_u(i,j)
1409 END DO
1410 END DO
1411 DO j=jstr,jend+1
1412 DO i=istr,iend
1413 cff1=max(va(i,j,k),0.0_r8)
1414 cff2=min(va(i,j,k),0.0_r8)
1415 fe(i,j)=(cff1*ta(i,j-1,k,itrc)+ &
1416 & cff2*ta(i,j ,k,itrc))* &
1417 & 0.5_r8*(hz(i,j,k)+hz(i,j-1,k))*om_v(i,j)
1418 END DO
1419 END DO
1420!
1421! Time-step corrected horizontal advection (Tunits m).
1422!
1423 DO j=jstr,jend
1424 DO i=istr,iend
1425 cff=dt(ng)*pm(i,j)*pn(i,j)
1426 cff1=cff*(fx(i+1,j)-fx(i,j))
1427 cff2=cff*(fe(i,j+1)-fe(i,j))
1428 cff3=cff1+cff2
1429 t(i,j,k,nnew,itrc)=ta(i,j,k,itrc)*hz(i,j,k)-cff3
1430# ifdef DIAGNOSTICS_TS
1431 diatwrk(i,j,k,itrc,itxadv)=diatwrk(i,j,k,itrc,itxadv)- &
1432 & cff1
1433 diatwrk(i,j,k,itrc,ityadv)=diatwrk(i,j,k,itrc,ityadv)- &
1434 & cff2
1435 diatwrk(i,j,k,itrc,ithadv)=diatwrk(i,j,k,itrc,ithadv)- &
1436 & cff3
1437# endif
1438 END DO
1439 END DO
1440 END DO
1441!
1442! Compute anti-diffusive corrected vertical advection flux.
1443!
1444 DO j=jstr,jend
1445 DO k=1,n(ng)-1
1446 DO i=istr,iend
1447 cff1=max(wa(i,j,k),0.0_r8)
1448 cff2=min(wa(i,j,k),0.0_r8)
1449 fc(i,k)=cff1*ta(i,j,k ,itrc)+ &
1450 & cff2*ta(i,j,k+1,itrc)
1451 END DO
1452 END DO
1453 DO i=istr,iend
1454 fc(i,0)=0.0_r8
1455 fc(i,n(ng))=0.0_r8
1456 END DO
1457!
1458! Time-step corrected vertical advection (Tunits).
1459# ifdef DIAGNOSTICS_TS
1460! Convert units of tracer diagnostic terms to Tunits.
1461# endif
1462!
1463 DO i=istr,iend
1464 cf(i,0)=dt(ng)*pm(i,j)*pn(i,j)
1465 END DO
1466 DO k=1,n(ng)
1467 DO i=istr,iend
1468 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1469 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
1470# ifdef DIAGNOSTICS_TS
1471 diatwrk(i,j,k,itrc,itvadv)=dvadv(i,j,k,itrc)- &
1472 & cff1
1473 DO idiag=1,ndt
1474 diatwrk(i,j,k,itrc,idiag)=diatwrk(i,j,k,itrc,idiag)* &
1475 & ohz(i,j,k)
1476 END DO
1477# endif
1478 END DO
1479 END DO
1480 END DO
1481 END IF mpdata
1482 END DO t_loop3
1483!
1484!-----------------------------------------------------------------------
1485! Add tracer divergence due to cell-centered (LwSrc) point sources.
1486!-----------------------------------------------------------------------
1487!
1488! When LTracerSrc is .true. the inflowing concentration is Tsrc.
1489! When LtracerSrc is .false. we add tracer mass to compensate for the
1490! added volume to keep the receiving cell concentration unchanged.
1491! J. Levin (Jupiter Intelligence Inc.) and J. Wilkin
1492!
1493! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
1494!
1495 IF (lwsrc(ng)) THEN
1496 DO itrc=1,nt(ng)
1497 IF (.not.((hadvection(itrc,ng)%MPDATA).and. &
1498 & (vadvection(itrc,ng)%MPDATA))) THEN
1499 DO is=1,nsrc(ng)
1500 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
1501 isrc=sources(ng)%Isrc(is)
1502 jsrc=sources(ng)%Jsrc(is)
1503 IF (((istr.le.isrc).and.(isrc.le.iend+1)).and. &
1504 & ((jstr.le.jsrc).and.(jsrc.le.jend+1))) THEN
1505 DO k=1,n(ng)
1506 cff=dt(ng)*pm(isrc,jsrc)*pn(isrc,jsrc)
1507# ifdef SPLINES_VDIFF
1508 cff=cff*ohz(isrc,jsrc,k)
1509# endif
1510 IF (ltracersrc(itrc,ng)) THEN
1511 cff3=sources(ng)%Tsrc(is,k,itrc)
1512 ELSE
1513 cff3=t(isrc,jsrc,k,3,itrc)
1514 END IF
1515 t(isrc,jsrc,k,nnew,itrc)=t(isrc,jsrc,k,nnew,itrc)+ &
1516 & cff*sources(ng)%Qsrc(is,k)*&
1517 & cff3
1518 END DO
1519 END IF
1520 END IF
1521 END DO
1522 END IF
1523 END DO
1524 END IF
1525
1526# if defined SEDIMENT && defined SED_MORPH
1527!
1528! If not MPDATA, add tracer mass to compensate for the added volume to
1529! keep the receiving cell concentration unchanged.
1530! Use 1/N to distribute the bed change in each water column layer.
1531!
1532 cff1=1.0_r8/n(ng)
1533 DO itrc=1,nt(ng)
1534 IF (.not.((hadvection(itrc,ng)%MPDATA).and. &
1535 & (vadvection(itrc,ng)%MPDATA))) THEN
1536 DO j=jstr,jend
1537 DO k=1,n(ng)
1538 DO i=istr,iend
1539 cff3=t(i,j,k,3,itrc)
1540# ifdef SPLINES_VDIFF
1541 cff3=cff3*ohz(i,j,k)
1542# endif
1543 cff=cff1*(bed_thick(i,j,nstp)-bed_thick(i,j,3))
1544 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)- &
1545 & cff*cff3
1546 END DO
1547 END DO
1548 END DO
1549 END IF
1550 END DO
1551# endif
1552
1553# ifdef OMEGA_IMPLICIT
1554!
1555!-----------------------------------------------------------------------
1556! Add adaptive, Courant-number based implicit vertical advection term.
1557!-----------------------------------------------------------------------
1558!
1559! Compute off-diagonal coefficients FC [dt*Wi*pm*pn] for the implicit
1560! vertical advection term located at horizontal RHO-points and vertical
1561! W-points. Also, set FC at the top and bottom levels.
1562!
1563! It needs to be before the Diffusion.
1564!
1565 DO j=jstr,jend
1566 DO itrc=1,nt(ng)
1567 IF (.not.vadvection(itrc,ng)%MPDATA) THEN
1568 cff=dt(ng)
1569 DO k=1,n(ng)-1
1570 DO i=istr,iend
1571 cff1=cff*pm(i,j)*pn(i,j)
1572 fcmax(i,k)=max(wi(i,j,k),0.0_r8)*cff1
1573 fcmin(i,k)=min(wi(i,j,k),0.0_r8)*cff1
1574 END DO
1575 END DO
1576 DO i=istr,iend
1577 fcmax(i,0)=0.0_r8
1578 fcmin(i,0)=0.0_r8
1579 fcmax(i,n(ng))=0.0_r8
1580 fcmin(i,n(ng))=0.0_r8
1581 END DO
1582!
1583! Compute diagonal matrix coefficients BC and load right-hand-side
1584! terms for the tracer equation into DC.
1585!
1586 DO k=1,n(ng)
1587 DO i=istr,iend
1588 bc(i,k)=hz(i,j,k)+fcmax(i,k)-fcmin(i,k-1)
1589# ifdef SPLINES_VDIFF
1590 dc(i,k)=t(i,j,k,nnew,itrc)*hz(i,j,k)
1591# else
1592 dc(i,k)=t(i,j,k,nnew,itrc)
1593# endif
1594 END DO
1595 END DO
1596!
1597! Solve the tridiagonal system.
1598!
1599 DO i=istr,iend
1600 cff=1.0_r8/bc(i,1)
1601 cf(i,1)=cff*fcmin(i,1)
1602 dc(i,1)=cff*dc(i,1)
1603 END DO
1604 DO k=2,n(ng)-1
1605 DO i=istr,iend
1606 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
1607 cf(i,k)=cff*fcmin(i,k)
1608 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
1609 END DO
1610 END DO
1611!
1612! Compute new solution by back substitution.
1613!
1614 DO i=istr,iend
1615# ifdef DIAGNOSTICS_TS
1616 cff1=t(i,j,n(ng),nnew,itrc)*ohz(i,j,n(ng))
1617# endif
1618 cff=1.0_r8/(bc(i,n(ng))+ &
1619 & fcmax(i,n(ng)-1)*cf(i,n(ng)-1))
1620 dc(i,n(ng))=cff*(dc(i,n(ng))+ &
1621 & fcmax(i,n(ng)-1)*dc(i,n(ng)-1))
1622# ifdef SPLINES_VDIFF
1623 t(i,j,n(ng),nnew,itrc)=dc(i,n(ng))
1624# else
1625 t(i,j,n(ng),nnew,itrc)=dc(i,n(ng))*hz(i,j,n(ng))
1626# endif
1627# ifdef DIAGNOSTICS_TS
1628 diatwrk(i,j,n(ng),itrc,itvadv)=diatwrk(i,j,n(ng),itrc, &
1629 & itvadv)+ &
1630 & dc(i,n(ng))-cff1
1631# endif
1632 END DO
1633!
1634 DO k=n(ng)-1,1,-1
1635 DO i=istr,iend
1636# ifdef DIAGNOSTICS_TS
1637 cff1=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1638# endif
1639 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1640# ifdef SPLINES_VDIFF
1641 t(i,j,k,nnew,itrc)=dc(i,k)
1642# else
1643 t(i,j,k,nnew,itrc)=dc(i,k)*hz(i,j,k)
1644# endif
1645# ifdef DIAGNOSTICS_TS
1646 diatwrk(i,j,k,itrc,itvadv)=diatwrk(i,j,k,itrc,itvadv)+ &
1647 & dc(i,k)-cff1
1648# endif
1649 END DO
1650 END DO
1651 END IF
1652 END DO
1653 END DO
1654# endif
1655!
1656!-----------------------------------------------------------------------
1657! Time-step vertical diffusion term.
1658!-----------------------------------------------------------------------
1659!
1660 j_loop2 : DO j=jstr,jend ! start pipelined J-loop
1661 DO itrc=1,nt(ng)
1662 ltrc=min(nat,itrc)
1663
1664# ifdef SPLINES_VDIFF
1665 IF (.not.((hadvection(itrc,ng)%MPDATA).and. &
1666 & (vadvection(itrc,ng)%MPDATA))) THEN
1667!
1668! Use conservative, parabolic spline reconstruction of vertical
1669! diffusion derivatives. Then, time step vertical diffusion term
1670! implicitly.
1671!
1672 cff1=1.0_r8/6.0_r8
1673 DO k=1,n(ng)-1
1674 DO i=istr,iend
1675 fc(i,k)=cff1*hz(i,j,k )- &
1676 & dt(ng)*akt(i,j,k-1,ltrc)*ohz(i,j,k )
1677 cf(i,k)=cff1*hz(i,j,k+1)- &
1678 & dt(ng)*akt(i,j,k+1,ltrc)*ohz(i,j,k+1)
1679 END DO
1680 END DO
1681 DO i=istr,iend
1682 cf(i,0)=0.0_r8
1683 dc(i,0)=0.0_r8
1684 END DO
1685!
1686! LU decomposition and forward substitution.
1687!
1688 cff1=1.0_r8/3.0_r8
1689 DO k=1,n(ng)-1
1690 DO i=istr,iend
1691 bc(i,k)=cff1*(hz(i,j,k)+hz(i,j,k+1))+ &
1692 & dt(ng)*akt(i,j,k,ltrc)*(ohz(i,j,k)+ohz(i,j,k+1))
1693 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1694 cf(i,k)=cff*cf(i,k)
1695 dc(i,k)=cff*(t(i,j,k+1,nnew,itrc)-t(i,j,k,nnew,itrc)- &
1696 & fc(i,k)*dc(i,k-1))
1697 END DO
1698 END DO
1699!
1700! Backward substitution.
1701!
1702 DO i=istr,iend
1703 dc(i,n(ng))=0.0_r8
1704 END DO
1705 DO k=n(ng)-1,1,-1
1706 DO i=istr,iend
1707 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1708 END DO
1709 END DO
1710!
1711 DO k=1,n(ng)
1712 DO i=istr,iend
1713 dc(i,k)=dc(i,k)*akt(i,j,k,ltrc)
1714 cff1=dt(ng)*ohz(i,j,k)*(dc(i,k)-dc(i,k-1))
1715 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+cff1
1716# ifdef DIAGNOSTICS_TS
1717 diatwrk(i,j,k,itrc,itvdif)=diatwrk(i,j,k,itrc,itvdif)+ &
1718 & cff1
1719# endif
1720 END DO
1721 END DO
1722 ELSE
1723# endif
1724!
1725! Compute off-diagonal coefficients FC [lambda*dt*Akt/Hz] for the
1726! implicit vertical diffusion terms at future time step, located
1727! at horizontal RHO-points and vertical W-points.
1728! Also set FC at the top and bottom levels.
1729!
1730 cff=-dt(ng)*lambda
1731 DO k=1,n(ng)-1
1732 DO i=istr,iend
1733 cff1=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1734 fc(i,k)=cff*cff1*akt(i,j,k,ltrc)
1735 END DO
1736 END DO
1737 DO i=istr,iend
1738 fc(i,0)=0.0_r8
1739 fc(i,n(ng))=0.0_r8
1740 END DO
1741!
1742! Compute diagonal matrix coefficients BC and load right-hand-side
1743! terms for the tracer equation into DC.
1744!
1745 DO k=1,n(ng)
1746 DO i=istr,iend
1747 bc(i,k)=hz(i,j,k)-fc(i,k)-fc(i,k-1)
1748 dc(i,k)=t(i,j,k,nnew,itrc)
1749 END DO
1750 END DO
1751!
1752! Solve the tridiagonal system.
1753!
1754 DO i=istr,iend
1755 cff=1.0_r8/bc(i,1)
1756 cf(i,1)=cff*fc(i,1)
1757 dc(i,1)=cff*dc(i,1)
1758 END DO
1759 DO k=2,n(ng)-1
1760 DO i=istr,iend
1761 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
1762 cf(i,k)=cff*fc(i,k)
1763 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
1764 END DO
1765 END DO
1766!
1767! Compute new solution by back substitution.
1768!
1769 DO i=istr,iend
1770# ifdef DIAGNOSTICS_TS
1771 cff1=t(i,j,n(ng),nnew,itrc)*ohz(i,j,n(ng))
1772# endif
1773 dc(i,n(ng))=(dc(i,n(ng))-fc(i,n(ng)-1)*dc(i,n(ng)-1))/ &
1774 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
1775 t(i,j,n(ng),nnew,itrc)=dc(i,n(ng))
1776# ifdef DIAGNOSTICS_TS
1777 diatwrk(i,j,n(ng),itrc,itvdif)= &
1778 & diatwrk(i,j,n(ng),itrc,itvdif)+ &
1779 & t(i,j,n(ng),nnew,itrc)-cff1
1780# endif
1781 END DO
1782 DO k=n(ng)-1,1,-1
1783 DO i=istr,iend
1784# ifdef DIAGNOSTICS_TS
1785 cff1=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1786# endif
1787 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1788 t(i,j,k,nnew,itrc)=dc(i,k)
1789# ifdef DIAGNOSTICS_TS
1790 diatwrk(i,j,k,itrc,itvdif)=diatwrk(i,j,k,itrc,itvdif)+ &
1791 & t(i,j,k,nnew,itrc)-cff1
1792# endif
1793 END DO
1794 END DO
1795# ifdef SPLINES_VDIFF
1796 END IF
1797# endif
1798 END DO
1799 END DO j_loop2
1800
1801# if defined AGE_MEAN && defined T_PASSIVE
1802!
1803!-----------------------------------------------------------------------
1804! If inert passive tracer and Mean Age, compute age concentration (even
1805! inert index) forced by the right-hand-side term that is concentration
1806! of an associated conservative passive tracer (odd inert index). Mean
1807! Age is age concentration divided by conservative passive tracer
1808! concentration. Code implements NPT/2 mean age tracer pairs.
1809!
1810! Implemented and tested by W.G. Zhang and J. Wilkin. See following
1811! reference for details.
1812!
1813! Zhang et al. (2010): Simulation of water age and residence time in
1814! the New York Bight, JPO, 40,965-982, doi:10.1175/2009JPO4249.1
1815!-----------------------------------------------------------------------
1816!
1817 DO itrc=1,npt,2
1818 iage=inert(itrc+1) ! even inert tracer index
1819 DO k=1,n(ng)
1820 DO j=jstr,jend
1821 DO i=istr,iend
1822 IF ((hadvection(inert(itrc),ng)%MPDATA).and. &
1823 & (vadvection(inert(itrc),ng)%MPDATA)) THEN
1824 t(i,j,k,nnew,iage)=t(i,j,k,nnew,iage)+ &
1825 & dt(ng)*t(i,j,k,nnew,inert(itrc))
1826 ELSE
1827 t(i,j,k,nnew,iage)=t(i,j,k,nnew,iage)+ &
1828 & dt(ng)*t(i,j,k,3,inert(itrc))
1829 END IF
1830 END DO
1831 END DO
1832 END DO
1833 END DO
1834# endif
1835!
1836!-----------------------------------------------------------------------
1837! Apply lateral boundary conditions and, if appropriate, nudge
1838! to tracer data and apply Land/Sea mask.
1839!-----------------------------------------------------------------------
1840!
1841! Initialize tracer counter index. The "tclm" array is only allocated
1842! to the NTCLM fields that need to be processed. This is done to
1843! reduce memory.
1844!
1845 ic=0
1846!
1847 DO itrc=1,nt(ng)
1848!
1849! Set compact reduced memory tracer index for nudging coefficients and
1850! climatology arrays.
1851!
1852 IF (ltracerclm(itrc,ng).and.lnudgetclm(itrc,ng)) THEN
1853 ic=ic+1
1854 END IF
1855!
1856! Set lateral boundary conditions.
1857!
1858 CALL t3dbc_tile (ng, tile, itrc, ic, &
1859 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
1860 & imins, imaxs, jmins, jmaxs, &
1861 & nstp, nnew, &
1862 & t)
1863!
1864! Nudge towards tracer climatology.
1865!
1866 IF (ltracerclm(itrc,ng).and.lnudgetclm(itrc,ng)) THEN
1867 DO k=1,n(ng)
1868 DO j=jstrr,jendr
1869 DO i=istrr,iendr
1870 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+ &
1871 & dt(ng)* &
1872 & clima(ng)%Tnudgcof(i,j,k,ic)* &
1873 & (clima(ng)%tclm(i,j,k,ic)- &
1874 & t(i,j,k,nnew,itrc))
1875 END DO
1876 END DO
1877 END DO
1878 END IF
1879
1880# ifdef MASKING
1881!
1882! Apply Land/Sea mask.
1883!
1884 DO k=1,n(ng)
1885 DO j=jstrr,jendr
1886 DO i=istrr,iendr
1887 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*rmask(i,j)
1888 END DO
1889 END DO
1890 END DO
1891# endif
1892# ifdef DIAGNOSTICS_TS
1893!
1894! Compute time-rate-of-change diagnostic term.
1895!
1896 DO k=1,n(ng)
1897 DO j=jstrr,jendr
1898 DO i=istrr,iendr
1899 diatwrk(i,j,k,itrc,itrate)=t(i,j,k,nnew,itrc)- &
1900 & diatwrk(i,j,k,itrc,itrate)
1901!! DiaTwrk(i,j,k,itrc,iTrate)=t(i,j,k,nnew,itrc)- &
1902!! & t(i,j,k,nstp,itrc)
1903 END DO
1904 END DO
1905 END DO
1906# endif
1907!
1908! Apply periodic boundary conditions.
1909!
1910 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1911 CALL exchange_r3d_tile (ng, tile, &
1912 & lbi, ubi, lbj, ubj, 1, n(ng), &
1913 & t(:,:,:,nnew,itrc))
1914 END IF
1915 END DO
1916# ifdef DISTRIBUTE
1917!
1918! Exchange boundary data.
1919!
1920 CALL mp_exchange4d (ng, tile, inlm, 1, &
1921 & lbi, ubi, lbj, ubj, 1, n(ng), 1, nt(ng), &
1922 & nghostpoints, &
1923 & ewperiodic(ng), nsperiodic(ng), &
1924 & t(:,:,:,nnew,:))
1925# endif
1926# if defined FLOATS && defined FLOAT_VWALK
1927!
1928!-----------------------------------------------------------------------
1929! Compute vertical gradient in vertical T-diffusion coefficient for
1930! floats random walk.
1931!-----------------------------------------------------------------------
1932!
1933 DO j=jstrr,jendr
1934 DO i=istrr,iendr
1935 DO k=1,n(ng)
1936 daktdz(i,j,k)=(akt(i,j,k,1)-akt(i,j,k-1,1))/hz(i,j,k)
1937 END DO
1938 END DO
1939 END DO
1940!
1941! Apply periodic boundary conditions.
1942!
1943 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1944 CALL exchange_r3d_tile (ng, tile, &
1945 & lbi, ubi, lbj, ubj, 1, n(ng), &
1946 & daktdz)
1947 END IF
1948
1949# ifdef DISTRIBUTE
1950 CALL mp_exchange3d (ng, tile, inlm, 1, &
1951 & lbi, ubi, lbj, ubj, 1, n(ng), &
1952 & nghostpoints, &
1953 & ewperiodic(ng), nsperiodic(ng), &
1954 & daktdz)
1955# endif
1956# endif
1957!
1958!-----------------------------------------------------------------------
1959! If applicable, deallocate local arrays.
1960!-----------------------------------------------------------------------
1961!
1962 IF (lmpdata) THEN
1963# ifdef DIAGNOSTICS_TS
1964 IF (allocated(dhadv)) deallocate (dhadv)
1965 IF (allocated(dvadv)) deallocate (dvadv)
1966# endif
1967 IF (allocated(ta)) deallocate (ta)
1968 IF (allocated(ua)) deallocate (ua)
1969 IF (allocated(va)) deallocate (va)
1970 IF (allocated(wa)) deallocate (wa)
1971 END IF
1972!
1973 RETURN
1974 END SUBROUTINE step3d_t_tile
1975#endif
1976 END MODULE step3d_t_mod
mpdata_adiff_tile
subroutine mpdata_adiff_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, pm, pn, omn, om_u, on_v, z_r, ohz, huon, hvom, w, t, ta, ua, va, wa)
Definition mpdata_adiff.F:58
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_r3d_tile
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:283
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_nesting
Definition mod_nesting.F:2
mod_nesting::bry_contact
type(t_bcp), dimension(:,:), allocatable bry_contact
Definition mod_nesting.F:335
mod_nesting::rcontact
type(t_ngc), dimension(:), allocatable rcontact
Definition mod_nesting.F:286
mod_nesting::ncontact
integer ncontact
Definition mod_nesting.F:256
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::hadvection
type(t_adv), dimension(:,:), allocatable hadvection
Definition mod_param.F:403
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::bounds
type(t_bounds), dimension(:), allocatable bounds
Definition mod_param.F:232
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::vadvection
type(t_adv), dimension(:,:), allocatable vadvection
Definition mod_param.F:404
mod_param::ndt
integer ndt
Definition mod_param.F:574
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_param::npt
integer npt
Definition mod_param.F:505
mod_scalars
Definition mod_scalars.F:2
mod_scalars::luvsrc
logical, dimension(:), allocatable luvsrc
Definition mod_scalars.F:542
mod_scalars::cc3
real(r8) cc3
Definition mod_scalars.F:395
mod_scalars::ityadv
integer ityadv
Definition mod_scalars.F:142
mod_scalars::ltracersrc
logical, dimension(:,:), allocatable ltracersrc
Definition mod_scalars.F:544
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::lambda
real(r8) lambda
Definition mod_scalars.F:750
mod_scalars::cc1
real(r8) cc1
Definition mod_scalars.F:393
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::ithadv
integer ithadv
Definition mod_scalars.F:140
mod_scalars::itvadv
integer itvadv
Definition mod_scalars.F:139
mod_scalars::lwsrc
logical, dimension(:), allocatable lwsrc
Definition mod_scalars.F:543
mod_scalars::itrate
integer itrate
Definition mod_scalars.F:138
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::cc2
real(r8) cc2
Definition mod_scalars.F:394
mod_scalars::itvdif
integer itvdif
Definition mod_scalars.F:143
mod_scalars::inert
integer, dimension(:), pointer inert
Definition mod_scalars.F:129
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::refinedgrid
logical, dimension(:), allocatable refinedgrid
Definition mod_scalars.F:494
mod_scalars::ltracerclm
logical, dimension(:,:), allocatable ltracerclm
Definition mod_scalars.F:528
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::lnudgetclm
logical, dimension(:,:), allocatable lnudgetclm
Definition mod_scalars.F:535
mod_scalars::itxadv
integer itxadv
Definition mod_scalars.F:141
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
mod_sources
Definition mod_sources.F:2
mod_sources::sources
type(t_sources), dimension(:), allocatable sources
Definition mod_sources.F:90
mod_sources::nsrc
integer, dimension(:), allocatable nsrc
Definition mod_sources.F:97
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange4d
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
Definition mp_exchange.F:3377
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
nesting_mod
Definition nesting.F:2
nesting_mod::bry_fluxes
subroutine, public bry_fluxes(dg, rg, cr, model, tile, imins, imaxs, jmins, jmaxs, ilb, iub, jlb, jub, scale, fx, fe, f_west, f_east, f_south, f_north)
Definition nesting.F:381
step3d_t_mod
Definition step3d_t.F:2
step3d_t_mod::step3d_t
subroutine, public step3d_t(ng, tile)
Definition step3d_t.F:41
step3d_t_mod::step3d_t_tile
subroutine step3d_t_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, omn, om_u, om_v, on_u, on_v, pm, pn, hz, huon, hvom, z_r, akt, w, bed_thick, daktdz, diatwrk, t)
Definition step3d_t.F:152
t3dbc_mod
Definition t3dbc_im.F:2
t3dbc_mod::t3dbc_tile
subroutine, public t3dbc_tile(ng, tile, itrc, ic, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nout, t)
Definition t3dbc_im.F:55
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3