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