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