ROMS
Loading...
Searching...
No Matches
t3dmix4_geo.h
Go to the documentation of this file.
1 MODULE t3dmix4_mod
2!
3!git $Id$
4!=======================================================================
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md Hernan G. Arango !
8!========================================== Alexander F. Shchepetkin ===
9! !
10! This subroutine computes horizontal biharmonic mixing of tracers !
11! along geopotential surfaces. !
12! !
13!=======================================================================
14!
15 implicit none
16!
17 PRIVATE
18 PUBLIC t3dmix4
19!
20 CONTAINS
21!
22!***********************************************************************
23 SUBROUTINE t3dmix4 (ng, tile)
24!***********************************************************************
25!
26 USE mod_param
27#ifdef TS_MIX_CLIMA
28 USE mod_clima
29#endif
30#ifdef DIAGNOSTICS_TS
31 USE mod_diags
32#endif
33 USE mod_grid
34 USE mod_mixing
35 USE mod_ocean
36 USE mod_stepping
37!
38! Imported variable declarations.
39!
40 integer, intent(in) :: ng, tile
41!
42! Local variable declarations.
43!
44 character (len=*), parameter :: myfile = &
45 & __FILE__
46!
47#include "tile.h"
48!
49#ifdef PROFILE
50 CALL wclock_on (ng, inlm, 28, __line__, myfile)
51#endif
52 CALL t3dmix4_geo_tile (ng, tile, &
53 & lbi, ubi, lbj, ubj, &
54 & imins, imaxs, jmins, jmaxs, &
55 & nrhs(ng), nstp(ng), nnew(ng), &
56#ifdef MASKING
57 & grid(ng) % umask, &
58 & grid(ng) % vmask, &
59#endif
60#ifdef WET_DRY
61 & grid(ng) % umask_wet, &
62 & grid(ng) % vmask_wet, &
63#endif
64 & grid(ng) % om_v, &
65 & grid(ng) % on_u, &
66 & grid(ng) % pm, &
67 & grid(ng) % pn, &
68 & grid(ng) % Hz, &
69 & grid(ng) % z_r, &
70#ifdef DIFF_3DCOEF
71# ifdef TS_U3ADV_SPLIT
72 & mixing(ng) % diff3d_u, &
73 & mixing(ng) % diff3d_v, &
74# else
75 & mixing(ng) % diff3d_r, &
76# endif
77#else
78 & mixing(ng) % diff4, &
79#endif
80#ifdef TS_MIX_CLIMA
81 & clima(ng) % tclm, &
82#endif
83#ifdef DIAGNOSTICS_TS
84 & diags(ng) % DiaTwrk, &
85#endif
86 & ocean(ng) % t)
87#ifdef PROFILE
88 CALL wclock_off (ng, inlm, 28, __line__, myfile)
89#endif
90!
91 RETURN
92 END SUBROUTINE t3dmix4
93!
94!***********************************************************************
95 SUBROUTINE t3dmix4_geo_tile (ng, tile, &
96 & LBi, UBi, LBj, UBj, &
97 & IminS, ImaxS, JminS, JmaxS, &
98 & nrhs, nstp, nnew, &
99#ifdef MASKING
100 & umask, vmask, &
101#endif
102#ifdef WET_DRY
103 & umask_wet, vmask_wet, &
104#endif
105 & om_v, on_u, pm, pn, &
106 & Hz, z_r, &
107#ifdef DIFF_3DCOEF
108# ifdef TS_U3ADV_SPLIT
109 & diff3d_u, diff3d_v, &
110# else
111 & diff3d_r, &
112# endif
113#else
114 & diff4, &
115#endif
116#ifdef TS_MIX_CLIMA
117 & tclm, &
118#endif
119#ifdef DIAGNOSTICS_TS
120 & DiaTwrk, &
121#endif
122 & t)
123!***********************************************************************
124!
125 USE mod_param
126 USE mod_ncparam
127 USE mod_scalars
128!
129! Imported variable declarations.
130!
131 integer, intent(in) :: ng, tile
132 integer, intent(in) :: LBi, UBi, LBj, UBj
133 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
134 integer, intent(in) :: nrhs, nstp, nnew
135
136#ifdef ASSUMED_SHAPE
137# ifdef MASKING
138 real(r8), intent(in) :: umask(LBi:,LBj:)
139 real(r8), intent(in) :: vmask(LBi:,LBj:)
140# endif
141# ifdef WET_DRY
142 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
143 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
144# endif
145# ifdef DIFF_3DCOEF
146# ifdef TS_U3ADV_SPLIT
147 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
148 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
149# else
150 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
151# endif
152# else
153 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
154# endif
155 real(r8), intent(in) :: om_v(LBi:,LBj:)
156 real(r8), intent(in) :: on_u(LBi:,LBj:)
157 real(r8), intent(in) :: pm(LBi:,LBj:)
158 real(r8), intent(in) :: pn(LBi:,LBj:)
159 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
160 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
161# ifdef TS_MIX_CLIMA
162 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
163# endif
164# ifdef DIAGNOSTICS_TS
165 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
166# endif
167 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
168#else
169# ifdef MASKING
170 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
171 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
172# endif
173# ifdef WET_DRY
174 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
175 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
176# endif
177# ifdef DIFF_3DCOEF
178# ifdef TS_U3ADV_SPLIT
179 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
180 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
181# else
182 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
183# endif
184# else
185 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
186# endif
187 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
188 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
189 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
190 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
191 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
192 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
193# ifdef TS_MIX_CLIMA
194 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
195# endif
196# ifdef DIAGNOSTICS_TS
197 real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
198 & NDT)
199# endif
200 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
201#endif
202!
203! Local variable declarations.
204!
205 integer :: Imin, Imax, Jmin, Jmax
206 integer :: i, itrc, j, k, k1, k2
207
208 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
209
210 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
211
212 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
213 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
214
215 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
216 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
217 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
218 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdz
219 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde
220 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx
221
222#include "set_bounds.h"
223!
224!-----------------------------------------------------------------------
225! Compute horizontal biharmonic diffusion along geopotential
226! surfaces. The biharmonic operator is computed by applying
227! the harmonic operator twice.
228!-----------------------------------------------------------------------
229!
230! Set local I- and J-ranges.
231!
232 IF (ewperiodic(ng)) THEN
233 imin=istr-1
234 imax=iend+1
235 ELSE
236 imin=max(istr-1,1)
237 imax=min(iend+1,lm(ng))
238 END IF
239 IF (nsperiodic(ng)) THEN
240 jmin=jstr-1
241 jmax=jend+1
242 ELSE
243 jmin=max(jstr-1,1)
244 jmax=min(jend+1,mm(ng))
245 END IF
246!
247! Compute horizontal and vertical gradients associated with the
248! first rotated harmonic operator. Notice the recursive blocking
249! sequence. The vertical placement of the gradients is:
250!
251! dTdx,dTde(:,:,k1) k rho-points
252! dTdx,dTde(:,:,k2) k+1 rho-points
253! FS,dTdz(:,:,k1) k-1/2 W-points
254! FS,dTdz(:,:,k2) k+1/2 W-points
255!
256#ifdef TS_MIX_STABILITY
257! In order to increase stability, the rotated biharmonic is applied
258! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
259!
260#endif
261
262 t_loop : DO itrc=1,nt(ng)
263 k2=1
264 k_loop1 : DO k=0,n(ng)
265 k1=k2
266 k2=3-k1
267 IF (k.lt.n(ng)) THEN
268 DO j=jmin,jmax
269 DO i=imin,imax+1
270 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
271#ifdef MASKING
272 cff=cff*umask(i,j)
273#endif
274#ifdef WET_DRY
275 cff=cff*umask_wet(i,j)
276#endif
277 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
278 & z_r(i-1,j,k+1))
279#if defined TS_MIX_STABILITY
280 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
281 & t(i-1,j,k+1,nrhs,itrc))+ &
282 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
283 & t(i-1,j,k+1,nstp,itrc)))
284#elif defined TS_MIX_CLIMA
285 IF (ltracerclm(itrc,ng)) THEN
286 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
287 & tclm(i ,j,k+1,itrc))- &
288 & (t(i-1,j,k+1,nrhs,itrc)- &
289 & tclm(i-1,j,k+1,itrc)))
290 ELSE
291 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
292 & t(i-1,j,k+1,nrhs,itrc))
293 END IF
294#else
295 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
296 & t(i-1,j,k+1,nrhs,itrc))
297#endif
298 END DO
299 END DO
300 DO j=jmin,jmax+1
301 DO i=imin,imax
302 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
303#ifdef MASKING
304 cff=cff*vmask(i,j)
305#endif
306#ifdef WET_DRY
307 cff=cff*vmask_wet(i,j)
308#endif
309 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
310 & z_r(i,j-1,k+1))
311#if defined TS_MIX_STABILITY
312 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
313 & t(i,j-1,k+1,nrhs,itrc))+ &
314 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
315 & t(i,j-1,k+1,nstp,itrc)))
316#elif defined TS_MIX_CLIMA
317 IF (ltracerclm(itrc,ng)) THEN
318 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
319 & tclm(i,j ,k+1,itrc))- &
320 & (t(i,j-1,k+1,nrhs,itrc)- &
321 & tclm(i,j-1,k+1,itrc)))
322 ELSE
323 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
324 & t(i,j-1,k+1,nrhs,itrc))
325 END IF
326#else
327 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
328 & t(i,j-1,k+1,nrhs,itrc))
329#endif
330 END DO
331 END DO
332 END IF
333 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
334 DO j=jmin-1,jmax+1
335 DO i=imin-1,imax+1
336 dtdz(i,j,k2)=0.0_r8
337 fs(i,j,k2)=0.0_r8
338 END DO
339 END DO
340 ELSE
341 DO j=jmin-1,jmax+1
342 DO i=imin-1,imax+1
343 cff=1.0_r8/(z_r(i,j,k+1)- &
344 & z_r(i,j,k ))
345#if defined TS_MIX_STABILITY
346 dtdz(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
347 & t(i,j,k ,nrhs,itrc))+ &
348 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
349 & t(i,j,k ,nstp,itrc)))
350#elif defined TS_MIX_CLIMA
351 IF (ltracerclm(itrc,ng)) THEN
352 dtdz(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
353 & tclm(i,j,k+1,itrc))- &
354 & (t(i,j,k ,nrhs,itrc)- &
355 & tclm(i,j,k ,itrc)))
356 ELSE
357 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
358 & t(i,j,k ,nrhs,itrc))
359 END IF
360#else
361 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
362 & t(i,j,k ,nrhs,itrc))
363#endif
364 END DO
365 END DO
366 END IF
367 IF (k.gt.0) THEN
368 DO j=jmin,jmax
369 DO i=imin,imax+1
370#ifdef DIFF_3DCOEF
371# ifdef TS_U3ADV_SPLIT
372 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
373# else
374 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
375 & on_u(i,j)
376# endif
377#else
378 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
379 & on_u(i,j)
380#endif
381 fx(i,j)=cff* &
382 & (hz(i,j,k)+hz(i-1,j,k))* &
383 & (dtdx(i,j,k1)- &
384 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
385 & (dtdz(i-1,j,k1)+ &
386 & dtdz(i ,j,k2))+ &
387 & max(dzdx(i,j,k1),0.0_r8)* &
388 & (dtdz(i-1,j,k2)+ &
389 & dtdz(i ,j,k1))))
390 END DO
391 END DO
392 DO j=jmin,jmax+1
393 DO i=imin,imax
394#ifdef DIFF_3DCOEF
395# ifdef TS_U3ADV_SPLIT
396 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
397# else
398 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
399 & om_v(i,j)
400# endif
401#else
402 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
403 & om_v(i,j)
404#endif
405 fe(i,j)=cff* &
406 & (hz(i,j,k)+hz(i,j-1,k))* &
407 & (dtde(i,j,k1)- &
408 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
409 & (dtdz(i,j-1,k1)+ &
410 & dtdz(i,j ,k2))+ &
411 & max(dzde(i,j,k1),0.0_r8)* &
412 & (dtdz(i,j-1,k2)+ &
413 & dtdz(i,j ,k1))))
414 END DO
415 END DO
416 IF (k.lt.n(ng)) THEN
417 DO j=jmin,jmax
418 DO i=imin,imax
419#ifdef DIFF_3DCOEF
420# ifdef TS_U3ADV_SPLIT
421 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
422 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
423 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
424 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
425# else
426 difx=0.5_r8*diff3d_r(i,j,k)
427 dife=difx
428# endif
429#else
430 difx=0.5_r8*diff4(i,j,itrc)
431 dife=difx
432#endif
433 cff1=min(dzdx(i ,j,k1),0.0_r8)
434 cff2=min(dzdx(i+1,j,k2),0.0_r8)
435 cff3=max(dzdx(i ,j,k2),0.0_r8)
436 cff4=max(dzdx(i+1,j,k1),0.0_r8)
437 fs(i,j,k2)=difx* &
438 & (cff1*(cff1*dtdz(i,j,k2)- &
439 & dtdx(i ,j,k1))+ &
440 & cff2*(cff2*dtdz(i,j,k2)- &
441 & dtdx(i+1,j,k2))+ &
442 & cff3*(cff3*dtdz(i,j,k2)- &
443 & dtdx(i ,j,k2))+ &
444 & cff4*(cff4*dtdz(i,j,k2)- &
445 & dtdx(i+1,j,k1)))
446!
447 cff1=min(dzde(i,j ,k1),0.0_r8)
448 cff2=min(dzde(i,j+1,k2),0.0_r8)
449 cff3=max(dzde(i,j ,k2),0.0_r8)
450 cff4=max(dzde(i,j+1,k1),0.0_r8)
451 fs(i,j,k2)=fs(i,j,k2)+ &
452 & dife* &
453 & (cff1*(cff1*dtdz(i,j,k2)- &
454 & dtde(i,j ,k1))+ &
455 & cff2*(cff2*dtdz(i,j,k2)- &
456 & dtde(i,j+1,k2))+ &
457 & cff3*(cff3*dtdz(i,j,k2)- &
458 & dtde(i,j ,k2))+ &
459 & cff4*(cff4*dtdz(i,j,k2)- &
460 & dtde(i,j+1,k1)))
461 END DO
462 END DO
463 END IF
464!
465! Compute first harmonic operator, without mixing coefficient.
466! Multiply by the metrics of the second harmonic operator. Save
467! into work array "LapT".
468!
469 DO j=jmin,jmax
470 DO i=imin,imax
471 cff=pm(i,j)*pn(i,j)
472 cff1=1.0_r8/hz(i,j,k)
473 lapt(i,j,k)=cff1*(cff* &
474 & (fx(i+1,j)-fx(i,j)+ &
475 & fe(i,j+1)-fe(i,j))+ &
476 & (fs(i,j,k2)-fs(i,j,k1)))
477 END DO
478 END DO
479 END IF
480 END DO k_loop1
481!
482! Apply boundary conditions (except periodic; closed or gradient)
483! to the first harmonic operator.
484!
485 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
486 IF (domain(ng)%Western_Edge(tile)) THEN
487 IF (lbc(iwest,istvar(itrc),ng)%closed) THEN
488 DO k=1,n(ng)
489 DO j=jmin,jmax
490 lapt(istr-1,j,k)=0.0_r8
491 END DO
492 END DO
493 ELSE
494 DO k=1,n(ng)
495 DO j=jmin,jmax
496 lapt(istr-1,j,k)=lapt(istr,j,k)
497 END DO
498 END DO
499 END IF
500 END IF
501 END IF
502!
503 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
504 IF (domain(ng)%Eastern_Edge(tile)) THEN
505 IF (lbc(ieast,istvar(itrc),ng)%closed) THEN
506 DO k=1,n(ng)
507 DO j=jmin,jmax
508 lapt(iend+1,j,k)=0.0_r8
509 END DO
510 END DO
511 ELSE
512 DO k=1,n(ng)
513 DO j=jmin,jmax
514 lapt(iend+1,j,k)=lapt(iend,j,k)
515 END DO
516 END DO
517 END IF
518 END IF
519 END IF
520!
521 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
522 IF (domain(ng)%Southern_Edge(tile)) THEN
523 IF (lbc(isouth,istvar(itrc),ng)%closed) THEN
524 DO k=1,n(ng)
525 DO i=imin,imax
526 lapt(i,jstr-1,k)=0.0_r8
527 END DO
528 END DO
529 ELSE
530 DO k=1,n(ng)
531 DO i=imin,imax
532 lapt(i,jstr-1,k)=lapt(i,jstr,k)
533 END DO
534 END DO
535 END IF
536 END IF
537 END IF
538!
539 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
540 IF (domain(ng)%Northern_Edge(tile)) THEN
541 IF (lbc(inorth,istvar(itrc),ng)%closed) THEN
542 DO k=1,n(ng)
543 DO i=imin,imax
544 lapt(i,jend+1,k)=0.0_r8
545 END DO
546 END DO
547 ELSE
548 DO k=1,n(ng)
549 DO i=imin,imax
550 lapt(i,jend+1,k)=lapt(i,jend,k)
551 END DO
552 END DO
553 END IF
554 END IF
555 END IF
556!
557 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
558 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
559 IF (domain(ng)%SouthWest_Corner(tile)) THEN
560 DO k=1,n(ng)
561 lapt(istr-1,jstr-1,k)=0.5_r8* &
562 & (lapt(istr ,jstr-1,k)+ &
563 & lapt(istr-1,jstr ,k))
564 END DO
565 END IF
566 END IF
567
568 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
569 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
570 IF (domain(ng)%SouthEast_Corner(tile)) THEN
571 DO k=1,n(ng)
572 lapt(iend+1,jstr-1,k)=0.5_r8* &
573 & (lapt(iend ,jstr-1,k)+ &
574 & lapt(iend+1,jstr ,k))
575 END DO
576 END IF
577 END IF
578
579 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
580 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
581 IF (domain(ng)%NorthWest_Corner(tile)) THEN
582 DO k=1,n(ng)
583 lapt(istr-1,jend+1,k)=0.5_r8* &
584 & (lapt(istr ,jend+1,k)+ &
585 & lapt(istr-1,jend ,k))
586 END DO
587 END IF
588 END IF
589
590 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
591 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
592 IF (domain(ng)%NorthEast_Corner(tile)) THEN
593 DO k=1,n(ng)
594 lapt(iend+1,jend+1,k)=0.5_r8* &
595 & (lapt(iend ,jend+1,k)+ &
596 & lapt(iend+1,jend ,k))
597 END DO
598 END IF
599 END IF
600!
601! Compute horizontal and vertical gradients associated with the
602! second rotated harmonic operator.
603!
604 k2=1
605 k_loop2: DO k=0,n(ng)
606 k1=k2
607 k2=3-k1
608 IF (k.lt.n(ng)) THEN
609 DO j=jstr,jend
610 DO i=istr,iend+1
611 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
612#ifdef MASKING
613 cff=cff*umask(i,j)
614#endif
615#ifdef WET_DRY
616 cff=cff*umask_wet(i,j)
617#endif
618 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
619 & z_r(i-1,j,k+1))
620 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
621 & lapt(i-1,j,k+1))
622 END DO
623 END DO
624 DO j=jstr,jend+1
625 DO i=istr,iend
626 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
627#ifdef MASKING
628 cff=cff*vmask(i,j)
629#endif
630#ifdef WET_DRY
631 cff=cff*vmask_wet(i,j)
632#endif
633 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
634 & z_r(i,j-1,k+1))
635 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
636 & lapt(i,j-1,k+1))
637 END DO
638 END DO
639 END IF
640 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
641 DO j=jstr-1,jend+1
642 DO i=istr-1,iend+1
643 dtdz(i,j,k2)=0.0_r8
644 fs(i,j,k2)=0.0_r8
645 END DO
646 END DO
647 ELSE
648 DO j=jstr-1,jend+1
649 DO i=istr-1,iend+1
650 cff=1.0_r8/(z_r(i,j,k+1)- &
651 & z_r(i,j,k ))
652 dtdz(i,j,k2)=cff*(lapt(i,j,k+1)- &
653 & lapt(i,j,k ))
654 END DO
655 END DO
656 END IF
657!
658! Compute components of the rotated tracer flux (T m4/s) along
659! geopotential surfaces.
660!
661 IF (k.gt.0) THEN
662 DO j=jstr,jend
663 DO i=istr,iend+1
664#ifdef DIFF_3DCOEF
665# ifdef TS_U3ADV_SPLIT
666 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
667# else
668 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
669 & on_u(i,j)
670# endif
671#else
672 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
673 & on_u(i,j)
674#endif
675 fx(i,j)=cff* &
676 & (hz(i,j,k)+hz(i-1,j,k))* &
677 & (dtdx(i ,j,k1)- &
678 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
679 & (dtdz(i-1,j,k1)+ &
680 & dtdz(i ,j,k2))+ &
681 & max(dzdx(i,j,k1),0.0_r8)* &
682 & (dtdz(i-1,j,k2)+ &
683 & dtdz(i ,j,k1))))
684 END DO
685 END DO
686 DO j=jstr,jend+1
687 DO i=istr,iend
688#ifdef DIFF_3DCOEF
689# ifdef TS_U3ADV_SPLIT
690 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
691# else
692 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
693 & om_v(i,j)
694# endif
695#else
696 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
697 & om_v(i,j)
698#endif
699 fe(i,j)=cff* &
700 & (hz(i,j,k)+hz(i,j-1,k))* &
701 & (dtde(i,j,k1)- &
702 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
703 & (dtdz(i,j-1,k1)+ &
704 & dtdz(i,j ,k2))+ &
705 & max(dzde(i,j,k1),0.0_r8)* &
706 & (dtdz(i,j-1,k2)+ &
707 & dtdz(i,j ,k1))))
708 END DO
709 END DO
710 IF (k.lt.n(ng)) THEN
711 DO j=jstr,jend
712 DO i=istr,iend
713#ifdef DIFF_3DCOEF
714# ifdef TS_U3ADV_SPLIT
715 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
716 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
717 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
718 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
719# else
720 difx=0.5_r8*diff3d_r(i,j,k)
721 dife=difx
722# endif
723#else
724 difx=0.5_r8*diff4(i,j,itrc)
725 dife=difx
726#endif
727 cff1=min(dzdx(i ,j,k1),0.0_r8)
728 cff2=min(dzdx(i+1,j,k2),0.0_r8)
729 cff3=max(dzdx(i ,j,k2),0.0_r8)
730 cff4=max(dzdx(i+1,j,k1),0.0_r8)
731 fs(i,j,k2)=difx* &
732 & (cff1*(cff1*dtdz(i,j,k2)- &
733 & dtdx(i ,j,k1))+ &
734 & cff2*(cff2*dtdz(i,j,k2)- &
735 & dtdx(i+1,j,k2))+ &
736 & cff3*(cff3*dtdz(i,j,k2)- &
737 & dtdx(i ,j,k2))+ &
738 & cff4*(cff4*dtdz(i,j,k2)- &
739 & dtdx(i+1,j,k1)))
740!
741 cff1=min(dzde(i,j ,k1),0.0_r8)
742 cff2=min(dzde(i,j+1,k2),0.0_r8)
743 cff3=max(dzde(i,j ,k2),0.0_r8)
744 cff4=max(dzde(i,j+1,k1),0.0_r8)
745 fs(i,j,k2)=fs(i,j,k2)+ &
746 & dife* &
747 & (cff1*(cff1*dtdz(i,j,k2)- &
748 & dtde(i,j ,k1))+ &
749 & cff2*(cff2*dtdz(i,j,k2)- &
750 & dtde(i,j+1,k2))+ &
751 & cff3*(cff3*dtdz(i,j,k2)- &
752 & dtde(i,j ,k2))+ &
753 & cff4*(cff4*dtdz(i,j,k2)- &
754 & dtde(i,j+1,k1)))
755 END DO
756 END DO
757 END IF
758!
759! Time-step biharmonic, geopotential diffusion term (m Tunits).
760!
761 DO j=jstr,jend
762 DO i=istr,iend
763 cff=dt(ng)*pm(i,j)*pn(i,j)
764 cff1=cff*(fx(i+1,j )-fx(i,j))
765 cff2=cff*(fe(i ,j+1)-fe(i,j))
766 cff3=dt(ng)*(fs(i,j,k2)-fs(i,j,k1))
767 cff4=cff1+cff2+cff3
768 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff4
769#ifdef DIAGNOSTICS_TS
770 diatwrk(i,j,k,itrc,itxdif)=-cff1
771 diatwrk(i,j,k,itrc,itydif)=-cff2
772 diatwrk(i,j,k,itrc,itsdif)=-cff3
773 diatwrk(i,j,k,itrc,ithdif)=-cff4
774#endif
775 END DO
776 END DO
777 END IF
778 END DO k_loop2
779 END DO t_loop
780!
781 RETURN
782 END SUBROUTINE t3dmix4_geo_tile
783
784 END MODULE t3dmix4_mod
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
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::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::lbc
type(t_lbc), dimension(:,:,:), allocatable lbc
Definition mod_param.F:375
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::itxdif
integer itxdif
Definition mod_scalars.F:145
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::itsdif
integer itsdif
Definition mod_scalars.F:147
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::itydif
integer itydif
Definition mod_scalars.F:146
mod_scalars::ltracerclm
logical, dimension(:,:), allocatable ltracerclm
Definition mod_scalars.F:528
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::ithdif
integer ithdif
Definition mod_scalars.F:144
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
t3dmix4_mod
Definition t3dmix4_geo.h:1
t3dmix4_mod::t3dmix4
subroutine, public t3dmix4(ng, tile)
Definition t3dmix4_geo.h:24
t3dmix4_mod::t3dmix4_geo_tile
subroutine 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, z_r, diff3d_u, diff3d_v, diff3d_r, diff4, tclm, diatwrk, t)
Definition t3dmix4_geo.h:123
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