ROMS
Loading...
Searching...
No Matches
conv_3d.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3 MODULE conv_3d_mod
4
5#if defined NONLINEAR && 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! A 3D state variable to diffuse. !
48! !
49! On Output: !
50! !
51! A Convolved 3D state variable. !
52! !
53! Routines: !
54! !
55! conv_r3d_tile Nonlinear 3D convolution at RHO-points !
56! conv_u3d_tile Nonlinear 3D convolution at U-points !
57! conv_v3d_tile Nonlinear 3D convolution at V-points !
58! !
59!=======================================================================
60!
61 implicit none
62!
63 PUBLIC
64!
65 CONTAINS
66!
67!***********************************************************************
68 SUBROUTINE 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 & 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 real(r8), intent(inout) :: A(LBi:,LBj:,LBk:)
125# else
126 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
127 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
128# ifdef GEOPOTENTIAL_HCONV
129 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
130 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
131# else
132 real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
133 real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
134# endif
135# ifdef MASKING
136 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
137 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
138 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
139# endif
140 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
141 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
142
143 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
144 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
145 real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj,LBk:UBk)
146# endif
147!
148! Local variable declarations.
149!
150 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
151
152 real(r8) :: cff, cff1, cff2, cff3, cff4
153
154 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: Awrk
155
156 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
157 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
158 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
159
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 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
171# ifdef SPLINES_VCONV
172 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
173# endif
174# else
175 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FS
176# endif
177# endif
178# ifdef GEOPOTENTIAL_HCONV
179 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FZ
180 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdz
181 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdx
182 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAde
183 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx
184 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde
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# ifdef DISTRIBUTE
230 CALL mp_exchange3d (ng, tile, model, 1, &
231 & lbi, ubi, lbj, ubj, lbk, ubk, &
232 & nghost, &
233 & ewperiodic(ng), nsperiodic(ng), &
234 & a)
235# endif
236 DO k=1,n(ng)
237 DO j=jstr-1,jend+1
238 DO i=istr-1,iend+1
239 awrk(i,j,k,nold)=a(i,j,k)
240 END DO
241 END DO
242 END DO
243!
244!-----------------------------------------------------------------------
245! Integrate horizontal diffusion equation.
246!-----------------------------------------------------------------------
247!
248 DO step=1,nhsteps
249
250# ifdef GEOPOTENTIAL_HCONV
251!
252! Diffusion along geopotential surfaces: Compute horizontal and
253! vertical gradients. Notice the recursive blocking sequence. The
254! vertical placement of the gradients is:
255!
256! dAdx,dAde(:,:,k1) k rho-points
257! dAdx,dAde(:,:,k2) k+1 rho-points
258! FZ,dAdz(:,:,k1) k-1/2 W-points
259! FZ,dAdz(:,:,k2) k+1/2 W-points
260!
261 k2=1
262 k_loop : DO k=0,n(ng)
263 k1=k2
264 k2=3-k1
265 IF (k.lt.n(ng)) THEN
266 DO j=jstr,jend
267 DO i=istr,iend+1
268 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
269# ifdef MASKING
270 cff=cff*umask(i,j)
271# endif
272 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
273 & z_r(i-1,j,k+1))
274# ifdef MASKING
275 dadx(i,j,k2)=cff*(awrk(i ,j,k+1,nold)*rmask(i ,j)- &
276 & awrk(i-1,j,k+1,nold)*rmask(i-1,j))
277# else
278 dadx(i,j,k2)=cff*(awrk(i ,j,k+1,nold)- &
279 & awrk(i-1,j,k+1,nold))
280# endif
281 END DO
282 END DO
283 DO j=jstr,jend+1
284 DO i=istr,iend
285 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
286# ifdef MASKING
287 cff=cff*vmask(i,j)
288# endif
289 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
290 & z_r(i,j-1,k+1))
291# ifdef MASKING
292 dade(i,j,k2)=cff*(awrk(i,j ,k+1,nold)*rmask(i,j )- &
293 & awrk(i,j-1,k+1,nold)*rmask(i,j-1))
294# else
295 dade(i,j,k2)=cff*(awrk(i,j ,k+1,nold)- &
296 & awrk(i,j-1,k+1,nold))
297# endif
298 END DO
299 END DO
300 END IF
301 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
302 DO j=jstr-1,jend+1
303 DO i=istr-1,iend+1
304 dadz(i,j,k2)=0.0_r8
305 fz(i,j,k2)=0.0_r8
306 END DO
307 END DO
308 ELSE
309 DO j=jstr-1,jend+1
310 DO i=istr-1,iend+1
311 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
312 dadz(i,j,k2)=cff*(awrk(i,j,k+1,nold)- &
313 & awrk(i,j,k ,nold))
314# ifdef MASKING
315 dadz(i,j,k2)=dadz(i,j,k2)*rmask(i,j)
316# endif
317 END DO
318 END DO
319 END IF
320!
321! Compute components of the rotated A flux (A m3/s) along geopotential
322! surfaces.
323!
324 IF (k.gt.0) THEN
325 DO j=jstr,jend
326 DO i=istr,iend+1
327 cff=0.25_r8*(kh(i-1,j)+kh(i-1,j))*on_u(i,j)
328 cff1=min(dzdx(i,j,k1),0.0_r8)
329 cff2=max(dzdx(i,j,k1),0.0_r8)
330 fx(i,j)=cff* &
331 & (hz(i,j,k)+hz(i-1,j,k))* &
332 & (dadx(i,j,k1)- &
333 & 0.5_r8*(cff1*(dadz(i-1,j,k1)+ &
334 & dadz(i ,j,k2))+ &
335 & cff2*(dadz(i-1,j,k2)+ &
336 & dadz(i ,j,k1))))
337 END DO
338 END DO
339 DO j=jstr,jend+1
340 DO i=istr,iend
341 cff=0.25_r8*(kh(i,j-1)+kh(i,j))*om_v(i,j)
342 cff1=min(dzde(i,j,k1),0.0_r8)
343 cff2=max(dzde(i,j,k1),0.0_r8)
344 fe(i,j)=cff* &
345 & (hz(i,j,k)+hz(i,j-1,k))* &
346 & (dade(i,j,k1)- &
347 & 0.5_r8*(cff1*(dadz(i,j-1,k1)+ &
348 & dadz(i,j ,k2))+ &
349 & cff2*(dadz(i,j-1,k2)+ &
350 & dadz(i,j ,k1))))
351 END DO
352 END DO
353 IF (k.lt.n(ng)) THEN
354 DO j=jstr,jend
355 DO i=istr,iend
356 cff=0.5_r8*kh(i,j)
357 cff1=min(dzdx(i ,j,k1),0.0_r8)
358 cff2=min(dzdx(i+1,j,k2),0.0_r8)
359 cff3=max(dzdx(i ,j,k2),0.0_r8)
360 cff4=max(dzdx(i+1,j,k1),0.0_r8)
361 fz(i,j,k2)=cff* &
362 & (cff1*(cff1*dadz(i,j,k2)-dadx(i ,j,k1))+ &
363 & cff2*(cff2*dadz(i,j,k2)-dadx(i+1,j,k2))+ &
364 & cff3*(cff3*dadz(i,j,k2)-dadx(i ,j,k2))+ &
365 & cff4*(cff4*dadz(i,j,k2)-dadx(i+1,j,k1)))
366 cff1=min(dzde(i,j ,k1),0.0_r8)
367 cff2=min(dzde(i,j+1,k2),0.0_r8)
368 cff3=max(dzde(i,j ,k2),0.0_r8)
369 cff4=max(dzde(i,j+1,k1),0.0_r8)
370 fz(i,j,k2)=fz(i,j,k2)+ &
371 & cff* &
372 & (cff1*(cff1*dadz(i,j,k2)-dade(i,j ,k1))+ &
373 & cff2*(cff2*dadz(i,j,k2)-dade(i,j+1,k2))+ &
374 & cff3*(cff3*dadz(i,j,k2)-dade(i,j ,k2))+ &
375 & cff4*(cff4*dadz(i,j,k2)-dade(i,j+1,k1)))
376 END DO
377 END DO
378 END IF
379!
380! Time-step harmonic, geopotential diffusion term (m Tunits).
381!
382 DO j=jstr,jend
383 DO i=istr,iend
384 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
385 & hfac(i,j)* &
386 & (fx(i+1,j )-fx(i,j)+ &
387 & fe(i ,j+1)-fe(i,j))+ &
388 & dtsizeh* &
389 & (fz(i,j,k2)-fz(i,j,k1))
390 END DO
391 END DO
392 END IF
393 END DO k_loop
394
395# else
396
397!
398! Diffusion along S-coordinates: compute XI- and ETA-components of
399! diffusive flux.
400!
401 DO k=1,n(ng)
402 DO j=jstr,jend
403 DO i=istr,iend+1
404 fx(i,j)=pmon_u(i,j)*0.5_r8*(kh(i-1,j)+kh(i,j))* &
405 & (awrk(i,j,k,nold)-awrk(i-1,j,k,nold))
406# ifdef MASKING
407 fx(i,j)=fx(i,j)*umask(i,j)
408# endif
409 END DO
410 END DO
411 DO j=jstr,jend+1
412 DO i=istr,iend
413 fe(i,j)=pnom_v(i,j)*0.5_r8*(kh(i,j-1)+kh(i,j))* &
414 & (awrk(i,j,k,nold)-awrk(i,j-1,k,nold))
415# ifdef MASKING
416 fe(i,j)=fe(i,j)*vmask(i,j)
417# endif
418 END DO
419 END DO
420!
421! Time-step horizontal diffusion equation.
422!
423 DO j=jstr,jend
424 DO i=istr,iend
425 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
426 & hfac(i,j)* &
427 & (fx(i+1,j)-fx(i,j)+ &
428 & fe(i,j+1)-fe(i,j))
429 END DO
430 END DO
431 END DO
432# endif
433!
434! Apply boundary conditions.
435!
436 CALL dabc_r3d_tile (ng, tile, &
437 & lbi, ubi, lbj, ubj, lbk, ubk, &
438 & awrk(:,:,:,nnew))
439# ifdef DISTRIBUTE
440 CALL mp_exchange3d (ng, tile, model, 1, &
441 & lbi, ubi, lbj, ubj, lbk, ubk, &
442 & nghost, &
443 & ewperiodic(ng), nsperiodic(ng), &
444 & awrk(:,:,:,nnew))
445# endif
446!
447! Update integration indices.
448!
449 nsav=nold
450 nold=nnew
451 nnew=nsav
452 END DO
453
454# ifdef VCONVOLUTION
455# ifdef IMPLICIT_VCONV
456# ifdef SPLINES_VCONV
457!
458!-----------------------------------------------------------------------
459! Integrate vertical diffusion equation implicitly using parabolic
460! splines.
461!-----------------------------------------------------------------------
462!
463 DO step=1,nvsteps
464!
465! Use conservative, parabolic spline reconstruction of vertical
466! diffusion derivatives. Then, time step vertical diffusion term
467! implicitly.
468!
469 DO j=jstr,jend
470 cff1=1.0_r8/6.0_r8
471 DO k=1,n(ng)-1
472 DO i=istr,iend
473 fc(i,k)=cff1*hz(i,j,k )- &
474 & dtsizev*kv(i,j,k-1)*ohz(i,j,k )
475 cf(i,k)=cff1*hz(i,j,k+1)- &
476 & dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
477 END DO
478 END DO
479 DO i=istr,iend
480 cf(i,0)=0.0_r8
481 dc(i,0)=0.0_r8
482 END DO
483!
484! LU decomposition and forward substitution.
485!
486 cff1=1.0_r8/3.0_r8
487 DO k=1,n(ng)-1
488 DO i=istr,iend
489 bc(i,k)=cff1*(hz(i,j,k)+hz(i,j,k+1))+ &
490 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
491 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
492 cf(i,k)=cff*cf(i,k)
493 dc(i,k)=cff*(awrk(i,j,k+1,nold)-awrk(i,j,k,nold)- &
494 & fc(i,k)*dc(i,k-1))
495 END DO
496 END DO
497!
498! Backward substitution.
499!
500 DO i=istr,iend
501 dc(i,n(ng))=0.0_r8
502 END DO
503 DO k=n(ng)-1,1,-1
504 DO i=istr,iend
505 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
506 END DO
507 END DO
508!
509 DO k=1,n(ng)
510 DO i=istr,iend
511 dc(i,k)=dc(i,k)*kv(i,j,k)
512 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
513 & dtsizev*ohz(i,j,k)* &
514 & (dc(i,k)-dc(i,k-1))
515 END DO
516 END DO
517 END DO
518!
519! Update integration indices.
520!
521 nsav=nold
522 nold=nnew
523 nnew=nsav
524 END DO
525# else
526!
527!-----------------------------------------------------------------------
528! Integrate vertical diffusion equation implicitly.
529!-----------------------------------------------------------------------
530!
531 DO step=1,nvsteps
532!
533! Compute diagonal matrix coefficients BC and load right-hand-side
534! terms for the diffusion equation into DC.
535!
536 DO j=jstr,jend
537 DO k=1,n(ng)
538 DO i=istr,iend
539 bc(i,k)=hz(i,j,k)-fc(i,j,k)-fc(i,j,k-1)
540 dc(i,k)=awrk(i,j,k,nold)*hz(i,j,k)
541 END DO
542 END DO
543!
544! Solve the tridiagonal system.
545!
546 DO i=istr,iend
547 cff=1.0_r8/bc(i,1)
548 cf(i,1)=cff*fc(i,j,1)
549 dc(i,1)=cff*dc(i,1)
550 END DO
551 DO k=2,n(ng)-1
552 DO i=istr,iend
553 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
554 cf(i,k)=cff*fc(i,j,k)
555 dc(i,k)=cff*(dc(i,k)-fc(i,j,k-1)*dc(i,k-1))
556 END DO
557 END DO
558!
559! Compute new solution by back substitution.
560!
561 DO i=istr,iend
562 dc(i,n(ng))=(dc(i,n(ng))- &
563 & fc(i,j,n(ng)-1)*dc(i,n(ng)-1))/ &
564 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
565 awrk(i,j,n(ng),nnew)=dc(i,n(ng))
566# ifdef MASKING
567 awrk(i,j,n(ng),nnew)=awrk(i,j,n(ng),nnew)*rmask(i,j)
568# endif
569 END DO
570 DO k=n(ng)-1,1,-1
571 DO i=istr,iend
572 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
573 awrk(i,j,k,nnew)=dc(i,k)
574# ifdef MASKING
575 awrk(i,j,k,nnew)=awrk(i,j,k,nnew)*rmask(i,j)
576# endif
577 END DO
578 END DO
579 END DO
580!
581! Update integration indices.
582!
583 nsav=nold
584 nold=nnew
585 nnew=nsav
586 END DO
587# endif
588# else
589!
590!-----------------------------------------------------------------------
591! Integrate vertical diffusion equation explicitly.
592!-----------------------------------------------------------------------
593!
594 DO step=1,nvsteps
595!
596! Compute vertical diffusive flux. Notice that "FC" is assumed to
597! be time invariant.
598!
599 DO j=jstr,jend
600 DO k=1,n(ng)-1
601 DO i=istr,iend
602 fs(i,k)=fc(i,j,k)*(awrk(i,j,k+1,nold)- &
603 & awrk(i,j,k ,nold))
604# ifdef MASKING
605 fs(i,k)=fs(i,k)*rmask(i,j)
606# endif
607 END DO
608 END DO
609 DO i=istr,iend
610 fs(i,0)=0.0_r8
611 fs(i,n(ng))=0.0_r8
612 END DO
613!
614! Time-step vertical diffusive term. Notice that "oHz" is assumed to
615! be time invariant.
616!
617 DO k=1,n(ng)
618 DO i=istr,iend
619 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
620 & ohz(i,j,k)*(fs(i,k )- &
621 & fs(i,k-1))
622 END DO
623 END DO
624 END DO
625!
626! Update integration indices.
627!
628 nsav=nold
629 nold=nnew
630 nnew=nsav
631 END DO
632# endif
633# endif
634!
635!-----------------------------------------------------------------------
636! Load convolved solution.
637!-----------------------------------------------------------------------
638!
639 DO k=1,n(ng)
640 DO j=jstr,jend
641 DO i=istr,iend
642 a(i,j,k)=awrk(i,j,k,nold)
643 END DO
644 END DO
645 END DO
646 CALL dabc_r3d_tile (ng, tile, &
647 & lbi, ubi, lbj, ubj, lbk, ubk, &
648 & a)
649# ifdef DISTRIBUTE
650 CALL mp_exchange3d (ng, tile, model, 1, &
651 & lbi, ubi, lbj, ubj, lbk, ubk, &
652 & nghost, &
653 & ewperiodic(ng), nsperiodic(ng), &
654 & a)
655# endif
656
657 RETURN
658 END SUBROUTINE conv_r3d_tile
659!
660!***********************************************************************
661 SUBROUTINE conv_u3d_tile (ng, tile, model, &
662 & LBi, UBi, LBj, UBj, LBk, UBk, &
663 & IminS, ImaxS, JminS, JmaxS, &
664 & Nghost, NHsteps, NVsteps, &
665 & DTsizeH, DTsizeV, &
666 & Kh, Kv, &
667 & pm, pn, &
668# ifdef GEOPOTENTIAL_HCONV
669 & on_r, om_p, &
670# else
671 & pmon_r, pnom_p, &
672# endif
673# ifdef MASKING
674# ifdef GEOPOTENTIAL_HCONV
675 & pmask, rmask, umask, vmask, &
676# else
677 & umask, pmask, &
678# endif
679# endif
680 & Hz, z_r, &
681 & A)
682!***********************************************************************
683!
684 USE mod_param
685 USE mod_scalars
686!
687 USE bc_3d_mod, ONLY: dabc_u3d_tile
688# ifdef DISTRIBUTE
689 USE mp_exchange_mod, ONLY : mp_exchange3d
690# endif
691!
692! Imported variable declarations.
693!
694 integer, intent(in) :: ng, tile, model
695 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
696 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
697 integer, intent(in) :: Nghost, NHsteps, NVsteps
698
699 real(r8), intent(in) :: DTsizeH, DTsizeV
700!
701# ifdef ASSUMED_SHAPE
702 real(r8), intent(in) :: pm(LBi:,LBj:)
703 real(r8), intent(in) :: pn(LBi:,LBj:)
704# ifdef GEOPOTENTIAL_HCONV
705 real(r8), intent(in) :: on_r(LBi:,LBj:)
706 real(r8), intent(in) :: om_p(LBi:,LBj:)
707# else
708 real(r8), intent(in) :: pmon_r(LBi:,LBj:)
709 real(r8), intent(in) :: pnom_p(LBi:,LBj:)
710# endif
711# ifdef MASKING
712# ifdef GEOPOTENTIAL_HCONV
713 real(r8), intent(in) :: pmask(LBi:,LBj:)
714 real(r8), intent(in) :: rmask(LBi:,LBj:)
715 real(r8), intent(in) :: umask(LBi:,LBj:)
716 real(r8), intent(in) :: vmask(LBi:,LBj:)
717# else
718 real(r8), intent(in) :: umask(LBi:,LBj:)
719 real(r8), intent(in) :: pmask(LBi:,LBj:)
720# endif
721# endif
722 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
723 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
724
725 real(r8), intent(in) :: Kh(LBi:,LBj:)
726 real(r8), intent(in) :: Kv(LBi:,LBj:,0:)
727 real(r8), intent(inout) :: A(LBi:,LBj:,LBk:)
728# else
729 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
730 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
731# ifdef GEOPOTENTIAL_HCONV
732 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
733 real(r8), intent(in) :: om_p(LBi:UBi,LBj:UBj)
734# else
735 real(r8), intent(in) :: pmon_r(LBi:UBi,LBj:UBj)
736 real(r8), intent(in) :: pnom_p(LBi:UBi,LBj:UBj)
737# endif
738# ifdef MASKING
739# ifdef GEOPOTENTIAL_HCONV
740 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
741 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
742 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
743 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
744# else
745 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
746 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
747# endif
748# endif
749 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
750 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
751
752 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
753 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
754 real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj,LBk:UBk)
755# endif
756!
757! Local variable declarations.
758!
759 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
760
761 real(r8) :: cff, cff1, cff2, cff3, cff4
762
763 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: Awrk
764
765 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
766 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
767 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
768# ifdef VCONVOLUTION
769# ifndef SPLINES_VCONV
770 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FC
771# endif
772# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
773 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
774# endif
775# if defined IMPLICIT_VCONV || defined SPLINES_VCONV
776 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
777 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
778 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
779# ifdef SPLINES_VCONV
780 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
781 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
782# endif
783# else
784 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FS
785# endif
786# endif
787# ifdef GEOPOTENTIAL_HCONV
788 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZdx
789 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZde
790
791 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FZ
792 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdz
793 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdx
794 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAde
795 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
796 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
797# endif
798
799# include "set_bounds.h"
800!
801!-----------------------------------------------------------------------
802! Space convolution of the diffusion equation for a 3D state variable
803! at U-points.
804!-----------------------------------------------------------------------
805!
806! Compute metrics factors. Notice that "z_r" and "Hz" are assumed to
807! be time invariant in the vertical convolution. Scratch array are
808! used for efficiency.
809!
810 cff=dtsizeh*0.25_r8
811 DO j=jstr-1,jend+1
812 DO i=istru-1,iend+1
813 hfac(i,j)=cff*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
814# ifdef VCONVOLUTION
815# ifndef SPLINES_VCONV
816 fc(i,j,n(ng))=0.0_r8
817 DO k=1,n(ng)-1
818# ifdef IMPLICIT_VCONV
819 fc(i,j,k)=-dtsizev*(kv(i-1,j,k)+kv(i,j,k))/ &
820 & (z_r(i-1,j,k+1)+z_r(i,j,k+1)- &
821 & z_r(i-1,j,k )-z_r(i,j,k ))
822# else
823 fc(i,j,k)=dtsizev*(kv(i-1,j,k)+kv(i,j,k))/ &
824 & (z_r(i-1,j,k+1)+z_r(i,j,k+1)- &
825 & z_r(i-1,j,k )-z_r(i,j,k ))
826# endif
827 END DO
828 fc(i,j,0)=0.0_r8
829# endif
830# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
831 DO k=1,n(ng)
832 ohz(i,j,k)=2.0_r8/(hz(i-1,j,k)+hz(i,j,k))
833 END DO
834# endif
835# endif
836 END DO
837 END DO
838!
839! Set integration indices and initial conditions.
840!
841 nold=1
842 nnew=2
843 CALL dabc_u3d_tile (ng, tile, &
844 & lbi, ubi, lbj, ubj, lbk, ubk, &
845 & a)
846# ifdef DISTRIBUTE
847 CALL mp_exchange3d (ng, tile, model, 1, &
848 & lbi, ubi, lbj, ubj, lbk, ubk, &
849 & nghost, &
850 & ewperiodic(ng), nsperiodic(ng), &
851 & a)
852# endif
853 DO k=1,n(ng)
854 DO j=jstr-1,jend+1
855 DO i=istru-1,iend+1
856 awrk(i,j,k,nold)=a(i,j,k)
857 END DO
858 END DO
859 END DO
860!
861!-----------------------------------------------------------------------
862! Integrate horizontal diffusion equation.
863!-----------------------------------------------------------------------
864!
865 DO step=1,nhsteps
866
867# ifdef GEOPOTENTIAL_HCONV
868!
869! Diffusion along geopotential surfaces: Compute horizontal and
870! vertical gradients. Notice the recursive blocking sequence. The
871! vertical placement of the gradients is:
872!
873! dAdx,dAde(:,:,k1) k rho-points
874! dAdx,dAde(:,:,k2) k+1 rho-points
875! FZ,dAdz(:,:,k1) k-1/2 W-points
876! FZ,dAdz(:,:,k2) k+1/2 W-points
877!
878 k2=1
879 k_loop : DO k=0,n(ng)
880 k1=k2
881 k2=3-k1
882 IF (k.lt.n(ng)) THEN
883 DO j=jstr,jend
884 DO i=istru-1,iend+1
885 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
886# ifdef MASKING
887 cff=cff*umask(i,j)
888# endif
889 dzdx(i,j)=cff*(z_r(i ,j,k+1)- &
890 & z_r(i-1,j,k+1))
891 END DO
892 END DO
893!
894 DO j=jstr,jend+1
895 DO i=istru-1,iend
896 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
897# ifdef MASKING
898 cff=cff*vmask(i,j)
899# endif
900 dzde(i,j)=cff*(z_r(i,j ,k+1)- &
901 & z_r(i,j-1,k+1))
902 END DO
903 END DO
904!
905 DO j=jstr,jend
906 DO i=istru-1,iend
907# ifdef MASKING
908 dadx(i,j,k2)=pm(i,j)* &
909 & (awrk(i+1,j,k+1,nold)*umask(i+1,j)- &
910 & awrk(i ,j,k+1,nold)*umask(i ,j))
911 dadx(i,j,k2)=dadx(i,j,k2)*rmask(i,j)
912# else
913 dadx(i,j,k2)=pm(i,j)*(awrk(i+1,j,k+1,nold)- &
914 & awrk(i ,j,k+1,nold))
915# endif
916 dzdx_r(i,j,k2)=0.5_r8*(dzdx(i ,j)+ &
917 & dzdx(i+1,j))
918 END DO
919 END DO
920!
921 DO j=jstr,jend+1
922 DO i=istru,iend
923 cff=0.25_r8*(pn(i-1,j )+pn(i,j )+ &
924 & pn(i-1,j-1)+pn(i,j-1))
925# ifdef MASKING
926 dade(i,j,k2)=cff* &
927 & (awrk(i,j ,k+1,nold)*umask(i,j )- &
928 & awrk(i,j-1,k+1,nold)*umask(i,j-1))
929 dade(i,j,k2)=dade(i,j,k2)*pmask(i,j)
930# else
931 dade(i,j,k2)=cff*(awrk(i,j ,k+1,nold)- &
932 & awrk(i,j-1,k+1,nold))
933# endif
934 dzde_p(i,j,k2)=0.5_r8*(dzde(i-1,j)+ &
935 & dzde(i ,j))
936 END DO
937 END DO
938 END IF
939!
940 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
941 DO j=jstr-1,jend+1
942 DO i=istru-1,iend+1
943 dadz(i,j,k2)=0.0_r8
944 fz(i,j,k2)=0.0_r8
945 END DO
946 END DO
947 ELSE
948 DO j=jstr-1,jend+1
949 DO i=istru-1,iend+1
950 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
951 dadz(i,j,k2)=cff*(awrk(i,j,k+1,nold)- &
952 & awrk(i,j,k ,nold))
953# ifdef MASKING
954 dadz(i,j,k2)=dadz(i,j,k2)*umask(i,j)
955# endif
956 END DO
957 END DO
958 END IF
959!
960! Compute components of the rotated A flux (A m3/s) along geopotential
961! surfaces.
962!
963 IF (k.gt.0) THEN
964 DO j=jstr,jend
965 DO i=istru-1,iend
966 cff=kh(i,j)*on_r(i,j)
967 cff1=min(dzdx_r(i,j,k1),0.0_r8)
968 cff2=max(dzdx_r(i,j,k1),0.0_r8)
969 fx(i,j)=cff* &
970 & hz(i,j,k)* &
971 & (dadx(i,j,k1)- &
972 & 0.5_r8*(cff1*(dadz(i ,j,k1)+ &
973 & dadz(i+1,j,k2))+ &
974 & cff2*(dadz(i ,j,k2)+ &
975 & dadz(i+1,j,k1))))
976 END DO
977 END DO
978 DO j=jstr,jend+1
979 DO i=istru,iend
980 cff=0.0625_r8*(kh(i-1,j-1)+kh(i-1,j)+ &
981 & kh(i ,j-1)+kh(i ,j))*om_p(i,j)
982 cff1=min(dzde_p(i,j,k1),0.0_r8)
983 cff2=max(dzde_p(i,j,k1),0.0_r8)
984 fe(i,j)=cff* &
985 & (hz(i-1,j-1,k)+hz(i-1,j,k)+ &
986 & hz(i ,j-1,k)+hz(i ,j,k))* &
987 & (dade(i,j,k1)- &
988 & 0.5_r8*(cff1*(dadz(i,j-1,k1)+ &
989 & dadz(i,j ,k2))+ &
990 & cff2*(dadz(i,j-1,k2)+ &
991 & dadz(i,j ,k1))))
992 END DO
993 END DO
994 IF (k.lt.n(ng)) THEN
995 DO j=jstr,jend
996 DO i=istru,iend
997 cff=0.25_r8*(kh(i-1,j)+kh(i,j))
998 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
999 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
1000 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
1001 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
1002 fz(i,j,k2)=cff* &
1003 & (cff1*(cff1*dadz(i,j,k2)-dadx(i-1,j,k1))+ &
1004 & cff2*(cff2*dadz(i,j,k2)-dadx(i ,j,k2))+ &
1005 & cff3*(cff3*dadz(i,j,k2)-dadx(i-1,j,k2))+ &
1006 & cff4*(cff4*dadz(i,j,k2)-dadx(i ,j,k1)))
1007 cff1=min(dzde_p(i,j ,k1),0.0_r8)
1008 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
1009 cff3=max(dzde_p(i,j ,k2),0.0_r8)
1010 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
1011 fz(i,j,k2)=fz(i,j,k2)+ &
1012 & cff* &
1013 & (cff1*(cff1*dadz(i,j,k2)-dade(i,j ,k1))+ &
1014 & cff2*(cff2*dadz(i,j,k2)-dade(i,j+1,k2))+ &
1015 & cff3*(cff3*dadz(i,j,k2)-dade(i,j ,k2))+ &
1016 & cff4*(cff4*dadz(i,j,k2)-dade(i,j+1,k1)))
1017 END DO
1018 END DO
1019 END IF
1020!
1021! Time-step harmonic, geopotential diffusion term (m Tunits).
1022!
1023 DO j=jstr,jend
1024 DO i=istru,iend
1025 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1026 & hfac(i,j)* &
1027 & (fx(i,j )-fx(i-1,j)+ &
1028 & fe(i,j+1)-fe(i ,j))+ &
1029 & dtsizeh* &
1030 & (fz(i,j,k2)-fz(i,j,k1))
1031 END DO
1032 END DO
1033 END IF
1034 END DO k_loop
1035
1036# else
1037
1038!
1039! Diffusion along S-coordinates: compute XI- and ETA-components of
1040! diffusive flux.
1041!
1042 DO k=1,n(ng)
1043 DO j=jstr,jend
1044 DO i=istru-1,iend
1045 fx(i,j)=pmon_r(i,j)*kh(i,j)* &
1046 & (awrk(i+1,j,k,nold)-awrk(i,j,k,nold))
1047 END DO
1048 END DO
1049 DO j=jstr,jend+1
1050 DO i=istru,iend
1051 fe(i,j)=pnom_p(i,j)*0.25_r8*(kh(i-1,j )+kh(i,j )+ &
1052 & kh(i-1,j-1)+kh(i,j-1))* &
1053 & (awrk(i,j,k,nold)-awrk(i,j-1,k,nold))
1054# ifdef MASKING
1055 fe(i,j)=fe(i,j)*pmask(i,j)
1056# endif
1057 END DO
1058 END DO
1059!
1060! Time-step horizontal diffusion equation.
1061!
1062 DO j=jstr,jend
1063 DO i=istru,iend
1064 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1065 & hfac(i,j)* &
1066 & (fx(i,j)-fx(i-1,j)+ &
1067 & fe(i,j+1)-fe(i,j))
1068 END DO
1069 END DO
1070 END DO
1071# endif
1072!
1073! Apply boundary conditions.
1074!
1075 CALL dabc_u3d_tile (ng, tile, &
1076 & lbi, ubi, lbj, ubj, lbk, ubk, &
1077 & awrk(:,:,:,nnew))
1078# ifdef DISTRIBUTE
1079 CALL mp_exchange3d (ng, tile, model, 1, &
1080 & lbi, ubi, lbj, ubj, lbk, ubk, &
1081 & nghost, &
1082 & ewperiodic(ng), nsperiodic(ng), &
1083 & awrk(:,:,:,nnew))
1084# endif
1085!
1086! Update integration indices.
1087!
1088 nsav=nold
1089 nold=nnew
1090 nnew=nsav
1091 END DO
1092
1093# ifdef VCONVOLUTION
1094# ifdef IMPLICIT_VCONV
1095# ifdef SPLINES_VCONV
1096!
1097!-----------------------------------------------------------------------
1098! Integrate vertical diffusion equation implicitly using parabolic
1099! splines.
1100!-----------------------------------------------------------------------
1101!
1102 DO step=1,nvsteps
1103!
1104! Use conservative, parabolic spline reconstruction of vertical
1105! diffusion derivatives. Then, time step vertical diffusion term
1106! implicitly.
1107!
1108 DO j=jstr,jend
1109 DO k=1,n(ng)
1110 DO i=istru,iend
1111 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
1112 & hz(i ,j,k))
1113 END DO
1114 END DO
1115 cff1=1.0_r8/6.0_r8
1116 DO k=1,n(ng)-1
1117 DO i=istru,iend
1118 fc(i,k)=cff1*hzk(i,k )-dtsizev*kv(i,j,k-1)*ohz(i,j,k )
1119 cf(i,k)=cff1*hzk(i,k+1)-dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
1120 END DO
1121 END DO
1122 DO i=istru,iend
1123 cf(i,0)=0.0_r8
1124 dc(i,0)=0.0_r8
1125 END DO
1126!
1127! LU decomposition and forward substitution.
1128!
1129 cff1=1.0_r8/3.0_r8
1130 DO k=1,n(ng)-1
1131 DO i=istru,iend
1132 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1133 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
1134 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1135 cf(i,k)=cff*cf(i,k)
1136 dc(i,k)=cff*(awrk(i,j,k+1,nold)-awrk(i,j,k,nold)- &
1137 & fc(i,k)*dc(i,k-1))
1138 END DO
1139 END DO
1140!
1141! Backward substitution.
1142!
1143 DO i=istru,iend
1144 dc(i,n(ng))=0.0_r8
1145 END DO
1146 DO k=n(ng)-1,1,-1
1147 DO i=istru,iend
1148 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1149 END DO
1150 END DO
1151!
1152 DO k=1,n(ng)
1153 DO i=istru,iend
1154 dc(i,k)=dc(i,k)*kv(i,j,k)
1155 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1156 & dtsizev*ohz(i,j,k)* &
1157 & (dc(i,k)-dc(i,k-1))
1158 END DO
1159 END DO
1160 END DO
1161!
1162! Update integration indices.
1163!
1164 nsav=nold
1165 nold=nnew
1166 nnew=nsav
1167 END DO
1168# else
1169!
1170!-----------------------------------------------------------------------
1171! Integrate vertical diffusion equation implicitly.
1172!-----------------------------------------------------------------------
1173!
1174 DO step=1,nvsteps
1175!
1176! Compute diagonal matrix coefficients BC and load right-hand-side
1177! terms for the diffusion equation into DC.
1178!
1179 DO j=jstr,jend
1180 DO k=1,n(ng)
1181 DO i=istru,iend
1182 cff=0.5_r8*(hz(i-1,j,k)+hz(i,j,k))
1183 bc(i,k)=cff-fc(i,j,k)-fc(i,j,k-1)
1184 dc(i,k)=awrk(i,j,k,nold)*cff
1185 END DO
1186 END DO
1187!
1188! Solve the tridiagonal system.
1189!
1190 DO i=istru,iend
1191 cff=1.0_r8/bc(i,1)
1192 cf(i,1)=cff*fc(i,j,1)
1193 dc(i,1)=cff*dc(i,1)
1194 END DO
1195 DO k=2,n(ng)-1
1196 DO i=istru,iend
1197 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
1198 cf(i,k)=cff*fc(i,j,k)
1199 dc(i,k)=cff*(dc(i,k)-fc(i,j,k-1)*dc(i,k-1))
1200 END DO
1201 END DO
1202!
1203! Compute new solution by back substitution.
1204!
1205 DO i=istru,iend
1206 dc(i,n(ng))=(dc(i,n(ng))- &
1207 & fc(i,j,n(ng)-1)*dc(i,n(ng)-1))/ &
1208 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
1209 awrk(i,j,n(ng),nnew)=dc(i,n(ng))
1210# ifdef MASKING
1211 awrk(i,j,n(ng),nnew)=awrk(i,j,n(ng),nnew)*umask(i,j)
1212# endif
1213 END DO
1214 DO k=n(ng)-1,1,-1
1215 DO i=istru,iend
1216 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1217 awrk(i,j,k,nnew)=dc(i,k)
1218# ifdef MASKING
1219 awrk(i,j,k,nnew)=awrk(i,j,k,nnew)*umask(i,j)
1220# endif
1221 END DO
1222 END DO
1223 END DO
1224!
1225! Update integration indices.
1226!
1227 nsav=nold
1228 nold=nnew
1229 nnew=nsav
1230 END DO
1231# endif
1232# else
1233!
1234!-----------------------------------------------------------------------
1235! Integrate vertical diffusion equation explicitly.
1236!-----------------------------------------------------------------------
1237!
1238 DO step=1,nvsteps
1239!
1240! Compute vertical diffusive flux. Notice that "FC" is assumed to
1241! be time invariant.
1242!
1243 DO j=jstr,jend
1244 DO k=1,n(ng)-1
1245 DO i=istru,iend
1246 fs(i,k)=fc(i,j,k)*(awrk(i,j,k+1,nold)- &
1247 & awrk(i,j,k ,nold))
1248# ifdef MASKING
1249 fs(i,k)=fs(i,k)*umask(i,j)
1250# endif
1251 END DO
1252 END DO
1253 DO i=istru,iend
1254 fs(i,0)=0.0_r8
1255 fs(i,n(ng))=0.0_r8
1256 END DO
1257!
1258! Time-step vertical diffusive term. Notice that "oHz" is assumed to
1259! be time invariant.
1260!
1261 DO k=1,n(ng)
1262 DO i=istru,iend
1263 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1264 & ohz(i,j,k)*(fs(i,k )- &
1265 & fs(i,k-1))
1266 END DO
1267 END DO
1268 END DO
1269!
1270! Update integration indices.
1271!
1272 nsav=nold
1273 nold=nnew
1274 nnew=nsav
1275 END DO
1276# endif
1277# endif
1278!
1279!-----------------------------------------------------------------------
1280! Load convolved solution.
1281!-----------------------------------------------------------------------
1282!
1283 DO k=1,n(ng)
1284 DO j=jstr,jend
1285 DO i=istru,iend
1286 a(i,j,k)=awrk(i,j,k,nold)
1287 END DO
1288 END DO
1289 END DO
1290 CALL dabc_u3d_tile (ng, tile, &
1291 & lbi, ubi, lbj, ubj, lbk, ubk, &
1292 & a)
1293# ifdef DISTRIBUTE
1294 CALL mp_exchange3d (ng, tile, model, 1, &
1295 & lbi, ubi, lbj, ubj, lbk, ubk, &
1296 & nghost, &
1297 & ewperiodic(ng), nsperiodic(ng), &
1298 & a)
1299# endif
1300
1301 RETURN
1302 END SUBROUTINE conv_u3d_tile
1303!
1304!***********************************************************************
1305 SUBROUTINE conv_v3d_tile (ng, tile, model, &
1306 & LBi, UBi, LBj, UBj, LBk, UBk, &
1307 & IminS, ImaxS, JminS, JmaxS, &
1308 & Nghost, NHsteps, NVsteps, &
1309 & DTsizeH, DTsizeV, &
1310 & Kh, Kv, &
1311 & pm, pn, &
1312# ifdef GEOPOTENTIAL_HCONV
1313 & on_p, om_r, &
1314# else
1315 & pmon_p, pnom_r, &
1316# endif
1317# ifdef MASKING
1318# ifdef GEOPOTENTIAL_HCONV
1319 & pmask, rmask, umask, vmask, &
1320# else
1321 & vmask, pmask, &
1322# endif
1323# endif
1324 & Hz, z_r, &
1325 & A)
1326!***********************************************************************
1327!
1328 USE mod_param
1329 USE mod_scalars
1330!
1331 USE bc_3d_mod, ONLY: dabc_v3d_tile
1332# ifdef DISTRIBUTE
1333 USE mp_exchange_mod, ONLY : mp_exchange3d
1334# endif
1335!
1336! Imported variable declarations.
1337!
1338 integer, intent(in) :: ng, tile, model
1339 integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk
1340 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1341 integer, intent(in) :: Nghost, NHsteps, NVsteps
1342
1343 real(r8), intent(in) :: DTsizeH, DTsizeV
1344!
1345# ifdef ASSUMED_SHAPE
1346 real(r8), intent(in) :: pm(LBi:,LBj:)
1347 real(r8), intent(in) :: pn(LBi:,LBj:)
1348# ifdef GEOPOTENTIAL_HCONV
1349 real(r8), intent(in) :: on_p(LBi:,LBj:)
1350 real(r8), intent(in) :: om_r(LBi:,LBj:)
1351# else
1352 real(r8), intent(in) :: pmon_p(LBi:,LBj:)
1353 real(r8), intent(in) :: pnom_r(LBi:,LBj:)
1354# endif
1355# ifdef MASKING
1356# ifdef GEOPOTENTIAL_HCONV
1357 real(r8), intent(in) :: pmask(LBi:,LBj:)
1358 real(r8), intent(in) :: rmask(LBi:,LBj:)
1359 real(r8), intent(in) :: umask(LBi:,LBj:)
1360 real(r8), intent(in) :: vmask(LBi:,LBj:)
1361# else
1362 real(r8), intent(in) :: vmask(LBi:,LBj:)
1363 real(r8), intent(in) :: pmask(LBi:,LBj:)
1364# endif
1365# endif
1366 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
1367 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
1368
1369 real(r8), intent(in) :: Kh(LBi:,LBj:)
1370 real(r8), intent(in) :: Kv(LBi:,LBj:,0:)
1371 real(r8), intent(inout) :: A(LBi:,LBj:,LBk:)
1372# else
1373 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
1374 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
1375# ifdef GEOPOTENTIAL_HCONV
1376 real(r8), intent(in) :: on_p(LBi:UBi,LBj:UBj)
1377 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
1378# else
1379 real(r8), intent(in) :: pmon_p(LBi:UBi,LBj:UBj)
1380 real(r8), intent(in) :: pnom_r(LBi:UBi,LBj:UBj)
1381# endif
1382# ifdef MASKING
1383# ifdef GEOPOTENTIAL_HCONV
1384 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
1385 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1386 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1387 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1388# else
1389 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1390 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
1391# endif
1392# endif
1393 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
1394 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
1395
1396 real(r8), intent(in) :: Kh(LBi:UBi,LBj:UBj)
1397 real(r8), intent(in) :: Kv(LBi:UBi,LBj:UBj,0:UBk)
1398 real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj,LBk:UBk)
1399# endif
1400!
1401! Local variable declarations.
1402!
1403 integer :: Nnew, Nold, Nsav, i, j, k, k1, k2, step
1404
1405 real(r8) :: cff, cff1, cff2, cff3, cff4
1406
1407 real(r8), dimension(LBi:UBi,LBj:UBj,LBk:UBk,2) :: Awrk
1408
1409 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
1410 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
1411 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hfac
1412# ifdef VCONVOLUTION
1413# ifndef SPLINES_VCONV
1414 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: FC
1415# endif
1416# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
1417 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
1418# endif
1419# if defined IMPLICIT_VCONV || defined SPLINES_VCONV
1420 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
1421 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
1422 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
1423# ifdef SPLINES_VCONV
1424 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
1425 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
1426# endif
1427# else
1428 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FS
1429# endif
1430# endif
1431# ifdef GEOPOTENTIAL_HCONV
1432 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZdx
1433 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: dZde
1434
1435 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FZ
1436 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdz
1437 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAdx
1438 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dAde
1439 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
1440 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
1441# endif
1442
1443# include "set_bounds.h"
1444!
1445!-----------------------------------------------------------------------
1446! Space convolution of the diffusion equation for a 3D state variable
1447! at V-points.
1448!-----------------------------------------------------------------------
1449!
1450! Compute metrics factors. Notice that "z_r" and "Hz" are assumed to
1451! be time invariant in the vertical convolution. Scratch array are
1452! used for efficiency.
1453!
1454 cff=dtsizeh*0.25_r8
1455 DO j=jstrv-1,jend+1
1456 DO i=istr-1,iend+1
1457 hfac(i,j)=cff*(pm(i,j-1)+pm(i,j))*(pn(i,j-1)+pn(i,j))
1458# ifdef VCONVOLUTION
1459# ifndef SPLINES_VCONV
1460 fc(i,j,n(ng))=0.0_r8
1461 DO k=1,n(ng)-1
1462# ifdef IMPLICIT_VCONV
1463 fc(i,j,k)=-dtsizev*(kv(i,j-1,k)+kv(i,j,k))/ &
1464 & (z_r(i,j-1,k+1)+z_r(i,j,k+1)- &
1465 & z_r(i,j-1,k )-z_r(i,j,k ))
1466# else
1467 fc(i,j,k)=dtsizev*(kv(i,j-1,k)+kv(i,j,k))/ &
1468 & (z_r(i,j-1,k+1)+z_r(i,j,k+1)- &
1469 & z_r(i,j-1,k )-z_r(i,j,k ))
1470# endif
1471 END DO
1472 fc(i,j,0)=0.0_r8
1473# endif
1474# if !defined IMPLICIT_VCONV || defined SPLINES_VCONV
1475 DO k=1,n(ng)
1476 ohz(i,j,k)=2.0_r8/(hz(i,j-1,k)+hz(i,j,k))
1477 END DO
1478# endif
1479# endif
1480 END DO
1481 END DO
1482!
1483! Set integration indices and initial conditions.
1484!
1485 nold=1
1486 nnew=2
1487 CALL dabc_v3d_tile (ng, tile, &
1488 & lbi, ubi, lbj, ubj, lbk, ubk, &
1489 & a)
1490# ifdef DISTRIBUTE
1491 CALL mp_exchange3d (ng, tile, model, 1, &
1492 & lbi, ubi, lbj, ubj, lbk, ubk, &
1493 & nghost, &
1494 & ewperiodic(ng), nsperiodic(ng), &
1495 & a)
1496# endif
1497 DO k=1,n(ng)
1498 DO j=jstrv-1,jend+1
1499 DO i=istr-1,iend+1
1500 awrk(i,j,k,nold)=a(i,j,k)
1501 END DO
1502 END DO
1503 END DO
1504!
1505!-----------------------------------------------------------------------
1506! Integrate horizontal diffusion equation.
1507!-----------------------------------------------------------------------
1508!
1509 DO step=1,nhsteps
1510
1511# ifdef GEOPOTENTIAL_HCONV
1512!
1513! Diffusion along geopotential surfaces: Compute horizontal and
1514! vertical gradients. Notice the recursive blocking sequence. The
1515! vertical placement of the gradients is:
1516!
1517! dAdx,dAde(:,:,k1) k rho-points
1518! dAdx,dAde(:,:,k2) k+1 rho-points
1519! FZ,dAdz(:,:,k1) k-1/2 W-points
1520! FZ,dAdz(:,:,k2) k+1/2 W-points
1521!
1522 k2=1
1523 k_loop : DO k=0,n(ng)
1524 k1=k2
1525 k2=3-k1
1526 IF (k.lt.n(ng)) THEN
1527 DO j=jstrv-1,jend
1528 DO i=istr,iend+1
1529 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
1530# ifdef MASKING
1531 cff=cff*umask(i,j)
1532# endif
1533 dzdx(i,j)=cff*(z_r(i ,j,k+1)- &
1534 & z_r(i-1,j,k+1))
1535 END DO
1536 END DO
1537!
1538 DO j=jstrv-1,jend+1
1539 DO i=istr,iend
1540 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
1541# ifdef MASKING
1542 cff=cff*vmask(i,j)
1543# endif
1544 dzde(i,j)=cff*(z_r(i,j ,k+1)- &
1545 & z_r(i,j-1,k+1))
1546 END DO
1547 END DO
1548!
1549 DO j=jstrv,jend
1550 DO i=istr,iend+1
1551 cff=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+ &
1552 & pm(i ,j-1)+pm(i ,j))
1553# ifdef MASKING
1554 dadx(i,j,k2)=cff* &
1555 & (awrk(i ,j,k+1,nold)*vmask(i ,j)- &
1556 & awrk(i-1,j,k+1,nold)*vmask(i-1,j))
1557 dadx(i,j,k2)=dadx(i,j,k2)*pmask(i,j)
1558# else
1559 dadx(i,j,k2)=cff*(awrk(i ,j,k+1,nold)- &
1560 & awrk(i-1,j,k+1,nold))
1561# endif
1562 dzdx_p(i,j,k2)=0.5_r8*(dzdx(i,j-1)+ &
1563 & dzdx(i,j ))
1564 END DO
1565 END DO
1566!
1567 DO j=jstrv-1,jend
1568 DO i=istr,iend
1569# ifdef MASKING
1570 dade(i,j,k2)=pn(i,j)* &
1571 & (awrk(i,j+1,k+1,nold)*vmask(i,j+1)- &
1572 & awrk(i,j ,k+1,nold)*vmask(i,j ))
1573 dade(i,j,k2)=dade(i,j,k2)*rmask(i,j)
1574# else
1575 dade(i,j,k2)=pn(i,j)*(awrk(i,j+1,k+1,nold)- &
1576 & awrk(i,j ,k+1,nold))
1577# endif
1578 dzde_r(i,j,k2)=0.5_r8*(dzde(i,j )+ &
1579 & dzde(i,j+1))
1580 END DO
1581 END DO
1582 END IF
1583!
1584 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1585 DO j=jstrv-1,jend+1
1586 DO i=istr-1,iend+1
1587 dadz(i,j,k2)=0.0_r8
1588 fz(i,j,k2)=0.0_r8
1589 END DO
1590 END DO
1591 ELSE
1592 DO j=jstrv-1,jend+1
1593 DO i=istr-1,iend+1
1594 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1595 dadz(i,j,k2)=cff*(awrk(i,j,k+1,nold)- &
1596 & awrk(i,j,k ,nold))
1597# ifdef MASKING
1598 dadz(i,j,k2)=dadz(i,j,k2)*vmask(i,j)
1599# endif
1600 END DO
1601 END DO
1602 END IF
1603!
1604! Compute components of the rotated A flux (A m3/s) along geopotential
1605! surfaces.
1606!
1607 IF (k.gt.0) THEN
1608 DO j=jstrv,jend
1609 DO i=istr,iend+1
1610 cff=0.0625_r8*(kh(i-1,j-1)+kh(i-1,j)+ &
1611 & kh(i ,j-1)+kh(i ,j))*on_p(i,j)
1612 cff1=min(dzdx_p(i,j,k1),0.0_r8)
1613 cff2=max(dzdx_p(i,j,k1),0.0_r8)
1614 fx(i,j)=cff* &
1615 & (hz(i-1,j-1,k)+hz(i-1,j,k)+ &
1616 & hz(i ,j-1,k)+hz(i ,j,k))* &
1617 & (dadx(i,j,k1)- &
1618 & 0.5_r8*(cff1*(dadz(i-1,j,k1)+ &
1619 & dadz(i ,j,k2))+ &
1620 & cff2*(dadz(i-1,j,k2)+ &
1621 & dadz(i ,j,k1))))
1622 END DO
1623 END DO
1624 DO j=jstrv-1,jend
1625 DO i=istr,iend
1626 cff=kh(i,j)*om_r(i,j)
1627 cff1=min(dzde_r(i,j,k1),0.0_r8)
1628 cff2=max(dzde_r(i,j,k1),0.0_r8)
1629 fe(i,j)=cff* &
1630 & hz(i,j,k)* &
1631 & (dade(i,j,k1)- &
1632 & 0.5_r8*(cff1*(dadz(i,j ,k1)+ &
1633 & dadz(i,j+1,k2))+ &
1634 & cff2*(dadz(i,j ,k2)+ &
1635 & dadz(i,j+1,k1))))
1636 END DO
1637 END DO
1638 IF (k.lt.n(ng)) THEN
1639 DO j=jstrv,jend
1640 DO i=istr,iend
1641 cff=0.5_r8*(kh(i,j-1)+kh(i,j))
1642 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
1643 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
1644 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
1645 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
1646 fz(i,j,k2)=cff* &
1647 & (cff1*(cff1*dadz(i,j,k2)-dadx(i ,j,k1))+ &
1648 & cff2*(cff2*dadz(i,j,k2)-dadx(i+1,j,k2))+ &
1649 & cff3*(cff3*dadz(i,j,k2)-dadx(i ,j,k2))+ &
1650 & cff4*(cff4*dadz(i,j,k2)-dadx(i+1,j,k1)))
1651 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
1652 cff2=min(dzde_r(i,j ,k2),0.0_r8)
1653 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
1654 cff4=max(dzde_r(i,j ,k1),0.0_r8)
1655 fz(i,j,k2)=fz(i,j,k2)+ &
1656 & cff* &
1657 & (cff1*(cff1*dadz(i,j,k2)-dade(i,j-1,k1))+ &
1658 & cff2*(cff2*dadz(i,j,k2)-dade(i,j ,k2))+ &
1659 & cff3*(cff3*dadz(i,j,k2)-dade(i,j-1,k2))+ &
1660 & cff4*(cff4*dadz(i,j,k2)-dade(i,j ,k1)))
1661 END DO
1662 END DO
1663 END IF
1664!
1665! Time-step harmonic, geopotential diffusion term (m Tunits).
1666!
1667 DO j=jstrv,jend
1668 DO i=istr,iend
1669 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1670 & hfac(i,j)* &
1671 & (fx(i+1,j)-fx(i,j )+ &
1672 & fe(i ,j)-fe(i,j-1))+ &
1673 & dtsizeh* &
1674 & (fz(i,j,k2)-fz(i,j,k1))
1675 END DO
1676 END DO
1677 END IF
1678 END DO k_loop
1679
1680# else
1681
1682!
1683! Compute XI- and ETA-components of diffusive flux.
1684!
1685 DO k=1,n(ng)
1686 DO j=jstrv,jend
1687 DO i=istr,iend+1
1688 fx(i,j)=pmon_p(i,j)*0.25_r8*(kh(i-1,j )+kh(i,j )+ &
1689 & kh(i-1,j-1)+kh(i,j-1))* &
1690 & (awrk(i,j,k,nold)-awrk(i-1,j,k,nold))
1691# ifdef MASKING
1692 fx(i,j)=fx(i,j)*pmask(i,j)
1693# endif
1694 END DO
1695 END DO
1696 DO j=jstrv-1,jend
1697 DO i=istr,iend
1698 fe(i,j)=pnom_r(i,j)*kh(i,j)* &
1699 & (awrk(i,j+1,k,nold)-awrk(i,j,k,nold))
1700 END DO
1701 END DO
1702!
1703! Time-step horizontal diffusion equation.
1704!
1705 DO j=jstrv,jend
1706 DO i=istr,iend
1707 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1708 & hfac(i,j)* &
1709 & (fx(i+1,j)-fx(i,j)+ &
1710 & fe(i,j)-fe(i,j-1))
1711 END DO
1712 END DO
1713 END DO
1714# endif
1715!
1716! Apply boundary conditions.
1717!
1718 CALL dabc_v3d_tile (ng, tile, &
1719 & lbi, ubi, lbj, ubj, lbk, ubk, &
1720 & awrk(:,:,:,nnew))
1721# ifdef DISTRIBUTE
1722 CALL mp_exchange3d (ng, tile, model, 1, &
1723 & lbi, ubi, lbj, ubj, lbk, ubk, &
1724 & nghost, &
1725 & ewperiodic(ng), nsperiodic(ng), &
1726 & awrk(:,:,:,nnew))
1727# endif
1728!
1729! Update integration indices.
1730!
1731 nsav=nold
1732 nold=nnew
1733 nnew=nsav
1734 END DO
1735
1736# ifdef VCONVOLUTION
1737# ifdef IMPLICIT_VCONV
1738# ifdef SPLINES_VCONV
1739!
1740!-----------------------------------------------------------------------
1741! Integrate vertical diffusion equation implicitly using parabolic
1742! splines.
1743!-----------------------------------------------------------------------
1744!
1745 DO step=1,nvsteps
1746!
1747! Use conservative, parabolic spline reconstruction of vertical
1748! diffusion derivatives. Then, time step vertical diffusion term
1749! implicitly.
1750!
1751 DO j=jstrv,jend
1752 DO k=1,n(ng)
1753 DO i=istr,iend
1754 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
1755 & hz(i,j ,k))
1756 END DO
1757 END DO
1758 cff1=1.0_r8/6.0_r8
1759 DO k=1,n(ng)-1
1760 DO i=istr,iend
1761 fc(i,k)=cff1*hzk(i,k )-dtsizev*kv(i,j,k-1)*ohz(i,j,k )
1762 cf(i,k)=cff1*hzk(i,k+1)-dtsizev*kv(i,j,k+1)*ohz(i,j,k+1)
1763 END DO
1764 END DO
1765 DO i=istr,iend
1766 cf(i,0)=0.0_r8
1767 dc(i,0)=0.0_r8
1768 END DO
1769!
1770! LU decomposition and forward substitution.
1771!
1772 cff1=1.0_r8/3.0_r8
1773 DO k=1,n(ng)-1
1774 DO i=istr,iend
1775 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1776 & dtsizev*kv(i,j,k)*(ohz(i,j,k)+ohz(i,j,k+1))
1777 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1778 cf(i,k)=cff*cf(i,k)
1779 dc(i,k)=cff*(awrk(i,j,k+1,nold)-awrk(i,j,k,nold)- &
1780 & fc(i,k)*dc(i,k-1))
1781 END DO
1782 END DO
1783!
1784! Backward substitution.
1785!
1786 DO i=istr,iend
1787 dc(i,n(ng))=0.0_r8
1788 END DO
1789 DO k=n(ng)-1,1,-1
1790 DO i=istr,iend
1791 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1792 END DO
1793 END DO
1794!
1795 DO k=1,n(ng)
1796 DO i=istr,iend
1797 dc(i,k)=dc(i,k)*kv(i,j,k)
1798 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1799 & dtsizev*ohz(i,j,k)* &
1800 & (dc(i,k)-dc(i,k-1))
1801 END DO
1802 END DO
1803 END DO
1804!
1805! Update integration indices.
1806!
1807 nsav=nold
1808 nold=nnew
1809 nnew=nsav
1810 END DO
1811# else
1812!
1813!-----------------------------------------------------------------------
1814! Integerate vertical diffusion equation implicitly.
1815!-----------------------------------------------------------------------
1816!
1817 DO step=1,nvsteps
1818!
1819! Compute diagonal matrix coefficients BC and load right-hand-side
1820! terms for the diffusion equation into DC.
1821!
1822 DO j=jstrv,jend
1823 DO k=1,n(ng)
1824 DO i=istr,iend
1825 cff=0.5_r8*(hz(i,j-1,k)+hz(i,j,k))
1826 bc(i,k)=cff-fc(i,j,k)-fc(i,j,k-1)
1827 dc(i,k)=awrk(i,j,k,nold)*cff
1828 END DO
1829 END DO
1830!
1831! Solve the tridiagonal system.
1832!
1833 DO i=istr,iend
1834 cff=1.0_r8/bc(i,1)
1835 cf(i,1)=cff*fc(i,j,1)
1836 dc(i,1)=cff*dc(i,1)
1837 END DO
1838 DO k=2,n(ng)-1
1839 DO i=istr,iend
1840 cff=1.0_r8/(bc(i,k)-fc(i,j,k-1)*cf(i,k-1))
1841 cf(i,k)=cff*fc(i,j,k)
1842 dc(i,k)=cff*(dc(i,k)-fc(i,j,k-1)*dc(i,k-1))
1843 END DO
1844 END DO
1845!
1846! Compute new solution by back substitution.
1847!
1848 DO i=istr,iend
1849 dc(i,n(ng))=(dc(i,n(ng))- &
1850 & fc(i,j,n(ng)-1)*dc(i,n(ng)-1))/ &
1851 & (bc(i,n(ng))-fc(i,j,n(ng)-1)*cf(i,n(ng)-1))
1852 awrk(i,j,n(ng),nnew)=dc(i,n(ng))
1853# ifdef MASKING
1854 awrk(i,j,n(ng),nnew)=awrk(i,j,n(ng),nnew)*vmask(i,j)
1855# endif
1856 END DO
1857 DO k=n(ng)-1,1,-1
1858 DO i=istr,iend
1859 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1860 awrk(i,j,k,nnew)=dc(i,k)
1861# ifdef MASKING
1862 awrk(i,j,k,nnew)=awrk(i,j,k,nnew)*vmask(i,j)
1863# endif
1864 END DO
1865 END DO
1866 END DO
1867!
1868! Update integration indices.
1869!
1870 nsav=nold
1871 nold=nnew
1872 nnew=nsav
1873 END DO
1874# endif
1875# else
1876!
1877!-----------------------------------------------------------------------
1878! Integerate vertical diffusion equation explicitly.
1879!-----------------------------------------------------------------------
1880!
1881 DO step=1,nvsteps
1882!
1883! Compute vertical diffusive flux. Notice that "FC" is assumed to
1884! be time invariant.
1885!
1886 DO j=jstrv,jend
1887 DO k=1,n(ng)-1
1888 DO i=istr,iend
1889 fs(i,k)=fc(i,j,k)*(awrk(i,j,k+1,nold)- &
1890 & awrk(i,j,k ,nold))
1891# ifdef MASKING
1892 fs(i,k)=fs(i,k)*vmask(i,j)
1893# endif
1894 END DO
1895 END DO
1896 DO i=istr,iend
1897 fs(i,0)=0.0_r8
1898 fs(i,n(ng))=0.0_r8
1899 END DO
1900!
1901! Time-step vertical diffusive term. Notice that "oHz" is assumed to
1902! be time invariant.
1903!
1904 DO k=1,n(ng)
1905 DO i=istr,iend
1906 awrk(i,j,k,nnew)=awrk(i,j,k,nold)+ &
1907 & ohz(i,j,k)*(fs(i,k )- &
1908 & fs(i,k-1))
1909 END DO
1910 END DO
1911 END DO
1912!
1913! Update integration indices.
1914!
1915 nsav=nold
1916 nold=nnew
1917 nnew=nsav
1918 END DO
1919# endif
1920# endif
1921!
1922!-----------------------------------------------------------------------
1923! Load convolved solution.
1924!-----------------------------------------------------------------------
1925!
1926 DO k=1,n(ng)
1927 DO j=jstrv,jend
1928 DO i=istr,iend
1929 a(i,j,k)=awrk(i,j,k,nold)
1930 END DO
1931 END DO
1932 END DO
1933 CALL dabc_v3d_tile (ng, tile, &
1934 & lbi, ubi, lbj, ubj, lbk, ubk, &
1935 & a)
1936# ifdef DISTRIBUTE
1937 CALL mp_exchange3d (ng, tile, model, 1, &
1938 & lbi, ubi, lbj, ubj, lbk, ubk, &
1939 & nghost, &
1940 & ewperiodic(ng), nsperiodic(ng), &
1941 & a)
1942# endif
1943
1944 RETURN
1945 END SUBROUTINE conv_v3d_tile
1946#endif
1947 END MODULE 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
conv_3d_mod
Definition conv_3d.F:3
conv_3d_mod::conv_u3d_tile
subroutine 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, a)
Definition conv_3d.F:682
conv_3d_mod::conv_v3d_tile
subroutine 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, a)
Definition conv_3d.F:1326
conv_3d_mod::conv_r3d_tile
subroutine 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, a)
Definition conv_3d.F:85
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