ROMS
Loading...
Searching...
No Matches
rp_t3dmix4_geo.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 geopotential surfaces. !
12! !
13! BASIC STATE variables needed: diff4, Hz, t, 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, 28, __line__, myfile)
53#endif
54 CALL rp_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, irpm, 28, __line__, myfile)
94#endif
95!
96 RETURN
97 END SUBROUTINE rp_t3dmix4
98!
99!***********************************************************************
100 SUBROUTINE rp_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#ifdef TL_IOMS
397 & 2.0_r8*cff
398#endif
399#if defined TS_MIX_STABILITY
400 dtdz(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
401 & t(i,j,k ,nrhs,itrc))+ &
402 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
403 & t(i,j,k ,nstp,itrc)))
404 tl_dtdz(i,j,k2)=tl_cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
405 & t(i,j,k ,nrhs,itrc))+ &
406 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
407 & t(i,j,k ,nstp,itrc)))+&
408 & cff*(0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
409 & tl_t(i,j,k ,nrhs,itrc))+ &
410 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
411 & tl_t(i,j,k ,nstp,itrc)))-&
412# ifdef TL_IOMS
413 & dtdz(i,j,k2)
414# endif
415#elif defined TS_MIX_CLIMA
416 IF (ltracerclm(itrc,ng)) THEN
417 dtdz(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
418 & tclm(i,j,k+1,itrc))- &
419 & (t(i,j,k ,nrhs,itrc)- &
420 & tclm(i,j,k ,itrc)))
421 tl_dtdz(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
422 & tclm(i,j,k+1,itrc))- &
423 & (t(i,j,k ,nrhs,itrc)- &
424 & tclm(i,j,k ,itrc)))+ &
425 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
426 & tl_t(i,j,k ,nrhs,itrc))- &
427# ifdef TL_IOMS
428 & dtdz(i,j,k2)
429# endif
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# ifdef TL_IOMS
438 & dtdz(i,j,k2)
439# endif
440 END IF
441#else
442 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
443 & t(i,j,k ,nrhs,itrc))
444 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
445 & t(i,j,k ,nrhs,itrc))+ &
446 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
447 & tl_t(i,j,k ,nrhs,itrc))- &
448# ifdef TL_IOMS
449 & dtdz(i,j,k2)
450# endif
451#endif
452 END DO
453 END DO
454 END IF
455 IF (k.gt.0) THEN
456 DO j=jmin,jmax
457 DO i=imin,imax+1
458#ifdef DIFF_3DCOEF
459# ifdef TS_U3ADV_SPLIT
460 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
461# else
462 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
463 & on_u(i,j)
464# endif
465#else
466 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
467 & on_u(i,j)
468#endif
469 fx(i,j)=cff* &
470 & (hz(i,j,k)+hz(i-1,j,k))* &
471 & (dtdx(i,j,k1)- &
472 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
473 & (dtdz(i-1,j,k1)+ &
474 & dtdz(i ,j,k2))+ &
475 & max(dzdx(i,j,k1),0.0_r8)* &
476 & (dtdz(i-1,j,k2)+ &
477 & dtdz(i ,j,k1))))
478 tl_fx(i,j)=cff* &
479 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
480 & (dtdx(i,j,k1)- &
481 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
482 & (dtdz(i-1,j,k1)+ &
483 & dtdz(i ,j,k2))+ &
484 & max(dzdx(i,j,k1),0.0_r8)* &
485 & (dtdz(i-1,j,k2)+ &
486 & dtdz(i ,j,k1))))+ &
487 & (hz(i,j,k)+hz(i-1,j,k))* &
488 & (tl_dtdx(i,j,k1)- &
489 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
490 & (tl_dtdz(i-1,j,k1)+ &
491 & tl_dtdz(i ,j,k2))+ &
492 & max(dzdx(i,j,k1),0.0_r8)* &
493 & (tl_dtdz(i-1,j,k2)+ &
494 & tl_dtdz(i ,j,k1)))- &
495 & 0.5_r8*((0.5_r8+ &
496 & sign(0.5_r8,-dzdx(i,j,k1)))* &
497 & tl_dzdx(i,j,k1)* &
498 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
499 & (0.5_r8+ &
500 & sign(0.5_r8, dzdx(i,j,k1)))* &
501 & tl_dzdx(i,j,k1)* &
502 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))- &
503#ifdef TL_IOMS
504 & cff* &
505 & (hz(i,j,k)+hz(i-1,j,k))* &
506 & (dtdx(i,j,k1)- &
507 & (min(dzdx(i,j,k1),0.0_r8)* &
508 & (dtdz(i-1,j,k1)+ &
509 & dtdz(i ,j,k2))+ &
510 & max(dzdx(i,j,k1),0.0_r8)* &
511 & (dtdz(i-1,j,k2)+ &
512 & dtdz(i ,j,k1))))
513#endif
514 END DO
515 END DO
516 DO j=jmin,jmax+1
517 DO i=imin,imax
518#ifdef DIFF_3DCOEF
519# ifdef TS_U3ADV_SPLIT
520 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
521# else
522 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
523 & om_v(i,j)
524# endif
525#else
526 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
527 & om_v(i,j)
528#endif
529 fe(i,j)=cff* &
530 & (hz(i,j,k)+hz(i,j-1,k))* &
531 & (dtde(i,j,k1)- &
532 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
533 & (dtdz(i,j-1,k1)+ &
534 & dtdz(i,j ,k2))+ &
535 & max(dzde(i,j,k1),0.0_r8)* &
536 & (dtdz(i,j-1,k2)+ &
537 & dtdz(i,j ,k1))))
538 tl_fe(i,j)=cff* &
539 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
540 & (dtde(i,j,k1)- &
541 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
542 & (dtdz(i,j-1,k1)+ &
543 & dtdz(i,j ,k2))+ &
544 & max(dzde(i,j,k1),0.0_r8)* &
545 & (dtdz(i,j-1,k2)+ &
546 & dtdz(i,j ,k1))))+ &
547 & (hz(i,j,k)+hz(i,j-1,k))* &
548 & (tl_dtde(i,j,k1)- &
549 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
550 & (tl_dtdz(i,j-1,k1)+ &
551 & tl_dtdz(i,j ,k2))+ &
552 & max(dzde(i,j,k1),0.0_r8)* &
553 & (tl_dtdz(i,j-1,k2)+ &
554 & tl_dtdz(i,j ,k1)))- &
555 & 0.5_r8*((0.5_r8+ &
556 & sign(0.5_r8,-dzde(i,j,k1)))* &
557 & tl_dzde(i,j,k1)* &
558 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
559 & (0.5_r8+ &
560 & sign(0.5_r8, dzde(i,j,k1)))* &
561 & tl_dzde(i,j,k1)* &
562 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))- &
563#ifdef TL_IOMS
564 & cff* &
565 & (hz(i,j,k)+hz(i,j-1,k))* &
566 & (dtde(i,j,k1)- &
567 & (min(dzde(i,j,k1),0.0_r8)* &
568 & (dtdz(i,j-1,k1)+ &
569 & dtdz(i,j ,k2))+ &
570 & max(dzde(i,j,k1),0.0_r8)* &
571 & (dtdz(i,j-1,k2)+ &
572 & dtdz(i,j ,k1))))
573#endif
574 END DO
575 END DO
576 IF (k.lt.n(ng)) THEN
577 DO j=jmin,jmax
578 DO i=imin,imax
579#ifdef DIFF_3DCOEF
580# ifdef TS_U3ADV_SPLIT
581 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
582 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
583 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
584 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
585# else
586 difx=0.5_r8*diff3d_r(i,j,k)
587 dife=difx
588# endif
589#else
590 difx=0.5_r8*diff4(i,j,itrc)
591 dife=difx
592#endif
593 cff1=min(dzdx(i ,j,k1),0.0_r8)
594 cff2=min(dzdx(i+1,j,k2),0.0_r8)
595 cff3=max(dzdx(i ,j,k2),0.0_r8)
596 cff4=max(dzdx(i+1,j,k1),0.0_r8)
597 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
598 & tl_dzdx(i ,j,k1)
599 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
600 & tl_dzdx(i+1,j,k2)
601 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
602 & tl_dzdx(i ,j,k2)
603 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
604 & tl_dzdx(i+1,j,k1)
605 fs(i,j,k2)=difx* &
606 & (cff1*(cff1*dtdz(i,j,k2)- &
607 & dtdx(i ,j,k1))+ &
608 & cff2*(cff2*dtdz(i,j,k2)- &
609 & dtdx(i+1,j,k2))+ &
610 & cff3*(cff3*dtdz(i,j,k2)- &
611 & dtdx(i ,j,k2))+ &
612 & cff4*(cff4*dtdz(i,j,k2)- &
613 & dtdx(i+1,j,k1)))
614 tl_fs(i,j,k2)=difx* &
615 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
616 & dtdx(i ,j,k1))+ &
617 & tl_cff2*(cff2*dtdz(i,j,k2)- &
618 & dtdx(i+1,j,k2))+ &
619 & tl_cff3*(cff3*dtdz(i,j,k2)- &
620 & dtdx(i ,j,k2))+ &
621 & tl_cff4*(cff4*dtdz(i,j,k2)- &
622 & dtdx(i+1,j,k1))+ &
623 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
624 & cff1*tl_dtdz(i,j,k2)- &
625 & tl_dtdx(i ,j,k1))+ &
626 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
627 & cff2*tl_dtdz(i,j,k2)- &
628 & tl_dtdx(i+1,j,k2))+ &
629 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
630 & cff3*tl_dtdz(i,j,k2)- &
631 & tl_dtdx(i ,j,k2))+ &
632 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
633 & cff4*tl_dtdz(i,j,k2)- &
634 & tl_dtdx(i+1,j,k1)))- &
635#ifdef TL_IOMS
636 & difx* &
637 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
638 & dtdx(i ,j,k1))+ &
639 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
640 & dtdx(i+1,j,k2))+ &
641 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
642 & dtdx(i ,j,k2))+ &
643 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
644 & dtdx(i+1,j,k1)))
645#endif
646!
647 cff1=min(dzde(i,j ,k1),0.0_r8)
648 cff2=min(dzde(i,j+1,k2),0.0_r8)
649 cff3=max(dzde(i,j ,k2),0.0_r8)
650 cff4=max(dzde(i,j+1,k1),0.0_r8)
651 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
652 & tl_dzde(i,j ,k1)
653 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
654 & tl_dzde(i,j+1,k2)
655 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
656 & tl_dzde(i,j ,k2)
657 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
658 & tl_dzde(i,j+1,k1)
659 fs(i,j,k2)=fs(i,j,k2)+ &
660 & dife* &
661 & (cff1*(cff1*dtdz(i,j,k2)- &
662 & dtde(i,j ,k1))+ &
663 & cff2*(cff2*dtdz(i,j,k2)- &
664 & dtde(i,j+1,k2))+ &
665 & cff3*(cff3*dtdz(i,j,k2)- &
666 & dtde(i,j ,k2))+ &
667 & cff4*(cff4*dtdz(i,j,k2)- &
668 & dtde(i,j+1,k1)))
669 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
670 & dife* &
671 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
672 & dtde(i,j ,k1))+ &
673 & tl_cff2*(cff2*dtdz(i,j,k2)- &
674 & dtde(i,j+1,k2))+ &
675 & tl_cff3*(cff3*dtdz(i,j,k2)- &
676 & dtde(i,j ,k2))+ &
677 & tl_cff4*(cff4*dtdz(i,j,k2)- &
678 & dtde(i,j+1,k1))+ &
679 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
680 & cff1*tl_dtdz(i,j,k2)- &
681 & tl_dtde(i,j ,k1))+ &
682 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
683 & cff2*tl_dtdz(i,j,k2)- &
684 & tl_dtde(i,j+1,k2))+ &
685 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
686 & cff3*tl_dtdz(i,j,k2)- &
687 & tl_dtde(i,j ,k2))+ &
688 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
689 & cff4*tl_dtdz(i,j,k2)- &
690 & tl_dtde(i,j+1,k1)))- &
691#ifdef TL_IOMS
692 & dife* &
693 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
694 & dtde(i,j ,k1))+ &
695 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
696 & dtde(i,j+1,k2))+ &
697 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
698 & dtde(i,j ,k2))+ &
699 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
700 & dtde(i,j+1,k1)))
701#endif
702 END DO
703 END DO
704 END IF
705!
706! Compute first harmonic operator, without mixing coefficient.
707! Multiply by the metrics of the second harmonic operator. Save
708! into work array "LapT".
709!
710 DO j=jmin,jmax
711 DO i=imin,imax
712 cff=pm(i,j)*pn(i,j)
713 cff1=1.0_r8/hz(i,j,k)
714 tl_cff1=-cff1*cff1*tl_hz(i,j,k)+ &
715#ifdef TL_IOMS
716 & 2.0_r8*cff1
717#endif
718 lapt(i,j,k)=cff1*(cff* &
719 & (fx(i+1,j)-fx(i,j)+ &
720 & fe(i,j+1)-fe(i,j))+ &
721 & (fs(i,j,k2)-fs(i,j,k1)))
722 tl_lapt(i,j,k)=tl_cff1*(cff* &
723 & (fx(i+1,j)-fx(i,j)+ &
724 & fe(i,j+1)-fe(i,j))+ &
725 & (fs(i,j,k2)-fs(i,j,k1)))+ &
726 & cff1*(cff* &
727 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
728 & tl_fe(i,j+1)-tl_fe(i,j))+ &
729 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))- &
730#ifdef TL_IOMS
731 & lapt(i,j,k)
732#endif
733 END DO
734 END DO
735 END IF
736 END DO k_loop1
737!
738! Apply boundary conditions (except periodic; closed or gradient)
739! to the first harmonic operator.
740!
741 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
742 IF (domain(ng)%Western_Edge(tile)) THEN
743 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
744 DO k=1,n(ng)
745 DO j=jmin,jmax
746 lapt(istr-1,j,k)=0.0_r8
747 tl_lapt(istr-1,j,k)=0.0_r8
748 END DO
749 END DO
750 ELSE
751 DO k=1,n(ng)
752 DO j=jmin,jmax
753 lapt(istr-1,j,k)=lapt(istr,j,k)
754 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
755 END DO
756 END DO
757 END IF
758 END IF
759 END IF
760!
761 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
762 IF (domain(ng)%Eastern_Edge(tile)) THEN
763 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
764 DO k=1,n(ng)
765 DO j=jmin,jmax
766 lapt(iend+1,j,k)=0.0_r8
767 tl_lapt(iend+1,j,k)=0.0_r8
768 END DO
769 END DO
770 ELSE
771 DO k=1,n(ng)
772 DO j=jmin,jmax
773 lapt(iend+1,j,k)=lapt(iend,j,k)
774 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
775 END DO
776 END DO
777 END IF
778 END IF
779 END IF
780!
781 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
782 IF (domain(ng)%Southern_Edge(tile)) THEN
783 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
784 DO k=1,n(ng)
785 DO i=imin,imax
786 lapt(i,jstr-1,k)=0.0_r8
787 tl_lapt(i,jstr-1,k)=0.0_r8
788 END DO
789 END DO
790 ELSE
791 DO k=1,n(ng)
792 DO i=imin,imax
793 lapt(i,jstr-1,k)=lapt(i,jstr,k)
794 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
795 END DO
796 END DO
797 END IF
798 END IF
799 END IF
800!
801 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
802 IF (domain(ng)%Northern_Edge(tile)) THEN
803 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
804 DO k=1,n(ng)
805 DO i=imin,imax
806 lapt(i,jend+1,k)=0.0_r8
807 tl_lapt(i,jend+1,k)=0.0_r8
808 END DO
809 END DO
810 ELSE
811 DO k=1,n(ng)
812 DO i=imin,imax
813 lapt(i,jend+1,k)=lapt(i,jend,k)
814 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
815 END DO
816 END DO
817 END IF
818 END IF
819 END IF
820!
821 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
822 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
823 IF (domain(ng)%SouthWest_Corner(tile)) THEN
824 DO k=1,n(ng)
825 lapt(istr-1,jstr-1,k)=0.5_r8* &
826 & (lapt(istr ,jstr-1,k)+ &
827 & lapt(istr-1,jstr ,k))
828 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
829 & (tl_lapt(istr ,jstr-1,k)+ &
830 & tl_lapt(istr-1,jstr ,k))
831 END DO
832 END IF
833 END IF
834
835 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
836 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
837 IF (domain(ng)%SouthEast_Corner(tile)) THEN
838 DO k=1,n(ng)
839 lapt(iend+1,jstr-1,k)=0.5_r8* &
840 & (lapt(iend ,jstr-1,k)+ &
841 & lapt(iend+1,jstr ,k))
842 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
843 & (tl_lapt(iend ,jstr-1,k)+ &
844 & tl_lapt(iend+1,jstr ,k))
845 END DO
846 END IF
847 END IF
848
849 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
850 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
851 IF (domain(ng)%NorthWest_Corner(tile)) THEN
852 DO k=1,n(ng)
853 lapt(istr-1,jend+1,k)=0.5_r8* &
854 & (lapt(istr ,jend+1,k)+ &
855 & lapt(istr-1,jend ,k))
856 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
857 & (tl_lapt(istr ,jend+1,k)+ &
858 & tl_lapt(istr-1,jend ,k))
859 END DO
860 END IF
861 END IF
862
863 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
864 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
865 IF (domain(ng)%NorthEast_Corner(tile)) THEN
866 DO k=1,n(ng)
867 lapt(iend+1,jend+1,k)=0.5_r8* &
868 & (lapt(iend ,jend+1,k)+ &
869 & lapt(iend+1,jend ,k))
870 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
871 & (tl_lapt(iend ,jend+1,k)+ &
872 & tl_lapt(iend+1,jend ,k))
873 END DO
874 END IF
875 END IF
876!
877! Compute tangent linear horizontal and vertical gradients associated
878! with the second rotated harmonic operator.
879!
880 k2=1
881 k_loop2: DO k=0,n(ng)
882 k1=k2
883 k2=3-k1
884 IF (k.lt.n(ng)) THEN
885 DO j=jstr,jend
886 DO i=istr,iend+1
887 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
888#ifdef MASKING
889 cff=cff*umask(i,j)
890#endif
891#ifdef WET_DRY_NOT_YET
892 cff=cff*umask_wet(i,j)
893#endif
894 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
895 & z_r(i-1,j,k+1))
896 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
897 & tl_z_r(i-1,j,k+1))
898 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
899 & lapt(i-1,j,k+1))
900 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
901 & tl_lapt(i-1,j,k+1))
902 END DO
903 END DO
904 DO j=jstr,jend+1
905 DO i=istr,iend
906 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
907#ifdef MASKING
908 cff=cff*vmask(i,j)
909#endif
910#ifdef WET_DRY_NOT_YET
911 cff=cff*vmask_wet(i,j)
912#endif
913 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
914 & z_r(i,j-1,k+1))
915 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
916 & tl_z_r(i,j-1,k+1))
917 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
918 & lapt(i,j-1,k+1))
919 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
920 & tl_lapt(i,j-1,k+1))
921 END DO
922 END DO
923 END IF
924 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
925 DO j=jstr-1,jend+1
926 DO i=istr-1,iend+1
927 dtdz(i,j,k2)=0.0_r8
928 tl_dtdz(i,j,k2)=0.0_r8
929 fs(i,j,k2)=0.0_r8
930 tl_fs(i,j,k2)=0.0_r8
931 END DO
932 END DO
933 ELSE
934 DO j=jstr-1,jend+1
935 DO i=istr-1,iend+1
936 cff=1.0_r8/(z_r(i,j,k+1)- &
937 & z_r(i,j,k ))
938 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
939 & tl_z_r(i,j,k ))+ &
940#ifdef TL_IOMS
941 & 2.0_r8*cff
942#endif
943 dtdz(i,j,k2)=cff*(lapt(i,j,k+1)- &
944 & lapt(i,j,k ))
945 tl_dtdz(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
946 & lapt(i,j,k ))+ &
947 & cff*(tl_lapt(i,j,k+1)- &
948 & tl_lapt(i,j,k ))- &
949#ifdef TL_IOMS
950 & dtdz(i,j,k2)
951#endif
952 END DO
953 END DO
954 END IF
955!
956! Compute components of the rotated tracer flux (T m4/s) along
957! geopotential surfaces.
958!
959 IF (k.gt.0) THEN
960 DO j=jstr,jend
961 DO i=istr,iend+1
962#ifdef DIFF_3DCOEF
963# ifdef TS_U3ADV_SPLIT
964 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
965# else
966 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
967 & on_u(i,j)
968# endif
969#else
970 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
971 & on_u(i,j)
972#endif
973!^ FX(i,j)=cff*
974!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
975!^ & (dTdx(i ,j,k1)- &
976!^ & 0.5_r8*(MIN(dZdx(i,j,k1),0.0_r8)* &
977!^ & (dTdz(i-1,j,k1)+ &
978!^ & dTdz(i ,j,k2))+ &
979!^ & MAX(dZdx(i,j,k1),0.0_r8)* &
980!^ & (dTdz(i-1,j,k2)+ &
981!^ & dTdz(i ,j,k1))))
982!^
983 tl_fx(i,j)=cff* &
984 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
985 & (dtdx(i,j,k1)- &
986 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
987 & (dtdz(i-1,j,k1)+ &
988 & dtdz(i ,j,k2))+ &
989 & max(dzdx(i,j,k1),0.0_r8)* &
990 & (dtdz(i-1,j,k2)+ &
991 & dtdz(i ,j,k1))))+ &
992 & (hz(i,j,k)+hz(i-1,j,k))* &
993 & (tl_dtdx(i,j,k1)- &
994 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
995 & (tl_dtdz(i-1,j,k1)+ &
996 & tl_dtdz(i ,j,k2))+ &
997 & max(dzdx(i,j,k1),0.0_r8)* &
998 & (tl_dtdz(i-1,j,k2)+ &
999 & tl_dtdz(i ,j,k1)))- &
1000 & 0.5_r8*((0.5_r8+ &
1001 & sign(0.5_r8,-dzdx(i,j,k1)))* &
1002 & tl_dzdx(i,j,k1)* &
1003 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
1004 & (0.5_r8+ &
1005 & sign(0.5_r8, dzdx(i,j,k1)))* &
1006 & tl_dzdx(i,j,k1)* &
1007 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))- &
1008#ifdef TL_IOMS
1009 & cff* &
1010 & (hz(i,j,k)+hz(i-1,j,k))* &
1011 & (dtdx(i ,j,k1)- &
1012 & (min(dzdx(i,j,k1),0.0_r8)* &
1013 & (dtdz(i-1,j,k1)+ &
1014 & dtdz(i ,j,k2))+ &
1015 & max(dzdx(i,j,k1),0.0_r8)* &
1016 & (dtdz(i-1,j,k2)+ &
1017 & dtdz(i ,j,k1))))
1018#endif
1019 END DO
1020 END DO
1021 DO j=jstr,jend+1
1022 DO i=istr,iend
1023#ifdef DIFF_3DCOEF
1024# ifdef TS_U3ADV_SPLIT
1025 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1026# else
1027 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1028 & om_v(i,j)
1029# endif
1030#else
1031 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1032 & om_v(i,j)
1033#endif
1034!^ FE(i,j)=cff* &
1035!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1036!^ & (dTde(i,j,k1)- &
1037!^ & 0.5_r8*(MIN(dZde(i,j,k1),0.0_r8)* &
1038!^ & (dTdz(i,j-1,k1)+ &
1039!^ & dTdz(i,j ,k2))+ &
1040!^ & MAX(dZde(i,j,k1),0.0_r8)* &
1041!^ & (dTdz(i,j-1,k2)+ &
1042!^ & dTdz(i,j ,k1))))
1043!^
1044 tl_fe(i,j)=cff* &
1045 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
1046 & (dtde(i,j,k1)- &
1047 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
1048 & (dtdz(i,j-1,k1)+ &
1049 & dtdz(i,j ,k2))+ &
1050 & max(dzde(i,j,k1),0.0_r8)* &
1051 & (dtdz(i,j-1,k2)+ &
1052 & dtdz(i,j ,k1))))+ &
1053 & (hz(i,j,k)+hz(i,j-1,k))* &
1054 & (tl_dtde(i,j,k1)- &
1055 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
1056 & (tl_dtdz(i,j-1,k1)+ &
1057 & tl_dtdz(i,j ,k2))+ &
1058 & max(dzde(i,j,k1),0.0_r8)* &
1059 & (tl_dtdz(i,j-1,k2)+ &
1060 & tl_dtdz(i,j ,k1)))- &
1061 & 0.5_r8*((0.5_r8+ &
1062 & sign(0.5_r8,-dzde(i,j,k1)))* &
1063 & tl_dzde(i,j,k1)* &
1064 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
1065 & (0.5_r8+ &
1066 & sign(0.5_r8, dzde(i,j,k1)))* &
1067 & tl_dzde(i,j,k1)* &
1068 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))- &
1069#ifdef TL_IOMS
1070 & cff* &
1071 & (hz(i,j,k)+hz(i,j-1,k))* &
1072 & (dtde(i,j,k1)- &
1073 & (min(dzde(i,j,k1),0.0_r8)* &
1074 & (dtdz(i,j-1,k1)+ &
1075 & dtdz(i,j ,k2))+ &
1076 & max(dzde(i,j,k1),0.0_r8)* &
1077 & (dtdz(i,j-1,k2)+ &
1078 & dtdz(i,j ,k1))))
1079#endif
1080 END DO
1081 END DO
1082 IF (k.lt.n(ng)) THEN
1083 DO j=jstr,jend
1084 DO i=istr,iend
1085#ifdef DIFF_3DCOEF
1086# ifdef TS_U3ADV_SPLIT
1087 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1088 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1089 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1090 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1091# else
1092 difx=0.5_r8*diff3d_r(i,j,k)
1093 dife=difx
1094# endif
1095#else
1096 difx=0.5_r8*diff4(i,j,itrc)
1097 dife=difx
1098#endif
1099 cff1=min(dzdx(i ,j,k1),0.0_r8)
1100 cff2=min(dzdx(i+1,j,k2),0.0_r8)
1101 cff3=max(dzdx(i ,j,k2),0.0_r8)
1102 cff4=max(dzdx(i+1,j,k1),0.0_r8)
1103 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
1104 & tl_dzdx(i ,j,k1)
1105 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
1106 & tl_dzdx(i+1,j,k2)
1107 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
1108 & tl_dzdx(i ,j,k2)
1109 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
1110 & tl_dzdx(i+1,j,k1)
1111!^
1112!^ FS(i,j,k2)=difx* &
1113!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1114!^ & dTdx(i ,j,k1))+ &
1115!^ & cff2*(cff2*dTdz(i,j,k2)- &
1116!^ & dTdx(i+1,j,k2))+ &
1117!^ & cff3*(cff3*dTdz(i,j,k2)- &
1118!^ & dTdx(i ,j,k2))+ &
1119!^ & cff4*(cff4*dTdz(i,j,k2)- &
1120!^ & dTdx(i+1,j,k1)))
1121!^
1122 tl_fs(i,j,k2)=difx* &
1123 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1124 & dtdx(i ,j,k1))+ &
1125 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1126 & dtdx(i+1,j,k2))+ &
1127 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1128 & dtdx(i ,j,k2))+ &
1129 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1130 & dtdx(i+1,j,k1))+ &
1131 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1132 & cff1*tl_dtdz(i,j,k2)- &
1133 & tl_dtdx(i ,j,k1))+ &
1134 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1135 & cff2*tl_dtdz(i,j,k2)- &
1136 & tl_dtdx(i+1,j,k2))+ &
1137 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1138 & cff3*tl_dtdz(i,j,k2)- &
1139 & tl_dtdx(i ,j,k2))+ &
1140 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1141 & cff4*tl_dtdz(i,j,k2)- &
1142 & tl_dtdx(i+1,j,k1)))- &
1143#ifdef TL_IOMS
1144 & difx* &
1145 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
1146 & dtdx(i ,j,k1))+ &
1147 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
1148 & dtdx(i+1,j,k2))+ &
1149 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
1150 & dtdx(i ,j,k2))+ &
1151 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
1152 & dtdx(i+1,j,k1)))
1153#endif
1154!
1155 cff1=min(dzde(i,j ,k1),0.0_r8)
1156 cff2=min(dzde(i,j+1,k2),0.0_r8)
1157 cff3=max(dzde(i,j ,k2),0.0_r8)
1158 cff4=max(dzde(i,j+1,k1),0.0_r8)
1159 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
1160 & tl_dzde(i,j ,k1)
1161 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
1162 & tl_dzde(i,j+1,k2)
1163 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
1164 & tl_dzde(i,j ,k2)
1165 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
1166 & tl_dzde(i,j+1,k1)
1167!^ FS(i,j,k2)=FS(i,j,k2)+ &
1168!^ & dife* &
1169!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1170!^ & dTde(i,j ,k1))+ &
1171!^ & cff2*(cff2*dTdz(i,j,k2)- &
1172!^ & dTde(i,j+1,k2))+ &
1173!^ & cff3*(cff3*dTdz(i,j,k2)- &
1174!^ & dTde(i,j ,k2))+ &
1175!^ & cff4*(cff4*dTdz(i,j,k2)- &
1176!^ & dTde(i,j+1,k1)))
1177!^
1178 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
1179 & dife* &
1180 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1181 & dtde(i,j ,k1))+ &
1182 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1183 & dtde(i,j+1,k2))+ &
1184 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1185 & dtde(i,j ,k2))+ &
1186 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1187 & dtde(i,j+1,k1))+ &
1188 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1189 & cff1*tl_dtdz(i,j,k2)- &
1190 & tl_dtde(i,j ,k1))+ &
1191 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1192 & cff2*tl_dtdz(i,j,k2)- &
1193 & tl_dtde(i,j+1,k2))+ &
1194 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1195 & cff3*tl_dtdz(i,j,k2)- &
1196 & tl_dtde(i,j ,k2))+ &
1197 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1198 & cff4*tl_dtdz(i,j,k2)- &
1199 & tl_dtde(i,j+1,k1)))- &
1200#ifdef TL_IOMS
1201 & dife* &
1202 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
1203 & dtde(i,j ,k1))+ &
1204 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
1205 & dtde(i,j+1,k2))+ &
1206 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
1207 & dtde(i,j ,k2))+ &
1208 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
1209 & dtde(i,j+1,k1)))
1210#endif
1211 END DO
1212 END DO
1213 END IF
1214!
1215! Time-step biharmonic, geopotential diffusion term (m Tunits).
1216!
1217 DO j=jstr,jend
1218 DO i=istr,iend
1219!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1220!^ & (FX(i+1,j )-FX(i,j)+ &
1221!^ & FE(i ,j+1)-FE(i,j))+ &
1222!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1223!^
1224 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1225 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1226 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1227 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1228!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1229!^
1230 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1231#ifdef DIAGNOSTICS_TS
1232!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1233#endif
1234 END DO
1235 END DO
1236 END IF
1237 END DO k_loop2
1238 END DO t_loop
1239!
1240 RETURN
1241 END SUBROUTINE rp_t3dmix4_geo_tile
1242
1243 END MODULE rp_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::irpm
integer, parameter irpm
Definition mod_param.F:664
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::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
rp_t3dmix4_mod
Definition rp_t3dmix4_geo.h:1
rp_t3dmix4_mod::rp_t3dmix4_geo_tile
subroutine rp_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 rp_t3dmix4_geo.h:129
rp_t3dmix4_mod::rp_t3dmix4
subroutine, public rp_t3dmix4(ng, tile)
Definition rp_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