ROMS
Loading...
Searching...
No Matches
tl_conv_3d.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3 MODULE tl_conv_3d_mod
4
5#if defined TANGENT && defined FOUR_DVAR && defined SOLVE3D
6!
7!git $Id$
8!================================================== Hernan G. Arango ===
9! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
10! Licensed under a MIT/X style license !
11! See License_ROMS.md !
12!=======================================================================
13! !
14! These routines applies the background error covariance to data !
15! assimilation fields via the space convolution of the diffusion !
16! equation (filter) for 3D state variables. The diffusion filter !
17! is solved using an implicit or explicit vertical algorithm. !
18! !
19! For Gaussian (bell-shaped) correlations, the space convolution !
20! of the diffusion operator is an efficient way to estimate the !
21! finite domain error covariances. !
22! !
23! Notice that "z_r" and "Hz" are assumed to be time invariant in !
24! the spatial convolution. That it, they are not adjointable. !
25! !
26! On Input: !
27! !
28! ng Nested grid number. !
29! model Calling model identifier. !
30! Istr Starting tile index in the I-direction. !
31! Iend Ending tile index in the I-direction. !
32! Jstr Starting tile index in the J-direction. !
33! Jend Ending tile index in the J-direction. !
34! LBi I-dimension lower bound. !
35! UBi I-dimension upper bound. !
36! LBj J-dimension lower bound. !
37! UBj J-dimension upper bound. !
38! LBk K-dimension lower bound. !
39! UBk K-dimension upper bound. !
40! Nghost Number of ghost points. !
41! NHsteps Number of horizontal diffusion integration steps. !
42! NVsteps Number of vertical diffusion integration steps. !
43! DTsizeH Horizontal diffusion pseudo time-step size. !
44! DTsizeV Vertical diffusion pseudo time-step size. !
45! Kh Horizontal diffusion coefficients. !
46! Kv Vertical diffusion coefficients. !
47! tl_A 3D tangent linear state variable to diffuse. !
48! !
49! On Output: !
50! !
51! tl_A Convolved 3D tangent linear state variable. !
52! !
53! Routines: !
54! !
55! tl_conv_r3d_tile Tangent linear 3D convolution at RHO-points !
56! tl_conv_u3d_tile Tangent linear 3D convolution at U-points !
57! tl_conv_v3d_tile Tangent linear 3D convolution at V-points !
58! !
59!=======================================================================
60!
61 implicit none
62!
63 PUBLIC
64!
65 CONTAINS
66!
67!***********************************************************************
68 SUBROUTINE tl_conv_r3d_tile (ng, tile, model, &
69 & LBi, UBi, LBj, UBj, LBk, UBk, &
70 & IminS, ImaxS, JminS, JmaxS, &
71 & Nghost, NHsteps, NVsteps, &
72 & DTsizeH, DTsizeV, &
73 & Kh, Kv, &
74 & pm, pn, &
75# ifdef GEOPOTENTIAL_HCONV
76 & on_u, om_v, &
77# else
78 & pmon_u, pnom_v, &
79# endif
80# ifdef MASKING
81 & rmask, umask, vmask, &
82# endif
83 & Hz, z_r, &
84 & tl_A)
85!***********************************************************************
86!
87 USE mod_param
88 USE mod_scalars
89!
90 USE bc_3d_mod, ONLY: dabc_r3d_tile
91# ifdef DISTRIBUTE
92 USE mp_exchange_mod, ONLY : mp_exchange3d
93# endif
94!
95! Imported variable declarations.
96!
97 integer, intent(in) :: ng, tile, model
98 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
99 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
100 integer, intent(in) :: Nghost, NHsteps, NVsteps
101
102 real(r8), intent(in) :: DTsizeH, DTsizeV
103!
104# ifdef ASSUMED_SHAPE
105 real(r8), intent(in) :: pm(LBi:,LBj:)
106 real(r8), intent(in) :: pn(LBi:,LBj:)
107# ifdef GEOPOTENTIAL_HCONV
108 real(r8), intent(in) :: on_u(LBi:,LBj:)
109 real(r8), intent(in) :: om_v(LBi:,LBj:)
110# else
111 real(r8), intent(in) :: pmon_u(LBi:,LBj:)
112 real(r8), intent(in) :: pnom_v(LBi:,LBj:)
113# endif
114# ifdef MASKING
115 real(r8), intent(in) :: rmask(LBi:,LBj:)
116 real(r8), intent(in) :: umask(LBi:,LBj:)
117 real(r8), intent(in) :: vmask(LBi:,LBj:)
118# endif
119 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
120 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
121
122 real(r8), intent(in) :: Kh(LBi:,LBj:)
123 real(r8), intent(in) :: Kv(LBi:,LBj:,0:)
124
125 real(r8), intent(inout) :: tl_A(LBi:,LBj:,LBk:)
126# else
127 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
128 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
129# ifdef GEOPOTENTIAL_HCONV
130 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
131 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
132# else
133 real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
134 real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
135# endif
136# ifdef MASKING
137 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
138 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
139 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
140# endif
141 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
142 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
143
144 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
145 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
146 real(r8), intent(inout) :: tl_A(LBi:UBi,LBj:UBj,LBk:UBk)
147# endif
148!
149! Local variable declarations.
150!
151 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
152
153 real(r8) :: cff, cff1, cff2, cff3, cff4
154
155 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: tl_Awrk
156
157 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
158 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
159 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
160# ifdef VCONVOLUTION
161# ifndef SPLINES_VCONV
162 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FC
163# endif
164# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
165 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
166# endif
167# if defined IMPLICIT_VCONV || defined SPLINES_VCONV
168 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
169 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
170# ifdef SPLINES_VCONV
171 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
172# endif
173 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
174# else
175 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FS
176# endif
177# endif
178# ifdef GEOPOTENTIAL_HCONV
179 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx
180 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde
181 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FZ
182 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdz
183 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdx
184 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAde
185# endif
186
187# include "set_bounds.h"
188!
189!-----------------------------------------------------------------------
190! Space convolution of the diffusion equation for a 3D state variable
191! at RHO-points.
192!-----------------------------------------------------------------------
193!
194! Compute metrics factors. Notice that "z_r" and "Hz" are assumed to
195! be time invariant in the vertical convolution. Scratch array are
196! used for efficiency.
197!
198 DO j=jstr-1,jend+1
199 DO i=istr-1,iend+1
200 hfac(i,j)=dtsizeh*pm(i,j)*pn(i,j)
201# ifdef VCONVOLUTION
202# ifndef SPLINES_VCONV
203 fc(i,j,n(ng))=0.0_r8
204 DO k=1,n(ng)-1
205# ifdef IMPLICIT_VCONV
206 fc(i,j,k)=-dtsizev*kv(i,j,k)/(z_r(i,j,k+1)-z_r(i,j,k))
207# else
208 fc(i,j,k)=dtsizev*kv(i,j,k)/(z_r(i,j,k+1)-z_r(i,j,k))
209# endif
210 END DO
211 fc(i,j,0)=0.0_r8
212# endif
213# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
214 DO k=1,n(ng)
215 ohz(i,j,k)=1.0_r8/hz(i,j,k)
216 END DO
217# endif
218# endif
219 END DO
220 END DO
221!
222! Set integration indices and initial conditions.
223!
224 nold=1
225 nnew=2
226!^ CALL dabc_r3d_tile (ng, tile, &
227!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
228!^ & A)
229!^
230 CALL dabc_r3d_tile (ng, tile, &
231 & lbi, ubi, lbj, ubj, lbk, ubk, &
232 & tl_a)
233# ifdef DISTRIBUTE
234!^ CALL mp_exchange3d (ng, tile, model, 1, &
235!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
236!^ & Nghost, &
237!^ & EWperiodic(ng), NSperiodic(ng), &
238!^ & A)
239!^
240 CALL mp_exchange3d (ng, tile, model, 1, &
241 & lbi, ubi, lbj, ubj, lbk, ubk, &
242 & nghost, &
243 & ewperiodic(ng), nsperiodic(ng), &
244 & tl_a)
245# endif
246 DO k=1,n(ng)
247 DO j=jstr-1,jend+1
248 DO i=istr-1,iend+1
249!^ Awrk(i,j,k,Nold)=A(i,j,k)
250!^
251 tl_awrk(i,j,k,nold)=tl_a(i,j,k)
252 END DO
253 END DO
254 END DO
255!
256!-----------------------------------------------------------------------
257! Integrate horizontal diffusion equation.
258!-----------------------------------------------------------------------
259!
260 DO step=1,nhsteps
261
262# ifdef GEOPOTENTIAL_HCONV
263!
264! Diffusion along geopotential surfaces: Compute horizontal and
265! vertical gradients. Notice the recursive blocking sequence. The
266! vertical placement of the gradients is:
267!
268! dAdx,dAde(:,:,k1) k rho-points
269! dAdx,dAde(:,:,k2) k+1 rho-points
270! FZ,dAdz(:,:,k1) k-1/2 W-points
271! FZ,dAdz(:,:,k2) k+1/2 W-points
272!
273 k2=1
274 k_loop : DO k=0,n(ng)
275 k1=k2
276 k2=3-k1
277 IF (k.lt.n(ng)) THEN
278 DO j=jstr,jend
279 DO i=istr,iend+1
280 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
281# ifdef MASKING
282 cff=cff*umask(i,j)
283# endif
284 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
285 & z_r(i-1,j,k+1))
286# ifdef MASKING
287!^ dAdx(i,j,k2)=cff*(Awrk(i ,j,k+1,Nold)*rmask(i ,j)- &
288!^ & Awrk(i-1,j,k+1,Nold)*rmask(i-1,j))
289!^
290 tl_dadx(i,j,k2)=cff* &
291 & (tl_awrk(i ,j,k+1,nold)*rmask(i ,j)- &
292 & tl_awrk(i-1,j,k+1,nold)*rmask(i-1,j))
293# else
294!^ dAdx(i,j,k2)=cff*(Awrk(i ,j,k+1,Nold)- &
295!^ & Awrk(i-1,j,k+1,Nold))
296!^
297 tl_dadx(i,j,k2)=cff*(tl_awrk(i ,j,k+1,nold)- &
298 & tl_awrk(i-1,j,k+1,nold))
299# endif
300 END DO
301 END DO
302 DO j=jstr,jend+1
303 DO i=istr,iend
304 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
305# ifdef MASKING
306 cff=cff*vmask(i,j)
307# endif
308 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
309 & z_r(i,j-1,k+1))
310# ifdef MASKING
311!^ dAde(i,j,k2)=cff*(Awrk(i,j ,k+1,Nold)*rmask(i,j )- &
312!^ & Awrk(i,j-1,k+1,Nold)*rmask(i,j-1))
313!^
314 tl_dade(i,j,k2)=cff* &
315 & (tl_awrk(i,j ,k+1,nold)*rmask(i,j )- &
316 & tl_awrk(i,j-1,k+1,nold)*rmask(i,j-1))
317# else
318!^ dAde(i,j,k2)=cff*(Awrk(i,j ,k+1,Nold)- &
319!^ & Awrk(i,j-1,k+1,Nold))
320!^
321 tl_dade(i,j,k2)=cff*(tl_awrk(i,j ,k+1,nold)- &
322 & tl_awrk(i,j-1,k+1,nold))
323# endif
324 END DO
325 END DO
326 END IF
327 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
328 DO j=jstr-1,jend+1
329 DO i=istr-1,iend+1
330!^ dAdz(i,j,k2)=0.0_r8
331!^
332 tl_dadz(i,j,k2)=0.0_r8
333!^ FZ(i,j,k2)=0.0_r8
334!^
335 tl_fz(i,j,k2)=0.0_r8
336 END DO
337 END DO
338 ELSE
339 DO j=jstr-1,jend+1
340 DO i=istr-1,iend+1
341 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
342!^ dAdz(i,j,k2)=cff*(Awrk(i,j,k+1,Nold)- &
343!^ & Awrk(i,j,k ,Nold))
344!^
345 tl_dadz(i,j,k2)=cff*(tl_awrk(i,j,k+1,nold)- &
346 & tl_awrk(i,j,k ,nold))
347# ifdef MASKING
348!^ dAdz(i,j,k2)=dAdz(i,j,k2)*rmask(i,j)
349!^
350 tl_dadz(i,j,k2)=tl_dadz(i,j,k2)*rmask(i,j)
351# endif
352 END DO
353 END DO
354 END IF
355!
356! Compute components of the rotated A flux (A m3/s) along geopotential
357! surfaces.
358!
359 IF (k.gt.0) THEN
360 DO j=jstr,jend
361 DO i=istr,iend+1
362 cff=0.25_r8*(kh(i-1,j)+kh(i-1,j))*on_u(i,j)
363 cff1=min(dzdx(i,j,k1),0.0_r8)
364 cff2=max(dzdx(i,j,k1),0.0_r8)
365!^ FX(i,j)=cff* &
366!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
367!^ & (dAdx(i,j,k1)- &
368!^ & 0.5_r8*(cff1*(dAdz(i-1,j,k1)+ &
369!^ & dAdz(i ,j,k2))+ &
370!^ & cff2*(dAdz(i-1,j,k2)+ &
371!^ & dAdz(i ,j,k1))))
372!^
373 tl_fx(i,j)=cff* &
374 & (hz(i,j,k)+hz(i-1,j,k))* &
375 & (tl_dadx(i,j,k1)- &
376 & 0.5_r8*(cff1*(tl_dadz(i-1,j,k1)+ &
377 & tl_dadz(i ,j,k2))+ &
378 & cff2*(tl_dadz(i-1,j,k2)+ &
379 & tl_dadz(i ,j,k1))))
380 END DO
381 END DO
382 DO j=jstr,jend+1
383 DO i=istr,iend
384 cff=0.25_r8*(kh(i,j-1)+kh(i,j))*om_v(i,j)
385 cff1=min(dzde(i,j,k1),0.0_r8)
386 cff2=max(dzde(i,j,k1),0.0_r8)
387!^ FE(i,j)=cff* &
388!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
389!^ & (dAde(i,j,k1)- &
390!^ & 0.5_r8*(cff1*(dAdz(i,j-1,k1)+ &
391!^ & dAdz(i,j ,k2))+ &
392!^ & cff2*(dAdz(i,j-1,k2)+ &
393!^ & dAdz(i,j ,k1))))
394!^
395 tl_fe(i,j)=cff* &
396 & (hz(i,j,k)+hz(i,j-1,k))* &
397 & (tl_dade(i,j,k1)- &
398 & 0.5_r8*(cff1*(tl_dadz(i,j-1,k1)+ &
399 & tl_dadz(i,j ,k2))+ &
400 & cff2*(tl_dadz(i,j-1,k2)+ &
401 & tl_dadz(i,j ,k1))))
402 END DO
403 END DO
404 IF (k.lt.n(ng)) THEN
405 DO j=jstr,jend
406 DO i=istr,iend
407 cff=0.5_r8*kh(i,j)
408 cff1=min(dzdx(i ,j,k1),0.0_r8)
409 cff2=min(dzdx(i+1,j,k2),0.0_r8)
410 cff3=max(dzdx(i ,j,k2),0.0_r8)
411 cff4=max(dzdx(i+1,j,k1),0.0_r8)
412!^ FZ(i,j,k2)=cff* &
413!^ & (cff1*(cff1*dAdz(i,j,k2)-dAdx(i ,j,k1))+ &
414!^ & cff2*(cff2*dAdz(i,j,k2)-dAdx(i+1,j,k2))+ &
415!^ & cff3*(cff3*dAdz(i,j,k2)-dAdx(i ,j,k2))+ &
416!^ & cff4*(cff4*dAdz(i,j,k2)-dAdx(i+1,j,k1)))
417!^
418 tl_fz(i,j,k2)=cff* &
419 & (cff1*(cff1*tl_dadz(i,j,k2)- &
420 & tl_dadx(i ,j,k1))+ &
421 & cff2*(cff2*tl_dadz(i,j,k2)- &
422 & tl_dadx(i+1,j,k2))+ &
423 & cff3*(cff3*tl_dadz(i,j,k2)- &
424 & tl_dadx(i ,j,k2))+ &
425 & cff4*(cff4*tl_dadz(i,j,k2)- &
426 & tl_dadx(i+1,j,k1)))
427 cff1=min(dzde(i,j ,k1),0.0_r8)
428 cff2=min(dzde(i,j+1,k2),0.0_r8)
429 cff3=max(dzde(i,j ,k2),0.0_r8)
430 cff4=max(dzde(i,j+1,k1),0.0_r8)
431!^ FZ(i,j,k2)=FZ(i,j,k2)+ &
432!^ & cff* &
433!^ & (cff1*(cff1*dAdz(i,j,k2)-dAde(i,j ,k1))+ &
434!^ & cff2*(cff2*dAdz(i,j,k2)-dAde(i,j+1,k2))+ &
435!^ & cff3*(cff3*dAdz(i,j,k2)-dAde(i,j ,k2))+ &
436!^ & cff4*(cff4*dAdz(i,j,k2)-dAde(i,j+1,k1)))
437!^
438 tl_fz(i,j,k2)=tl_fz(i,j,k2)+ &
439 & cff* &
440 & (cff1*(cff1*tl_dadz(i,j,k2)- &
441 & tl_dade(i,j ,k1))+ &
442 & cff2*(cff2*tl_dadz(i,j,k2)- &
443 & tl_dade(i,j+1,k2))+ &
444 & cff3*(cff3*tl_dadz(i,j,k2)- &
445 & tl_dade(i,j ,k2))+ &
446 & cff4*(cff4*tl_dadz(i,j,k2)- &
447 & tl_dade(i,j+1,k1)))
448 END DO
449 END DO
450 END IF
451!
452! Time-step harmonic, geopotential diffusion term (m Tunits).
453!
454 DO j=jstr,jend
455 DO i=istr,iend
456!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
457!^ & Hfac(i,j)* &
458!^ & (FX(i+1,j )-FX(i,j)+ &
459!^ & FE(i ,j+1)-FE(i,j))+ &
460!^ & DTsizeH* &
461!^ & (FZ(i,j,k2)-FZ(i,j,k1))
462!^
463 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
464 & hfac(i,j)* &
465 & (tl_fx(i+1,j )-tl_fx(i,j)+ &
466 & tl_fe(i ,j+1)-tl_fe(i,j))+ &
467 & dtsizeh* &
468 & (tl_fz(i,j,k2)-tl_fz(i,j,k1))
469 END DO
470 END DO
471 END IF
472 END DO k_loop
473
474# else
475
476!
477! Diffusion along S-coordinates: compute XI- and ETA-components of
478! diffusive flux.
479!
480 DO k=1,n(ng)
481 DO j=jstr,jend
482 DO i=istr,iend+1
483!^ FX(i,j)=pmon_u(i,j)*0.5_r8*(Kh(i-1,j)+Kh(i,j))* &
484!^ & (Awrk(i,j,k,Nold)-Awrk(i-1,j,k,Nold))
485!^
486 tl_fx(i,j)=pmon_u(i,j)*0.5_r8*(kh(i-1,j)+kh(i,j))* &
487 & (tl_awrk(i,j,k,nold)-tl_awrk(i-1,j,k,nold))
488# ifdef MASKING
489!^ FX(i,j)=FX(i,j)*umask(i,j)
490!^
491 tl_fx(i,j)=tl_fx(i,j)*umask(i,j)
492# endif
493 END DO
494 END DO
495 DO j=jstr,jend+1
496 DO i=istr,iend
497!^ FE(i,j)=pnom_v(i,j)*0.5_r8*(Kh(i,j-1)+Kh(i,j))* &
498!^ & (Awrk(i,j,k,Nold)-Awrk(i,j-1,k,Nold))
499!^
500 tl_fe(i,j)=pnom_v(i,j)*0.5_r8*(kh(i,j-1)+kh(i,j))* &
501 & (tl_awrk(i,j,k,nold)-tl_awrk(i,j-1,k,nold))
502# ifdef MASKING
503!^ FE(i,j)=FE(i,j)*vmask(i,j)
504!^
505 tl_fe(i,j)=tl_fe(i,j)*vmask(i,j)
506# endif
507 END DO
508 END DO
509!
510! Time-step horizontal diffusion equation.
511!
512 DO j=jstr,jend
513 DO i=istr,iend
514!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
515!^ & Hfac(i,j)* &
516!^ & (FX(i+1,j)-FX(i,j)+ &
517!^ & FE(i,j+1)-FE(i,j))
518!^
519 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
520 & hfac(i,j)* &
521 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
522 & tl_fe(i,j+1)-tl_fe(i,j))
523 END DO
524 END DO
525 END DO
526# endif
527!
528! Apply boundary conditions.
529!
530!^ CALL dabc_r3d_tile (ng, tile, &
531!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
532!^ & Awrk(:,:,:,Nnew))
533!^
534 CALL dabc_r3d_tile (ng, tile, &
535 & lbi, ubi, lbj, ubj, lbk, ubk, &
536 & tl_awrk(:,:,:,nnew))
537# ifdef DISTRIBUTE
538!^ CALL mp_exchange3d (ng, tile, model, 1, &
539!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
540!^ & Nghost, &
541!^ & EWperiodic(ng), NSperiodic(ng), &
542!^ & Awrk(:,:,:,Nnew))
543!^
544 CALL mp_exchange3d (ng, tile, model, 1, &
545 & lbi, ubi, lbj, ubj, lbk, ubk, &
546 & nghost, &
547 & ewperiodic(ng), nsperiodic(ng), &
548 & tl_awrk(:,:,:,nnew))
549# endif
550!
551! Update integration indices.
552!
553 nsav=nold
554 nold=nnew
555 nnew=nsav
556 END DO
557
558# ifdef VCONVOLUTION
559# ifdef IMPLICIT_VCONV
560# ifdef SPLINES_VCONV
561!
562!-----------------------------------------------------------------------
563! Integrate vertical diffusion equation implicitly using parabolic
564! splines.
565!-----------------------------------------------------------------------
566!
567 DO step=1,nvsteps
568!
569! Use conservative, parabolic spline reconstruction of vertical
570! diffusion derivatives. Then, time step vertical diffusion term
571! implicitly.
572!
573 DO j=jstr,jend
574 cff1=1.0_r8/6.0_r8
575 DO k=1,n(ng)-1
576 DO i=istr,iend
577 fc(i,k)=cff1*hz(i,j,k )- &
578 & dtsizev*kv(i,j,k-1)*ohz(i,j,k )
579 cf(i,k)=cff1*hz(i,j,k+1)- &
580 & dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
581 END DO
582 END DO
583 DO i=istr,iend
584 cf(i,0)=0.0_r8
585!^ DC(i,0)=0.0_r8
586!^
587 tl_dc(i,0)=0.0_r8
588 END DO
589!
590! LU decomposition and forward substitution.
591!
592 cff1=1.0_r8/3.0_r8
593 DO k=1,n(ng)-1
594 DO i=istr,iend
595 bc(i,k)=cff1*(hz(i,j,k)+hz(i,j,k+1))+ &
596 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
597 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
598 cf(i,k)=cff*cf(i,k)
599!^ DC(i,k)=cff*(Awrk(i,j,k+1,Nold)-Awrk(i,j,k,Nold)- &
600!^ & FC(i,k)*DC(i,k-1))
601!^
602 tl_dc(i,k)=cff*(tl_awrk(i,j,k+1,nold)- &
603 & tl_awrk(i,j,k ,nold)- &
604 & fc(i,k)*tl_dc(i,k-1))
605 END DO
606 END DO
607!
608! Backward substitution.
609!
610 DO i=istr,iend
611!^ DC(i,N(ng))=0.0_r8
612!^
613 tl_dc(i,n(ng))=0.0_r8
614 END DO
615 DO k=n(ng)-1,1,-1
616 DO i=istr,iend
617!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
618!^
619 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
620 END DO
621 END DO
622!
623 DO k=1,n(ng)
624 DO i=istr,iend
625!^ DC(i,k)=DC(i,k)*Kv(i,j,k)
626!^
627 tl_dc(i,k)=tl_dc(i,k)*kv(i,j,k)
628!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
629!^ & DTsizeV*oHz(i,j,k)* &
630!^ & (DC(i,k)-DC(i,k-1))
631!^
632 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
633 & dtsizev*ohz(i,j,k)* &
634 & (tl_dc(i,k)-tl_dc(i,k-1))
635 END DO
636 END DO
637 END DO
638!
639! Update integration indices.
640!
641 nsav=nold
642 nold=nnew
643 nnew=nsav
644 END DO
645# else
646!
647!-----------------------------------------------------------------------
648! Integrate vertical diffusion equation implicitly.
649!-----------------------------------------------------------------------
650!
651 DO step=1,nvsteps
652!
653! Compute diagonal matrix coefficients BC and load right-hand-side
654! terms for the tangent diffusion equation into tl_DC.
655!
656 DO j=jstr,jend
657 DO k=1,n(ng)
658 DO i=istr,iend
659 bc(i,k)=hz(i,j,k)-fc(i,j,k)-fc(i,j,k-1)
660!^ DC(i,k)=Awrk(i,j,k,Nold)*Hz(i,j,k)
661!^
662 tl_dc(i,k)=tl_awrk(i,j,k,nold)*hz(i,j,k)
663 END DO
664 END DO
665!
666! Solve the tangent linear tridiagonal system.
667!
668 DO i=istr,iend
669 cff=1.0_r8/bc(i,1)
670 cf(i,1)=cff*fc(i,j,1)
671!^ DC(i,1)=cff*DC(i,1)
672!^
673 tl_dc(i,1)=cff*tl_dc(i,1)
674 END DO
675 DO k=2,n(ng)-1
676 DO i=istr,iend
677 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
678 cf(i,k)=cff*fc(i,j,k)
679!^ DC(i,k)=cff*(DC(i,k)-FC(i,j,k-1)*DC(i,k-1))
680!^
681 tl_dc(i,k)=cff*(tl_dc(i,k)-fc(i,j,k-1)*tl_dc(i,k-1))
682 END DO
683 END DO
684!
685! Compute new solution by back substitution.
686!
687 DO i=istr,iend
688!^ DC(i,N(ng))=(DC(i,N(ng))- &
689!^ & FC(i,j,N(ng)-1)*DC(i,N(ng)-1))/ &
690!^ & (BC(i,N(ng))-FC(i,j,N(ng)-1)*CF(i,N(ng)-1))
691!^
692 tl_dc(i,n(ng))=(tl_dc(i,n(ng))- &
693 & fc(i,j,n(ng)-1)*tl_dc(i,n(ng)-1))/ &
694 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
695!^ Awrk(i,j,N(ng),Nnew)=DC(i,N(ng))
696!^
697 tl_awrk(i,j,n(ng),nnew)=tl_dc(i,n(ng))
698# ifdef MASKING
699!^ Awrk(i,j,N(ng),Nnew)=Awrk(i,j,N(ng),Nnew)*rmask(i,j)
700!^
701 tl_awrk(i,j,n(ng),nnew)=tl_awrk(i,j,n(ng),nnew)*rmask(i,j)
702# endif
703 END DO
704 DO k=n(ng)-1,1,-1
705 DO i=istr,iend
706!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
707!^
708 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
709!^ Awrk(i,j,k,Nnew)=DC(i,k)
710!^
711 tl_awrk(i,j,k,nnew)=tl_dc(i,k)
712# ifdef MASKING
713!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nnew)*rmask(i,j)
714!^
715 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nnew)*rmask(i,j)
716# endif
717 END DO
718 END DO
719 END DO
720!
721! Update integration indices.
722!
723 nsav=nold
724 nold=nnew
725 nnew=nsav
726 END DO
727# endif
728# else
729!
730!-----------------------------------------------------------------------
731! Integrate vertical diffusion equation explicitly.
732!-----------------------------------------------------------------------
733!
734 DO step=1,nvsteps
735!
736! Compute vertical diffusive flux. Notice that "FC" is assumed to
737! be time invariant.
738!
739 DO j=jstr,jend
740 DO k=1,n(ng)-1
741 DO i=istr,iend
742!^ FS(i,k)=FC(i,j,k)*(Awrk(i,j,k+1,Nold)- &
743!^ & Awrk(i,j,k ,Nold))
744!^
745 tl_fs(i,k)=fc(i,j,k)*(tl_awrk(i,j,k+1,nold)- &
746 & tl_awrk(i,j,k ,nold))
747# ifdef MASKING
748!^ FS(i,k)=FS(i,k)*rmask(i,j)
749!^
750 tl_fs(i,k)=tl_fs(i,k)*rmask(i,j)
751# endif
752 END DO
753 END DO
754 DO i=istr,iend
755!^ FS(i,0)=0.0_r8
756!^
757 tl_fs(i,0)=0.0_r8
758!^ FS(i,N(ng))=0.0_r8
759!^
760 tl_fs(i,n(ng))=0.0_r8
761 END DO
762!
763! Time-step vertical diffusive term. Notice that "oHz" is assumed to
764! be time invariant.
765!
766 DO k=1,n(ng)
767 DO i=istr,iend
768!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
769!^ & oHz(i,j,k)*(FS(i,k )- &
770!^ & FS(i,k-1))
771!^
772 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
773 & ohz(i,j,k)*(tl_fs(i,k )- &
774 & tl_fs(i,k-1))
775 END DO
776 END DO
777 END DO
778!
779! Update integration indices.
780!
781 nsav=nold
782 nold=nnew
783 nnew=nsav
784 END DO
785# endif
786# endif
787!
788!-----------------------------------------------------------------------
789! Load convolved solution.
790!-----------------------------------------------------------------------
791!
792 DO k=1,n(ng)
793 DO j=jstr,jend
794 DO i=istr,iend
795!^ A(i,j,k)=Awrk(i,j,k,Nold)
796!^
797 tl_a(i,j,k)=tl_awrk(i,j,k,nold)
798 END DO
799 END DO
800 END DO
801!^ CALL dabc_r3d_tile (ng, tile, &
802!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
803!^ & A)
804!^
805 CALL dabc_r3d_tile (ng, tile, &
806 & lbi, ubi, lbj, ubj, lbk, ubk, &
807 & tl_a)
808# ifdef DISTRIBUTE
809!^ CALL mp_exchange3d (ng, tile, model, 1, &
810!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
811!^ & Nghost, &
812!^ & EWperiodic(ng), NSperiodic(ng), &
813!^ & A)
814!^
815 CALL mp_exchange3d (ng, tile, model, 1, &
816 & lbi, ubi, lbj, ubj, lbk, ubk, &
817 & nghost, &
818 & ewperiodic(ng), nsperiodic(ng), &
819 & tl_a)
820# endif
821
822 RETURN
823 END SUBROUTINE tl_conv_r3d_tile
824!
825!***********************************************************************
826 SUBROUTINE tl_conv_u3d_tile (ng, tile, model, &
827 & LBi, UBi, LBj, UBj, LBk, UBk, &
828 & IminS, ImaxS, JminS, JmaxS, &
829 & Nghost, NHsteps, NVsteps, &
830 & DTsizeH, DTsizeV, &
831 & Kh, Kv, &
832 & pm, pn, &
833# ifdef GEOPOTENTIAL_HCONV
834 & on_r, om_p, &
835# else
836 & pmon_r, pnom_p, &
837# endif
838# ifdef MASKING
839# ifdef GEOPOTENTIAL_HCONV
840 & pmask, rmask, umask, vmask, &
841# else
842 & umask, pmask, &
843# endif
844# endif
845 & Hz, z_r, &
846 & tl_A)
847!***********************************************************************
848!
849 USE mod_param
850 USE mod_scalars
851!
852 USE bc_3d_mod, ONLY: dabc_u3d_tile
853# ifdef DISTRIBUTE
854 USE mp_exchange_mod, ONLY : mp_exchange3d
855# endif
856!
857! Imported variable declarations.
858!
859 integer, intent(in) :: ng, tile, model
860 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
861 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
862 integer, intent(in) :: Nghost, NHsteps, NVsteps
863
864 real(r8), intent(in) :: DTsizeH, DTsizeV
865!
866# ifdef ASSUMED_SHAPE
867 real(r8), intent(in) :: pm(LBi:,LBj:)
868 real(r8), intent(in) :: pn(LBi:,LBj:)
869# ifdef GEOPOTENTIAL_HCONV
870 real(r8), intent(in) :: on_r(LBi:,LBj:)
871 real(r8), intent(in) :: om_p(LBi:,LBj:)
872# else
873 real(r8), intent(in) :: pmon_r(LBi:,LBj:)
874 real(r8), intent(in) :: pnom_p(LBi:,LBj:)
875# endif
876# ifdef MASKING
877# ifdef GEOPOTENTIAL_HCONV
878 real(r8), intent(in) :: pmask(LBi:,LBj:)
879 real(r8), intent(in) :: rmask(LBi:,LBj:)
880 real(r8), intent(in) :: umask(LBi:,LBj:)
881 real(r8), intent(in) :: vmask(LBi:,LBj:)
882# else
883 real(r8), intent(in) :: umask(LBi:,LBj:)
884 real(r8), intent(in) :: pmask(LBi:,LBj:)
885# endif
886# endif
887 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
888 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
889
890 real(r8), intent(in) :: Kh(LBi:,LBj:)
891 real(r8), intent(in) :: Kv(LBi:,LBj:,0:)
892
893 real(r8), intent(inout) :: tl_A(LBi:,LBj:,LBk:)
894# else
895 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
896 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
897# ifdef GEOPOTENTIAL_HCONV
898 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
899 real(r8), intent(in) :: om_p(LBi:UBi,LBj:UBj)
900# else
901 real(r8), intent(in) :: pmon_r(LBi:UBi,LBj:UBj)
902 real(r8), intent(in) :: pnom_p(LBi:UBi,LBj:UBj)
903# endif
904# ifdef MASKING
905# ifdef GEOPOTENTIAL_HCONV
906 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
907 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
908 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
909 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
910# else
911 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
912 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
913# endif
914# endif
915 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
916 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
917
918 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
919 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
920 real(r8), intent(inout) :: tl_A(LBi:UBi,LBj:UBj,LBk:UBk)
921# endif
922!
923! Local variable declarations.
924!
925 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
926
927 real(r8) :: cff, cff1, cff2, cff3, cff4
928
929 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: tl_Awrk
930
931 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
932 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
933 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
934# ifdef VCONVOLUTION
935# ifndef SPLINES_VCONV
936 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FC
937# endif
938# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
939 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
940# endif
941# if defined IMPLICIT_VCONV || defined SPLINES_VCONV
942 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
943 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
944# ifdef SPLINES_VCONV
945 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
946 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
947# endif
948 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
949# else
950 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FS
951# endif
952# endif
953# ifdef GEOPOTENTIAL_HCONV
954 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZdx
955 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZde
956
957 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
958 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
959 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FZ
960 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdz
961 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdx
962 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAde
963# endif
964
965# include "set_bounds.h"
966!
967!-----------------------------------------------------------------------
968! Space convolution of the diffusion equation for a 3D state variable
969! at U-points.
970!-----------------------------------------------------------------------
971!
972! Compute metrics factors. Notice that "z_r" and "Hz" are assumed to
973! be time invariant in the vertical convolution. Scratch array are
974! used for efficiency.
975!
976 cff=dtsizeh*0.25_r8
977 DO j=jstr-1,jend+1
978 DO i=istru-1,iend+1
979 hfac(i,j)=cff*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
980# ifdef VCONVOLUTION
981# ifndef SPLINES_VCONV
982 fc(i,j,n(ng))=0.0_r8
983 DO k=1,n(ng)-1
984# ifdef IMPLICIT_VCONV
985 fc(i,j,k)=-dtsizev*(kv(i-1,j,k)+kv(i,j,k))/ &
986 & (z_r(i-1,j,k+1)+z_r(i,j,k+1)- &
987 & z_r(i-1,j,k )-z_r(i,j,k ))
988# else
989 fc(i,j,k)=dtsizev*(kv(i-1,j,k)+kv(i,j,k))/ &
990 & (z_r(i-1,j,k+1)+z_r(i,j,k+1)- &
991 & z_r(i-1,j,k )-z_r(i,j,k ))
992# endif
993 END DO
994 fc(i,j,0)=0.0_r8
995# endif
996# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
997 DO k=1,n(ng)
998 ohz(i,j,k)=2.0_r8/(hz(i-1,j,k)+hz(i,j,k))
999 END DO
1000# endif
1001# endif
1002 END DO
1003 END DO
1004!
1005! Set integration indices and initial conditions.
1006!
1007 nold=1
1008 nnew=2
1009!^ CALL dabc_u3d_tile (ng, tile, &
1010!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1011!^ & A)
1012!^
1013 CALL dabc_u3d_tile (ng, tile, &
1014 & lbi, ubi, lbj, ubj, lbk, ubk, &
1015 & tl_a)
1016# ifdef DISTRIBUTE
1017!^ CALL mp_exchange3d (ng, tile, model, 1, &
1018!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1019!^ & Nghost, &
1020!^ & EWperiodic(ng), NSperiodic(ng), &
1021!^ & A)
1022!^
1023 CALL mp_exchange3d (ng, tile, model, 1, &
1024 & lbi, ubi, lbj, ubj, lbk, ubk, &
1025 & nghost, &
1026 & ewperiodic(ng), nsperiodic(ng), &
1027 & tl_a)
1028# endif
1029 DO k=1,n(ng)
1030 DO j=jstr-1,jend+1
1031 DO i=istru-1,iend+1
1032!^ Awrk(i,j,k,Nold)=A(i,j,k)
1033!^
1034 tl_awrk(i,j,k,nold)=tl_a(i,j,k)
1035 END DO
1036 END DO
1037 END DO
1038!
1039!-----------------------------------------------------------------------
1040! Integrate horizontal diffusion equation.
1041!-----------------------------------------------------------------------
1042!
1043 DO step=1,nhsteps
1044
1045# ifdef GEOPOTENTIAL_HCONV
1046!
1047! Diffusion along geopotential surfaces: Compute horizontal and
1048! vertical gradients. Notice the recursive blocking sequence. The
1049! vertical placement of the gradients is:
1050!
1051! dAdx,dAde(:,:,k1) k rho-points
1052! dAdx,dAde(:,:,k2) k+1 rho-points
1053! FZ,dAdz(:,:,k1) k-1/2 W-points
1054! FZ,dAdz(:,:,k2) k+1/2 W-points
1055!
1056 k2=1
1057 k_loop : DO k=0,n(ng)
1058 k1=k2
1059 k2=3-k1
1060 IF (k.lt.n(ng)) THEN
1061 DO j=jstr,jend
1062 DO i=istru-1,iend+1
1063 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
1064# ifdef MASKING
1065 cff=cff*umask(i,j)
1066# endif
1067 dzdx(i,j)=cff*(z_r(i ,j,k+1)- &
1068 & z_r(i-1,j,k+1))
1069 END DO
1070 END DO
1071!
1072 DO j=jstr,jend+1
1073 DO i=istru-1,iend
1074 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
1075# ifdef MASKING
1076 cff=cff*vmask(i,j)
1077# endif
1078 dzde(i,j)=cff*(z_r(i,j ,k+1)- &
1079 & z_r(i,j-1,k+1))
1080 END DO
1081 END DO
1082!
1083 DO j=jstr,jend
1084 DO i=istru-1,iend
1085# ifdef MASKING
1086!^ dAdx(i,j,k2)=pm(i,j)* &
1087!^ & (Awrk(i+1,j,k+1,Nold)*umask(i+1,j)- &
1088!^ & Awrk(i ,j,k+1,Nold)*umask(i ,j))
1089!^
1090 tl_dadx(i,j,k2)=pm(i,j)* &
1091 & (tl_awrk(i+1,j,k+1,nold)*umask(i+1,j)- &
1092 & tl_awrk(i ,j,k+1,nold)*umask(i ,j))
1093!^ dAdx(i,j,k2)=dAdx(i,j,k2)*rmask(i,j)
1094!^
1095 tl_dadx(i,j,k2)=tl_dadx(i,j,k2)*rmask(i,j)
1096# else
1097!^ dAdx(i,j,k2)=pm(i,j)*(Awrk(i+1,j,k+1,Nold)- &
1098!^ & Awrk(i ,j,k+1,Nold))
1099!^
1100 tl_dadx(i,j,k2)=pm(i,j)*(tl_awrk(i+1,j,k+1,nold)- &
1101 & tl_awrk(i ,j,k+1,nold))
1102# endif
1103 dzdx_r(i,j,k2)=0.5_r8*(dzdx(i ,j)+ &
1104 & dzdx(i+1,j))
1105 END DO
1106 END DO
1107!
1108 DO j=jstr,jend+1
1109 DO i=istru,iend
1110 cff=0.25_r8*(pn(i-1,j )+pn(i,j )+ &
1111 & pn(i-1,j-1)+pn(i,j-1))
1112# ifdef MASKING
1113!^ dAde(i,j,k2)=cff* &
1114!^ & (Awrk(i,j ,k+1,Nold)*umask(i,j )- &
1115!^ & Awrk(i,j-1,k+1,Nold)*umask(i,j-1))
1116!^
1117 tl_dade(i,j,k2)=cff* &
1118 & (tl_awrk(i,j ,k+1,nold)*umask(i,j )- &
1119 & tl_awrk(i,j-1,k+1,nold)*umask(i,j-1))
1120
1121!^ dAde(i,j,k2)=dAde(i,j,k2)*pmask(i,j)
1122!^
1123 tl_dade(i,j,k2)=tl_dade(i,j,k2)*pmask(i,j)
1124# else
1125!^ dAde(i,j,k2)=cff*(Awrk(i,j ,k+1,Nold)- &
1126!^ & Awrk(i,j-1,k+1,Nold))
1127!^
1128 tl_dade(i,j,k2)=cff*(tl_awrk(i,j ,k+1,nold)- &
1129 & tl_awrk(i,j-1,k+1,nold))
1130# endif
1131 dzde_p(i,j,k2)=0.5_r8*(dzde(i-1,j)+ &
1132 & dzde(i ,j))
1133 END DO
1134 END DO
1135 END IF
1136!
1137 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1138 DO j=jstr-1,jend+1
1139 DO i=istru-1,iend+1
1140!^ dAdz(i,j,k2)=0.0_r8
1141!^
1142 tl_dadz(i,j,k2)=0.0_r8
1143!^ FZ(i,j,k2)=0.0_r8
1144!^
1145 tl_fz(i,j,k2)=0.0_r8
1146 END DO
1147 END DO
1148 ELSE
1149 DO j=jstr-1,jend+1
1150 DO i=istru-1,iend+1
1151 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1152!^ dAdz(i,j,k2)=cff*(Awrk(i,j,k+1,Nold)- &
1153!^ & Awrk(i,j,k ,Nold))
1154!^
1155 tl_dadz(i,j,k2)=cff*(tl_awrk(i,j,k+1,nold)- &
1156 & tl_awrk(i,j,k ,nold))
1157# ifdef MASKING
1158!^ dAdz(i,j,k2)=dAdz(i,j,k2)*umask(i,j)
1159!^
1160 tl_dadz(i,j,k2)=tl_dadz(i,j,k2)*umask(i,j)
1161# endif
1162 END DO
1163 END DO
1164 END IF
1165!
1166! Compute components of the rotated A flux (A m3/s) along geopotential
1167! surfaces.
1168!
1169 IF (k.gt.0) THEN
1170 DO j=jstr,jend
1171 DO i=istru-1,iend
1172 cff=kh(i,j)*on_r(i,j)
1173 cff1=min(dzdx_r(i,j,k1),0.0_r8)
1174 cff2=max(dzdx_r(i,j,k1),0.0_r8)
1175!^ FX(i,j)=cff* &
1176!^ & Hz(i,j,k)* &
1177!^ & (dAdx(i,j,k1)- &
1178!^ & 0.5_r8*(cff1*(dAdz(i ,j,k1)+ &
1179!^ & dAdz(i+1,j,k2))+ &
1180!^ & cff2*(dAdz(i ,j,k2)+ &
1181!^ & dAdz(i+1,j,k1))))
1182!^
1183 tl_fx(i,j)=cff* &
1184 & hz(i,j,k)* &
1185 & (tl_dadx(i,j,k1)- &
1186 & 0.5_r8*(cff1*(tl_dadz(i ,j,k1)+ &
1187 & tl_dadz(i+1,j,k2))+ &
1188 & cff2*(tl_dadz(i ,j,k2)+ &
1189 & tl_dadz(i+1,j,k1))))
1190 END DO
1191 END DO
1192 DO j=jstr,jend+1
1193 DO i=istru,iend
1194 cff=0.0625_r8*(kh(i-1,j-1)+kh(i-1,j)+ &
1195 & kh(i ,j-1)+kh(i ,j))*om_p(i,j)
1196 cff1=min(dzde_p(i,j,k1),0.0_r8)
1197 cff2=max(dzde_p(i,j,k1),0.0_r8)
1198!^ FE(i,j)=cff* &
1199!^ & (Hz(i-1,j-1,k)+Hz(i-1,j,k)+ &
1200!^ & Hz(i ,j-1,k)+Hz(i ,j,k))* &
1201!^ & (dAde(i,j,k1)- &
1202!^ & 0.5_r8*(cff1*(dAdz(i,j-1,k1)+ &
1203!^ & dAdz(i,j ,k2))+ &
1204!^ & cff2*(dAdz(i,j-1,k2)+ &
1205!^ & dAdz(i,j ,k1))))
1206!^
1207 tl_fe(i,j)=cff* &
1208 & (hz(i-1,j-1,k)+hz(i-1,j,k)+ &
1209 & hz(i ,j-1,k)+hz(i ,j,k))* &
1210 & (tl_dade(i,j,k1)- &
1211 & 0.5_r8*(cff1*(tl_dadz(i,j-1,k1)+ &
1212 & tl_dadz(i,j ,k2))+ &
1213 & cff2*(tl_dadz(i,j-1,k2)+ &
1214 & tl_dadz(i,j ,k1))))
1215 END DO
1216 END DO
1217 IF (k.lt.n(ng)) THEN
1218 DO j=jstr,jend
1219 DO i=istru,iend
1220 cff=0.25_r8*(kh(i-1,j)+kh(i,j))
1221 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
1222 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
1223 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
1224 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
1225!^ FZ(i,j,k2)=cff* &
1226!^ & (cff1*(cff1*dAdz(i,j,k2)-dAdx(i-1,j,k1))+ &
1227!^ & cff2*(cff2*dAdz(i,j,k2)-dAdx(i ,j,k2))+ &
1228!^ & cff3*(cff3*dAdz(i,j,k2)-dAdx(i-1,j,k2))+ &
1229!^ & cff4*(cff4*dAdz(i,j,k2)-dAdx(i ,j,k1)))
1230!^
1231 tl_fz(i,j,k2)=cff* &
1232 & (cff1*(cff1*tl_dadz(i,j,k2)- &
1233 & tl_dadx(i-1,j,k1))+ &
1234 & cff2*(cff2*tl_dadz(i,j,k2)- &
1235 & tl_dadx(i ,j,k2))+ &
1236 & cff3*(cff3*tl_dadz(i,j,k2)- &
1237 & tl_dadx(i-1,j,k2))+ &
1238 & cff4*(cff4*tl_dadz(i,j,k2)- &
1239 & tl_dadx(i ,j,k1)))
1240 cff1=min(dzde_p(i,j ,k1),0.0_r8)
1241 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
1242 cff3=max(dzde_p(i,j ,k2),0.0_r8)
1243 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
1244!^ FZ(i,j,k2)=FZ(i,j,k2)+ &
1245!^ & cff* &
1246!^ & (cff1*(cff1*dAdz(i,j,k2)-dAde(i,j ,k1))+ &
1247!^ & cff2*(cff2*dAdz(i,j,k2)-dAde(i,j+1,k2))+ &
1248!^ & cff3*(cff3*dAdz(i,j,k2)-dAde(i,j ,k2))+ &
1249!^ & cff4*(cff4*dAdz(i,j,k2)-dAde(i,j+1,k1)))
1250!^
1251 tl_fz(i,j,k2)=tl_fz(i,j,k2)+ &
1252 & cff* &
1253 & (cff1*(cff1*tl_dadz(i,j,k2)- &
1254 & tl_dade(i,j ,k1))+ &
1255 & cff2*(cff2*tl_dadz(i,j,k2)- &
1256 & tl_dade(i,j+1,k2))+ &
1257 & cff3*(cff3*tl_dadz(i,j,k2)- &
1258 & tl_dade(i,j ,k2))+ &
1259 & cff4*(cff4*tl_dadz(i,j,k2)- &
1260 & tl_dade(i,j+1,k1)))
1261 END DO
1262 END DO
1263 END IF
1264!
1265! Time-step harmonic, geopotential diffusion term (m Tunits).
1266!
1267 DO j=jstr,jend
1268 DO i=istru,iend
1269!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
1270!^ & Hfac(i,j)* &
1271!^ & (FX(i,j )-FX(i-1,j)+ &
1272!^ & FE(i,j+1)-FE(i ,j))+ &
1273!^ & DTsizeH* &
1274!^ & (FZ(i,j,k2)-FZ(i,j,k1))
1275!^
1276 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
1277 & hfac(i,j)* &
1278 & (tl_fx(i,j )-tl_fx(i-1,j)+ &
1279 & tl_fe(i,j+1)-tl_fe(i ,j))+ &
1280 & dtsizeh* &
1281 & (tl_fz(i,j,k2)-tl_fz(i,j,k1))
1282 END DO
1283 END DO
1284 END IF
1285 END DO k_loop
1286
1287# else
1288
1289!
1290! Diffusion along S-coordinates: compute XI- and ETA-components of
1291! diffusive flux.
1292!
1293 DO k=1,n(ng)
1294 DO j=jstr,jend
1295 DO i=istru-1,iend
1296!^ FX(i,j)=pmon_r(i,j)*Kh(i,j)* &
1297!^ & (Awrk(i+1,j,k,Nold)-Awrk(i,j,k,Nold))
1298!^
1299 tl_fx(i,j)=pmon_r(i,j)*kh(i,j)* &
1300 & (tl_awrk(i+1,j,k,nold)-tl_awrk(i,j,k,nold))
1301 END DO
1302 END DO
1303 DO j=jstr,jend+1
1304 DO i=istru,iend
1305!^ FE(i,j)=pnom_p(i,j)*0.25_r8*(Kh(i-1,j )+Kh(i,j )+ &
1306!^ & Kh(i-1,j-1)+Kh(i,j-1))* &
1307!^ & (Awrk(i,j,k,Nold)-Awrk(i,j-1,k,Nold))
1308!^
1309 tl_fe(i,j)=pnom_p(i,j)*0.25_r8*(kh(i-1,j )+kh(i,j )+ &
1310 & kh(i-1,j-1)+kh(i,j-1))* &
1311 & (tl_awrk(i,j,k,nold)-tl_awrk(i,j-1,k,nold))
1312# ifdef MASKING
1313!^ FE(i,j)=FE(i,j)*pmask(i,j)
1314!^
1315 tl_fe(i,j)=tl_fe(i,j)*pmask(i,j)
1316# endif
1317 END DO
1318 END DO
1319!
1320! Time-step horizontal diffusion equation.
1321!
1322 DO j=jstr,jend
1323 DO i=istru,iend
1324!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
1325!^ & Hfac(i,j)* &
1326!^ & (FX(i,j)-FX(i-1,j)+ &
1327!^ & FE(i,j+1)-FE(i,j))
1328!^
1329 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
1330 & hfac(i,j)* &
1331 & (tl_fx(i,j)-tl_fx(i-1,j)+ &
1332 & tl_fe(i,j+1)-tl_fe(i,j))
1333 END DO
1334 END DO
1335 END DO
1336# endif
1337!
1338! Apply boundary conditions.
1339!
1340!^ CALL dabc_u3d_tile (ng, tile, &
1341!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1342!^ & Awrk(:,:,:,Nnew))
1343!^
1344 CALL dabc_u3d_tile (ng, tile, &
1345 & lbi, ubi, lbj, ubj, lbk, ubk, &
1346 & tl_awrk(:,:,:,nnew))
1347# ifdef DISTRIBUTE
1348!^ CALL mp_exchange3d (ng, tile, model, 1, &
1349!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1350!^ & Nghost, &
1351!^ & EWperiodic(ng), NSperiodic(ng), &
1352!^ & Awrk(:,:,:,Nnew))
1353!^
1354 CALL mp_exchange3d (ng, tile, model, 1, &
1355 & lbi, ubi, lbj, ubj, lbk, ubk, &
1356 & nghost, &
1357 & ewperiodic(ng), nsperiodic(ng), &
1358 & tl_awrk(:,:,:,nnew))
1359# endif
1360!
1361! Update integration indices.
1362!
1363 nsav=nold
1364 nold=nnew
1365 nnew=nsav
1366 END DO
1367
1368# ifdef VCONVOLUTION
1369# ifdef IMPLICIT_VCONV
1370# ifdef SPLINES_VCONV
1371!
1372!-----------------------------------------------------------------------
1373! Integrate vertical diffusion equation implicitly using parabolic
1374! splines.
1375!-----------------------------------------------------------------------
1376!
1377 DO step=1,nvsteps
1378!
1379! Use conservative, parabolic spline reconstruction of vertical
1380! diffusion derivatives. Then, time step vertical diffusion term
1381! implicitly.
1382!
1383 DO j=jstr,jend
1384 DO k=1,n(ng)
1385 DO i=istru,iend
1386 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
1387 & hz(i ,j,k))
1388 END DO
1389 END DO
1390 cff1=1.0_r8/6.0_r8
1391 DO k=1,n(ng)-1
1392 DO i=istru,iend
1393 fc(i,k)=cff1*hzk(i,k )-dtsizev*kv(i,j,k-1)*ohz(i,j,k )
1394 cf(i,k)=cff1*hzk(i,k+1)-dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
1395 END DO
1396 END DO
1397 DO i=istru,iend
1398 cf(i,0)=0.0_r8
1399!^ DC(i,0)=0.0_r8
1400!^
1401 tl_dc(i,0)=0.0_r8
1402 END DO
1403!
1404! LU decomposition and forward substitution.
1405!
1406 cff1=1.0_r8/3.0_r8
1407 DO k=1,n(ng)-1
1408 DO i=istru,iend
1409 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1410 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
1411 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1412 cf(i,k)=cff*cf(i,k)
1413!^ DC(i,k)=cff*(Awrk(i,j,k+1,Nold)-Awrk(i,j,k,Nold)- &
1414!^ & FC(i,k)*DC(i,k-1))
1415!^
1416 tl_dc(i,k)=cff*(tl_awrk(i,j,k+1,nold)- &
1417 & tl_awrk(i,j,k ,nold)- &
1418 & fc(i,k)*tl_dc(i,k-1))
1419 END DO
1420 END DO
1421!
1422! Backward substitution.
1423!
1424 DO i=istru,iend
1425!^ DC(i,N(ng))=0.0_r8
1426!^
1427 tl_dc(i,n(ng))=0.0_r8
1428 END DO
1429 DO k=n(ng)-1,1,-1
1430 DO i=istru,iend
1431!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
1432!^
1433 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
1434 END DO
1435 END DO
1436!
1437 DO k=1,n(ng)
1438 DO i=istru,iend
1439!^ DC(i,k)=DC(i,k)*Kv(i,j,k)
1440!^
1441 tl_dc(i,k)=tl_dc(i,k)*kv(i,j,k)
1442!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
1443!^ & DTsizeV*oHz(i,j,k)* &
1444!^ & (DC(i,k)-DC(i,k-1))
1445!^
1446 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
1447 & dtsizev*ohz(i,j,k)* &
1448 & (tl_dc(i,k)-tl_dc(i,k-1))
1449 END DO
1450 END DO
1451 END DO
1452!
1453! Update integration indices.
1454!
1455 nsav=nold
1456 nold=nnew
1457 nnew=nsav
1458 END DO
1459# else
1460!
1461!-----------------------------------------------------------------------
1462! Integrate vertical diffusion equation implicitly.
1463!-----------------------------------------------------------------------
1464!
1465 DO step=1,nvsteps
1466!
1467! Compute diagonal matrix coefficients BC and load right-hand-side
1468! terms for the tangent linear diffusion equation into tl_DC.
1469!
1470 DO j=jstr,jend
1471 DO k=1,n(ng)
1472 DO i=istru,iend
1473 cff=0.5_r8*(hz(i-1,j,k)+hz(i,j,k))
1474 bc(i,k)=cff-fc(i,j,k)-fc(i,j,k-1)
1475!^ DC(i,k)=Awrk(i,j,k,Nold)*cff
1476!^
1477 tl_dc(i,k)=tl_awrk(i,j,k,nold)*cff
1478 END DO
1479 END DO
1480!
1481! Solve the tridiagonal system.
1482!
1483 DO i=istru,iend
1484 cff=1.0_r8/bc(i,1)
1485 cf(i,1)=cff*fc(i,j,1)
1486!^ DC(i,1)=cff*DC(i,1)
1487!^
1488 tl_dc(i,1)=cff*tl_dc(i,1)
1489 END DO
1490 DO k=2,n(ng)-1
1491 DO i=istru,iend
1492 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
1493 cf(i,k)=cff*fc(i,j,k)
1494!^ DC(i,k)=cff*(DC(i,k)-FC(i,j,k-1)*DC(i,k-1))
1495!^
1496 tl_dc(i,k)=cff*(tl_dc(i,k)-fc(i,j,k-1)*tl_dc(i,k-1))
1497 END DO
1498 END DO
1499!
1500! Compute new solution by back substitution.
1501!
1502 DO i=istru,iend
1503!^ DC(i,N(ng))=(DC(i,N(ng))- &
1504!^ & FC(i,j,N(ng)-1)*DC(i,N(ng)-1))/ &
1505!^ & (BC(i,N(ng))-FC(i,j,N(ng)-1)*CF(i,N(ng)-1))
1506!^
1507 tl_dc(i,n(ng))=(tl_dc(i,n(ng))- &
1508 & fc(i,j,n(ng)-1)*tl_dc(i,n(ng)-1))/ &
1509 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
1510!^ Awrk(i,j,N(ng),Nnew)=DC(i,N(ng))
1511!^
1512 tl_awrk(i,j,n(ng),nnew)=tl_dc(i,n(ng))
1513# ifdef MASKING
1514!^ Awrk(i,j,N(ng),Nnew)=Awrk(i,j,N(ng),Nnew)*umask(i,j)
1515!^
1516 tl_awrk(i,j,n(ng),nnew)=tl_awrk(i,j,n(ng),nnew)*umask(i,j)
1517# endif
1518 END DO
1519 DO k=n(ng)-1,1,-1
1520 DO i=istru,iend
1521!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
1522!^
1523 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
1524!^ Awrk(i,j,k,Nnew)=DC(i,k)
1525!^
1526 tl_awrk(i,j,k,nnew)=tl_dc(i,k)
1527# ifdef MASKING
1528!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nnew)*umask(i,j)
1529!^
1530 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nnew)*umask(i,j)
1531# endif
1532 END DO
1533 END DO
1534 END DO
1535!
1536! Update integration indices.
1537!
1538 nsav=nold
1539 nold=nnew
1540 nnew=nsav
1541 END DO
1542# endif
1543# else
1544!
1545!-----------------------------------------------------------------------
1546! Integrate vertical diffusion equation explicitly.
1547!-----------------------------------------------------------------------
1548!
1549 DO step=1,nvsteps
1550!
1551! Compute vertical diffusive flux. Notice that "FC" is assumed to
1552! be time invariant.
1553!
1554 DO j=jstr,jend
1555 DO k=1,n(ng)-1
1556 DO i=istru,iend
1557!^ FS(i,k)=FC(i,j,k)*(Awrk(i,j,k+1,Nold)- &
1558!^ & Awrk(i,j,k ,Nold))
1559!^
1560 tl_fs(i,k)=fc(i,j,k)*(tl_awrk(i,j,k+1,nold)- &
1561 & tl_awrk(i,j,k ,nold))
1562# ifdef MASKING
1563!^ FS(i,k)=FS(i,k)*umask(i,j)
1564!^
1565 tl_fs(i,k)=tl_fs(i,k)*umask(i,j)
1566# endif
1567 END DO
1568 END DO
1569 DO i=istru,iend
1570!^ FS(i,0)=0.0_r8
1571!^
1572 tl_fs(i,0)=0.0_r8
1573!^ FS(i,N(ng))=0.0_r8
1574!^
1575 tl_fs(i,n(ng))=0.0_r8
1576 END DO
1577!
1578! Time-step vertical diffusive term. Notice that "oHz" is assumed to
1579! be time invariant.
1580!
1581 DO k=1,n(ng)
1582 DO i=istru,iend
1583!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
1584!^ & oHz(i,j,k)*(FS(i,k )- &
1585!^ & FS(i,k-1))
1586!^
1587 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
1588 & ohz(i,j,k)*(tl_fs(i,k )- &
1589 & tl_fs(i,k-1))
1590 END DO
1591 END DO
1592 END DO
1593!
1594! Update integration indices.
1595!
1596 nsav=nold
1597 nold=nnew
1598 nnew=nsav
1599 END DO
1600# endif
1601# endif
1602!
1603!-----------------------------------------------------------------------
1604! Load convolved solution.
1605!-----------------------------------------------------------------------
1606!
1607 DO k=1,n(ng)
1608 DO j=jstr,jend
1609 DO i=istru,iend
1610!^ A(i,j,k)=Awrk(i,j,k,Nold)
1611!^
1612 tl_a(i,j,k)=tl_awrk(i,j,k,nold)
1613 END DO
1614 END DO
1615 END DO
1616!^ CALL dabc_u3d_tile (ng, tile, &
1617!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1618!^ & A)
1619!^
1620 CALL dabc_u3d_tile (ng, tile, &
1621 & lbi, ubi, lbj, ubj, lbk, ubk, &
1622 & tl_a)
1623# ifdef DISTRIBUTE
1624!^ CALL mp_exchange3d (ng, tile, model, 1, &
1625!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1626!^ & Nghost, &
1627!^ & EWperiodic(ng), NSperiodic(ng), &
1628!^ & A)
1629!^
1630 CALL mp_exchange3d (ng, tile, model, 1, &
1631 & lbi, ubi, lbj, ubj, lbk, ubk, &
1632 & nghost, &
1633 & ewperiodic(ng), nsperiodic(ng), &
1634 & tl_a)
1635# endif
1636
1637 RETURN
1638 END SUBROUTINE tl_conv_u3d_tile
1639!
1640!***********************************************************************
1641 SUBROUTINE tl_conv_v3d_tile (ng, tile, model, &
1642 & LBi, UBi, LBj, UBj, LBk, UBk, &
1643 & IminS, ImaxS, JminS, JmaxS, &
1644 & Nghost, NHsteps, NVsteps, &
1645 & DTsizeH, DTsizeV, &
1646 & Kh, Kv, &
1647 & pm, pn, &
1648# ifdef GEOPOTENTIAL_HCONV
1649 & on_p, om_r, &
1650# else
1651 & pmon_p, pnom_r, &
1652# endif
1653# ifdef MASKING
1654# ifdef GEOPOTENTIAL_HCONV
1655 & pmask, rmask, umask, vmask, &
1656# else
1657 & vmask, pmask, &
1658# endif
1659# endif
1660 & Hz, z_r, &
1661 & tl_A)
1662!***********************************************************************
1663!
1664 USE mod_param
1665 USE mod_scalars
1666!
1667 USE bc_3d_mod, ONLY: dabc_v3d_tile
1668# ifdef DISTRIBUTE
1669 USE mp_exchange_mod, ONLY : mp_exchange3d
1670# endif
1671!
1672! Imported variable declarations.
1673!
1674 integer, intent(in) :: ng, tile, model
1675 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
1676 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1677 integer, intent(in) :: Nghost, NHsteps, NVsteps
1678
1679 real(r8), intent(in) :: DTsizeH, DTsizeV
1680!
1681# ifdef ASSUMED_SHAPE
1682 real(r8), intent(in) :: pm(LBi:,LBj:)
1683 real(r8), intent(in) :: pn(LBi:,LBj:)
1684# ifdef GEOPOTENTIAL_HCONV
1685 real(r8), intent(in) :: on_p(LBi:,LBj:)
1686 real(r8), intent(in) :: om_r(LBi:,LBj:)
1687# else
1688 real(r8), intent(in) :: pmon_p(LBi:,LBj:)
1689 real(r8), intent(in) :: pnom_r(LBi:,LBj:)
1690# endif
1691# ifdef MASKING
1692# ifdef GEOPOTENTIAL_HCONV
1693 real(r8), intent(in) :: pmask(LBi:,LBj:)
1694 real(r8), intent(in) :: rmask(LBi:,LBj:)
1695 real(r8), intent(in) :: umask(LBi:,LBj:)
1696 real(r8), intent(in) :: vmask(LBi:,LBj:)
1697# else
1698 real(r8), intent(in) :: vmask(LBi:,LBj:)
1699 real(r8), intent(in) :: pmask(LBi:,LBj:)
1700# endif
1701# endif
1702 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
1703 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
1704
1705 real(r8), intent(in) :: Kh(LBi:,LBj:)
1706 real(r8), intent(in) :: Kv(LBi:,LBj:,0:)
1707
1708 real(r8), intent(inout) :: tl_A(LBi:,LBj:,LBk:)
1709# else
1710 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
1711 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
1712# ifdef GEOPOTENTIAL_HCONV
1713 real(r8), intent(in) :: on_p(LBi:UBi,LBj:UBj)
1714 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
1715# else
1716 real(r8), intent(in) :: pmon_p(LBi:UBi,LBj:UBj)
1717 real(r8), intent(in) :: pnom_r(LBi:UBi,LBj:UBj)
1718# endif
1719# ifdef MASKING
1720# ifdef GEOPOTENTIAL_HCONV
1721 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
1722 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1723 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1724 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1725# else
1726 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1727 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
1728# endif
1729# endif
1730 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
1731 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
1732
1733 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
1734 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
1735
1736 real(r8), intent(inout) :: tl_A(LBi:UBi,LBj:UBj,LBk:UBk)
1737# endif
1738!
1739! Local variable declarations.
1740!
1741 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
1742
1743 real(r8) :: cff, cff1, cff2, cff3, cff4
1744
1745 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: tl_Awrk
1746
1747 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
1748 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
1749 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
1750# ifdef VCONVOLUTION
1751# ifndef SPLINES_VCONV
1752 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FC
1753# endif
1754# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
1755 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
1756# endif
1757# if defined IMPLICIT_VCONV || defined SPLINES_VCONV
1758 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
1759 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
1760 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
1761# ifdef SPLINES_VCONV
1762 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
1763 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
1764# endif
1765# else
1766 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FS
1767# endif
1768# endif
1769# ifdef GEOPOTENTIAL_HCONV
1770 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZdx
1771 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZde
1772
1773 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
1774 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
1775 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FZ
1776 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdz
1777 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAdx
1778 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dAde
1779# endif
1780
1781# include "set_bounds.h"
1782!
1783!-----------------------------------------------------------------------
1784! Space convolution of the diffusion equation for a 3D state variable
1785! at V-points.
1786!-----------------------------------------------------------------------
1787!
1788! Compute metrics factors. Notice that "z_r" and "Hz" are assumed to
1789! be time invariant in the vertical convolution. Scratch array are
1790! used for efficiency.
1791!
1792 cff=dtsizeh*0.25_r8
1793 DO j=jstrv-1,jend+1
1794 DO i=istr-1,iend+1
1795 hfac(i,j)=cff*(pm(i,j-1)+pm(i,j))*(pn(i,j-1)+pn(i,j))
1796# ifdef VCONVOLUTION
1797# ifndef SPLINES_VCONV
1798 fc(i,j,n(ng))=0.0_r8
1799 DO k=1,n(ng)-1
1800# ifdef IMPLICIT_VCONV
1801 fc(i,j,k)=-dtsizev*(kv(i,j-1,k)+kv(i,j,k))/ &
1802 & (z_r(i,j-1,k+1)+z_r(i,j,k+1)- &
1803 & z_r(i,j-1,k )-z_r(i,j,k ))
1804# else
1805 fc(i,j,k)=dtsizev*(kv(i,j-1,k)+kv(i,j,k))/ &
1806 & (z_r(i,j-1,k+1)+z_r(i,j,k+1)- &
1807 & z_r(i,j-1,k )-z_r(i,j,k ))
1808# endif
1809 END DO
1810 fc(i,j,0)=0.0_r8
1811# endif
1812# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
1813 DO k=1,n(ng)
1814 ohz(i,j,k)=2.0_r8/(hz(i,j-1,k)+hz(i,j,k))
1815 END DO
1816# endif
1817# endif
1818 END DO
1819 END DO
1820!
1821! Set integration indices and initial conditions.
1822!
1823 nold=1
1824 nnew=2
1825!^ CALL dabc_v3d_tile (ng, tile, &
1826!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1827!^ & A)
1828!^
1829 CALL dabc_v3d_tile (ng, tile, &
1830 & lbi, ubi, lbj, ubj, lbk, ubk, &
1831 & tl_a)
1832# ifdef DISTRIBUTE
1833!^ CALL mp_exchange3d (ng, tile, model, 1, &
1834!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
1835!^ & Nghost, &
1836!^ & EWperiodic(ng), NSperiodic(ng), &
1837!^ & A)
1838!^
1839 CALL mp_exchange3d (ng, tile, model, 1, &
1840 & lbi, ubi, lbj, ubj, lbk, ubk, &
1841 & nghost, &
1842 & ewperiodic(ng), nsperiodic(ng), &
1843 & tl_a)
1844# endif
1845 DO k=1,n(ng)
1846 DO j=jstrv-1,jend+1
1847 DO i=istr-1,iend+1
1848!^ Awrk(i,j,k,Nold)=A(i,j,k)
1849!^
1850 tl_awrk(i,j,k,nold)=tl_a(i,j,k)
1851 END DO
1852 END DO
1853 END DO
1854!
1855!-----------------------------------------------------------------------
1856! Integrate horizontal diffusion equation.
1857!-----------------------------------------------------------------------
1858!
1859 DO step=1,nhsteps
1860
1861# ifdef GEOPOTENTIAL_HCONV
1862!
1863! Diffusion along geopotential surfaces: Compute horizontal and
1864! vertical gradients. Notice the recursive blocking sequence. The
1865! vertical placement of the gradients is:
1866!
1867! dAdx,dAde(:,:,k1) k rho-points
1868! dAdx,dAde(:,:,k2) k+1 rho-points
1869! FZ,dAdz(:,:,k1) k-1/2 W-points
1870! FZ,dAdz(:,:,k2) k+1/2 W-points
1871!
1872 k2=1
1873 k_loop : DO k=0,n(ng)
1874 k1=k2
1875 k2=3-k1
1876 IF (k.lt.n(ng)) THEN
1877 DO j=jstrv-1,jend
1878 DO i=istr,iend+1
1879 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
1880# ifdef MASKING
1881 cff=cff*umask(i,j)
1882# endif
1883 dzdx(i,j)=cff*(z_r(i ,j,k+1)- &
1884 & z_r(i-1,j,k+1))
1885 END DO
1886 END DO
1887!
1888 DO j=jstrv-1,jend+1
1889 DO i=istr,iend
1890 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
1891# ifdef MASKING
1892 cff=cff*vmask(i,j)
1893# endif
1894 dzde(i,j)=cff*(z_r(i,j ,k+1)- &
1895 & z_r(i,j-1,k+1))
1896 END DO
1897 END DO
1898!
1899 DO j=jstrv,jend
1900 DO i=istr,iend+1
1901 cff=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+ &
1902 & pm(i ,j-1)+pm(i ,j))
1903# ifdef MASKING
1904!^ dAdx(i,j,k2)=cff* &
1905!^ & (Awrk(i ,j,k+1,Nold)*vmask(i ,j)- &
1906!^ & Awrk(i-1,j,k+1,Nold)*vmask(i-1,j))
1907!^
1908 tl_dadx(i,j,k2)=cff* &
1909 & (tl_awrk(i ,j,k+1,nold)*vmask(i ,j)- &
1910 & tl_awrk(i-1,j,k+1,nold)*vmask(i-1,j))
1911!^ dAdx(i,j,k2)=dAdx(i,j,k2)*pmask(i,j)
1912!^
1913 tl_dadx(i,j,k2)=tl_dadx(i,j,k2)*pmask(i,j)
1914# else
1915!^ dAdx(i,j,k2)=cff*(Awrk(i ,j,k+1,Nold)- &
1916!^ & Awrk(i-1,j,k+1,Nold))
1917!^
1918 tl_dadx(i,j,k2)=cff*(tl_awrk(i ,j,k+1,nold)- &
1919 & tl_awrk(i-1,j,k+1,nold))
1920# endif
1921 dzdx_p(i,j,k2)=0.5_r8*(dzdx(i,j-1)+ &
1922 & dzdx(i,j ))
1923 END DO
1924 END DO
1925!
1926 DO j=jstrv-1,jend
1927 DO i=istr,iend
1928# ifdef MASKING
1929!^ dAde(i,j,k2)=pn(i,j)* &
1930!^ & (Awrk(i,j+1,k+1,Nold)*vmask(i,j+1)- &
1931!^ & Awrk(i,j ,k+1,Nold)*vmask(i,j ))
1932!^
1933 tl_dade(i,j,k2)=pn(i,j)* &
1934 & (tl_awrk(i,j+1,k+1,nold)*vmask(i,j+1)- &
1935 & tl_awrk(i,j ,k+1,nold)*vmask(i,j ))
1936!^ dAde(i,j,k2)=dAde(i,j,k2)*rmask(i,j)
1937!^
1938 tl_dade(i,j,k2)=tl_dade(i,j,k2)*rmask(i,j)
1939# else
1940!^ dAde(i,j,k2)=pn(i,j)*(Awrk(i,j+1,k+1,Nold)- &
1941!^ & Awrk(i,j ,k+1,Nold))
1942!^
1943 tl_dade(i,j,k2)=pn(i,j)*(tl_awrk(i,j+1,k+1,nold)- &
1944 & tl_awrk(i,j ,k+1,nold))
1945# endif
1946 dzde_r(i,j,k2)=0.5_r8*(dzde(i,j )+ &
1947 & dzde(i,j+1))
1948 END DO
1949 END DO
1950 END IF
1951!
1952 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1953 DO j=jstrv-1,jend+1
1954 DO i=istr-1,iend+1
1955!^ dAdz(i,j,k2)=0.0_r8
1956!^
1957 tl_dadz(i,j,k2)=0.0_r8
1958!^ FZ(i,j,k2)=0.0_r8
1959!^
1960 tl_fz(i,j,k2)=0.0_r8
1961 END DO
1962 END DO
1963 ELSE
1964 DO j=jstrv-1,jend+1
1965 DO i=istr-1,iend+1
1966 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1967!^ dAdz(i,j,k2)=cff*(Awrk(i,j,k+1,Nold)- &
1968!^ & Awrk(i,j,k ,Nold))
1969!^
1970 tl_dadz(i,j,k2)=cff*(tl_awrk(i,j,k+1,nold)- &
1971 & tl_awrk(i,j,k ,nold))
1972# ifdef MASKING
1973!^ dAdz(i,j,k2)=dAdz(i,j,k2)*vmask(i,j)
1974!^
1975 tl_dadz(i,j,k2)=tl_dadz(i,j,k2)*vmask(i,j)
1976# endif
1977 END DO
1978 END DO
1979 END IF
1980!
1981! Compute components of the rotated A flux (A m3/s) along geopotential
1982! surfaces.
1983!
1984 IF (k.gt.0) THEN
1985 DO j=jstrv,jend
1986 DO i=istr,iend+1
1987 cff=0.0625_r8*(kh(i-1,j-1)+kh(i-1,j)+ &
1988 & kh(i ,j-1)+kh(i ,j))*on_p(i,j)
1989 cff1=min(dzdx_p(i,j,k1),0.0_r8)
1990 cff2=max(dzdx_p(i,j,k1),0.0_r8)
1991!^ FX(i,j)=cff* &
1992!^ & (Hz(i-1,j-1,k)+Hz(i-1,j,k)+ &
1993!^ & Hz(i ,j-1,k)+Hz(i ,j,k))* &
1994!^ & (dAdx(i,j,k1)- &
1995!^ & 0.5_r8*(cff1*(dAdz(i-1,j,k1)+ &
1996!^ & dAdz(i ,j,k2))+ &
1997!^ & cff2*(dAdz(i-1,j,k2)+ &
1998!^ & dAdz(i ,j,k1))))
1999!^
2000 tl_fx(i,j)=cff* &
2001 & (hz(i-1,j-1,k)+hz(i-1,j,k)+ &
2002 & hz(i ,j-1,k)+hz(i ,j,k))* &
2003 & (tl_dadx(i,j,k1)- &
2004 & 0.5_r8*(cff1*(tl_dadz(i-1,j,k1)+ &
2005 & tl_dadz(i ,j,k2))+ &
2006 & cff2*(tl_dadz(i-1,j,k2)+ &
2007 & tl_dadz(i ,j,k1))))
2008 END DO
2009 END DO
2010 DO j=jstrv-1,jend
2011 DO i=istr,iend
2012 cff=kh(i,j)*om_r(i,j)
2013 cff1=min(dzde_r(i,j,k1),0.0_r8)
2014 cff2=max(dzde_r(i,j,k1),0.0_r8)
2015!^ FE(i,j)=cff* &
2016!^ & Hz(i,j,k)* &
2017!^ & (dAde(i,j,k1)- &
2018!^ & 0.5_r8*(cff1*(dAdz(i,j ,k1)+ &
2019!^ & dAdz(i,j+1,k2))+ &
2020!^ & cff2*(dAdz(i,j ,k2)+ &
2021!^ & dAdz(i,j+1,k1))))
2022!^
2023 tl_fe(i,j)=cff* &
2024 & hz(i,j,k)* &
2025 & (tl_dade(i,j,k1)- &
2026 & 0.5_r8*(cff1*(tl_dadz(i,j ,k1)+ &
2027 & tl_dadz(i,j+1,k2))+ &
2028 & cff2*(tl_dadz(i,j ,k2)+ &
2029 & tl_dadz(i,j+1,k1))))
2030 END DO
2031 END DO
2032 IF (k.lt.n(ng)) THEN
2033 DO j=jstrv,jend
2034 DO i=istr,iend
2035 cff=0.5_r8*(kh(i,j-1)+kh(i,j))
2036 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
2037 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
2038 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
2039 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
2040!^ FZ(i,j,k2)=cff* &
2041!^ & (cff1*(cff1*dAdz(i,j,k2)-dAdx(i ,j,k1))+ &
2042!^ & cff2*(cff2*dAdz(i,j,k2)-dAdx(i+1,j,k2))+ &
2043!^ & cff3*(cff3*dAdz(i,j,k2)-dAdx(i ,j,k2))+ &
2044!^ & cff4*(cff4*dAdz(i,j,k2)-dAdx(i+1,j,k1)))
2045!^
2046 tl_fz(i,j,k2)=cff* &
2047 & (cff1*(cff1*tl_dadz(i,j,k2)- &
2048 & tl_dadx(i ,j,k1))+ &
2049 & cff2*(cff2*tl_dadz(i,j,k2)- &
2050 & tl_dadx(i+1,j,k2))+ &
2051 & cff3*(cff3*tl_dadz(i,j,k2)- &
2052 & tl_dadx(i ,j,k2))+ &
2053 & cff4*(cff4*tl_dadz(i,j,k2)- &
2054 & tl_dadx(i+1,j,k1)))
2055 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
2056 cff2=min(dzde_r(i,j ,k2),0.0_r8)
2057 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
2058 cff4=max(dzde_r(i,j ,k1),0.0_r8)
2059!^ FZ(i,j,k2)=FZ(i,j,k2)+ &
2060!^ & cff* &
2061!^ & (cff1*(cff1*dAdz(i,j,k2)-dAde(i,j-1,k1))+ &
2062!^ & cff2*(cff2*dAdz(i,j,k2)-dAde(i,j ,k2))+ &
2063!^ & cff3*(cff3*dAdz(i,j,k2)-dAde(i,j-1,k2))+ &
2064!^ & cff4*(cff4*dAdz(i,j,k2)-dAde(i,j ,k1)))
2065!^
2066 tl_fz(i,j,k2)=tl_fz(i,j,k2)+ &
2067 & cff* &
2068 & (cff1*(cff1*tl_dadz(i,j,k2)- &
2069 & tl_dade(i,j-1,k1))+ &
2070 & cff2*(cff2*tl_dadz(i,j,k2)- &
2071 & tl_dade(i,j ,k2))+ &
2072 & cff3*(cff3*tl_dadz(i,j,k2)- &
2073 & tl_dade(i,j-1,k2))+ &
2074 & cff4*(cff4*tl_dadz(i,j,k2)- &
2075 & tl_dade(i,j ,k1)))
2076 END DO
2077 END DO
2078 END IF
2079!
2080! Time-step harmonic, geopotential diffusion term (m Tunits).
2081!
2082 DO j=jstrv,jend
2083 DO i=istr,iend
2084!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
2085!^ & Hfac(i,j)* &
2086!^ & (FX(i+1,j)-FX(i,j )+ &
2087!^ & FE(i ,j)-FE(i,j-1))+ &
2088!^ & DTsizeH* &
2089!^ & (FZ(i,j,k2)-FZ(i,j,k1))
2090!^
2091 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
2092 & hfac(i,j)* &
2093 & (tl_fx(i+1,j)-tl_fx(i,j )+ &
2094 & tl_fe(i ,j)-tl_fe(i,j-1))+ &
2095 & dtsizeh* &
2096 & (tl_fz(i,j,k2)-tl_fz(i,j,k1))
2097 END DO
2098 END DO
2099 END IF
2100 END DO k_loop
2101
2102# else
2103
2104!
2105! Compute XI- and ETA-components of diffusive flux.
2106!
2107 DO k=1,n(ng)
2108 DO j=jstrv,jend
2109 DO i=istr,iend+1
2110!^ FX(i,j)=pmon_p(i,j)*0.25_r8*(Kh(i-1,j )+Kh(i,j )+ &
2111!^ & Kh(i-1,j-1)+Kh(i,j-1))* &
2112!^ & (Awrk(i,j,k,Nold)-Awrk(i-1,j,k,Nold))
2113!^
2114 tl_fx(i,j)=pmon_p(i,j)*0.25_r8*(kh(i-1,j )+kh(i,j )+ &
2115 & kh(i-1,j-1)+kh(i,j-1))* &
2116 & (tl_awrk(i,j,k,nold)-tl_awrk(i-1,j,k,nold))
2117# ifdef MASKING
2118!^ FX(i,j)=FX(i,j)*pmask(i,j)
2119!^
2120 tl_fx(i,j)=tl_fx(i,j)*pmask(i,j)
2121# endif
2122 END DO
2123 END DO
2124 DO j=jstrv-1,jend
2125 DO i=istr,iend
2126!^ FE(i,j)=pnom_r(i,j)*Kh(i,j)* &
2127!^ & (Awrk(i,j+1,k,Nold)-Awrk(i,j,k,Nold))
2128!^
2129 tl_fe(i,j)=pnom_r(i,j)*kh(i,j)* &
2130 & (tl_awrk(i,j+1,k,nold)-tl_awrk(i,j,k,nold))
2131 END DO
2132 END DO
2133!
2134! Time-step horizontal diffusion equation.
2135!
2136 DO j=jstrv,jend
2137 DO i=istr,iend
2138!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
2139!^ & Hfac(i,j)* &
2140!^ & (FX(i+1,j)-FX(i,j)+ &
2141!^ & FE(i,j)-FE(i,j-1))
2142!^
2143 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
2144 & hfac(i,j)* &
2145 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
2146 & tl_fe(i,j)-tl_fe(i,j-1))
2147 END DO
2148 END DO
2149 END DO
2150# endif
2151!
2152! Apply boundary conditions.
2153!
2154!^ CALL dabc_v3d_tile (ng, tile, &
2155!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
2156!^ & Awrk(:,:,:,Nnew))
2157!^
2158 CALL dabc_v3d_tile (ng, tile, &
2159 & lbi, ubi, lbj, ubj, lbk, ubk, &
2160 & tl_awrk(:,:,:,nnew))
2161# ifdef DISTRIBUTE
2162!^ CALL mp_exchange3d (ng, tile, model, 1, &
2163!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
2164!^ & Nghost, &
2165!^ & EWperiodic(ng), NSperiodic(ng), &
2166!^ & Awrk(:,:,:,Nnew))
2167!^
2168 CALL mp_exchange3d (ng, tile, model, 1, &
2169 & lbi, ubi, lbj, ubj, lbk, ubk, &
2170 & nghost, &
2171 & ewperiodic(ng), nsperiodic(ng), &
2172 & tl_awrk(:,:,:,nnew))
2173# endif
2174!
2175! Update integration indices.
2176!
2177 nsav=nold
2178 nold=nnew
2179 nnew=nsav
2180 END DO
2181
2182# ifdef VCONVOLUTION
2183# ifdef IMPLICIT_VCONV
2184# ifdef SPLINES_VCONV
2185!
2186!-----------------------------------------------------------------------
2187! Integrate vertical diffusion equation implicitly using parabolic
2188! splines.
2189!-----------------------------------------------------------------------
2190!
2191 DO step=1,nvsteps
2192!
2193! Use conservative, parabolic spline reconstruction of vertical
2194! diffusion derivatives. Then, time step vertical diffusion term
2195! implicitly.
2196!
2197 DO j=jstrv,jend
2198 DO k=1,n(ng)
2199 DO i=istr,iend
2200 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
2201 & hz(i,j ,k))
2202 END DO
2203 END DO
2204 cff1=1.0_r8/6.0_r8
2205 DO k=1,n(ng)-1
2206 DO i=istr,iend
2207 fc(i,k)=cff1*hzk(i,k )-dtsizev*kv(i,j,k-1)*ohz(i,j,k )
2208 cf(i,k)=cff1*hzk(i,k+1)-dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
2209 END DO
2210 END DO
2211 DO i=istr,iend
2212 cf(i,0)=0.0_r8
2213!^ DC(i,0)=0.0_r8
2214!^
2215 tl_dc(i,0)=0.0_r8
2216 END DO
2217!
2218! LU decomposition and forward substitution.
2219!
2220 cff1=1.0_r8/3.0_r8
2221 DO k=1,n(ng)-1
2222 DO i=istr,iend
2223 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
2224 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
2225 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
2226 cf(i,k)=cff*cf(i,k)
2227!^ DC(i,k)=cff*(Awrk(i,j,k+1,Nold)- &
2228!^ & Awrk(i,j,k ,Nold)- &
2229!^ & FC(i,k)*DC(i,k-1))
2230!^
2231 tl_dc(i,k)=cff*(tl_awrk(i,j,k+1,nold)- &
2232 & tl_awrk(i,j,k ,nold)- &
2233 & fc(i,k)*tl_dc(i,k-1))
2234 END DO
2235 END DO
2236!
2237! Backward substitution.
2238!
2239 DO i=istr,iend
2240!^ DC(i,N(ng))=0.0_r8
2241!^
2242 tl_dc(i,n(ng))=0.0_r8
2243 END DO
2244 DO k=n(ng)-1,1,-1
2245 DO i=istr,iend
2246!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
2247!^
2248 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
2249 END DO
2250 END DO
2251!
2252 DO k=1,n(ng)
2253 DO i=istr,iend
2254!^ DC(i,k)=DC(i,k)*Kv(i,j,k)
2255!^
2256 tl_dc(i,k)=tl_dc(i,k)*kv(i,j,k)
2257!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
2258!^ & DTsizeV*oHz(i,j,k)* &
2259!^ & (DC(i,k)-DC(i,k-1))
2260!^
2261 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
2262 & dtsizev*ohz(i,j,k)* &
2263 & (tl_dc(i,k)-tl_dc(i,k-1))
2264 END DO
2265 END DO
2266 END DO
2267!
2268! Update integration indices.
2269!
2270 nsav=nold
2271 nold=nnew
2272 nnew=nsav
2273 END DO
2274# else
2275!
2276!-----------------------------------------------------------------------
2277! Integerate vertical diffusion equation implicitly.
2278!-----------------------------------------------------------------------
2279!
2280 DO step=1,nvsteps
2281!
2282! Compute diagonal matrix coefficients BC and load right-hand-side
2283! terms for the tangent linear diffusion equation into tl_DC.
2284!
2285 DO j=jstrv,jend
2286 DO k=1,n(ng)
2287 DO i=istr,iend
2288 cff=0.5_r8*(hz(i,j-1,k)+hz(i,j,k))
2289 bc(i,k)=cff-fc(i,j,k)-fc(i,j,k-1)
2290!^ DC(i,k)=Awrk(i,j,k,Nold)*cff
2291!^
2292 tl_dc(i,k)=tl_awrk(i,j,k,nold)*cff
2293 END DO
2294 END DO
2295!
2296! Solve the tridiagonal system.
2297!
2298 DO i=istr,iend
2299 cff=1.0_r8/bc(i,1)
2300 cf(i,1)=cff*fc(i,j,1)
2301!^ DC(i,1)=cff*DC(i,1)
2302!^
2303 tl_dc(i,1)=cff*tl_dc(i,1)
2304 END DO
2305 DO k=2,n(ng)-1
2306 DO i=istr,iend
2307 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
2308 cf(i,k)=cff*fc(i,j,k)
2309!^ DC(i,k)=cff*(DC(i,k)-FC(i,j,k-1)*DC(i,k-1))
2310!^
2311 tl_dc(i,k)=cff*(tl_dc(i,k)-fc(i,j,k-1)*tl_dc(i,k-1))
2312 END DO
2313 END DO
2314!
2315! Compute new solution by back substitution.
2316!
2317 DO i=istr,iend
2318!^ DC(i,N(ng))=(DC(i,N(ng))- &
2319!^ & FC(i,j,N(ng)-1)*DC(i,N(ng)-1))/ &
2320!^ & (BC(i,N(ng))-FC(i,j,N(ng)-1)*CF(i,N(ng)-1))
2321!^
2322 tl_dc(i,n(ng))=(tl_dc(i,n(ng))- &
2323 & fc(i,j,n(ng)-1)*tl_dc(i,n(ng)-1))/ &
2324 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
2325!^ Awrk(i,j,N(ng),Nnew)=DC(i,N(ng))
2326!^
2327 tl_awrk(i,j,n(ng),nnew)=tl_dc(i,n(ng))
2328# ifdef MASKING
2329!^ Awrk(i,j,N(ng),Nnew)=Awrk(i,j,N(ng),Nnew)*vmask(i,j)
2330!^
2331 tl_awrk(i,j,n(ng),nnew)=tl_awrk(i,j,n(ng),nnew)*vmask(i,j)
2332# endif
2333 END DO
2334 DO k=n(ng)-1,1,-1
2335 DO i=istr,iend
2336!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
2337!^
2338 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
2339!^ Awrk(i,j,k,Nnew)=DC(i,k)
2340!^
2341 tl_awrk(i,j,k,nnew)=tl_dc(i,k)
2342# ifdef MASKING
2343!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nnew)*vmask(i,j)
2344!^
2345 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nnew)*vmask(i,j)
2346# endif
2347 END DO
2348 END DO
2349 END DO
2350!
2351! Update integration indices.
2352!
2353 nsav=nold
2354 nold=nnew
2355 nnew=nsav
2356 END DO
2357# endif
2358# else
2359!
2360!-----------------------------------------------------------------------
2361! Integerate vertical diffusion equation explicitly.
2362!-----------------------------------------------------------------------
2363!
2364 DO step=1,nvsteps
2365!
2366! Compute vertical diffusive flux. Notice that "FC" is assumed to
2367! be time invariant.
2368!
2369 DO j=jstrv,jend
2370 DO k=1,n(ng)-1
2371 DO i=istr,iend
2372!^ FS(i,k)=FC(i,j,k)*(Awrk(i,j,k+1,Nold)- &
2373!^ & Awrk(i,j,k ,Nold))
2374!^
2375 tl_fs(i,k)=fc(i,j,k)*(tl_awrk(i,j,k+1,nold)- &
2376 & tl_awrk(i,j,k ,nold))
2377# ifdef MASKING
2378!^ FS(i,k)=FS(i,k)*vmask(i,j)
2379!^
2380 tl_fs(i,k)=tl_fs(i,k)*vmask(i,j)
2381# endif
2382 END DO
2383 END DO
2384 DO i=istr,iend
2385!^ FS(i,0)=0.0_r8
2386!^
2387 tl_fs(i,0)=0.0_r8
2388!^ FS(i,N(ng))=0.0_r8
2389!^
2390 tl_fs(i,n(ng))=0.0_r8
2391 END DO
2392!
2393! Time-step vertical diffusive term. Notice that "oHz" is assumed to
2394! be time invariant.
2395!
2396 DO k=1,n(ng)
2397 DO i=istr,iend
2398!^ Awrk(i,j,k,Nnew)=Awrk(i,j,k,Nold)+ &
2399!^ & oHz(i,j,k)*(FS(i,k )- &
2400!^ & FS(i,k-1))
2401!^
2402 tl_awrk(i,j,k,nnew)=tl_awrk(i,j,k,nold)+ &
2403 & ohz(i,j,k)*(tl_fs(i,k )- &
2404 & tl_fs(i,k-1))
2405 END DO
2406 END DO
2407 END DO
2408!
2409! Update integration indices.
2410!
2411 nsav=nold
2412 nold=nnew
2413 nnew=nsav
2414 END DO
2415# endif
2416# endif
2417!
2418!-----------------------------------------------------------------------
2419! Load convolved solution.
2420!-----------------------------------------------------------------------
2421!
2422 DO k=1,n(ng)
2423 DO j=jstrv,jend
2424 DO i=istr,iend
2425!^ A(i,j,k)=Awrk(i,j,k,Nold)
2426!^
2427 tl_a(i,j,k)=tl_awrk(i,j,k,nold)
2428 END DO
2429 END DO
2430 END DO
2431!^ CALL dabc_v3d_tile (ng, tile, &
2432!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
2433!^ & A)
2434!^
2435 CALL dabc_v3d_tile (ng, tile, &
2436 & lbi, ubi, lbj, ubj, lbk, ubk, &
2437 & tl_a)
2438# ifdef DISTRIBUTE
2439!^ CALL mp_exchange3d (ng, tile, model, 1, &
2440!^ & LBi, UBi, LBj, UBj, LBk, UBk, &
2441!^ & Nghost, &
2442!^ & EWperiodic(ng), NSperiodic(ng), &
2443!^ & A)
2444!^
2445 CALL mp_exchange3d (ng, tile, model, 1, &
2446 & lbi, ubi, lbj, ubj, lbk, ubk, &
2447 & nghost, &
2448 & ewperiodic(ng), nsperiodic(ng), &
2449 & tl_a)
2450# endif
2451
2452 RETURN
2453 END SUBROUTINE tl_conv_v3d_tile
2454#endif
2455 END MODULE tl_conv_3d_mod
bc_3d_mod
Definition bc_3d.F:2
bc_3d_mod::dabc_u3d_tile
subroutine dabc_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition bc_3d.F:869
bc_3d_mod::dabc_r3d_tile
subroutine dabc_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition bc_3d.F:730
bc_3d_mod::dabc_v3d_tile
subroutine dabc_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition bc_3d.F:1010
mod_param
Definition mod_param.F:2
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
tl_conv_3d_mod
Definition tl_conv_3d.F:3
tl_conv_3d_mod::tl_conv_u3d_tile
subroutine tl_conv_u3d_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbk, ubk, imins, imaxs, jmins, jmaxs, nghost, nhsteps, nvsteps, dtsizeh, dtsizev, kh, kv, pm, pn, on_r, om_p, pmask, rmask, umask, vmask, hz, z_r, tl_a)
Definition tl_conv_3d.F:847
tl_conv_3d_mod::tl_conv_v3d_tile
subroutine tl_conv_v3d_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbk, ubk, imins, imaxs, jmins, jmaxs, nghost, nhsteps, nvsteps, dtsizeh, dtsizev, kh, kv, pm, pn, on_p, om_r, pmask, rmask, umask, vmask, hz, z_r, tl_a)
Definition tl_conv_3d.F:1662
tl_conv_3d_mod::tl_conv_r3d_tile
subroutine tl_conv_r3d_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbk, ubk, imins, imaxs, jmins, jmaxs, nghost, nhsteps, nvsteps, dtsizeh, dtsizev, kh, kv, pm, pn, on_u, om_v, rmask, umask, vmask, hz, z_r, tl_a)
Definition tl_conv_3d.F:85