ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
tl_step3d_t_mod Module Reference

Functions/Subroutines

subroutine, public tl_step3d_t (ng, tile)
 
subroutine tl_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)
 

Function/Subroutine Documentation

◆ tl_step3d_t()

subroutine, public tl_step3d_t_mod::tl_step3d_t ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 46 of file tl_step3d_t.F.

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, itlm, 35, __line__, myfile)
71# endif
72 CALL tl_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, itlm, 35, __line__, myfile)
115# endif
116!
117 RETURN
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_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::itlm
integer, parameter itlm
Definition mod_param.F:663
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
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

References mod_grid::grid, mod_param::itlm, mod_mixing::mixing, mod_stepping::nnew, mod_stepping::nrhs, mod_stepping::nstp, mod_ocean::ocean, tl_step3d_t_tile(), wclock_off(), and wclock_on().

Referenced by tl_main3d().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ tl_step3d_t_tile()

subroutine tl_step3d_t_mod::tl_step3d_t_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
integer, intent(in) nstp,
integer, intent(in) nnew,
real(r8), dimension(lbi:,lbj:), intent(in) rmask,
real(r8), dimension(lbi:,lbj:), intent(in) umask,
real(r8), dimension(lbi:,lbj:), intent(in) vmask,
real(r8), dimension(lbi:,lbj:), intent(in) rmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) umask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) vmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) omn,
real(r8), dimension(lbi:,lbj:), intent(in) om_u,
real(r8), dimension(lbi:,lbj:), intent(in) om_v,
real(r8), dimension(lbi:,lbj:), intent(in) on_u,
real(r8), dimension(lbi:,lbj:), intent(in) on_v,
real(r8), dimension(lbi:,lbj:), intent(in) pm,
real(r8), dimension(lbi:,lbj:), intent(in) pn,
real(r8), dimension(lbi:,lbj:,:), intent(in) hz,
real(r8), dimension(lbi:,lbj:,:), intent(in) tl_hz,
real(r8), dimension(lbi:,lbj:,:), intent(in) huon,
real(r8), dimension(lbi:,lbj:,:), intent(in) tl_huon,
real(r8), dimension(lbi:,lbj:,:), intent(in) hvom,
real(r8), dimension(lbi:,lbj:,:), intent(in) tl_hvom,
real(r8), dimension(lbi:,lbj:,:), intent(in) z_r,
real(r8), dimension(lbi:,lbj:,:), intent(in) tl_z_r,
real(r8), dimension(lbi:,lbj:,0:,:), intent(in) akt,
real(r8), dimension(lbi:,lbj:,0:,:), intent(in) tl_akt,
real(r8), dimension(lbi:,lbj:,0:), intent(in) w,
real(r8), dimension(lbi:,lbj:,0:), intent(in) tl_w,
real(r8), dimension(lbi:,lbj:,:,:,:), intent(in) t,
real(r8), dimension(lbi:,lbj:,:,:,:), intent(inout) tl_t )
private

Definition at line 121 of file tl_step3d_t.F.

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

References mod_clima::clima, mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, exchange_3d_mod::exchange_r3d_tile(), mod_scalars::ieast, mod_scalars::inert, mod_param::inlm, mod_scalars::inorth, mod_scalars::isouth, mod_param::itlm, mod_scalars::iwest, mod_scalars::lambda, mod_scalars::lnudgetclm, mod_scalars::ltracerclm, mod_scalars::ltracersrc, mod_scalars::luvsrc, mod_scalars::lwsrc, mp_exchange_mod::mp_exchange3d(), mp_exchange_mod::mp_exchange4d(), mod_param::nghostpoints, mod_param::npt, mod_scalars::nsperiodic, mod_sources::nsrc, mod_sources::sources, mod_param::tl_hadvection, tl_t3dbc_mod::tl_t3dbc_tile(), and mod_param::tl_vadvection.

Referenced by tl_step3d_t().

Here is the call graph for this function:
Here is the caller graph for this function: