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