ROMS
Loading...
Searching...
No Matches
tl_t3dmix4_geo.h
Go to the documentation of this file.
1 MODULE tl_t3dmix4_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! This subroutine computes tangent linear horizontal biharmonic !
11! mixing of tracers along geopotential surfaces. !
12! !
13! BASIC STATE variables needed: diff4, Hz, t, z_r !
14! !
15!=======================================================================
16!
17 implicit none
18!
19 PRIVATE
20 PUBLIC tl_t3dmix4
21!
22 CONTAINS
23!
24!***********************************************************************
25 SUBROUTINE tl_t3dmix4 (ng, tile)
26!***********************************************************************
27!
28 USE mod_param
29#ifdef TS_MIX_CLIMA
30 USE mod_clima
31#endif
32#ifdef DIAGNOSTICS_TS
33!! USE mod_diags
34#endif
35 USE mod_grid
36 USE mod_mixing
37 USE mod_ocean
38 USE mod_stepping
39!
40! Imported variable declarations.
41!
42 integer, intent(in) :: ng, tile
43!
44! Local variable declarations.
45!
46 character (len=*), parameter :: myfile = &
47 & __FILE__
48!
49#include "tile.h"
50!
51#ifdef PROFILE
52 CALL wclock_on (ng, itlm, 28, __line__, myfile)
53#endif
54 CALL tl_t3dmix4_geo_tile (ng, tile, &
55 & lbi, ubi, lbj, ubj, &
56 & imins, imaxs, jmins, jmaxs, &
57 & nrhs(ng), nstp(ng), nnew(ng), &
58#ifdef MASKING
59 & grid(ng) % umask, &
60 & grid(ng) % vmask, &
61#endif
62#ifdef WET_DRY_NOT_YET
63 & grid(ng) % umask_wet, &
64 & grid(ng) % vmask_wet, &
65#endif
66 & grid(ng) % om_v, &
67 & grid(ng) % on_u, &
68 & grid(ng) % pm, &
69 & grid(ng) % pn, &
70 & grid(ng) % Hz, &
71 & grid(ng) % tl_Hz, &
72 & grid(ng) % z_r, &
73 & grid(ng) % tl_z_r, &
74#ifdef DIFF_3DCOEF
75# ifdef TS_U3ADV_SPLIT
76 & mixing(ng) % diff3d_u, &
77 & mixing(ng) % diff3d_v, &
78# else
79 & mixing(ng) % diff3d_r, &
80# endif
81#else
82 & mixing(ng) % diff4, &
83#endif
84#ifdef TS_MIX_CLIMA
85 & clima(ng) % tclm, &
86#endif
87#ifdef DIAGNOSTICS_TS
88!! & DIAGS(ng) % DiaTwrk, &
89#endif
90 & ocean(ng) % t, &
91 & ocean(ng) % tl_t)
92#ifdef PROFILE
93 CALL wclock_off (ng, itlm, 28, __line__, myfile)
94#endif
95!
96 RETURN
97 END SUBROUTINE tl_t3dmix4
98!
99!***********************************************************************
100 SUBROUTINE tl_t3dmix4_geo_tile (ng, tile, &
101 & LBi, UBi, LBj, UBj, &
102 & IminS, ImaxS, JminS, JmaxS, &
103 & nrhs, nstp, nnew, &
104#ifdef MASKING
105 & umask, vmask, &
106#endif
107#ifdef WET_DRY_NOT_YET
108 & umask_wet, vmask_wet, &
109#endif
110 & om_v, on_u, pm, pn, &
111 & Hz, tl_Hz, &
112 & z_r, tl_z_r, &
113#ifdef DIFF_3DCOEF
114# ifdef TS_U3ADV_SPLIT
115 & diff3d_u, diff3d_v, &
116# else
117 & diff3d_r, &
118# endif
119#else
120 & diff4, &
121#endif
122#ifdef TS_MIX_CLIMA
123 & tclm, &
124#endif
125#ifdef DIAGNOSTICS_TS
126!! & DiaTwrk, &
127#endif
128 & t, tl_t)
129!***********************************************************************
130!
131 USE mod_param
132 USE mod_ncparam
133 USE mod_scalars
134!
135! Imported variable declarations.
136!
137 integer, intent(in) :: ng, tile
138 integer, intent(in) :: LBi, UBi, LBj, UBj
139 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
140 integer, intent(in) :: nrhs, nstp, nnew
141
142#ifdef ASSUMED_SHAPE
143# ifdef MASKING
144 real(r8), intent(in) :: umask(LBi:,LBj:)
145 real(r8), intent(in) :: vmask(LBi:,LBj:)
146# endif
147# ifdef WET_DRY_NOT_YET
148 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
149 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
150# endif
151# ifdef DIFF_3DCOEF
152# ifdef TS_U3ADV_SPLIT
153 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
154 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
155# else
156 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
157# endif
158# else
159 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
160# endif
161 real(r8), intent(in) :: om_v(LBi:,LBj:)
162 real(r8), intent(in) :: on_u(LBi:,LBj:)
163 real(r8), intent(in) :: pm(LBi:,LBj:)
164 real(r8), intent(in) :: pn(LBi:,LBj:)
165 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
166 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
167 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
168# ifdef TS_MIX_CLIMA
169 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
170# endif
171 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
172 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
173# ifdef DIAGNOSTICS_TS
174 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
175# endif
176 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
177#else
178# ifdef MASKING
179 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
180 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
181# endif
182# ifdef WET_DRY_NOT_YET
183 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
184 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
185# endif
186# ifdef DIFF_3DCOEF
187# ifdef TS_U3ADV_SPLIT
188 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
189 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
190# else
191 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
192# endif
193# else
194 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
195# endif
196 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
197 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
198 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
199 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
200 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
201 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
202 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
203# ifdef TS_MIX_CLIMA
204 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
205# endif
206 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
207 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
208# ifdef DIAGNOSTICS_TS
209!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
210!! & NDT)
211# endif
212 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
213#endif
214!
215! Local variable declarations.
216!
217 integer :: Imin, Imax, Jmin, Jmax
218 integer :: i, itrc, j, k, k1, k2
219
220 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
221 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3, tl_cff4
222
223 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
224
225 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_LapT
226
227 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
228 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
229
230 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
231 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
232
233 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
234 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
235 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdz
237 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx
239
240 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FS
241 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTde
242 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdx
243 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdz
244 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dZde
245 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dZdx
246
247#include "set_bounds.h"
248!
249!-----------------------------------------------------------------------
250! Compute tangent linear horizontal biharmonic diffusion along
251! geopotential surfaces. The biharmonic operator is computed by
252! applying the harmonic operator twice.
253!-----------------------------------------------------------------------
254!
255! Set local I- and J-ranges.
256!
257 IF (ewperiodic(ng)) THEN
258 imin=istr-1
259 imax=iend+1
260 ELSE
261 imin=max(istr-1,1)
262 imax=min(iend+1,lm(ng))
263 END IF
264 IF (nsperiodic(ng)) THEN
265 jmin=jstr-1
266 jmax=jend+1
267 ELSE
268 jmin=max(jstr-1,1)
269 jmax=min(jend+1,mm(ng))
270 END IF
271!
272! Compute horizontal and vertical gradients associated with the
273! first rotated harmonic operator. Notice the recursive blocking
274! sequence. The vertical placement of the gradients is:
275!
276! dTdx,dTde(:,:,k1) k rho-points
277! dTdx,dTde(:,:,k2) k+1 rho-points
278! FS,dTdz(:,:,k1) k-1/2 W-points
279! FS,dTdz(:,:,k2) k+1/2 W-points
280!
281#ifdef TS_MIX_STABILITY
282! In order to increase stability, the rotated biharmonic is applied
283! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
284!
285#endif
286
287 t_loop : DO itrc=1,nt(ng)
288 k2=1
289 k_loop1 : DO k=0,n(ng)
290 k1=k2
291 k2=3-k1
292 IF (k.lt.n(ng)) THEN
293 DO j=jmin,jmax
294 DO i=imin,imax+1
295 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
296#ifdef MASKING
297 cff=cff*umask(i,j)
298#endif
299#ifdef WET_DRY_NOT_YET
300 cff=cff*umask_wet(i,j)
301#endif
302 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
303 & z_r(i-1,j,k+1))
304 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
305 & tl_z_r(i-1,j,k+1))
306#if defined TS_MIX_STABILITY
307 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
308 & t(i-1,j,k+1,nrhs,itrc))+ &
309 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
310 & t(i-1,j,k+1,nstp,itrc)))
311 tl_dtdx(i,j,k2)=cff*
312 & (0.75_r8*(tl_t(i ,j,k+1,nrhs,itrc)- &
313 & tl_t(i-1,j,k+1,nrhs,itrc))+ &
314 & 0.25_r8*(tl_t(i ,j,k+1,nstp,itrc)- &
315 & tl_t(i-1,j,k+1,nstp,itrc)))
316#elif defined TS_MIX_CLIMA
317 IF (ltracerclm(itrc,ng)) THEN
318 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
319 & tclm(i ,j,k+1,itrc))- &
320 & (t(i-1,j,k+1,nrhs,itrc)- &
321 & tclm(i-1,j,k+1,itrc)))
322 ELSE
323 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
324 & t(i-1,j,k+1,nrhs,itrc))
325 END IF
326 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
327 & tl_t(i-1,j,k+1,nrhs,itrc))
328#else
329 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
330 & t(i-1,j,k+1,nrhs,itrc))
331 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
332 & tl_t(i-1,j,k+1,nrhs,itrc))
333#endif
334 END DO
335 END DO
336 DO j=jmin,jmax+1
337 DO i=imin,imax
338 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
339#ifdef MASKING
340 cff=cff*vmask(i,j)
341#endif
342#ifdef WET_DRY_NOT_YET
343 cff=cff*vmask_wet(i,j)
344#endif
345 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
346 & z_r(i,j-1,k+1))
347 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
348 & tl_z_r(i,j-1,k+1))
349#if defined TS_MIX_STABILITY
350 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
351 & t(i,j-1,k+1,nrhs,itrc))+ &
352 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
353 & t(i,j-1,k+1,nstp,itrc)))
354 tl_dtde(i,j,k2)=cff*
355 & (0.75_r8*(tl_t(i,j ,k+1,nrhs,itrc)- &
356 & tl_t(i,j-1,k+1,nrhs,itrc))+ &
357 & 0.25_r8*(tl_t(i,j ,k+1,nstp,itrc)- &
358 & tl_t(i,j-1,k+1,nstp,itrc)))
359#elif defined TS_MIX_CLIMA
360 IF (ltracerclm(itrc,ng)) THEN
361 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
362 & tclm(i,j ,k+1,itrc))- &
363 & (t(i,j-1,k+1,nrhs,itrc)- &
364 & tclm(i,j-1,k+1,itrc)))
365 ELSE
366 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
367 & t(i,j-1,k+1,nrhs,itrc))
368 END IF
369 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
370 & tl_t(i,j-1,k+1,nrhs,itrc))
371#else
372 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
373 & t(i,j-1,k+1,nrhs,itrc))
374 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
375 & tl_t(i,j-1,k+1,nrhs,itrc))
376#endif
377 END DO
378 END DO
379 END IF
380 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
381 DO j=jmin-1,jmax+1
382 DO i=imin-1,imax+1
383 dtdz(i,j,k2)=0.0_r8
384 tl_dtdz(i,j,k2)=0.0_r8
385 fs(i,j,k2)=0.0_r8
386 tl_fs(i,j,k2)=0.0_r8
387 END DO
388 END DO
389 ELSE
390 DO j=jmin-1,jmax+1
391 DO i=imin-1,imax+1
392 cff=1.0_r8/(z_r(i,j,k+1)- &
393 & z_r(i,j,k ))
394 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
395 & tl_z_r(i,j,k ))
396#if defined TS_MIX_STABILITY
397 dtdz(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
398 & t(i,j,k ,nrhs,itrc))+ &
399 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
400 & t(i,j,k ,nstp,itrc)))
401 tl_dtdz(i,j,k2)=tl_cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
402 & t(i,j,k ,nrhs,itrc))+ &
403 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
404 & t(i,j,k ,nstp,itrc)))+&
405 & cff*(0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
406 & tl_t(i,j,k ,nrhs,itrc))+ &
407 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
408 & tl_t(i,j,k ,nstp,itrc)))
409#elif defined TS_MIX_CLIMA
410 IF (ltracerclm(itrc,ng)) THEN
411 dtdz(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
412 & tclm(i,j,k+1,itrc))- &
413 & (t(i,j,k ,nrhs,itrc)- &
414 & tclm(i,j,k ,itrc)))
415 tl_dtdz(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
416 & tclm(i,j,k+1,itrc))- &
417 & (t(i,j,k ,nrhs,itrc)- &
418 & tclm(i,j,k ,itrc)))+ &
419 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
420 & tl_t(i,j,k ,nrhs,itrc))
421
422 ELSE
423 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
424 & t(i,j,k ,nrhs,itrc))
425 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
426 & t(i,j,k ,nrhs,itrc))+ &
427 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
428 & tl_t(i,j,k ,nrhs,itrc))
429 END IF
430#else
431 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
432 & t(i,j,k ,nrhs,itrc))
433 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
434 & t(i,j,k ,nrhs,itrc))+ &
435 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
436 & tl_t(i,j,k ,nrhs,itrc))
437#endif
438 END DO
439 END DO
440 END IF
441 IF (k.gt.0) THEN
442 DO j=jmin,jmax
443 DO i=imin,imax+1
444#ifdef DIFF_3DCOEF
445# ifdef TS_U3ADV_SPLIT
446 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
447# else
448 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
449 & on_u(i,j)
450# endif
451#else
452 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
453 & on_u(i,j)
454#endif
455 fx(i,j)=cff* &
456 & (hz(i,j,k)+hz(i-1,j,k))* &
457 & (dtdx(i,j,k1)- &
458 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
459 & (dtdz(i-1,j,k1)+ &
460 & dtdz(i ,j,k2))+ &
461 & max(dzdx(i,j,k1),0.0_r8)* &
462 & (dtdz(i-1,j,k2)+ &
463 & dtdz(i ,j,k1))))
464 tl_fx(i,j)=cff* &
465 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
466 & (dtdx(i,j,k1)- &
467 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
468 & (dtdz(i-1,j,k1)+ &
469 & dtdz(i ,j,k2))+ &
470 & max(dzdx(i,j,k1),0.0_r8)* &
471 & (dtdz(i-1,j,k2)+ &
472 & dtdz(i ,j,k1))))+ &
473 & (hz(i,j,k)+hz(i-1,j,k))* &
474 & (tl_dtdx(i,j,k1)- &
475 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
476 & (tl_dtdz(i-1,j,k1)+ &
477 & tl_dtdz(i ,j,k2))+ &
478 & max(dzdx(i,j,k1),0.0_r8)* &
479 & (tl_dtdz(i-1,j,k2)+ &
480 & tl_dtdz(i ,j,k1)))- &
481 & 0.5_r8*((0.5_r8+ &
482 & sign(0.5_r8,-dzdx(i,j,k1)))* &
483 & tl_dzdx(i,j,k1)* &
484 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
485 & (0.5_r8+ &
486 & sign(0.5_r8, dzdx(i,j,k1)))* &
487 & tl_dzdx(i,j,k1)* &
488 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))
489 END DO
490 END DO
491 DO j=jmin,jmax+1
492 DO i=imin,imax
493#ifdef DIFF_3DCOEF
494# ifdef TS_U3ADV_SPLIT
495 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
496# else
497 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
498 & om_v(i,j)
499# endif
500#else
501 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
502 & om_v(i,j)
503#endif
504 fe(i,j)=cff* &
505 & (hz(i,j,k)+hz(i,j-1,k))* &
506 & (dtde(i,j,k1)- &
507 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
508 & (dtdz(i,j-1,k1)+ &
509 & dtdz(i,j ,k2))+ &
510 & max(dzde(i,j,k1),0.0_r8)* &
511 & (dtdz(i,j-1,k2)+ &
512 & dtdz(i,j ,k1))))
513 tl_fe(i,j)=cff* &
514 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
515 & (dtde(i,j,k1)- &
516 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
517 & (dtdz(i,j-1,k1)+ &
518 & dtdz(i,j ,k2))+ &
519 & max(dzde(i,j,k1),0.0_r8)* &
520 & (dtdz(i,j-1,k2)+ &
521 & dtdz(i,j ,k1))))+ &
522 & (hz(i,j,k)+hz(i,j-1,k))* &
523 & (tl_dtde(i,j,k1)- &
524 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
525 & (tl_dtdz(i,j-1,k1)+ &
526 & tl_dtdz(i,j ,k2))+ &
527 & max(dzde(i,j,k1),0.0_r8)* &
528 & (tl_dtdz(i,j-1,k2)+ &
529 & tl_dtdz(i,j ,k1)))- &
530 & 0.5_r8*((0.5_r8+ &
531 & sign(0.5_r8,-dzde(i,j,k1)))* &
532 & tl_dzde(i,j,k1)* &
533 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
534 & (0.5_r8+ &
535 & sign(0.5_r8, dzde(i,j,k1)))* &
536 & tl_dzde(i,j,k1)* &
537 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))
538 END DO
539 END DO
540 IF (k.lt.n(ng)) THEN
541 DO j=jmin,jmax
542 DO i=imin,imax
543#ifdef DIFF_3DCOEF
544# ifdef TS_U3ADV_SPLIT
545 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
546 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
547 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
548 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
549# else
550 difx=0.5_r8*diff3d_r(i,j,k)
551 dife=difx
552# endif
553#else
554 difx=0.5_r8*diff4(i,j,itrc)
555 dife=difx
556#endif
557 cff1=min(dzdx(i ,j,k1),0.0_r8)
558 cff2=min(dzdx(i+1,j,k2),0.0_r8)
559 cff3=max(dzdx(i ,j,k2),0.0_r8)
560 cff4=max(dzdx(i+1,j,k1),0.0_r8)
561 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
562 & tl_dzdx(i ,j,k1)
563 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
564 & tl_dzdx(i+1,j,k2)
565 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
566 & tl_dzdx(i ,j,k2)
567 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
568 & tl_dzdx(i+1,j,k1)
569 fs(i,j,k2)=difx* &
570 & (cff1*(cff1*dtdz(i,j,k2)- &
571 & dtdx(i ,j,k1))+ &
572 & cff2*(cff2*dtdz(i,j,k2)- &
573 & dtdx(i+1,j,k2))+ &
574 & cff3*(cff3*dtdz(i,j,k2)- &
575 & dtdx(i ,j,k2))+ &
576 & cff4*(cff4*dtdz(i,j,k2)- &
577 & dtdx(i+1,j,k1)))
578 tl_fs(i,j,k2)=difx* &
579 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
580 & dtdx(i ,j,k1))+ &
581 & tl_cff2*(cff2*dtdz(i,j,k2)- &
582 & dtdx(i+1,j,k2))+ &
583 & tl_cff3*(cff3*dtdz(i,j,k2)- &
584 & dtdx(i ,j,k2))+ &
585 & tl_cff4*(cff4*dtdz(i,j,k2)- &
586 & dtdx(i+1,j,k1))+ &
587 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
588 & cff1*tl_dtdz(i,j,k2)- &
589 & tl_dtdx(i ,j,k1))+ &
590 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
591 & cff2*tl_dtdz(i,j,k2)- &
592 & tl_dtdx(i+1,j,k2))+ &
593 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
594 & cff3*tl_dtdz(i,j,k2)- &
595 & tl_dtdx(i ,j,k2))+ &
596 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
597 & cff4*tl_dtdz(i,j,k2)- &
598 & tl_dtdx(i+1,j,k1)))
599!
600 cff1=min(dzde(i,j ,k1),0.0_r8)
601 cff2=min(dzde(i,j+1,k2),0.0_r8)
602 cff3=max(dzde(i,j ,k2),0.0_r8)
603 cff4=max(dzde(i,j+1,k1),0.0_r8)
604 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
605 & tl_dzde(i,j ,k1)
606 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
607 & tl_dzde(i,j+1,k2)
608 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
609 & tl_dzde(i,j ,k2)
610 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
611 & tl_dzde(i,j+1,k1)
612 fs(i,j,k2)=fs(i,j,k2)+ &
613 & dife* &
614 & (cff1*(cff1*dtdz(i,j,k2)- &
615 & dtde(i,j ,k1))+ &
616 & cff2*(cff2*dtdz(i,j,k2)- &
617 & dtde(i,j+1,k2))+ &
618 & cff3*(cff3*dtdz(i,j,k2)- &
619 & dtde(i,j ,k2))+ &
620 & cff4*(cff4*dtdz(i,j,k2)- &
621 & dtde(i,j+1,k1)))
622 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
623 & dife* &
624 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
625 & dtde(i,j ,k1))+ &
626 & tl_cff2*(cff2*dtdz(i,j,k2)- &
627 & dtde(i,j+1,k2))+ &
628 & tl_cff3*(cff3*dtdz(i,j,k2)- &
629 & dtde(i,j ,k2))+ &
630 & tl_cff4*(cff4*dtdz(i,j,k2)- &
631 & dtde(i,j+1,k1))+ &
632 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
633 & cff1*tl_dtdz(i,j,k2)- &
634 & tl_dtde(i,j ,k1))+ &
635 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
636 & cff2*tl_dtdz(i,j,k2)- &
637 & tl_dtde(i,j+1,k2))+ &
638 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
639 & cff3*tl_dtdz(i,j,k2)- &
640 & tl_dtde(i,j ,k2))+ &
641 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
642 & cff4*tl_dtdz(i,j,k2)- &
643 & tl_dtde(i,j+1,k1)))
644 END DO
645 END DO
646 END IF
647!
648! Compute first harmonic operator, without mixing coefficient.
649! Multiply by the metrics of the second harmonic operator. Save
650! into work array "LapT".
651!
652 DO j=jmin,jmax
653 DO i=imin,imax
654 cff=pm(i,j)*pn(i,j)
655 cff1=1.0_r8/hz(i,j,k)
656 tl_cff1=-cff1*cff1*tl_hz(i,j,k)
657 lapt(i,j,k)=cff1*(cff* &
658 & (fx(i+1,j)-fx(i,j)+ &
659 & fe(i,j+1)-fe(i,j))+ &
660 & (fs(i,j,k2)-fs(i,j,k1)))
661 tl_lapt(i,j,k)=tl_cff1*(cff* &
662 & (fx(i+1,j)-fx(i,j)+ &
663 & fe(i,j+1)-fe(i,j))+ &
664 & (fs(i,j,k2)-fs(i,j,k1)))+ &
665 & cff1*(cff* &
666 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
667 & tl_fe(i,j+1)-tl_fe(i,j))+ &
668 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))
669 END DO
670 END DO
671 END IF
672 END DO k_loop1
673!
674! Apply boundary conditions (except periodic; closed or gradient)
675! to the first harmonic operator.
676!
677 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
678 IF (domain(ng)%Western_Edge(tile)) THEN
679 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
680 DO k=1,n(ng)
681 DO j=jmin,jmax
682 lapt(istr-1,j,k)=0.0_r8
683 tl_lapt(istr-1,j,k)=0.0_r8
684 END DO
685 END DO
686 ELSE
687 DO k=1,n(ng)
688 DO j=jmin,jmax
689 lapt(istr-1,j,k)=lapt(istr,j,k)
690 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
691 END DO
692 END DO
693 END IF
694 END IF
695 END IF
696!
697 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
698 IF (domain(ng)%Eastern_Edge(tile)) THEN
699 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
700 DO k=1,n(ng)
701 DO j=jmin,jmax
702 lapt(iend+1,j,k)=0.0_r8
703 tl_lapt(iend+1,j,k)=0.0_r8
704 END DO
705 END DO
706 ELSE
707 DO k=1,n(ng)
708 DO j=jmin,jmax
709 lapt(iend+1,j,k)=lapt(iend,j,k)
710 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
711 END DO
712 END DO
713 END IF
714 END IF
715 END IF
716!
717 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
718 IF (domain(ng)%Southern_Edge(tile)) THEN
719 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
720 DO k=1,n(ng)
721 DO i=imin,imax
722 lapt(i,jstr-1,k)=0.0_r8
723 tl_lapt(i,jstr-1,k)=0.0_r8
724 END DO
725 END DO
726 ELSE
727 DO k=1,n(ng)
728 DO i=imin,imax
729 lapt(i,jstr-1,k)=lapt(i,jstr,k)
730 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
731 END DO
732 END DO
733 END IF
734 END IF
735 END IF
736!
737 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
738 IF (domain(ng)%Northern_Edge(tile)) THEN
739 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
740 DO k=1,n(ng)
741 DO i=imin,imax
742 lapt(i,jend+1,k)=0.0_r8
743 tl_lapt(i,jend+1,k)=0.0_r8
744 END DO
745 END DO
746 ELSE
747 DO k=1,n(ng)
748 DO i=imin,imax
749 lapt(i,jend+1,k)=lapt(i,jend,k)
750 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
751 END DO
752 END DO
753 END IF
754 END IF
755 END IF
756!
757 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
758 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
759 IF (domain(ng)%SouthWest_Corner(tile)) THEN
760 DO k=1,n(ng)
761 lapt(istr-1,jstr-1,k)=0.5_r8* &
762 & (lapt(istr ,jstr-1,k)+ &
763 & lapt(istr-1,jstr ,k))
764 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
765 & (tl_lapt(istr ,jstr-1,k)+ &
766 & tl_lapt(istr-1,jstr ,k))
767 END DO
768 END IF
769 END IF
770
771 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
772 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
773 IF (domain(ng)%SouthEast_Corner(tile)) THEN
774 DO k=1,n(ng)
775 lapt(iend+1,jstr-1,k)=0.5_r8* &
776 & (lapt(iend ,jstr-1,k)+ &
777 & lapt(iend+1,jstr ,k))
778 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
779 & (tl_lapt(iend ,jstr-1,k)+ &
780 & tl_lapt(iend+1,jstr ,k))
781 END DO
782 END IF
783 END IF
784
785 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
786 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
787 IF (domain(ng)%NorthWest_Corner(tile)) THEN
788 DO k=1,n(ng)
789 lapt(istr-1,jend+1,k)=0.5_r8* &
790 & (lapt(istr ,jend+1,k)+ &
791 & lapt(istr-1,jend ,k))
792 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
793 & (tl_lapt(istr ,jend+1,k)+ &
794 & tl_lapt(istr-1,jend ,k))
795 END DO
796 END IF
797 END IF
798
799 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
800 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
801 IF (domain(ng)%NorthEast_Corner(tile)) THEN
802 DO k=1,n(ng)
803 lapt(iend+1,jend+1,k)=0.5_r8* &
804 & (lapt(iend ,jend+1,k)+ &
805 & lapt(iend+1,jend ,k))
806 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
807 & (tl_lapt(iend ,jend+1,k)+ &
808 & tl_lapt(iend+1,jend ,k))
809 END DO
810 END IF
811 END IF
812!
813! Compute tangent linear horizontal and vertical gradients associated
814! with the second rotated harmonic operator.
815!
816 k2=1
817 k_loop2: DO k=0,n(ng)
818 k1=k2
819 k2=3-k1
820 IF (k.lt.n(ng)) THEN
821 DO j=jstr,jend
822 DO i=istr,iend+1
823 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
824#ifdef MASKING
825 cff=cff*umask(i,j)
826#endif
827#ifdef WET_DRY_NOT_YET
828 cff=cff*umask_wet(i,j)
829#endif
830 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
831 & z_r(i-1,j,k+1))
832 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
833 & tl_z_r(i-1,j,k+1))
834 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
835 & lapt(i-1,j,k+1))
836 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
837 & tl_lapt(i-1,j,k+1))
838 END DO
839 END DO
840 DO j=jstr,jend+1
841 DO i=istr,iend
842 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
843#ifdef MASKING
844 cff=cff*vmask(i,j)
845#endif
846#ifdef WET_DRY_NOT_YET
847 cff=cff*vmask_wet(i,j)
848#endif
849 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
850 & z_r(i,j-1,k+1))
851 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
852 & tl_z_r(i,j-1,k+1))
853 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
854 & lapt(i,j-1,k+1))
855 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
856 & tl_lapt(i,j-1,k+1))
857 END DO
858 END DO
859 END IF
860 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
861 DO j=jstr-1,jend+1
862 DO i=istr-1,iend+1
863 dtdz(i,j,k2)=0.0_r8
864 tl_dtdz(i,j,k2)=0.0_r8
865 fs(i,j,k2)=0.0_r8
866 tl_fs(i,j,k2)=0.0_r8
867 END DO
868 END DO
869 ELSE
870 DO j=jstr-1,jend+1
871 DO i=istr-1,iend+1
872 cff=1.0_r8/(z_r(i,j,k+1)- &
873 & z_r(i,j,k ))
874 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
875 & tl_z_r(i,j,k ))
876 dtdz(i,j,k2)=cff*(lapt(i,j,k+1)- &
877 & lapt(i,j,k ))
878 tl_dtdz(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
879 & lapt(i,j,k ))+ &
880 & cff*(tl_lapt(i,j,k+1)- &
881 & tl_lapt(i,j,k ))
882 END DO
883 END DO
884 END IF
885!
886! Compute components of the rotated tracer flux (T m4/s) along
887! geopotential surfaces.
888!
889 IF (k.gt.0) THEN
890 DO j=jstr,jend
891 DO i=istr,iend+1
892#ifdef DIFF_3DCOEF
893# ifdef TS_U3ADV_SPLIT
894 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
895# else
896 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
897 & on_u(i,j)
898# endif
899#else
900 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
901 & on_u(i,j)
902#endif
903!^ FX(i,j)=cff*
904!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
905!^ & (dTdx(i ,j,k1)- &
906!^ & 0.5_r8*(MIN(dZdx(i,j,k1),0.0_r8)* &
907!^ & (dTdz(i-1,j,k1)+ &
908!^ & dTdz(i ,j,k2))+ &
909!^ & MAX(dZdx(i,j,k1),0.0_r8)* &
910!^ & (dTdz(i-1,j,k2)+ &
911!^ & dTdz(i ,j,k1))))
912!^
913 tl_fx(i,j)=cff* &
914 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
915 & (dtdx(i,j,k1)- &
916 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
917 & (dtdz(i-1,j,k1)+ &
918 & dtdz(i ,j,k2))+ &
919 & max(dzdx(i,j,k1),0.0_r8)* &
920 & (dtdz(i-1,j,k2)+ &
921 & dtdz(i ,j,k1))))+ &
922 & (hz(i,j,k)+hz(i-1,j,k))* &
923 & (tl_dtdx(i,j,k1)- &
924 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
925 & (tl_dtdz(i-1,j,k1)+ &
926 & tl_dtdz(i ,j,k2))+ &
927 & max(dzdx(i,j,k1),0.0_r8)* &
928 & (tl_dtdz(i-1,j,k2)+ &
929 & tl_dtdz(i ,j,k1)))- &
930 & 0.5_r8*((0.5_r8+ &
931 & sign(0.5_r8,-dzdx(i,j,k1)))* &
932 & tl_dzdx(i,j,k1)* &
933 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
934 & (0.5_r8+ &
935 & sign(0.5_r8, dzdx(i,j,k1)))* &
936 & tl_dzdx(i,j,k1)* &
937 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))
938 END DO
939 END DO
940 DO j=jstr,jend+1
941 DO i=istr,iend
942#ifdef DIFF_3DCOEF
943# ifdef TS_U3ADV_SPLIT
944 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
945# else
946 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
947 & om_v(i,j)
948# endif
949#else
950 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
951 & om_v(i,j)
952#endif
953!^ FE(i,j)=cff* &
954!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
955!^ & (dTde(i,j,k1)- &
956!^ & 0.5_r8*(MIN(dZde(i,j,k1),0.0_r8)* &
957!^ & (dTdz(i,j-1,k1)+ &
958!^ & dTdz(i,j ,k2))+ &
959!^ & MAX(dZde(i,j,k1),0.0_r8)* &
960!^ & (dTdz(i,j-1,k2)+ &
961!^ & dTdz(i,j ,k1))))
962!^
963 tl_fe(i,j)=cff* &
964 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
965 & (dtde(i,j,k1)- &
966 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
967 & (dtdz(i,j-1,k1)+ &
968 & dtdz(i,j ,k2))+ &
969 & max(dzde(i,j,k1),0.0_r8)* &
970 & (dtdz(i,j-1,k2)+ &
971 & dtdz(i,j ,k1))))+ &
972 & (hz(i,j,k)+hz(i,j-1,k))* &
973 & (tl_dtde(i,j,k1)- &
974 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
975 & (tl_dtdz(i,j-1,k1)+ &
976 & tl_dtdz(i,j ,k2))+ &
977 & max(dzde(i,j,k1),0.0_r8)* &
978 & (tl_dtdz(i,j-1,k2)+ &
979 & tl_dtdz(i,j ,k1)))- &
980 & 0.5_r8*((0.5_r8+ &
981 & sign(0.5_r8,-dzde(i,j,k1)))* &
982 & tl_dzde(i,j,k1)* &
983 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
984 & (0.5_r8+ &
985 & sign(0.5_r8, dzde(i,j,k1)))* &
986 & tl_dzde(i,j,k1)* &
987 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))
988 END DO
989 END DO
990 IF (k.lt.n(ng)) THEN
991 DO j=jstr,jend
992 DO i=istr,iend
993#ifdef DIFF_3DCOEF
994# ifdef TS_U3ADV_SPLIT
995 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
996 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
997 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
998 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
999# else
1000 difx=0.5_r8*diff3d_r(i,j,k)
1001 dife=difx
1002# endif
1003#else
1004 difx=0.5_r8*diff4(i,j,itrc)
1005 dife=difx
1006#endif
1007 cff1=min(dzdx(i ,j,k1),0.0_r8)
1008 cff2=min(dzdx(i+1,j,k2),0.0_r8)
1009 cff3=max(dzdx(i ,j,k2),0.0_r8)
1010 cff4=max(dzdx(i+1,j,k1),0.0_r8)
1011 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
1012 & tl_dzdx(i ,j,k1)
1013 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
1014 & tl_dzdx(i+1,j,k2)
1015 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
1016 & tl_dzdx(i ,j,k2)
1017 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
1018 & tl_dzdx(i+1,j,k1)
1019!^
1020!^ FS(i,j,k2)=difx* &
1021!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1022!^ & dTdx(i ,j,k1))+ &
1023!^ & cff2*(cff2*dTdz(i,j,k2)- &
1024!^ & dTdx(i+1,j,k2))+ &
1025!^ & cff3*(cff3*dTdz(i,j,k2)- &
1026!^ & dTdx(i ,j,k2))+ &
1027!^ & cff4*(cff4*dTdz(i,j,k2)- &
1028!^ & dTdx(i+1,j,k1)))
1029!^
1030 tl_fs(i,j,k2)=difx* &
1031 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1032 & dtdx(i ,j,k1))+ &
1033 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1034 & dtdx(i+1,j,k2))+ &
1035 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1036 & dtdx(i ,j,k2))+ &
1037 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1038 & dtdx(i+1,j,k1))+ &
1039 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1040 & cff1*tl_dtdz(i,j,k2)- &
1041 & tl_dtdx(i ,j,k1))+ &
1042 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1043 & cff2*tl_dtdz(i,j,k2)- &
1044 & tl_dtdx(i+1,j,k2))+ &
1045 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1046 & cff3*tl_dtdz(i,j,k2)- &
1047 & tl_dtdx(i ,j,k2))+ &
1048 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1049 & cff4*tl_dtdz(i,j,k2)- &
1050 & tl_dtdx(i+1,j,k1)))
1051!
1052 cff1=min(dzde(i,j ,k1),0.0_r8)
1053 cff2=min(dzde(i,j+1,k2),0.0_r8)
1054 cff3=max(dzde(i,j ,k2),0.0_r8)
1055 cff4=max(dzde(i,j+1,k1),0.0_r8)
1056 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
1057 & tl_dzde(i,j ,k1)
1058 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
1059 & tl_dzde(i,j+1,k2)
1060 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
1061 & tl_dzde(i,j ,k2)
1062 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
1063 & tl_dzde(i,j+1,k1)
1064!^ FS(i,j,k2)=FS(i,j,k2)+ &
1065!^ & dife* &
1066!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1067!^ & dTde(i,j ,k1))+ &
1068!^ & cff2*(cff2*dTdz(i,j,k2)- &
1069!^ & dTde(i,j+1,k2))+ &
1070!^ & cff3*(cff3*dTdz(i,j,k2)- &
1071!^ & dTde(i,j ,k2))+ &
1072!^ & cff4*(cff4*dTdz(i,j,k2)- &
1073!^ & dTde(i,j+1,k1)))
1074!^
1075 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
1076 & dife* &
1077 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1078 & dtde(i,j ,k1))+ &
1079 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1080 & dtde(i,j+1,k2))+ &
1081 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1082 & dtde(i,j ,k2))+ &
1083 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1084 & dtde(i,j+1,k1))+ &
1085 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1086 & cff1*tl_dtdz(i,j,k2)- &
1087 & tl_dtde(i,j ,k1))+ &
1088 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1089 & cff2*tl_dtdz(i,j,k2)- &
1090 & tl_dtde(i,j+1,k2))+ &
1091 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1092 & cff3*tl_dtdz(i,j,k2)- &
1093 & tl_dtde(i,j ,k2))+ &
1094 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1095 & cff4*tl_dtdz(i,j,k2)- &
1096 & tl_dtde(i,j+1,k1)))
1097 END DO
1098 END DO
1099 END IF
1100!
1101! Time-step biharmonic, geopotential diffusion term (m Tunits).
1102!
1103 DO j=jstr,jend
1104 DO i=istr,iend
1105!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1106!^ & (FX(i+1,j )-FX(i,j)+ &
1107!^ & FE(i ,j+1)-FE(i,j))+ &
1108!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1109!^
1110 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1111 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1112 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1113 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1114!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1115!^
1116 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1117#ifdef DIAGNOSTICS_TS
1118!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1119#endif
1120 END DO
1121 END DO
1122 END IF
1123 END DO k_loop2
1124 END DO t_loop
1125!
1126 RETURN
1127 END SUBROUTINE tl_t3dmix4_geo_tile
1128
1129 END MODULE tl_t3dmix4_mod
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::istvar
integer, dimension(:), allocatable istvar
Definition mod_ncparam.F:520
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_lbc
type(t_lbc), dimension(:,:,:), allocatable tl_lbc
Definition mod_param.F:379
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
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::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
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_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
tl_t3dmix4_mod
Definition tl_t3dmix4_geo.h:1
tl_t3dmix4_mod::tl_t3dmix4_geo_tile
subroutine tl_t3dmix4_geo_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, tl_hz, z_r, tl_z_r, diff3d_u, diff3d_v, diff3d_r, diff4, tclm, t, tl_t)
Definition tl_t3dmix4_geo.h:129
tl_t3dmix4_mod::tl_t3dmix4
subroutine, public tl_t3dmix4(ng, tile)
Definition tl_t3dmix4_geo.h:26
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