ROMS
Loading...
Searching...
No Matches
tl_pre_step3d.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE tl_pre_step3d_mod
3#if defined TANGENT && defined SOLVE3D
4# define NEUMANN
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! This subroutine initialize computations for new time step of the !
14! tangent linear 3D primitive variables. !
15! !
16! Both n-1 and n-2 time-step contributions of the Adams/Bashforth !
17! scheme are added here to u and v at time index "nnew", since the !
18! right-hand-side arrays ru and rv at n-2 will be overwriten in !
19! subsequent calls to routines within the 3D engine. !
20! !
21! It also computes the time "n" vertical viscosity and diffusion !
22! contributions of the Crank-Nicholson implicit scheme because the !
23! thicknesses "Hz" will be overwriten at the end of the 2D engine !
24! (barotropic mode) computations. !
25! !
26! The actual time step will be carried out in routines "step3d_uv" !
27! and "step3d_t". !
28! !
29! BASIC STATE variables needed: AKt, AKv, Huon, Hvom, Hz, Tsrc, W, !
30! bustr, bvstr, ghats, srflx, sustr, !
31! svstr, t, z_r, z_w !
32! !
33!=======================================================================
34!
35 implicit none
36!
37 PRIVATE
38 PUBLIC :: tl_pre_step3d
39!
40 CONTAINS
41!
42!***********************************************************************
43 SUBROUTINE tl_pre_step3d (ng, tile)
44!***********************************************************************
45!
46 USE mod_param
47# ifdef DIAGNOSTICS
48!! USE mod_diags
49# endif
50 USE mod_forces
51 USE mod_grid
52 USE mod_mixing
53 USE mod_ocean
54 USE mod_stepping
55!
56! Imported variable declarations.
57!
58 integer, intent(in) :: ng, tile
59!
60! Local variable declarations.
61!
62 character (len=*), parameter :: myfile = &
63 & __FILE__
64!
65# include "tile.h"
66!
67# ifdef PROFILE
68 CALL wclock_on (ng, itlm, 22, __line__, myfile)
69# endif
70 CALL tl_pre_step3d_tile (ng, tile, &
71 & lbi, ubi, lbj, ubj, &
72 & imins, imaxs, jmins, jmaxs, &
73 & nrhs(ng), nstp(ng), nnew(ng), &
74# ifdef MASKING
75 & grid(ng) % rmask, &
76 & grid(ng) % umask, &
77 & grid(ng) % vmask, &
78# if defined SOLAR_SOURCE && defined WET_DRY
79 & grid(ng) % rmask_wet, &
80# endif
81# endif
82 & grid(ng) % pm, &
83 & grid(ng) % pn, &
84 & grid(ng) % Hz, &
85 & grid(ng) % tl_Hz, &
86 & grid(ng) % Huon, &
87 & grid(ng) % tl_Huon, &
88 & grid(ng) % Hvom, &
89 & grid(ng) % tl_Hvom, &
90 & grid(ng) % z_r, &
91 & grid(ng) % tl_z_r, &
92 & grid(ng) % z_w, &
93 & grid(ng) % tl_z_w, &
94 & forces(ng) % tl_btflx, &
95 & forces(ng) % tl_bustr, &
96 & forces(ng) % tl_bvstr, &
97 & forces(ng) % tl_stflx, &
98 & forces(ng) % tl_sustr, &
99 & forces(ng) % tl_svstr, &
100# ifdef SOLAR_SOURCE
101 & forces(ng) % srflx, &
102# endif
103 & mixing(ng) % Akt, &
104 & mixing(ng) % tl_Akt, &
105 & mixing(ng) % Akv, &
106 & mixing(ng) % tl_Akv, &
107# ifdef LMD_NONLOCAL_NOT_YET
108 & mixing(ng) % tl_ghats, &
109# endif
110 & ocean(ng) % W, &
111 & ocean(ng) % tl_W, &
112 & ocean(ng) % tl_ru, &
113 & ocean(ng) % tl_rv, &
114# ifdef DIAGNOSTICS_TS
115!! & DIAGS(ng) % DiaTwrk, &
116# endif
117# ifdef DIAGNOSTICS_UV
118!! & DIAGS(ng) % DiaU3wrk, &
119!! & DIAGS(ng) % DiaV3wrk, &
120!! & DIAGS(ng) % DiaRU, &
121!! & DIAGS(ng) % DiaRV, &
122# endif
123 & ocean(ng) % t, &
124 & ocean(ng) % tl_t, &
125 & ocean(ng) % u, &
126 & ocean(ng) % tl_u, &
127 & ocean(ng) % v, &
128 & ocean(ng) % tl_v)
129# ifdef PROFILE
130 CALL wclock_off (ng, itlm, 22, __line__, myfile)
131# endif
132!
133 RETURN
134 END SUBROUTINE tl_pre_step3d
135!
136!***********************************************************************
137 SUBROUTINE tl_pre_step3d_tile (ng, tile, &
138 & LBi, UBi, LBj, UBj, &
139 & IminS, ImaxS, JminS, JmaxS, &
140 & nrhs, nstp, nnew, &
141# ifdef MASKING
142 & rmask, umask, vmask, &
143# if defined SOLAR_SOURCE && defined WET_DRY
144 & rmask_wet, &
145# endif
146# endif
147 & pm, pn, &
148 & Hz, tl_Hz, &
149 & Huon, tl_Huon, &
150 & Hvom, tl_Hvom, &
151 & z_r, tl_z_r, &
152 & z_w, tl_z_w, &
153 & tl_btflx, tl_bustr, tl_bvstr, &
154 & tl_stflx, tl_sustr, tl_svstr, &
155# ifdef SOLAR_SOURCE
156 & srflx, &
157# endif
158 & Akt, tl_Akt, &
159 & Akv, tl_Akv, &
160# ifdef LMD_NONLOCAL_NOT_YET
161 & tl_ghats, &
162# endif
163 & W, tl_W, &
164 & tl_ru, tl_rv, &
165# ifdef DIAGNOSTICS_TS
166!! & DiaTwrk, &
167# endif
168# ifdef DIAGNOSTICS_UV
169!! & DiaU3wrk, DiaV3wrk, &
170!! & DiaRU, DiaRV, &
171# endif
172 & t, tl_t, &
173 & u, tl_u, &
174 & v, tl_v)
175!***********************************************************************
176!
177 USE mod_param
178 USE mod_scalars
179 USE mod_sources
180!
181 USE exchange_3d_mod, ONLY : exchange_r3d_tile
182# ifdef DISTRIBUTE
183 USE mp_exchange_mod, ONLY : mp_exchange4d
184# endif
185 USE tl_t3dbc_mod, ONLY : tl_t3dbc_tile
186!
187! Imported variable declarations.
188!
189 integer, intent(in) :: ng, tile
190 integer, intent(in) :: LBi, UBi, LBj, UBj
191 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
192 integer, intent(in) :: nrhs, nstp, nnew
193!
194# ifdef ASSUMED_SHAPE
195# ifdef MASKING
196 real(r8), intent(in) :: rmask(LBi:,LBj:)
197 real(r8), intent(in) :: umask(LBi:,LBj:)
198 real(r8), intent(in) :: vmask(LBi:,LBj:)
199# if defined SOLAR_SOURCE && defined WET_DRY
200 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
201# endif
202# endif
203 real(r8), intent(in) :: pm(LBi:,LBj:)
204 real(r8), intent(in) :: pn(LBi:,LBj:)
205 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
206 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
207 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
208 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
209 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
210# ifdef SOLAR_SOURCE
211 real(r8), intent(in) :: srflx(LBi:,LBj:)
212# endif
213# ifdef SUN
214 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
215# else
216 real(r8), intent(in) :: Akt(LBi:,LBj:,0:,:)
217# endif
218 real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
219 real(r8), intent(in) :: W(LBi:,LBj:,0:)
220# ifdef SUN
221 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
222# else
223 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
224# endif
225 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
226 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
227
228 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
229 real(r8), intent(in) :: tl_Huon(LBi:,LBj:,:)
230 real(r8), intent(in) :: tl_Hvom(LBi:,LBj:,:)
231 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
232 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
233 real(r8), intent(in) :: tl_btflx(LBi:,LBj:,:)
234 real(r8), intent(in) :: tl_bustr(LBi:,LBj:)
235 real(r8), intent(in) :: tl_bvstr(LBi:,LBj:)
236 real(r8), intent(in) :: tl_stflx(LBi:,LBj:,:)
237 real(r8), intent(in) :: tl_sustr(LBi:,LBj:)
238 real(r8), intent(in) :: tl_svstr(LBi:,LBj:)
239 real(r8), intent(in) :: tl_ru(LBi:,LBj:,0:,:)
240 real(r8), intent(in) :: tl_rv(LBi:,LBj:,0:,:)
241# ifdef LMD_NONLOCAL_NOT_YET
242 real(r8), intent(in) :: tl_ghats(LBi:,LBj:,0:,:)
243# endif
244# ifdef SUN
245 real(r8), intent(in) :: tl_Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
246# else
247 real(r8), intent(in) :: tl_Akt(LBi:,LBj:,0:,:)
248# endif
249 real(r8), intent(in) :: tl_Akv(LBi:,LBj:,0:)
250 real(r8), intent(in) :: tl_W(LBi:,LBj:,0:)
251
252# ifdef DIAGNOSTICS_TS
253!! real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
254# endif
255# ifdef DIAGNOSTICS_UV
256!! real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
257!! real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
258!! real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
259!! real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
260# endif
261# ifdef SUN
262 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
263# else
264 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
265# endif
266 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
267 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
268
269# else
270
271# ifdef MASKING
272 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
273 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
274 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
275# if defined SOLAR_SOURCE && defined WET_DRY
276 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
277# endif
278# endif
279 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
280 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
281 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
282 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
283 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
284 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
285 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
286# ifdef SOLAR_SOURCE
287 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
288# endif
289 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
290 real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
291 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
292 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
293 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
294 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
295
296 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
297 real(r8), intent(in) :: tl_Huon(LBi:UBi,LBj:UBj,N(ng))
298 real(r8), intent(in) :: tl_Hvom(LBi:UBi,LBj:UBj,N(ng))
299 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
300 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
301 real(r8), intent(in) :: tl_btflx(LBi:UBi,LBj:UBj,NT(ng))
302 real(r8), intent(in) :: tl_bustr(LBi:UBi,LBj:UBj)
303 real(r8), intent(in) :: tl_bvstr(LBi:UBi,LBj:UBj)
304 real(r8), intent(in) :: tl_stflx(LBi:UBi,LBj:UBj,NT(ng))
305 real(r8), intent(in) :: tl_sustr(LBi:UBi,LBj:UBj)
306 real(r8), intent(in) :: tl_svstr(LBi:UBi,LBj:UBj)
307 real(r8), intent(in) :: tl_ru(LBi:UBi,LBj:UBj,0:N(ng),2)
308 real(r8), intent(in) :: tl_rv(LBi:UBi,LBj:UBj,0:N(ng),2)
309# ifdef LMD_NONLOCAL_NOT_YET
310 real(r8), intent(in) :: ghats(LBi:UBi,LBj:UBj,0:N(ng),NAT)
311# endif
312 real(r8), intent(in) :: tl_Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
313 real(r8), intent(in) :: tl_Akv(LBi:UBi,LBj:UBj,0:N(ng))
314 real(r8), intent(in) :: tl_W(LBi:UBi,LBj:UBj,0:N(ng))
315
316# ifdef DIAGNOSTICS_TS
317!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
318!! & NDT)
319# endif
320# ifdef DIAGNOSTICS_UV
321!! real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
322!! real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
323!! real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
324!! real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
325# endif
326 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
327 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
328 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
329# endif
330!
331! Local variable declarations.
332!
333 integer :: Isrc, Jsrc
334 integer :: i, ic, indx, is, itrc, j, k, ltrc
335# if defined AGE_MEAN && defined T_PASSIVE
336 integer :: iage
337# endif
338# if defined DIAGNOSTICS_TS || defined DIAGNOSTICS_UV
339 integer :: idiag
340# endif
341 real(r8), parameter :: eps = 1.0e-16_r8
342
343 real(r8) :: cff, cff1, cff2, cff3, cff4
344 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3, tl_cff4
345 real(r8) :: Gamma
346
347 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
348 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
349 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
350
351 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_CF
352 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
353 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
354
355# ifdef SOLAR_SOURCE
356 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: tl_swdk
357# endif
358
359 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
360 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
361 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curv
362 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
363
364 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
365 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
366 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_curv
367 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_grad
368
369# include "set_bounds.h"
370
371# ifndef TS_FIXED
372!
373!=======================================================================
374! Tracer equation(s).
375!=======================================================================
376
377# ifdef SOLAR_SOURCE
378!
379! Compute fraction of the solar shortwave radiation, "swdk"
380! (at vertical W-points) penetrating water column.
381!
382 DO k=1,n(ng)-1
383 DO j=jstr,jend
384 DO i=istr,iend
385 fx(i,j)=z_w(i,j,n(ng))-z_w(i,j,k)
386 tl_fx(i,j)=tl_z_w(i,j,n(ng))-tl_z_w(i,j,k)
387 END DO
388 END DO
389!^ CALL lmd_swfrac_tile (ng, tile, &
390!^ & LBi, UBi, LBj, UBj, &
391!^ & IminS, ImaxS, JminS, JmaxS, &
392!^ & -1.0_r8, FX, FE)
393!^
394 CALL tl_lmd_swfrac_tile (ng, tile, &
395 & lbi, ubi, lbj, ubj, &
396 & imins, imaxs, jmins, jmaxs, &
397 & -1.0_r8, fx, tl_fx, tl_fe)
398 DO j=jstr,jend
399 DO i=istr,iend
400!^ swdk(i,j,k)=FE(i,j)
401!^
402 tl_swdk(i,j,k)=tl_fe(i,j)
403 END DO
404 END DO
405 END DO
406# endif
407!
408!-----------------------------------------------------------------------
409! Compute intermediate tracer at n+1/2 time-step, t(i,j,k,3,itrc).
410!-----------------------------------------------------------------------
411!
412! Compute time rate of change of intermediate tracer due to
413! horizontal advection.
414!
415 t_loop1 : DO itrc=1,nt(ng)
416 k_loop : DO k=1,n(ng)
417!
418 hadv_flux : IF (tl_hadvection(itrc,ng)%CENTERED2) THEN
419!
420! Second-order, centered differences horizontal advective fluxes.
421!
422 DO j=jstr,jend
423 DO i=istr,iend+1
424!^ FX(i,j)=Huon(i,j,k)* &
425!^ & 0.5_r8*(t(i-1,j,k,nstp,itrc)+ &
426!^ & t(i ,j,k,nstp,itrc))
427!^
428 tl_fx(i,j)=0.5_r8* &
429 & (tl_huon(i,j,k)* &
430 & (t(i-1,j,k,nstp,itrc)+ &
431 & t(i ,j,k,nstp,itrc))+ &
432 & huon(i,j,k)* &
433 & (tl_t(i-1,j,k,nstp,itrc)+ &
434 & tl_t(i ,j,k,nstp,itrc)))
435 END DO
436 END DO
437 DO j=jstr,jend+1
438 DO i=istr,iend
439!^ FE(i,j)=Hvom(i,j,k)* &
440!^ & 0.5_r8*(t(i,j-1,k,nstp,itrc)+ &
441!^ & t(i,j ,k,nstp,itrc))
442!^
443 tl_fe(i,j)=0.5_r8* &
444 & (tl_hvom(i,j,k)* &
445 & (t(i,j-1,k,nstp,itrc)+ &
446 & t(i,j ,k,nstp,itrc))+ &
447 & hvom(i,j,k)* &
448 & (tl_t(i,j-1,k,nstp,itrc)+ &
449 & tl_t(i,j ,k,nstp,itrc)))
450 END DO
451 END DO
452!
453 ELSE IF ((tl_hadvection(itrc,ng)%MPDATA).or. &
454 & (tl_hadvection(itrc,ng)%HSIMT)) THEN
455!
456! First-order, upstream differences horizontal advective fluxes.
457!
458 DO j=jstr,jend
459 DO i=istr,iend+1
460 cff1=max(huon(i,j,k),0.0_r8)
461 cff2=min(huon(i,j,k),0.0_r8)
462 tl_cff1=(0.5_r8+sign(0.5_r8, huon(i,j,k)))* &
463 & tl_huon(i,j,k)
464 tl_cff2=(0.5_r8+sign(0.5_r8,-huon(i,j,k)))* &
465 & tl_huon(i,j,k)
466!^ FX(i,j)=cff1*t(i-1,j,k,nstp,itrc)+ &
467!^ & cff2*t(i ,j,k,nstp,itrc)
468!^
469 tl_fx(i,j)=tl_cff1*t(i-1,j,k,nstp,itrc)+ &
470 & cff1*tl_t(i-1,j,k,nstp,itrc)+ &
471 & tl_cff2*t(i ,j,k,nstp,itrc)+ &
472 & cff2*tl_t(i ,j,k,nstp,itrc)
473 END DO
474 END DO
475 DO j=jstr,jend+1
476 DO i=istr,iend
477 cff1=max(hvom(i,j,k),0.0_r8)
478 cff2=min(hvom(i,j,k),0.0_r8)
479 tl_cff1=(0.5_r8+sign(0.5_r8, hvom(i,j,k)))* &
480 & tl_hvom(i,j,k)
481 tl_cff2=(0.5_r8+sign(0.5_r8,-hvom(i,j,k)))* &
482 & tl_hvom(i,j,k)
483!^ FE(i,j)=cff1*t(i,j-1,k,nstp,itrc)+ &
484!^ & cff2*t(i,j ,k,nstp,itrc)
485!^
486 tl_fe(i,j)=tl_cff1*t(i,j-1,k,nstp,itrc)+ &
487 & cff1*tl_t(i,j-1,k,nstp,itrc)+ &
488 & tl_cff2*t(i,j ,k,nstp,itrc)+ &
489 & cff2*tl_t(i,j ,k,nstp,itrc)
490 END DO
491 END DO
492!
493 ELSE IF ((tl_hadvection(itrc,ng)%AKIMA4).or. &
494 & (tl_hadvection(itrc,ng)%CENTERED4).or. &
495 & (tl_hadvection(itrc,ng)%SPLIT_U3).or. &
496 & (tl_hadvection(itrc,ng)%UPSTREAM3)) THEN
497!
498! Fourth-order Akima, fourth-order centered differences, or third-order
499! upstream-biased horizontal advective fluxes.
500!
501 DO j=jstr,jend
502 DO i=istrm1,iendp2
503 fx(i,j)=t(i ,j,k,nstp,itrc)- &
504 & t(i-1,j,k,nstp,itrc)
505 tl_fx(i,j)=tl_t(i ,j,k,nstp,itrc)- &
506 & tl_t(i-1,j,k,nstp,itrc)
507# ifdef MASKING
508 fx(i,j)=fx(i,j)*umask(i,j)
509 tl_fx(i,j)=tl_fx(i,j)*umask(i,j)
510# endif
511 END DO
512 END DO
513 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
514 IF (domain(ng)%Western_Edge(tile)) THEN
515 DO j=jstr,jend
516 fx(istr-1,j)=fx(istr,j)
517 tl_fx(istr-1,j)=tl_fx(istr,j)
518 END DO
519 END IF
520 END IF
521 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
522 IF (domain(ng)%Eastern_Edge(tile)) THEN
523 DO j=jstr,jend
524 fx(iend+2,j)=fx(iend+1,j)
525 tl_fx(iend+2,j)=tl_fx(iend+1,j)
526 END DO
527 END IF
528 END IF
529!
530 DO j=jstr,jend
531 DO i=istr-1,iend+1
532 IF (tl_hadvection(itrc,ng)%UPSTREAM3) THEN
533 curv(i,j)=fx(i+1,j)-fx(i,j)
534 tl_curv(i,j)=tl_fx(i+1,j)-tl_fx(i,j)
535 ELSE IF (tl_hadvection(itrc,ng)%AKIMA4) THEN
536 cff=2.0_r8*fx(i+1,j)*fx(i,j)
537 tl_cff=2.0_r8*(tl_fx(i+1,j)*fx(i,j)+ &
538 & fx(i+1,j)*tl_fx(i,j))
539 IF (cff.gt.eps) THEN
540 grad(i,j)=cff/(fx(i+1,j)+fx(i,j))
541 tl_grad(i,j)=((fx(i+1,j)+fx(i,j))*tl_cff- &
542 & cff*(tl_fx(i+1,j)+tl_fx(i,j)))/ &
543 & ((fx(i+1,j)+fx(i,j))* &
544 & (fx(i+1,j)+fx(i,j)))
545 ELSE
546 grad(i,j)=0.0_r8
547 tl_grad(i,j)=0.0_r8
548 END IF
549 ELSE IF ((tl_hadvection(itrc,ng)%CENTERED4).or. &
550 & (tl_hadvection(itrc,ng)%SPLIT_U3)) THEN
551 grad(i,j)=0.5_r8*(fx(i+1,j)+fx(i,j))
552 tl_grad(i,j)=0.5_r8*(tl_fx(i+1,j)+tl_fx(i,j))
553 END IF
554 END DO
555 END DO
556!
557 cff1=1.0_r8/6.0_r8
558 cff2=1.0_r8/3.0_r8
559 DO j=jstr,jend
560 DO i=istr,iend+1
561 IF (tl_hadvection(itrc,ng)%UPSTREAM3) THEN
562!^ FX(i,j)=Huon(i,j,k)*0.5_r8* &
563!^ & (t(i-1,j,k,nstp,itrc)+ &
564!^ & t(i ,j,k,nstp,itrc))- &
565!^ & cff1*(curv(i-1,j)*MAX(Huon(i,j,k),0.0_r8)+ &
566!^ & curv(i ,j)*MIN(Huon(i,j,k),0.0_r8))
567!^
568 tl_fx(i,j)=0.5_r8* &
569 & (tl_huon(i,j,k)* &
570 & (t(i-1,j,k,nstp,itrc)+ &
571 & t(i ,j,k,nstp,itrc))+ &
572 & huon(i,j,k)* &
573 & (tl_t(i-1,j,k,nstp,itrc)+ &
574 & tl_t(i ,j,k,nstp,itrc)))- &
575 & cff1* &
576 & (tl_curv(i-1,j)*max(huon(i,j,k),0.0_r8)+ &
577 & curv(i-1,j)* &
578 & (0.5_r8+sign(0.5_r8, huon(i,j,k)))* &
579 & tl_huon(i,j,k)+ &
580 & tl_curv(i ,j)*min(huon(i,j,k),0.0_r8)+ &
581 & curv(i ,j)* &
582 & (0.5_r8+sign(0.5_r8,-huon(i,j,k)))* &
583 & tl_huon(i,j,k))
584 ELSE IF ((tl_hadvection(itrc,ng)%AKIMA4).or. &
585 & (tl_hadvection(itrc,ng)%CENTERED4).or. &
586 & (tl_hadvection(itrc,ng)%SPLIT_U3)) THEN
587!^ FX(i,j)=Huon(i,j,k)*0.5_r8* &
588!^ & (t(i-1,j,k,nstp,itrc)+ &
589!^ & t(i ,j,k,nstp,itrc)- &
590!^ & cff2*(grad(i ,j)- &
591!^ & grad(i-1,j)))
592!^
593 tl_fx(i,j)=0.5_r8* &
594 & (tl_huon(i,j,k)* &
595 & (t(i-1,j,k,nstp,itrc)+ &
596 & t(i ,j,k,nstp,itrc)- &
597 & cff2*(grad(i ,j)- &
598 & grad(i-1,j)))+ &
599 & huon(i,j,k)* &
600 & (tl_t(i-1,j,k,nstp,itrc)+ &
601 & tl_t(i ,j,k,nstp,itrc)- &
602 & cff2*(tl_grad(i ,j)- &
603 & tl_grad(i-1,j))))
604 END IF
605 END DO
606 END DO
607!
608 DO j=jstrm1,jendp2
609 DO i=istr,iend
610 fe(i,j)=t(i,j ,k,nstp,itrc)- &
611 & t(i,j-1,k,nstp,itrc)
612 tl_fe(i,j)=tl_t(i,j ,k,nstp,itrc)- &
613 & tl_t(i,j-1,k,nstp,itrc)
614# ifdef MASKING
615 fe(i,j)=fe(i,j)*vmask(i,j)
616 tl_fe(i,j)=tl_fe(i,j)*vmask(i,j)
617# endif
618 END DO
619 END DO
620 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
621 IF (domain(ng)%Southern_Edge(tile)) THEN
622 DO i=istr,iend
623 fe(i,jstr-1)=fe(i,jstr)
624 tl_fe(i,jstr-1)=tl_fe(i,jstr)
625 END DO
626 END IF
627 END IF
628 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
629 IF (domain(ng)%Northern_Edge(tile)) THEN
630 DO i=istr,iend
631 fe(i,jend+2)=fe(i,jend+1)
632 tl_fe(i,jend+2)=tl_fe(i,jend+1)
633 END DO
634 END IF
635 END IF
636!
637 DO j=jstr-1,jend+1
638 DO i=istr,iend
639 IF (tl_hadvection(itrc,ng)%UPSTREAM3) THEN
640 curv(i,j)=fe(i,j+1)-fe(i,j)
641 tl_curv(i,j)=tl_fe(i,j+1)-tl_fe(i,j)
642 ELSE IF (tl_hadvection(itrc,ng)%AKIMA4) THEN
643 cff=2.0_r8*fe(i,j+1)*fe(i,j)
644 tl_cff=2.0_r8*(tl_fe(i,j+1)*fe(i,j)+ &
645 & fe(i,j+1)*tl_fe(i,j))
646 IF (cff.gt.eps) THEN
647 grad(i,j)=cff/(fe(i,j+1)+fe(i,j))
648 tl_grad(i,j)=((fe(i,j+1)+fe(i,j))*tl_cff- &
649 & cff*(tl_fe(i,j+1)+tl_fe(i,j)))/ &
650 & ((fe(i,j+1)+fe(i,j))* &
651 & (fe(i,j+1)+fe(i,j)))
652 ELSE
653 grad(i,j)=0.0_r8
654 tl_grad(i,j)=0.0_r8
655 END IF
656 ELSE IF ((tl_hadvection(itrc,ng)%CENTERED4).or. &
657 & (tl_hadvection(itrc,ng)%SPLIT_U3)) THEN
658 grad(i,j)=0.5_r8*(fe(i,j+1)+fe(i,j))
659 tl_grad(i,j)=0.5_r8*(tl_fe(i,j+1)+tl_fe(i,j))
660 END IF
661 END DO
662 END DO
663!
664 cff1=1.0_r8/6.0_r8
665 cff2=1.0_r8/3.0_r8
666 DO j=jstr,jend+1
667 DO i=istr,iend
668 IF (tl_hadvection(itrc,ng)%UPSTREAM3) THEN
669!^ FE(i,j)=Hvom(i,j,k)*0.5_r8* &
670!^ & (t(i,j-1,k,nstp,itrc)+ &
671!^ & t(i,j ,k,nstp,itrc))- &
672!^ & cff1*(curv(i,j-1)*MAX(Hvom(i,j,k),0.0_r8)+ &
673!^ & curv(i,j )*MIN(Hvom(i,j,k),0.0_r8))
674!^
675 tl_fe(i,j)=0.5_r8* &
676 & (tl_hvom(i,j,k)* &
677 & (t(i,j-1,k,nstp,itrc)+ &
678 & t(i,j ,k,nstp,itrc))+ &
679 & hvom(i,j,k)* &
680 & (tl_t(i,j-1,k,nstp,itrc)+ &
681 & tl_t(i,j ,k,nstp,itrc)))- &
682 & cff1* &
683 & (tl_curv(i,j-1)*max(hvom(i,j,k),0.0_r8)+ &
684 & curv(i,j-1)* &
685 & (0.5_r8+sign(0.5_r8, hvom(i,j,k)))* &
686 & tl_hvom(i,j,k)+ &
687 & tl_curv(i,j )*min(hvom(i,j,k),0.0_r8)+ &
688 & curv(i,j )* &
689 & (0.5_r8+sign(0.5_r8,-hvom(i,j,k)))* &
690 & tl_hvom(i,j,k))
691 ELSE IF ((tl_hadvection(itrc,ng)%AKIMA4).or. &
692 & (tl_hadvection(itrc,ng)%CENTERED4).or. &
693 & (tl_hadvection(itrc,ng)%SPLIT_U3)) THEN
694!^ FE(i,j)=Hvom(i,j,k)*0.5_r8* &
695!^ & (t(i,j-1,k,nstp,itrc)+ &
696!^ & t(i,j ,k,nstp,itrc)- &
697!^ & cff2*(grad(i,j )- &
698!^ & grad(i,j-1)))
699!^
700 tl_fe(i,j)=0.5_r8* &
701 & (tl_hvom(i,j,k)* &
702 & (t(i,j-1,k,nstp,itrc)+ &
703 & t(i,j ,k,nstp,itrc)- &
704 & cff2*(grad(i,j )- &
705 & grad(i,j-1)))+ &
706 & hvom(i,j,k)* &
707 & (tl_t(i,j-1,k,nstp,itrc)+ &
708 & tl_t(i,j ,k,nstp,itrc)- &
709 & cff2*(tl_grad(i,j )- &
710 & tl_grad(i,j-1))))
711 END IF
712 END DO
713 END DO
714 END IF hadv_flux
715!
716! Apply tracers point sources to the horizontal advection terms,
717! if any.
718!
719! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
720! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
721!
722 IF (luvsrc(ng)) THEN
723 DO is=1,nsrc(ng)
724 isrc=sources(ng)%Isrc(is)
725 jsrc=sources(ng)%Jsrc(is)
726 IF (((istr.le.isrc).and.(isrc.le.iend+1)).and. &
727 & ((jstr.le.jsrc).and.(jsrc.le.jend+1))) THEN
728 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
729 IF (ltracersrc(itrc,ng)) THEN
730!^ FX(Isrc,Jsrc)=Huon(Isrc,Jsrc,k)* &
731!^ & SOURCES(ng)%Tsrc(is,k,itrc)
732!^
733 tl_fx(isrc,jsrc)=tl_huon(isrc,jsrc,k)* &
734 & sources(ng)%Tsrc(is,k,itrc)+ &
735 & huon(isrc,jsrc,k)* &
736 & sources(ng)%tl_Tsrc(is,k,itrc)
737 ELSE
738!^ FX(Isrc,Jsrc)=0.0_r8
739!^
740 tl_fx(isrc,jsrc)=0.0_r8
741 END IF
742 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
743 IF (ltracersrc(itrc,ng)) THEN
744!^ FE(Isrc,Jsrc)=Hvom(Isrc,Jsrc,k)* &
745!^ & SOURCES(ng)%Tsrc(is,k,itrc)
746!^
747 tl_fe(isrc,jsrc)=tl_hvom(isrc,jsrc,k)* &
748 & sources(ng)%Tsrc(is,k,itrc)+ &
749 & hvom(isrc,jsrc,k)* &
750 & sources(ng)%tl_Tsrc(is,k,itrc)
751 ELSE
752!^ FE(Isrc,Jsrc)=0.0_r8
753!^
754 tl_fe(isrc,jsrc)=0.0_r8
755 END IF
756 END IF
757 END IF
758 END DO
759 END IF
760!
761! Time-step horizontal advection (m Tunits).
762!
763 IF ((tl_hadvection(itrc,ng)%MPDATA).or. &
764 & (tl_hadvection(itrc,ng)%HSIMT)) THEN
765 gamma=0.5_r8
766 ELSE
767 gamma=1.0_r8/6.0_r8
768 END IF
769# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
770 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
771# else
772 IF (iic(ng).eq.ntfirst(ng)) THEN
773# endif
774 cff=0.5_r8*dt(ng)
775 cff1=1.0_r8
776 cff2=0.0_r8
777 ELSE
778 cff=(1.0_r8-gamma)*dt(ng)
779 cff1=0.5_r8+gamma
780 cff2=0.5_r8-gamma
781 END IF
782 DO j=jstr,jend
783 DO i=istr,iend
784!^ t(i,j,k,3,itrc)=Hz(i,j,k)*(cff1*t(i,j,k,nstp,itrc)+ &
785!^ & cff2*t(i,j,k,nnew,itrc))- &
786!^ & cff*pm(i,j)*pn(i,j)* &
787!^ & (FX(i+1,j)-FX(i,j)+ &
788!^ & FE(i,j+1)-FE(i,j))
789!^
790 tl_t(i,j,k,3,itrc)=tl_hz(i,j,k)* &
791 & (cff1*t(i,j,k,nstp,itrc)+ &
792 & cff2*t(i,j,k,nnew,itrc))+ &
793 & hz(i,j,k)* &
794 & (cff1*tl_t(i,j,k,nstp,itrc)+ &
795 & cff2*tl_t(i,j,k,nnew,itrc))- &
796 & cff*pm(i,j)*pn(i,j)* &
797 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
798 & tl_fe(i,j+1)-tl_fe(i,j))
799 END DO
800 END DO
801 END DO k_loop
802 END DO t_loop1
803
804# if defined AGE_MEAN && defined T_PASSIVE
805!
806! If inert passive tracer and Mean Age, compute age concentration (even
807! inert index) forced by the right-hand-side term that is concentration
808! of an associated conservative passive tracer (odd inert index).
809! Implemented and tested by W.G. Zhang and J. Wilkin. (m Tunits)
810!
811 DO itrc=1,npt,2
812 IF (.not.((tl_hadvection(inert(itrc),ng)%MPDATA).or. &
813 & (tl_hadvection(inert(itrc),ng)%HSIMT))) THEN
814 IF (iic(ng).eq.ntfirst(ng)) THEN
815 cff=0.5_r8*dt(ng)
816 ELSE
817 gamma=1.0_r8/6.0_r8
818 cff=(1.0_r8-gamma)*dt(ng)
819 END IF
820 iage=inert(itrc+1) ! even inert tracer index
821 DO k=1,n(ng)
822 DO j=jstr,jend
823 DO i=istr,iend
824!^ t(i,j,k,3,iage)=t(i,j,k,3,iage)+ &
825!^ & cff*Hz(i,j,k)* &
826!^ & t(i,j,k,nnew,inert(itrc))
827!^
828 tl_t(i,j,k,3,iage)=tl_t(i,j,k,3,iage)+ &
829 & cff* &
830 & (hz(i,j,k)* &
831 & tl_t(i,j,k,nnew,inert(itrc))+ &
832 & tl_hz(i,j,k)* &
833 & t(i,j,k,nnew,inert(itrc)))
834 END DO
835 END DO
836 END DO
837 END IF
838 END DO
839# endif
840!
841!-----------------------------------------------------------------------
842! Compute time rate of change of intermediate tracer due to vertical
843! advection. Impose artificial continuity equation.
844!-----------------------------------------------------------------------
845!
846 j_loop1 : DO j=jstr,jend
847 t_loop2 : DO itrc=1,nt(ng)
848!
849 vadv_flux : IF (tl_vadvection(itrc,ng)%SPLINES) THEN
850!
851! Build conservative parabolic splines for the BASIC STATE vertical
852! derivatives "FC" of the tracer.
853!
854 DO i=istr,iend
855# ifdef NEUMANN
856 fc(i,0)=1.5_r8*t(i,j,1,nstp,itrc)
857 cf(i,1)=0.5_r8
858# else
859 fc(i,0)=2.0_r8*t(i,j,1,nstp,itrc)
860 cf(i,1)=1.0_r8
861# endif
862 END DO
863 DO k=1,n(ng)-1
864 DO i=istr,iend
865 cff=1.0_r8/(2.0_r8*hz(i,j,k)+ &
866 & hz(i,j,k+1)*(2.0_r8-cf(i,k)))
867 cf(i,k+1)=cff*hz(i,j,k)
868 fc(i,k)=cff*(3.0_r8*(hz(i,j,k )*t(i,j,k+1,nstp,itrc)+ &
869 & hz(i,j,k+1)*t(i,j,k ,nstp,itrc))- &
870 & hz(i,j,k+1)*fc(i,k-1))
871 END DO
872 END DO
873 DO i=istr,iend
874# ifdef NEUMANN
875 fc(i,n(ng))=(3.0_r8*t(i,j,n(ng),nstp,itrc)- &
876 & fc(i,n(ng)-1))/(2.0_r8-cf(i,n(ng)))
877# else
878 fc(i,n(ng))=(2.0_r8*t(i,j,n(ng),nstp,itrc)- &
879 & fc(i,n(ng)-1))/(1.0_r8-cf(i,n(ng)))
880# endif
881 END DO
882 DO k=n(ng)-1,0,-1
883 DO i=istr,iend
884 fc(i,k)=fc(i,k)-cf(i,k+1)*fc(i,k+1)
885 END DO
886 END DO
887!
888! Now the tangent linear spline code.
889!
890 DO i=istr,iend
891# ifdef NEUMANN
892!^ FC(i,0)=1.5_r8*t(i,j,1,nstp,itrc)
893!^
894 tl_fc(i,0)=1.5_r8*tl_t(i,j,1,nstp,itrc)
895 cf(i,1)=0.5_r8
896# else
897!^ FC(i,0)=2.0_r8*t(i,j,1,nstp,itrc)
898!^
899 tl_fc(i,0)=2.0_r8*tl_t(i,j,1,nstp,itrc)
900 cf(i,1)=1.0_r8
901# endif
902 END DO
903 DO k=1,n(ng)-1
904 DO i=istr,iend
905 cff=1.0_r8/(2.0_r8*hz(i,j,k)+ &
906 & hz(i,j,k+1)*(2.0_r8-cf(i,k)))
907 cf(i,k+1)=cff*hz(i,j,k)
908 tl_fc(i,k)=cff* &
909 & (3.0_r8*(hz(i,j,k )* &
910 & tl_t(i,j,k+1,nstp,itrc)+ &
911 & hz(i,j,k+1)* &
912 & tl_t(i,j,k ,nstp,itrc)+ &
913 & tl_hz(i,j,k )* &
914 & t(i,j,k+1,nstp,itrc)+ &
915 & tl_hz(i,j,k+1)* &
916 & t(i,j,k ,nstp,itrc))- &
917 & (tl_hz(i,j,k+1)*fc(i,k-1)+ &
918 & 2.0_r8*(tl_hz(i,j,k )+ &
919 & tl_hz(i,j,k+1))*fc(i,k)+ &
920 & tl_hz(i,j,k )*fc(i,k+1))- &
921 & hz(i,j,k+1)*tl_fc(i,k-1))
922 END DO
923 END DO
924 DO i=istr,iend
925# ifdef NEUMANN
926!^ FC(i,N(ng))=(3.0_r8*t(i,j,N(ng),nstp,itrc)- &
927!^ & FC(i,N(ng)-1))/(2.0_r8-CF(i,N(ng)))
928!^
929 tl_fc(i,n(ng))=(3.0_r8*tl_t(i,j,n(ng),nstp,itrc)- &
930 & tl_fc(i,n(ng)-1))/ &
931 & (2.0_r8-cf(i,n(ng)))
932# else
933!^ FC(i,N(ng))=(2.0_r8*t(i,j,N(ng),nstp,itrc)- &
934!^ & FC(i,N(ng)-1))/(1.0_r8-CF(i,N(ng)))
935!^
936 tl_fc(i,n(ng))=(2.0_r8*tl_t(i,j,n(ng),nstp,itrc)- &
937 & tl_fc(i,n(ng)-1))/ &
938 & (1.0_r8-cf(i,n(ng)))
939# endif
940 END DO
941 DO k=n(ng)-1,0,-1
942 DO i=istr,iend
943!^ FC(i,k)=FC(i,k)-CF(i,k+1)*FC(i,k+1)
944!^
945 tl_fc(i,k)=tl_fc(i,k)-cf(i,k+1)*tl_fc(i,k+1)
946!^ FC(i,k+1)=W(i,j,k+1)*FC(i,k+1)
947!^
948 tl_fc(i,k+1)=tl_w(i,j,k+1)*fc(i,k+1)+ &
949 & w(i,j,k+1)*tl_fc(i,k+1)
950 END DO
951 END DO
952 DO i=istr,iend
953!^ FC(i,N(ng))=0.0_r8
954!^
955 tl_fc(i,n(ng))=0.0_r8
956!^ FC(i,0)=0.0_r8
957!^
958 tl_fc(i,0)=0.0_r8
959 END DO
960!
961! Now complete the computation of the flux array FC.
962!
963 DO k=n(ng)-1,0,-1
964 DO i=istr,iend
965 fc(i,k+1)=w(i,j,k+1)*fc(i,k+1)
966 END DO
967 END DO
968 DO i=istr,iend
969 fc(i,n(ng))=0.0_r8
970 fc(i,0)=0.0_r8
971 END DO
972!
973 ELSE IF (tl_vadvection(itrc,ng)%AKIMA4) THEN
974!
975! Fourth-order, Akima vertical advective flux.
976!
977 DO k=1,n(ng)-1
978 DO i=istr,iend
979 fc(i,k)=t(i,j,k+1,nstp,itrc)- &
980 & t(i,j,k ,nstp,itrc)
981 tl_fc(i,k)=tl_t(i,j,k+1,nstp,itrc)- &
982 & tl_t(i,j,k ,nstp,itrc)
983 END DO
984 END DO
985 DO i=istr,iend
986 fc(i,0)=fc(i,1)
987 tl_fc(i,0)=tl_fc(i,1)
988 fc(i,n(ng))=fc(i,n(ng)-1)
989 tl_fc(i,n(ng))=tl_fc(i,n(ng)-1)
990 END DO
991 DO k=1,n(ng)
992 DO i=istr,iend
993 cff=2.0_r8*fc(i,k)*fc(i,k-1)
994 tl_cff=2.0_r8*(tl_fc(i,k)*fc(i,k-1)+ &
995 & fc(i,k)*tl_fc(i,k-1))
996 IF (cff.gt.eps) THEN
997 cf(i,k)=cff/(fc(i,k)+fc(i,k-1))
998 tl_cf(i,k)=((fc(i,k)+fc(i,k-1))*tl_cff- &
999 & cff*(tl_fc(i,k)+tl_fc(i,k-1)))/ &
1000 & ((fc(i,k)+fc(i,k-1))*(fc(i,k)+fc(i,k-1)))
1001 ELSE
1002 cf(i,k)=0.0_r8
1003 tl_cf(i,k)=0.0_r8
1004 END IF
1005 END DO
1006 END DO
1007 cff1=1.0_r8/3.0_r8
1008 DO k=1,n(ng)-1
1009 DO i=istr,iend
1010 fc(i,k)=w(i,j,k)* &
1011 & 0.5_r8*(t(i,j,k ,nstp,itrc)+ &
1012 & t(i,j,k+1,nstp,itrc)- &
1013 & cff1*(cf(i,k+1)-cf(i,k)))
1014 tl_fc(i,k)=0.5_r8* &
1015 & (tl_w(i,j,k)* &
1016 & (t(i,j,k ,nstp,itrc)+ &
1017 & t(i,j,k+1,nstp,itrc)- &
1018 & cff1*(cf(i,k+1)-cf(i,k)))+ &
1019 & w(i,j,k)* &
1020 & (tl_t(i,j,k ,nstp,itrc)+ &
1021 & tl_t(i,j,k+1,nstp,itrc)- &
1022 & cff1*(tl_cf(i,k+1)-tl_cf(i,k))))
1023 END DO
1024 END DO
1025 DO i=istr,iend
1026 fc(i,0)=0.0_r8
1027 tl_fc(i,0)=0.0_r8
1028 fc(i,n(ng))=0.0_r8
1029 tl_fc(i,n(ng))=0.0_r8
1030 END DO
1031!
1032 ELSE IF (tl_vadvection(itrc,ng)%CENTERED2) THEN
1033!
1034! Second-order, central differences vertical advective flux.
1035!
1036 DO k=1,n(ng)-1
1037 DO i=istr,iend
1038 fc(i,k)=w(i,j,k)* &
1039 & 0.5_r8*(t(i,j,k ,nstp,itrc)+ &
1040 & t(i,j,k+1,nstp,itrc))
1041 tl_fc(i,k)=0.5_r8* &
1042 & (tl_w(i,j,k)* &
1043 & (t(i,j,k ,nstp,itrc)+ &
1044 & t(i,j,k+1,nstp,itrc))+ &
1045 & w(i,j,k)* &
1046 & (tl_t(i,j,k ,nstp,itrc)+ &
1047 & tl_t(i,j,k+1,nstp,itrc)))
1048 END DO
1049 END DO
1050 DO i=istr,iend
1051 fc(i,0)=0.0_r8
1052 tl_fc(i,0)=0.0_r8
1053 fc(i,n(ng))=0.0_r8
1054 tl_fc(i,n(ng))=0.0_r8
1055 END DO
1056!
1057 ELSE IF ((tl_vadvection(itrc,ng)%MPDATA).or. &
1058 & (tl_vadvection(itrc,ng)%HSIMT)) THEN
1059!
1060! First_order, upstream differences vertical advective flux.
1061!
1062 DO k=1,n(ng)-1
1063 DO i=istr,iend
1064 cff1=max(w(i,j,k),0.0_r8)
1065 cff2=min(w(i,j,k),0.0_r8)
1066 tl_cff1=(0.5_r8+sign(0.5_r8, w(i,j,k)))*tl_w(i,j,k)
1067 tl_cff2=(0.5_r8+sign(0.5_r8,-w(i,j,k)))*tl_w(i,j,k)
1068!^ FC(i,k)=cff1*t(i,j,k ,nstp,itrc)+ &
1069!^ & cff2*t(i,j,k+1,nstp,itrc)
1070!^
1071 tl_fc(i,k)=tl_cff1*t(i,j,k ,nstp,itrc)+ &
1072 & cff1*tl_t(i,j,k ,nstp,itrc)+ &
1073 & tl_cff2*t(i,j,k+1,nstp,itrc)+ &
1074 & cff2*tl_t(i,j,k+1,nstp,itrc)
1075 END DO
1076 END DO
1077 DO i=istr,iend
1078!^ FC(i,0)=0.0_r8
1079!^
1080 tl_fc(i,0)=0.0_r8
1081!^ FC(i,N(ng))=0.0_r8
1082!^
1083 tl_fc(i,n(ng))=0.0_r8
1084 END DO
1085!
1086 ELSE IF ((tl_vadvection(itrc,ng)%CENTERED4).or. &
1087 & (tl_vadvection(itrc,ng)%SPLIT_U3)) THEN
1088!
1089! Fourth-order, central differences vertical advective flux.
1090!
1091 cff1=0.5_r8
1092 cff2=7.0_r8/12.0_r8
1093 cff3=1.0_r8/12.0_r8
1094 DO k=2,n(ng)-2
1095 DO i=istr,iend
1096 fc(i,k)=w(i,j,k)* &
1097 & (cff2*(t(i,j,k ,nstp,itrc)+ &
1098 & t(i,j,k+1,nstp,itrc))- &
1099 & cff3*(t(i,j,k-1,nstp,itrc)+ &
1100 & t(i,j,k+2,nstp,itrc)))
1101 tl_fc(i,k)=tl_w(i,j,k)* &
1102 & (cff2*(t(i,j,k ,nstp,itrc)+ &
1103 & t(i,j,k+1,nstp,itrc))- &
1104 & cff3*(t(i,j,k-1,nstp,itrc)+ &
1105 & t(i,j,k+2,nstp,itrc)))+ &
1106 & w(i,j,k)* &
1107 & (cff2*(tl_t(i,j,k ,nstp,itrc)+ &
1108 & tl_t(i,j,k+1,nstp,itrc))- &
1109 & cff3*(tl_t(i,j,k-1,nstp,itrc)+ &
1110 & tl_t(i,j,k+2,nstp,itrc)))
1111 END DO
1112 END DO
1113 DO i=istr,iend
1114 fc(i,0)=0.0_r8
1115 tl_fc(i,0)=0.0_r8
1116 fc(i,1)=w(i,j,1)* &
1117 & (cff1*t(i,j,1,nstp,itrc)+ &
1118 & cff2*t(i,j,2,nstp,itrc)- &
1119 & cff3*t(i,j,3,nstp,itrc))
1120 tl_fc(i,1)=tl_w(i,j,1)* &
1121 & (cff1*t(i,j,1,nstp,itrc)+ &
1122 & cff2*t(i,j,2,nstp,itrc)- &
1123 & cff3*t(i,j,3,nstp,itrc))+ &
1124 & w(i,j,1)* &
1125 & (cff1*tl_t(i,j,1,nstp,itrc)+ &
1126 & cff2*tl_t(i,j,2,nstp,itrc)- &
1127 & cff3*tl_t(i,j,3,nstp,itrc))
1128 fc(i,n(ng)-1)=w(i,j,n(ng)-1)* &
1129 & (cff1*t(i,j,n(ng) ,nstp,itrc)+ &
1130 & cff2*t(i,j,n(ng)-1,nstp,itrc)- &
1131 & cff3*t(i,j,n(ng)-2,nstp,itrc))
1132 tl_fc(i,n(ng)-1)=tl_w(i,j,n(ng)-1)* &
1133 & (cff1*t(i,j,n(ng) ,nstp,itrc)+ &
1134 & cff2*t(i,j,n(ng)-1,nstp,itrc)- &
1135 & cff3*t(i,j,n(ng)-2,nstp,itrc))+ &
1136 & w(i,j,n(ng)-1)* &
1137 & (cff1*tl_t(i,j,n(ng) ,nstp,itrc)+ &
1138 & cff2*tl_t(i,j,n(ng)-1,nstp,itrc)- &
1139 & cff3*tl_t(i,j,n(ng)-2,nstp,itrc))
1140 fc(i,n(ng))=0.0_r8
1141 tl_fc(i,n(ng))=0.0_r8
1142 END DO
1143 END IF vadv_flux
1144!
1145! Compute artificial continuity equation and load it into private
1146! array DC (1/m). It is imposed to preserve tracer constancy. It is
1147! not the same for all the tracer advection schemes because of the
1148! Gamma value.
1149!
1150 IF ((vadvection(itrc,ng)%MPDATA).or. &
1151 & (vadvection(itrc,ng)%HSIMT)) THEN
1152 gamma=0.5_r8
1153 ELSE
1154 gamma=1.0_r8/6.0_r8
1155 END IF
1156# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
1157 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
1158# else
1159 IF (iic(ng).eq.ntfirst(ng)) THEN
1160# endif
1161 cff=0.5_r8*dt(ng)
1162 ELSE
1163 cff=(1.0_r8-gamma)*dt(ng)
1164 END IF
1165 DO k=1,n(ng)
1166 DO i=istr,iend
1167 dc(i,k)=1.0_r8/(hz(i,j,k)- &
1168 & cff*pm(i,j)*pn(i,j)* &
1169 & (huon(i+1,j,k)-huon(i,j,k)+ &
1170 & hvom(i,j+1,k)-hvom(i,j,k)+ &
1171 & (w(i,j,k)-w(i,j,k-1))))
1172 tl_dc(i,k)=-dc(i,k)*dc(i,k)* &
1173 & (tl_hz(i,j,k)- &
1174 & cff*pm(i,j)*pn(i,j)* &
1175 & (tl_huon(i+1,j,k)-tl_huon(i,j,k)+ &
1176 & tl_hvom(i,j+1,k)-tl_hvom(i,j,k)+ &
1177 & (tl_w(i,j,k)-tl_w(i,j,k-1))))
1178 END DO
1179 END DO
1180!
1181! Time-step vertical advection of tracers (Tunits). Impose artificial
1182! continuity equation.
1183!
1184! WARNING: t(:,:,:,3,itrc) at this point should be in units of
1185! ======= m Tunits. But, t(:,:,:,3,itrc) is read from a BASIC
1186! STATE file and is in Tunits, so we need to multiply
1187! by level thickness (Hz).
1188!
1189 DO k=1,n(ng)
1190 DO i=istr,iend
1191 cff1=cff*pm(i,j)*pn(i,j)
1192!^ t(i,j,k,3,itrc)=DC(i,k)* &
1193!^ & (t(i,j,k,3,itrc)- &
1194!^ & cff1*(FC(i,k)-FC(i,k-1)))
1195!^
1196 tl_t(i,j,k,3,itrc)=tl_dc(i,k)* &
1197 & (t(i,j,k,3,itrc)*hz(i,j,k)- &
1198 & cff1*(fc(i,k)-fc(i,k-1)))+ &
1199 & dc(i,k)* &
1200 & (tl_t(i,j,k,3,itrc)- &
1201 & cff1*(tl_fc(i,k)-tl_fc(i,k-1)))
1202 END DO
1203 END DO
1204 END DO t_loop2
1205 END DO j_loop1
1206!
1207!-----------------------------------------------------------------------
1208! Start computation of tracers at n+1 time-step, t(i,j,k,nnew,itrc).
1209!-----------------------------------------------------------------------
1210!
1211! Compute vertical diffusive fluxes "FC" of the tracer fields at
1212! current time step n, and at horizontal RHO-points and vertical
1213! W-points. Notice that the vertical diffusion coefficients for
1214! passive tracers is the same as that for salinity (ltrc=NAT).
1215!
1216 DO j=jstr,jend
1217 cff3=dt(ng)*(1.0_r8-lambda)
1218 DO itrc=1,nt(ng)
1219 ltrc=min(nat,itrc)
1220 DO k=1,n(ng)-1
1221 DO i=istr,iend
1222 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1223 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k))
1224!^ FC(i,k)=cff3*cff*Akt(i,j,k,ltrc)* &
1225!^ & (t(i,j,k+1,nstp,itrc)- &
1226!^ & t(i,j,k ,nstp,itrc))
1227!^
1228 tl_fc(i,k)=cff3* &
1229 & (cff*(tl_akt(i,j,k,ltrc)* &
1230 & (t(i,j,k+1,nstp,itrc)- &
1231 & t(i,j,k ,nstp,itrc))+ &
1232 & akt(i,j,k,ltrc)* &
1233 & (tl_t(i,j,k+1,nstp,itrc)- &
1234 & tl_t(i,j,k ,nstp,itrc)))+ &
1235 & tl_cff*(akt(i,j,k,ltrc)* &
1236 & (t(i,j,k+1,nstp,itrc)- &
1237 & t(i,j,k,nstp,itrc))))
1238 END DO
1239 END DO
1240
1241# ifdef LMD_NONLOCAL_NOT_YET
1242!
1243! Add in the nonlocal transport flux for unstable (convective)
1244! forcing conditions into matrix FC when using the Large et al.
1245! KPP scheme. The nonlocal transport is only applied to active
1246! tracers.
1247!
1248 IF (itrc.le.nat) THEN
1249 DO k=1,n(ng)-1
1250 DO i=istr,iend
1251!^ FC(i,k)=FC(i,k)- &
1252!^ & dt(ng)*Akt(i,j,k,itrc)*ghats(i,j,k,itrc)
1253!^
1254 tl_fc(i,k)=tl_fc(i,k)- &
1255 & dt(ng)*(tl_akt(i,j,k,itrc)* &
1256 & ghats(i,j,k,itrc)+ &
1257 & akt(i,j,k,itrc)* &
1258 & tl_ghats(i,j,k,itrc))
1259 END DO
1260 END DO
1261 END IF
1262# endif
1263# ifdef SOLAR_SOURCE
1264!
1265! Add in incoming solar radiation at interior W-points using decay
1266! decay penetration function based on Jerlov water type.
1267!
1268 IF (itrc.eq.itemp) THEN
1269 DO k=1,n(ng)-1
1270 DO i=istr,iend
1271!^ FC(i,k)=FC(i,k)+ &
1272!^ & dt(ng)*srflx(i,j)* &
1273# ifdef WET_DRY
1274!^ & rmask_wet(i,j)* &
1275# endif
1276!^ & swdk(i,j,k)
1277!^
1278 tl_fc(i,k)=tl_fc(i,k)+ &
1279 & dt(ng)*srflx(i,j)* &
1280# ifdef WET_DRY_NOT_YET
1281 & rmask_wet(i,j)* &
1282# endif
1283 & tl_swdk(i,j,k)
1284 END DO
1285 END DO
1286 END IF
1287# endif
1288!
1289! Apply bottom and surface tracer flux conditions.
1290!
1291 DO i=istr,iend
1292!^ FC(i,0)=dt(ng)*btflx(i,j,itrc)
1293!^
1294 tl_fc(i,0)=dt(ng)*tl_btflx(i,j,itrc)
1295!^ FC(i,N(ng))=dt(ng)*stflx(i,j,itrc)
1296!^
1297 tl_fc(i,n(ng))=dt(ng)*tl_stflx(i,j,itrc)
1298 END DO
1299!
1300! Compute new tracer field (m Tunits).
1301!
1302 DO k=1,n(ng)
1303 DO i=istr,iend
1304 cff1=hz(i,j,k)*t(i,j,k,nstp,itrc)
1305 tl_cff1=tl_hz(i,j,k)*t(i,j,k,nstp,itrc)+ &
1306 & hz(i,j,k)*tl_t(i,j,k,nstp,itrc)
1307 cff2=fc(i,k)-fc(i,k-1)
1308 tl_cff2=tl_fc(i,k)-tl_fc(i,k-1)
1309!^ t(i,j,k,nnew,itrc)=cff1+cff2
1310!^
1311 tl_t(i,j,k,nnew,itrc)=tl_cff1+tl_cff2
1312# ifdef DIAGNOSTICS_TS
1313!! DiaTwrk(i,j,k,itrc,iTrate)=cff1
1314!! DiaTwrk(i,j,k,itrc,iTvdif)=cff2
1315# endif
1316 END DO
1317 END DO
1318 END DO
1319 END DO
1320# endif /* !TS_FIXED */
1321!
1322!=======================================================================
1323! 3D momentum equation in the XI-direction.
1324!=======================================================================
1325!
1326! Compute U-component viscous vertical momentum fluxes "FC" at
1327! current time-step n, and at horizontal U-points and vertical
1328! W-points.
1329!
1330 j_loop2: DO j=jstr,jend
1331 cff3=dt(ng)*(1.0_r8-lambda)
1332 DO k=1,n(ng)-1
1333 DO i=istru,iend
1334 cff=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
1335 & z_r(i,j,k )-z_r(i-1,j,k ))
1336 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)+tl_z_r(i-1,j,k+1)- &
1337 & tl_z_r(i,j,k )-tl_z_r(i-1,j,k ))
1338!^ FC(i,k)=cff3*cff*(u(i,j,k+1,nstp)-u(i,j,k,nstp))* &
1339!^ & (Akv(i,j,k)+Akv(i-1,j,k))
1340!^
1341 tl_fc(i,k)=cff3* &
1342 & (cff*((tl_u(i,j,k+1,nstp)-tl_u(i,j,k,nstp))* &
1343 & (akv(i,j,k)+akv(i-1,j,k))+ &
1344 & (u(i,j,k+1,nstp)-u(i,j,k,nstp))* &
1345 & (tl_akv(i,j,k)+tl_akv(i-1,j,k)))+ &
1346 & tl_cff*(u(i,j,k+1,nstp)-u(i,j,k,nstp))* &
1347 & (akv(i,j,k)+akv(i-1,j,k)))
1348 END DO
1349 END DO
1350!
1351! Apply bottom and surface stresses, if so is prescribed.
1352!
1353 DO i=istru,iend
1354# ifdef BODYFORCE
1355!^ FC(i,0)=0.0_r8
1356!^
1357 tl_fc(i,0)=0.0_r8
1358!^ FC(i,N(ng))=0.0_r8
1359!^
1360 tl_fc(i,n(ng))=0.0_r8
1361# else
1362!^ FC(i,0)=dt(ng)*bustr(i,j)
1363!^
1364 tl_fc(i,0)=dt(ng)*tl_bustr(i,j)
1365!^ FC(i,N(ng))=dt(ng)*sustr(i,j)
1366!^
1367 tl_fc(i,n(ng))=dt(ng)*tl_sustr(i,j)
1368# endif
1369 END DO
1370!
1371! Compute new U-momentum (m m/s).
1372!
1373 cff=dt(ng)*0.25_r8
1374 DO i=istru,iend
1375 dc(i,0)=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
1376 END DO
1377 indx=3-nrhs
1378# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
1379 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
1380# else
1381 IF (iic(ng).eq.ntfirst(ng)) THEN
1382# endif
1383 DO k=1,n(ng)
1384 DO i=istru,iend
1385!^ cff1=u(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))
1386!^
1387 tl_cff1=0.5_r8*(tl_u(i,j,k,nstp)* &
1388 & (hz(i,j,k)+hz(i-1,j,k))+ &
1389 & u(i,j,k,nstp)* &
1390 & (tl_hz(i,j,k)+tl_hz(i-1,j,k)))
1391!^ cff2=FC(i,k)-FC(i,k-1)
1392!^
1393 tl_cff2=tl_fc(i,k)-tl_fc(i,k-1)
1394!^ u(i,j,k,nnew)=cff1+cff2
1395!^
1396 tl_u(i,j,k,nnew)=tl_cff1+tl_cff2
1397# ifdef DIAGNOSTICS_UV
1398!! DO idiag=1,M3pgrd
1399!! DiaU3wrk(i,j,k,idiag)=0.0_r8
1400!! END DO
1401!! DiaU3wrk(i,j,k,M3vvis)=cff2
1402!! DiaU3wrk(i,j,k,M3rate)=cff1
1403# endif
1404 END DO
1405 END DO
1406# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
1407 ELSE IF (iic(ng).eq.(ntfirst(ng)+1).and.soinitial(ng)) THEN
1408# else
1409 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
1410# endif
1411 DO k=1,n(ng)
1412 DO i=istru,iend
1413!^ cff1=u(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))
1414!^
1415 tl_cff1=0.5_r8*(tl_u(i,j,k,nstp)* &
1416 & (hz(i,j,k)+hz(i-1,j,k))+ &
1417 & u(i,j,k,nstp)* &
1418 & (tl_hz(i,j,k)+tl_hz(i-1,j,k)))
1419!^ cff2=FC(i,k)-FC(i,k-1)
1420!^
1421 tl_cff2=tl_fc(i,k)-tl_fc(i,k-1)
1422 cff3=0.5_r8*dc(i,0)
1423!^ u(i,j,k,nnew)=cff1- &
1424!^ & cff3*ru(i,j,k,indx)+ &
1425!^ & cff2
1426!^
1427 tl_u(i,j,k,nnew)=tl_cff1- &
1428 & cff3*tl_ru(i,j,k,indx)+ &
1429 & tl_cff2
1430# ifdef DIAGNOSTICS_UV
1431!! DO idiag=1,M3pgrd
1432!! DiaU3wrk(i,j,k,idiag)=-cff3*DiaRU(i,j,k,indx,idiag)
1433!! END DO
1434!! DiaU3wrk(i,j,k,M3vvis)=cff2
1435# ifdef BODYFORCE
1436!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)- &
1437!! & cff3*DiaRU(i,j,k,indx,M3vvis)
1438# endif
1439!! DiaU3wrk(i,j,k,M3rate)=cff1
1440# endif
1441 END DO
1442 END DO
1443 ELSE
1444 cff1= 5.0_r8/12.0_r8
1445 cff2=16.0_r8/12.0_r8
1446 DO k=1,n(ng)
1447 DO i=istru,iend
1448!^ cff3=u(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))
1449!^
1450 tl_cff3=0.5_r8*(tl_u(i,j,k,nstp)* &
1451 & (hz(i,j,k)+hz(i-1,j,k))+ &
1452 & u(i,j,k,nstp)* &
1453 & (tl_hz(i,j,k)+tl_hz(i-1,j,k)))
1454!^ cff4=FC(i,k)-FC(i,k-1)
1455!^
1456 tl_cff4=tl_fc(i,k)-tl_fc(i,k-1)
1457!^ u(i,j,k,nnew)=cff3+ &
1458!^ & DC(i,0)*(cff1*ru(i,j,k,nrhs)- &
1459!^ & cff2*ru(i,j,k,indx))+ &
1460!^ & cff4
1461!^
1462 tl_u(i,j,k,nnew)=tl_cff3+ &
1463 & dc(i,0)*(cff1*tl_ru(i,j,k,nrhs)- &
1464 & cff2*tl_ru(i,j,k,indx))+ &
1465 & tl_cff4
1466# ifdef DIAGNOSTICS_UV
1467!! DO idiag=1,M3pgrd
1468!! DiaU3wrk(i,j,k,idiag)=DC(i,0)* &
1469!! & (cff1*DiaRU(i,j,k,nrhs,idiag)- &
1470!! & cff2*DiaRU(i,j,k,indx,idiag))
1471!! END DO
1472!! DiaU3wrk(i,j,k,M3vvis)=cff4
1473# ifdef BODYFORCE
1474!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+ &
1475!! & DC(i,0)* &
1476!! & (cff1*DiaRU(i,j,k,nrhs,M3vvis)- &
1477!! & cff2*DiaRU(i,j,k,indx,M3vvis))
1478# endif
1479!! DiaU3wrk(i,j,k,M3rate)=cff3
1480# endif
1481 END DO
1482 END DO
1483 END IF
1484!
1485!=======================================================================
1486! 3D momentum equation in the ETA-direction.
1487!=======================================================================
1488!
1489! Compute V-component viscous vertical momentum fluxes "FC" at
1490! current time-step n, and at horizontal V-points and vertical
1491! W-points.
1492!
1493 IF (j.ge.jstrv) THEN
1494 cff3=dt(ng)*(1.0_r8-lambda)
1495 DO k=1,n(ng)-1
1496 DO i=istr,iend
1497 cff=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
1498 & z_r(i,j,k )-z_r(i,j-1,k ))
1499 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)+tl_z_r(i,j-1,k+1)- &
1500 & tl_z_r(i,j,k )-tl_z_r(i,j-1,k ))
1501!^ FC(i,k)=cff3*cff*(v(i,j,k+1,nstp)-v(i,j,k,nstp))* &
1502!^ & (Akv(i,j,k)+Akv(i,j-1,k))
1503!^
1504 tl_fc(i,k)=cff3* &
1505 & (cff*((tl_v(i,j,k+1,nstp)-tl_v(i,j,k,nstp))* &
1506 & (akv(i,j,k)+akv(i,j-1,k))+ &
1507 & (v(i,j,k+1,nstp)-v(i,j,k,nstp))* &
1508 & (tl_akv(i,j,k)+tl_akv(i,j-1,k)))+ &
1509 & tl_cff*(v(i,j,k+1,nstp)-v(i,j,k,nstp))* &
1510 & (akv(i,j,k)+akv(i,j-1,k)))
1511 END DO
1512 END DO
1513!
1514! Apply bottom and surface stresses, if so is prescribed.
1515!
1516 DO i=istr,iend
1517# ifdef BODYFORCE
1518!^ FC(i,0)=0.0_r8
1519!^
1520 tl_fc(i,0)=0.0_r8
1521!^ FC(i,N(ng))=0.0_r8
1522!^
1523 tl_fc(i,n(ng))=0.0_r8
1524# else
1525!^ FC(i,0)=dt(ng)*bvstr(i,j)
1526!^
1527 tl_fc(i,0)=dt(ng)*tl_bvstr(i,j)
1528!^ FC(i,N(ng))=dt(ng)*svstr(i,j)
1529!^
1530 tl_fc(i,n(ng))=dt(ng)*tl_svstr(i,j)
1531# endif
1532 END DO
1533!
1534! Compute new V-momentum (m m/s).
1535!
1536 cff=dt(ng)*0.25_r8
1537 DO i=istr,iend
1538 dc(i,0)=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
1539 END DO
1540# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
1541 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
1542# else
1543 IF (iic(ng).eq.ntfirst(ng)) THEN
1544# endif
1545 DO k=1,n(ng)
1546 DO i=istr,iend
1547!^ cff1=v(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k))
1548!^
1549 tl_cff1=0.5_r8*(tl_v(i,j,k,nstp)* &
1550 & (hz(i,j,k)+hz(i,j-1,k))+ &
1551 & v(i,j,k,nstp)* &
1552 & (tl_hz(i,j,k)+tl_hz(i,j-1,k)))
1553!^ cff2=FC(i,k)-FC(i,k-1)
1554!^
1555 tl_cff2=tl_fc(i,k)-tl_fc(i,k-1)
1556!^ v(i,j,k,nnew)=cff1+cff2
1557!^
1558 tl_v(i,j,k,nnew)=tl_cff1+tl_cff2
1559# ifdef DIAGNOSTICS_UV
1560!! DO idiag=1,M3pgrd
1561!! DiaV3wrk(i,j,k,idiag)=0.0_r8
1562!! END DO
1563!! DiaV3wrk(i,j,k,M3vvis)=cff2
1564!! DiaV3wrk(i,j,k,M3rate)=cff1
1565# endif
1566 END DO
1567 END DO
1568# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
1569 ELSE IF (iic(ng).eq.(ntfirst(ng)+1).and.soinitial(ng)) THEN
1570# else
1571 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
1572# endif
1573 DO k=1,n(ng)
1574 DO i=istr,iend
1575!^ cff1=v(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k))
1576!^
1577 tl_cff1=0.5_r8*(tl_v(i,j,k,nstp)* &
1578 & (hz(i,j,k)+hz(i,j-1,k))+ &
1579 & v(i,j,k,nstp)* &
1580 & (tl_hz(i,j,k)+tl_hz(i,j-1,k)))
1581!^ cff2=FC(i,k)-FC(i,k-1)
1582!^
1583 tl_cff2=tl_fc(i,k)-tl_fc(i,k-1)
1584 cff3=0.5_r8*dc(i,0)
1585!^ v(i,j,k,nnew)=cff1- &
1586!^ & cff3*rv(i,j,k,indx)+ &
1587!^ & cff2
1588!^
1589 tl_v(i,j,k,nnew)=tl_cff1- &
1590 & cff3*tl_rv(i,j,k,indx)+ &
1591 & tl_cff2
1592# ifdef DIAGNOSTICS_UV
1593!! DO idiag=1,M3pgrd
1594!! DiaV3wrk(i,j,k,idiag)=-cff3*DiaRV(i,j,k,indx,idiag)
1595!! END DO
1596!! DiaV3wrk(i,j,k,M3vvis)=cff2
1597# ifdef BODYFORCE
1598!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)- &
1599!! & cff3*DiaRV(i,j,k,indx,M3vvis)
1600# endif
1601!! DiaV3wrk(i,j,k,M3rate)=cff1
1602# endif
1603 END DO
1604 END DO
1605 ELSE
1606 cff1= 5.0_r8/12.0_r8
1607 cff2=16.0_r8/12.0_r8
1608 DO k=1,n(ng)
1609 DO i=istr,iend
1610!^ cff3=v(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k))
1611!^
1612 tl_cff3=0.5_r8*(tl_v(i,j,k,nstp)* &
1613 & (hz(i,j,k)+hz(i,j-1,k))+ &
1614 & v(i,j,k,nstp)* &
1615 & (tl_hz(i,j,k)+tl_hz(i,j-1,k)))
1616!^ cff4=FC(i,k)-FC(i,k-1)
1617!^
1618 tl_cff4=tl_fc(i,k)-tl_fc(i,k-1)
1619!^ v(i,j,k,nnew)=cff3+ &
1620!^ & DC(i,0)*(cff1*rv(i,j,k,nrhs)- &
1621!^ & cff2*rv(i,j,k,indx))+ &
1622!^ & cff4
1623!^
1624 tl_v(i,j,k,nnew)=tl_cff3+ &
1625 & dc(i,0)*(cff1*tl_rv(i,j,k,nrhs)- &
1626 & cff2*tl_rv(i,j,k,indx))+ &
1627 & tl_cff4
1628# ifdef DIAGNOSTICS_UV
1629!! DO idiag=1,M3pgrd
1630!! DiaV3wrk(i,j,k,idiag)=DC(i,0)* &
1631!! & (cff1*DiaRV(i,j,k,nrhs,idiag)- &
1632!! & cff2*DiaRV(i,j,k,indx,idiag))
1633!! END DO
1634!! DiaV3wrk(i,j,k,M3vvis)=cff4
1635# ifdef BODYFORCE
1636!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+ &
1637!! & DC(i,0)* &
1638!! & (cff1*DiaRV(i,j,k,nrhs,M3vvis)- &
1639!! & cff2*DiaRV(i,j,k,indx,M3vvis))
1640# endif
1641!! DiaV3wrk(i,j,k,M3rate)=cff3
1642# endif
1643 END DO
1644 END DO
1645 END IF
1646 END IF
1647 END DO j_loop2
1648
1649# ifndef TS_FIXED
1650!
1651!=======================================================================
1652! Apply intermediate tracers lateral boundary conditions.
1653!=======================================================================
1654!
1655 ic=0
1656 DO itrc=1,nt(ng)
1657 IF (ltracerclm(itrc,ng).and.lnudgetclm(itrc,ng)) THEN
1658 ic=ic+1
1659 END IF
1660!^ CALL t3dbc_tile (ng, tile, itrc, ic, &
1661!^ & LBi, UBi, LBj, UBj, N(ng), NT(ng), &
1662!^ & IminS, ImaxS, JminS, JmaxS, &
1663!^ & nstp, 3, &
1664!^ & t)
1665!^
1666 CALL tl_t3dbc_tile (ng, tile, itrc, ic, &
1667 & lbi, ubi, lbj, ubj, n(ng), nt(ng), &
1668 & imins, imaxs, jmins, jmaxs, &
1669 & nstp, 3, &
1670 & tl_t)
1671
1672 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1673!^ CALL exchange_r3d_tile (ng, tile, &
1674!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1675!^ & t(:,:,:,3,itrc))
1676!^
1677 CALL exchange_r3d_tile (ng, tile, &
1678 & lbi, ubi, lbj, ubj, 1, n(ng), &
1679 & tl_t(:,:,:,3,itrc))
1680 END IF
1681 END DO
1682
1683# ifdef DISTRIBUTE
1684!^ CALL mp_exchange4d (ng, tile, iNLM, 1, &
1685!^ & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), &
1686!^ & NghostPoints, &
1687!^ & EWperiodic(ng), NSperiodic(ng), &
1688!^ & t(:,:,:,3,:))
1689!^
1690 CALL mp_exchange4d (ng, tile, itlm, 1, &
1691 & lbi, ubi, lbj, ubj, 1, n(ng), 1, nt(ng), &
1692 & nghostpoints, &
1693 & ewperiodic(ng), nsperiodic(ng), &
1694 & tl_t(:,:,:,3,:))
1695# endif
1696# endif
1697!
1698 RETURN
1699 END SUBROUTINE tl_pre_step3d_tile
1700#endif
1701 END MODULE tl_pre_step3d_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_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
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::tl_hadvection
type(t_adv), dimension(:,:), allocatable tl_hadvection
Definition mod_param.F:411
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::vadvection
type(t_adv), dimension(:,:), allocatable vadvection
Definition mod_param.F:404
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
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::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
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::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
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::ntfirst
integer, dimension(:), allocatable ntfirst
Definition mod_scalars.F:377
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
tl_pre_step3d_mod
Definition tl_pre_step3d.F:2
tl_pre_step3d_mod::tl_pre_step3d_tile
subroutine tl_pre_step3d_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, rmask, umask, vmask, rmask_wet, pm, pn, hz, tl_hz, huon, tl_huon, hvom, tl_hvom, z_r, tl_z_r, z_w, tl_z_w, tl_btflx, tl_bustr, tl_bvstr, tl_stflx, tl_sustr, tl_svstr, srflx, akt, tl_akt, akv, tl_akv, w, tl_w, tl_ru, tl_rv, t, tl_t, u, tl_u, v, tl_v)
Definition tl_pre_step3d.F:175
tl_pre_step3d_mod::tl_pre_step3d
subroutine, public tl_pre_step3d(ng, tile)
Definition tl_pre_step3d.F:44
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
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
tl_lmd_swfrac_tile
subroutine tl_lmd_swfrac_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, zscale, z, tl_z, tl_swdk)
Definition tl_lmd_swfrac.F:10