ROMS
Loading...
Searching...
No Matches
rp_t3dmix4_iso.h
Go to the documentation of this file.
1 MODULE rp_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 representers tangent linear horizontal !
11! biharmonic 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 rp_t3dmix4
21!
22 CONTAINS
23!
24!***********************************************************************
25 SUBROUTINE rp_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, irpm, 29, __line__, myfile)
53#endif
54 CALL rp_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, irpm, 29, __line__, myfile)
96#endif
97!
98 RETURN
99 END SUBROUTINE rp_t3dmix4
100!
101!***********************************************************************
102 SUBROUTINE rp_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# ifdef TL_IOMS
427 & cff2
428# endif
429 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
430 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
431 & small))* &
432 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
433# ifdef TL_IOMS
434 & (0.5_r8-sign(0.5_r8, &
435 & pden(i,j,k)-pden(i,j,k+1)-small))* &
436 & small
437# endif
438 cff4=max(cff2,cff3)
439 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
440 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
441 cff=-1.0_r8/cff4
442 tl_cff=cff*cff*tl_cff4+ &
443# ifdef TL_IOMS
444 & 2.0_r8*cff
445# endif
446#elif defined TS_MIX_MIN_STRAT
447 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
448 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
449 tl_cff1=(0.5_r8+sign(0.5_r8, &
450 & pden(i,j,k)-pden(i,j,k+1)- &
451 & strat_min*(z_r(i,j,k+1)- &
452 & z_r(i,j,k ))))* &
453 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
454 & (0.5_r8-sign(0.5_r8, &
455 & pden(i,j,k)-pden(i,j,k+1)- &
456 & strat_min*(z_r(i,j,k+1)- &
457 & z_r(i,j,k ))))* &
458 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
459 cff=-1.0_r8/cff1
460 tl_cff=cff*cff*tl_cff1+ &
461# ifdef TL_IOMS
462 & 2.0_r8*cff
463# endif
464#else
465 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
466 tl_cff1=(0.5_r8+sign(0.5_r8, &
467 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
468 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
469# ifdef TL_IOMS
470 & (0.5_r8-sign(0.5_r8, &
471 & pden(i,j,k)-pden(i,j,k+1)-eps))*eps
472# endif
473 cff=-1.0_r8/cff1
474 tl_cff=cff*cff*tl_cff1+ &
475# ifdef TL_IOMS
476 & 2.0_r8*cff
477# endif
478#endif
479#if defined TS_MIX_STABILITY
480 dtdr(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
481 & t(i,j,k ,nrhs,itrc))+ &
482 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
483 & t(i,j,k ,nstp,itrc)))
484 tl_dtdr(i,j,k2)=tl_cff* &
485 & (0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
486 & t(i,j,k ,nrhs,itrc))+ &
487 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
488 & t(i,j,k ,nstp,itrc)))+ &
489 & cff*
490 & (0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
491 & tl_t(i,j,k ,nrhs,itrc))+ &
492 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
493 & tl_t(i,j,k ,nstp,itrc)))- &
494# ifdef TL_IOMS
495 & dtdr(i,j,k2)
496# endif
497#elif defined TS_MIX_CLIMA
498 IF (ltracerclm(itrc,ng)) THEN
499 dtdr(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
500 & tclm(i,j,k+1,itrc))- &
501 & (t(i,j,k ,nrhs,itrc)- &
502 & tclm(i,j,k ,itrc)))
503 tl_dtdr(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
504 & tclm(i,j,k+1,itrc))- &
505 & (t(i,j,k ,nrhs,itrc)- &
506 & tclm(i,j,k ,itrc)))+ &
507 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
508 & tl_t(i,j,k ,nrhs,itrc))- &
509# ifdef TL_IOMS
510 & dtdr(i,j,k2)
511# endif
512 ELSE
513 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
514 & t(i,j,k ,nrhs,itrc))
515 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
516 & t(i,j,k ,nrhs,itrc))+ &
517 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
518 & tl_t(i,j,k ,nrhs,itrc))- &
519# ifdef TL_IOMS
520 & dtdr(i,j,k2)
521# endif
522 END IF
523#else
524 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
525 & t(i,j,k ,nrhs,itrc))
526 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
527 & t(i,j,k ,nrhs,itrc))+ &
528 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
529 & tl_t(i,j,k ,nrhs,itrc))- &
530# ifdef TL_IOMS
531 & dtdr(i,j,k2)
532# endif
533#endif
534 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
535 & z_r(i,j,k ))
536 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
537 & z_r(i,j,k ))+ &
538 & cff*(tl_z_r(i,j,k+1)- &
539 & tl_z_r(i,j,k ))- &
540#ifdef TL_IOMS
541 & fs(i,j,k2)
542#endif
543 END DO
544 END DO
545 END IF
546 IF (k.gt.0) THEN
547 DO j=jmin,jmax
548 DO i=imin,imax+1
549#ifdef DIFF_3DCOEF
550# ifdef TS_U3ADV_SPLIT
551 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
552# else
553 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
554 & on_u(i,j)
555# endif
556#else
557 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
558 & on_u(i,j)
559#endif
560 fx(i,j)=cff* &
561 & (hz(i,j,k)+hz(i-1,j,k))* &
562 & (dtdx(i,j,k1)- &
563 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
564 & (dtdr(i-1,j,k1)+ &
565 & dtdr(i ,j,k2))+ &
566 & min(drdx(i,j,k1),0.0_r8)* &
567 & (dtdr(i-1,j,k2)+ &
568 & dtdr(i ,j,k1))))
569 tl_fx(i,j)=cff* &
570 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
571 & (dtdx(i,j,k1)- &
572 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
573 & (dtdr(i-1,j,k1)+ &
574 & dtdr(i ,j,k2))+ &
575 & min(drdx(i,j,k1),0.0_r8)* &
576 & (dtdr(i-1,j,k2)+ &
577 & dtdr(i ,j,k1))))+ &
578 & (hz(i,j,k)+hz(i-1,j,k))* &
579 & (tl_dtdx(i,j,k1)- &
580 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
581 & (tl_dtdr(i-1,j,k1)+ &
582 & tl_dtdr(i ,j,k2))+ &
583 & min(drdx(i,j,k1),0.0_r8)* &
584 & (tl_dtdr(i-1,j,k2)+ &
585 & tl_dtdr(i ,j,k1)))- &
586 & 0.5_r8*((0.5_r8+ &
587 & sign(0.5_r8, drdx(i,j,k1)))* &
588 & tl_drdx(i,j,k1)* &
589 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
590 & (0.5_r8+ &
591 & sign(0.5_r8,-drdx(i,j,k1)))* &
592 & tl_drdx(i,j,k1)* &
593 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))- &
594#ifdef TL_IOMS
595 & cff* &
596 & (hz(i,j,k)+hz(i-1,j,k))* &
597 & (dtdx(i,j,k1)- &
598 & (max(drdx(i,j,k1),0.0_r8)* &
599 & (dtdr(i-1,j,k1)+ &
600 & dtdr(i ,j,k2))+ &
601 & min(drdx(i,j,k1),0.0_r8)* &
602 & (dtdr(i-1,j,k2)+ &
603 & dtdr(i ,j,k1))))
604#endif
605 END DO
606 END DO
607 DO j=jmin,jmax+1
608 DO i=imin,imax
609#ifdef DIFF_3DCOEF
610# ifdef TS_U3ADV_SPLIT
611 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
612# else
613 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
614 & om_v(i,j)
615# endif
616#else
617 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
618 & om_v(i,j)
619#endif
620 fe(i,j)=cff* &
621 & (hz(i,j,k)+hz(i,j-1,k))* &
622 & (dtde(i,j,k1)- &
623 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
624 & (dtdr(i,j-1,k1)+ &
625 & dtdr(i,j ,k2))+ &
626 & min(drde(i,j,k1),0.0_r8)* &
627 & (dtdr(i,j-1,k2)+ &
628 & dtdr(i,j ,k1))))
629 tl_fe(i,j)=cff* &
630 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
631 & (dtde(i,j,k1)- &
632 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
633 & (dtdr(i,j-1,k1)+ &
634 & dtdr(i,j ,k2))+ &
635 & min(drde(i,j,k1),0.0_r8)* &
636 & (dtdr(i,j-1,k2)+ &
637 & dtdr(i,j ,k1))))+ &
638 & (hz(i,j,k)+hz(i,j-1,k))* &
639 & (tl_dtde(i,j,k1)- &
640 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
641 & (tl_dtdr(i,j-1,k1)+ &
642 & tl_dtdr(i,j ,k2))+ &
643 & min(drde(i,j,k1),0.0_r8)* &
644 & (tl_dtdr(i,j-1,k2)+ &
645 & tl_dtdr(i,j ,k1)))- &
646 & 0.5_r8*((0.5_r8+ &
647 & sign(0.5_r8, drde(i,j,k1)))* &
648 & tl_drde(i,j,k1)* &
649 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
650 & (0.5_r8+ &
651 & sign(0.5_r8,-drde(i,j,k1)))* &
652 & tl_drde(i,j,k1)* &
653 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))- &
654#ifdef TL_IOMS
655 & cff* &
656 & (hz(i,j,k)+hz(i,j-1,k))* &
657 & (dtde(i,j,k1)- &
658 & (max(drde(i,j,k1),0.0_r8)* &
659 & (dtdr(i,j-1,k1)+ &
660 & dtdr(i,j ,k2))+ &
661 & min(drde(i,j,k1),0.0_r8)* &
662 & (dtdr(i,j-1,k2)+ &
663 & dtdr(i,j ,k1))))
664#endif
665 END DO
666 END DO
667 IF (k.lt.n(ng)) THEN
668 DO j=jmin,jmax
669 DO i=imin,imax
670#ifdef DIFF_3DCOEF
671# ifdef TS_U3ADV_SPLIT
672 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
673 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
674 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
675 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
676# else
677 difx=0.5_r8*diff3d_r(i,j,k)
678 dife=difx
679# endif
680#else
681 difx=0.5_r8*diff4(i,j,itrc)
682 dife=difx
683#endif
684 cff1=max(drdx(i ,j,k1),0.0_r8)
685 cff2=max(drdx(i+1,j,k2),0.0_r8)
686 cff3=min(drdx(i ,j,k2),0.0_r8)
687 cff4=min(drdx(i+1,j,k1),0.0_r8)
688 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
689 & tl_drdx(i ,j,k1)
690 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
691 & tl_drdx(i+1,j,k2)
692 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
693 & tl_drdx(i ,j,k2)
694 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
695 & tl_drdx(i+1,j,k1)
696 cff=difx* &
697 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
698 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
699 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
700 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
701 tl_cff=difx* &
702 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
703 & dtdx(i ,j,k1))+ &
704 & tl_cff2*(cff2*dtdr(i,j,k2)- &
705 & dtdx(i+1,j,k2))+ &
706 & tl_cff3*(cff3*dtdr(i,j,k2)- &
707 & dtdx(i ,j,k2))+ &
708 & tl_cff4*(cff4*dtdr(i,j,k2)- &
709 & dtdx(i+1,j,k1))+ &
710 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
711 & cff1*tl_dtdr(i,j,k2)- &
712 & tl_dtdx(i ,j,k1))+ &
713 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
714 & cff2*tl_dtdr(i,j,k2)- &
715 & tl_dtdx(i+1,j,k2))+ &
716 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
717 & cff3*tl_dtdr(i,j,k2)- &
718 & tl_dtdx(i ,j,k2))+ &
719 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
720 & cff4*tl_dtdr(i,j,k2)- &
721 & tl_dtdx(i+1,j,k1)))- &
722#ifdef TL_IOMS
723 & difx* &
724 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
725 & dtdx(i ,j,k1))- &
726 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
727 & dtdx(i+1,j,k2))- &
728 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
729 & dtdx(i ,j,k2))- &
730 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
731 & dtdx(i+1,j,k1)))
732#endif
733!
734 cff1=max(drde(i,j ,k1),0.0_r8)
735 cff2=max(drde(i,j+1,k2),0.0_r8)
736 cff3=min(drde(i,j ,k2),0.0_r8)
737 cff4=min(drde(i,j+1,k1),0.0_r8)
738 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
739 & tl_drde(i,j ,k1)
740 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
741 & tl_drde(i,j+1,k2)
742 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
743 & tl_drde(i,j ,k2)
744 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
745 & tl_drde(i,j+1,k1)
746 cff=cff+ &
747 & dife* &
748 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
749 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
750 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
751 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
752 tl_cff=tl_cff+ &
753 & dife* &
754 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
755 & dtde(i,j ,k1))+ &
756 & tl_cff2*(cff2*dtdr(i,j,k2)- &
757 & dtde(i,j+1,k2))+ &
758 & tl_cff3*(cff3*dtdr(i,j,k2)- &
759 & dtde(i,j ,k2))+ &
760 & tl_cff4*(cff4*dtdr(i,j,k2)- &
761 & dtde(i,j+1,k1))+ &
762 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
763 & cff1*tl_dtdr(i,j,k2)- &
764 & tl_dtde(i,j ,k1))+ &
765 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
766 & cff2*tl_dtdr(i,j,k2)- &
767 & tl_dtde(i,j+1,k2))+ &
768 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
769 & cff3*tl_dtdr(i,j,k2)- &
770 & tl_dtde(i,j ,k2))+ &
771 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
772 & cff4*tl_dtdr(i,j,k2)- &
773 & tl_dtde(i,j+1,k1)))- &
774#ifdef TL_IOMS
775 & dife* &
776 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
777 & dtde(i,j ,k1))- &
778 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
779 & dtde(i,j+1,k2))- &
780 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
781 & dtde(i,j ,k2))- &
782 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)-
783 & dtde(i,j+1,k1)))
784#endif
785!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
786!^ ! compute after TLM
787!^
788 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
789 & cff*tl_fs(i,j,k2)- &
790#ifdef TL_IOMS
791 & cff*fs(i,j,k2)
792#endif
793 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
794 END DO
795 END DO
796 END IF
797!
798! Compute first harmonic operator, without mixing coefficient.
799! Multiply by the metrics of the second harmonic operator. Save
800! into work array "LapT".
801!
802 DO j=jmin,jmax
803 DO i=imin,imax
804 cff=pm(i,j)*pn(i,j)
805 cff1=1.0_r8/hz(i,j,k)
806 tl_cff1=-cff1*cff1*tl_hz(i,j,k)+ &
807#ifdef TL_IOMS
808 & 2.0_r8*cff1
809#endif
810 lapt(i,j,k)=cff1*(cff* &
811 & (fx(i+1,j)-fx(i,j)+ &
812 & fe(i,j+1)-fe(i,j))+ &
813 & (fs(i,j,k2)-fs(i,j,k1)))
814 tl_lapt(i,j,k)=tl_cff1*(cff* &
815 & (fx(i+1,j)-fx(i,j)+ &
816 & fe(i,j+1)-fe(i,j))+ &
817 & (fs(i,j,k2)-fs(i,j,k1)))+ &
818 & cff1*(cff* &
819 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
820 & tl_fe(i,j+1)-tl_fe(i,j))+ &
821 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))- &
822#ifdef TL_IOMS
823 & lapt(i,j,k)
824#endif
825 END DO
826 END DO
827 END IF
828 END DO k_loop1
829!
830! Apply boundary conditions (except periodic; closed or gradient)
831! to the first harmonic operator.
832!
833 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
834 IF (domain(ng)%Western_Edge(tile)) THEN
835 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
836 DO k=1,n(ng)
837 DO j=jmin,jmax
838 lapt(istr-1,j,k)=0.0_r8
839 tl_lapt(istr-1,j,k)=0.0_r8
840 END DO
841 END DO
842 ELSE
843 DO k=1,n(ng)
844 DO j=jmin,jmax
845 lapt(istr-1,j,k)=lapt(istr,j,k)
846 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
847 END DO
848 END DO
849 END IF
850 END IF
851 END IF
852!
853 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
854 IF (domain(ng)%Eastern_Edge(tile)) THEN
855 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
856 DO k=1,n(ng)
857 DO j=jmin,jmax
858 lapt(iend+1,j,k)=0.0_r8
859 tl_lapt(iend+1,j,k)=0.0_r8
860 END DO
861 END DO
862 ELSE
863 DO k=1,n(ng)
864 DO j=jmin,jmax
865 lapt(iend+1,j,k)=lapt(iend,j,k)
866 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
867 END DO
868 END DO
869 END IF
870 END IF
871 END IF
872!
873 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
874 IF (domain(ng)%Southern_Edge(tile)) THEN
875 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
876 DO k=1,n(ng)
877 DO i=imin,imax
878 lapt(i,jstr-1,k)=0.0_r8
879 tl_lapt(i,jstr-1,k)=0.0_r8
880 END DO
881 END DO
882 ELSE
883 DO k=1,n(ng)
884 DO i=imin,imax
885 lapt(i,jstr-1,k)=lapt(i,jstr,k)
886 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
887 END DO
888 END DO
889 END IF
890 END IF
891 END IF
892!
893 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
894 IF (domain(ng)%Northern_Edge(tile)) THEN
895 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
896 DO k=1,n(ng)
897 DO i=imin,imax
898 lapt(i,jend+1,k)=0.0_r8
899 tl_lapt(i,jend+1,k)=0.0_r8
900 END DO
901 END DO
902 ELSE
903 DO k=1,n(ng)
904 DO i=imin,imax
905 lapt(i,jend+1,k)=lapt(i,jend,k)
906 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
907 END DO
908 END DO
909 END IF
910 END IF
911 END IF
912!
913 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
914 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
915 IF (domain(ng)%SouthWest_Corner(tile)) THEN
916 DO k=1,n(ng)
917 lapt(istr-1,jstr-1,k)=0.5_r8* &
918 & (lapt(istr ,jstr-1,k)+ &
919 & lapt(istr-1,jstr ,k))
920 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
921 & (tl_lapt(istr ,jstr-1,k)+ &
922 tl_lapt(istr-1,jstr ,k))
923 END DO
924 END IF
925 END IF
926
927 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
928 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
929 IF (domain(ng)%SouthEast_Corner(tile)) THEN
930 DO k=1,n(ng)
931 lapt(iend+1,jstr-1,k)=0.5_r8* &
932 & (lapt(iend ,jstr-1,k)+ &
933 & lapt(iend+1,jstr ,k))
934 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
935 & (tl_lapt(iend ,jstr-1,k)+ &
936 & tl_lapt(iend+1,jstr ,k))
937 END DO
938 END IF
939 END IF
940
941 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
942 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
943 IF (domain(ng)%NorthWest_Corner(tile)) THEN
944 DO k=1,n(ng)
945 lapt(istr-1,jend+1,k)=0.5_r8* &
946 & (lapt(istr ,jend+1,k)+ &
947 & lapt(istr-1,jend ,k))
948 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
949 & (tl_lapt(istr ,jend+1,k)+ &
950 & tl_lapt(istr-1,jend ,k))
951 END DO
952 END IF
953 END IF
954
955 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
956 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
957 IF (domain(ng)%NorthEast_Corner(tile)) THEN
958 DO k=1,n(ng)
959 lapt(iend+1,jend+1,k)=0.5_r8* &
960 & (lapt(iend ,jend+1,k)+ &
961 & lapt(iend+1,jend ,k))
962 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
963 & (tl_lapt(iend ,jend+1,k)+ &
964 & tl_lapt(iend+1,jend ,k))
965 END DO
966 END IF
967 END IF
968!
969! Compute horizontal and density gradients associated with the
970! second rotated harmonic operator.
971!
972 k2=1
973 k_loop2: DO k=0,n(ng)
974 k1=k2
975 k2=3-k1
976 IF (k.lt.n(ng)) THEN
977 DO j=jstr,jend
978 DO i=istr,iend+1
979 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
980#ifdef MASKING
981 cff=cff*umask(i,j)
982#endif
983#ifdef WET_DRY_NOT_YET
984 cff=cff*umask_wet(i,j)
985#endif
986 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
987 & pden(i-1,j,k+1))
988 tl_drdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
989 & tl_pden(i-1,j,k+1))
990 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
991 & lapt(i-1,j,k+1))
992 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
993 & tl_lapt(i-1,j,k+1))
994 END DO
995 END DO
996 DO j=jstr,jend+1
997 DO i=istr,iend
998 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
999#ifdef MASKING
1000 cff=cff*vmask(i,j)
1001#endif
1002#ifdef WET_DRY_NOT_YET
1003 cff=cff*vmask_wet(i,j)
1004#endif
1005 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
1006 & pden(i,j-1,k+1))
1007 tl_drde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
1008 & tl_pden(i,j-1,k+1))
1009 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
1010 & lapt(i,j-1,k+1))
1011 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
1012 & tl_lapt(i,j-1,k+1))
1013 END DO
1014 END DO
1015 END IF
1016 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1017 DO j=jstr-1,jend+1
1018 DO i=istr-1,iend+1
1019 dtdr(i,j,k2)=0.0_r8
1020 tl_dtdr(i,j,k2)=0.0_r8
1021 fs(i,j,k2)=0.0_r8
1022 tl_fs(i,j,k2)=0.0_r8
1023 END DO
1024 END DO
1025 ELSE
1026 DO j=jstr-1,jend+1
1027 DO i=istr-1,iend+1
1028#if defined TS_MIX_MAX_SLOPE
1029 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
1030 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
1031 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
1032 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
1033 IF (cff1.ne.0.0_r8) THEN
1034 tl_cff1=(drdx(i ,j,k2)*tl_drdx(i ,j,k2)+ &
1035 & drdx(i+1,j,k2)*tl_drdx(i+1,j,k2)+ &
1036 & drdx(i ,j,k1)*tl_drdx(i ,j,k1)+ &
1037 & drdx(i+1,j,k1)*tl_drdx(i+1,j,k1)+ &
1038 & drde(i,j ,k2)*tl_drde(i,j ,k2)+ &
1039 & drde(i,j+1,k2)*tl_drde(i,j+1,k2)+ &
1040 & drde(i,j ,k1)*tl_drde(i,j ,k1)+ &
1041 & drde(i,j+1,k1)*tl_drde(i,j+1,k1))/cff1
1042 ELSE
1043 tl_cff1=0.0_r8
1044 END IF
1045 cff2=0.25_r8*slope_max* &
1046 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
1047 tl_cff2=0.25_r8*slope_max* &
1048 & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
1049 & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)- &
1050# ifdef TL_IOMS
1051 & cff2
1052# endif
1053 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
1054 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
1055 & small))* &
1056 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1057# ifdef TL_IOMS
1058 & (0.5_r8-sign(0.5_r8, &
1059 & pden(i,j,k)-pden(i,j,k+1)-small))* &
1060 & small
1061# endif
1062 cff4=max(cff2,cff3)
1063 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
1064 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
1065 cff=-1.0_r8/cff4
1066 tl_cff=cff*cff*tl_cff4+ &
1067# ifdef TL_IOMS
1068 & 2.0_r8*cff
1069# endif
1070#elif defined TS_MIX_MIN_STRAT
1071 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
1072 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
1073 tl_cff1=(0.5_r8+sign(0.5_r8, &
1074 & pden(i,j,k)-pden(i,j,k+1)- &
1075 & strat_min*(z_r(i,j,k+1)- &
1076 & z_r(i,j,k ))))* &
1077 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1078 & (0.5_r8-sign(0.5_r8, &
1079 & pden(i,j,k)-pden(i,j,k+1)- &
1080 & strat_min*(z_r(i,j,k+1)- &
1081 & z_r(i,j,k ))))* &
1082 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
1083 cff=-1.0_r8/cff1
1084 tl_cff=cff*cff*tl_cff1+ &
1085# ifdef TL_IOMS
1086 & 2.0_r8*cff
1087# endif
1088#else
1089 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
1090 tl_cff1=(0.5_r8+sign(0.5_r8, &
1091 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
1092 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1093# ifdef TL_IOMS
1094 & (0.5_r8-sign(0.5_r8, &
1095 & pden(i,j,k)-pden(i,j,k+1)-eps))*eps
1096# endif
1097 cff=-1.0_r8/cff1
1098 tl_cff=cff*cff*tl_cff1+ &
1099# ifdef TL_IOMS
1100 & 2.0_r8*cff
1101# endif
1102#endif
1103 dtdr(i,j,k2)=cff*(lapt(i,j,k+1)- &
1104 & lapt(i,j,k ))
1105 tl_dtdr(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
1106 & lapt(i,j,k ))+ &
1107 & cff*(tl_lapt(i,j,k+1)- &
1108 & tl_lapt(i,j,k ))- &
1109#ifdef TL_IOMS
1110 & dtdr(i,j,k2)
1111#endif
1112 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
1113 & z_r(i,j,k ))
1114 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
1115 & z_r(i,j,k ))+ &
1116 & cff*(tl_z_r(i,j,k+1)- &
1117 & tl_z_r(i,j,k ))- &
1118#ifdef TL_IOMS
1119 & fs(i,j,k2)
1120#endif
1121 END DO
1122 END DO
1123 END IF
1124!
1125! Compute components of the rotated tracer flux (T m4/s) along
1126! isopycnic surfaces.
1127!
1128 IF (k.gt.0) THEN
1129 DO j=jstr,jend
1130 DO i=istr,iend+1
1131#ifdef DIFF_3DCOEF
1132# ifdef TS_U3ADV_SPLIT
1133 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1134# else
1135 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1136 & on_u(i,j)
1137# endif
1138#else
1139 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1140 & on_u(i,j)
1141#endif
1142!^ FX(i,j)=cff* &
1143!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
1144!^ & (dTdx(i,j,k1)- &
1145!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
1146!^ & (dTdr(i-1,j,k1)+ &
1147!^ & dTdr(i ,j,k2))+ &
1148!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
1149!^ & (dTdr(i-1,j,k2)+ &
1150!^ & dTdr(i ,j,k1))))
1151!^
1152 tl_fx(i,j)=cff* &
1153 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
1154 & (dtdx(i,j,k1)- &
1155 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1156 & (dtdr(i-1,j,k1)+ &
1157 & dtdr(i ,j,k2))+ &
1158 & min(drdx(i,j,k1),0.0_r8)* &
1159 & (dtdr(i-1,j,k2)+ &
1160 & dtdr(i ,j,k1))))+ &
1161 & (hz(i,j,k)+hz(i-1,j,k))* &
1162 & (tl_dtdx(i,j,k1)- &
1163 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1164 & (tl_dtdr(i-1,j,k1)+ &
1165 & tl_dtdr(i ,j,k2))+ &
1166 & min(drdx(i,j,k1),0.0_r8)* &
1167 & (tl_dtdr(i-1,j,k2)+ &
1168 & tl_dtdr(i ,j,k1)))- &
1169 & 0.5_r8*((0.5_r8+ &
1170 & sign(0.5_r8, drdx(i,j,k1)))* &
1171 & tl_drdx(i,j,k1)* &
1172 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
1173 & (0.5_r8+ &
1174 & sign(0.5_r8,-drdx(i,j,k1)))* &
1175 & tl_drdx(i,j,k1)* &
1176 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))- &
1177#ifdef TL_IOMS
1178 cff* &
1179 & (hz(i,j,k)+hz(i-1,j,k))* &
1180 & (dtdx(i,j,k1)- &
1181 & (max(drdx(i,j,k1),0.0_r8)* &
1182 & (dtdr(i-1,j,k1)+ &
1183 & dtdr(i ,j,k2))+ &
1184 & min(drdx(i,j,k1),0.0_r8)* &
1185 & (dtdr(i-1,j,k2)+ &
1186 & dtdr(i ,j,k1))))
1187#endif
1188 END DO
1189 END DO
1190 DO j=jstr,jend+1
1191 DO i=istr,iend
1192#ifdef DIFF_3DCOEF
1193# ifdef TS_U3ADV_SPLIT
1194 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1195# else
1196 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1197 & om_v(i,j)
1198# endif
1199#else
1200 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1201 & om_v(i,j)
1202#endif
1203!^ FE(i,j)=cff* &
1204!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1205!^ & (dTde(i,j,k1)- &
1206!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
1207!^ & (dTdr(i,j-1,k1)+ &
1208!^ & dTdr(i,j ,k2))+ &
1209!^ & MIN(dRde(i,j,k1),0.0_r8)* &
1210!^ & (dTdr(i,j-1,k2)+ &
1211!^ & dTdr(i,j ,k1))))
1212!^
1213 tl_fe(i,j)=cff* &
1214 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
1215 & (dtde(i,j,k1)- &
1216 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1217 & (dtdr(i,j-1,k1)+ &
1218 & dtdr(i,j ,k2))+ &
1219 & min(drde(i,j,k1),0.0_r8)* &
1220 & (dtdr(i,j-1,k2)+ &
1221 & dtdr(i,j ,k1))))+ &
1222 & (hz(i,j,k)+hz(i,j-1,k))* &
1223 & (tl_dtde(i,j,k1)- &
1224 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1225 & (tl_dtdr(i,j-1,k1)+ &
1226 & tl_dtdr(i,j ,k2))+ &
1227 & min(drde(i,j,k1),0.0_r8)* &
1228 & (tl_dtdr(i,j-1,k2)+ &
1229 & tl_dtdr(i,j ,k1)))- &
1230 & 0.5_r8*((0.5_r8+ &
1231 & sign(0.5_r8, drde(i,j,k1)))* &
1232 & tl_drde(i,j,k1)* &
1233 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
1234 & (0.5_r8+ &
1235 & sign(0.5_r8,-drde(i,j,k1)))* &
1236 & tl_drde(i,j,k1)* &
1237 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))- &
1238#ifdef TL_IOMS
1239 & cff* &
1240 & (hz(i,j,k)+hz(i,j-1,k))* &
1241 & (dtde(i,j,k1)- &
1242 & (max(drde(i,j,k1),0.0_r8)* &
1243 & (dtdr(i,j-1,k1)+ &
1244 & dtdr(i,j ,k2))+ &
1245 & min(drde(i,j,k1),0.0_r8)* &
1246 & (dtdr(i,j-1,k2)+ &
1247 & dtdr(i,j ,k1))))
1248#endif
1249 END DO
1250 END DO
1251 IF (k.lt.n(ng)) THEN
1252 DO j=jstr,jend
1253 DO i=istr,iend
1254#ifdef DIFF_3DCOEF
1255# ifdef TS_U3ADV_SPLIT
1256 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1257 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1258 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1259 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1260# else
1261 difx=0.5_r8*diff3d_r(i,j,k)
1262 dife=difx
1263# endif
1264#else
1265 difx=0.5_r8*diff4(i,j,itrc)
1266 dife=difx
1267#endif
1268 cff1=max(drdx(i ,j,k1),0.0_r8)
1269 cff2=max(drdx(i+1,j,k2),0.0_r8)
1270 cff3=min(drdx(i ,j,k2),0.0_r8)
1271 cff4=min(drdx(i+1,j,k1),0.0_r8)
1272 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
1273 & tl_drdx(i ,j,k1)
1274 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
1275 & tl_drdx(i+1,j,k2)
1276 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
1277 & tl_drdx(i ,j,k2)
1278 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
1279 & tl_drdx(i+1,j,k1)
1280 cff=difx* &
1281 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
1282 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
1283 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
1284 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
1285 tl_cff=difx* &
1286 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
1287 & dtdx(i ,j,k1))+ &
1288 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1289 & dtdx(i+1,j,k2))+ &
1290 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1291 & dtdx(i ,j,k2))+ &
1292 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1293 & dtdx(i+1,j,k1))+ &
1294 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1295 & cff1*tl_dtdr(i,j,k2)- &
1296 & tl_dtdx(i ,j,k1))+ &
1297 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1298 & cff2*tl_dtdr(i,j,k2)- &
1299 & tl_dtdx(i+1,j,k2))+ &
1300 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1301 & cff3*tl_dtdr(i,j,k2)- &
1302 & tl_dtdx(i ,j,k2))+ &
1303 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1304 & cff4*tl_dtdr(i,j,k2)- &
1305 & tl_dtdx(i+1,j,k1)))- &
1306#ifdef TL_IOMS
1307 & difx* &
1308 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
1309 & dtdx(i ,j,k1))- &
1310 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
1311 & dtdx(i+1,j,k2))- &
1312 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
1313 & dtdx(i ,j,k2))- &
1314 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
1315 & dtdx(i+1,j,k1)))
1316#endif
1317!
1318 cff1=max(drde(i,j ,k1),0.0_r8)
1319 cff2=max(drde(i,j+1,k2),0.0_r8)
1320 cff3=min(drde(i,j ,k2),0.0_r8)
1321 cff4=min(drde(i,j+1,k1),0.0_r8)
1322 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
1323 & tl_drde(i,j ,k1)
1324 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
1325 & tl_drde(i,j+1,k2)
1326 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
1327 & tl_drde(i,j ,k2)
1328 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
1329 & tl_drde(i,j+1,k1)
1330 cff=cff+ &
1331 & dife* &
1332 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
1333 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
1334 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
1335 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
1336 tl_cff=tl_cff+ &
1337 & dife* &
1338 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
1339 & dtde(i,j ,k1))+ &
1340 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1341 & dtde(i,j+1,k2))+ &
1342 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1343 & dtde(i,j ,k2))+ &
1344 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1345 & dtde(i,j+1,k1))+ &
1346 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1347 & cff1*tl_dtdr(i,j,k2)- &
1348 & tl_dtde(i,j ,k1))+ &
1349 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1350 & cff2*tl_dtdr(i,j,k2)- &
1351 & tl_dtde(i,j+1,k2))+ &
1352 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1353 & cff3*tl_dtdr(i,j,k2)- &
1354 & tl_dtde(i,j ,k2))+ &
1355 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1356 & cff4*tl_dtdr(i,j,k2)- &
1357 & tl_dtde(i,j+1,k1)))- &
1358#ifdef TL_IOMS
1359 & dife* &
1360 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
1361 & dtde(i,j ,k1))- &
1362 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
1363 & dtde(i,j+1,k2))- &
1364 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
1365 & dtde(i,j ,k2))- &
1366 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
1367 & dtde(i,j+1,k1)))
1368#endif
1369!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
1370!^ ! compute after TLM
1371!^
1372 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
1373 & cff*tl_fs(i,j,k2)
1374 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
1375#ifdef TL_IOMS
1376 tl_fs(i,j,k2)=tl_fs(i,j,k2)-fs(i,j,k2)
1377#endif
1378 END DO
1379 END DO
1380 END IF
1381!
1382! Time-step biharmonic, isopycnal diffusion term (m Tunits).
1383!
1384 DO j=jstr,jend
1385 DO i=istr,iend
1386!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1387!^ & (FX(i+1,j )-FX(i,j)+ &
1388!^ & FE(i ,j+1)-FE(i,j))+ &
1389!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1390!^
1391 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1392 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1393 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1394 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1395!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1396!^
1397 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1398#ifdef DIAGNOSTICS_TS
1399!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1400#endif
1401 END DO
1402 END DO
1403 END IF
1404 END DO k_loop2
1405 END DO t_loop
1406!
1407 RETURN
1408 END SUBROUTINE rp_t3dmix4_iso_tile
1409
1410 END MODULE rp_t3dmix4_mod
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
integer, dimension(:), allocatable istvar
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer, parameter irpm
Definition mod_param.F:664
type(t_lbc), dimension(:,:,:), allocatable tl_lbc
Definition mod_param.F:379
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer, dimension(:), allocatable lm
Definition mod_param.F:455
integer, dimension(:), allocatable mm
Definition mod_param.F:456
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ewperiodic
integer, parameter iwest
logical, dimension(:), allocatable nsperiodic
logical, dimension(:,:), allocatable compositegrid
integer, parameter isouth
integer, parameter ieast
logical, dimension(:,:), allocatable ltracerclm
integer, parameter inorth
integer, dimension(:), allocatable nrhs
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable nstp
subroutine rp_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)
subroutine, public rp_t3dmix4(ng, tile)
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3