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