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