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