source: palm/trunk/SOURCE/tridia_solver.f90 @ 1320

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

  • Property svn:keywords set to Id
File size: 22.0 KB
Line 
1 MODULE tridia_solver
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! ONLY-attribute added to USE-statements,
23! kind-parameters added to all INTEGER and REAL declaration statements,
24! kinds are defined in new module kinds,
25! old module precision_kind is removed,
26! revision history before 2012 removed,
27! comment fields (!:) to be used for variable explanations added to
28! all variable declaration statements
29!
30! Former revisions:
31! -----------------
32! $Id: tridia_solver.f90 1320 2014-03-20 08:40:49Z raasch $
33!
34! 1257 2013-11-08 15:18:40Z raasch
35! openacc loop and loop vector clauses removed, declare create moved after
36! the FORTRAN declaration statement
37!
38! 1221 2013-09-10 08:59:13Z raasch
39! dummy argument tri in 1d-routines replaced by tri_for_1d because of name
40! conflict with arry tri in module arrays_3d
41!
42! 1216 2013-08-26 09:31:42Z raasch
43! +tridia_substi_overlap for handling overlapping fft / transposition
44!
45! 1212 2013-08-15 08:46:27Z raasch
46! Initial revision.
47! Routines have been moved to seperate module from former file poisfft to here.
48! The tridiagonal matrix coefficients of array tri are calculated only once at
49! the beginning, i.e. routine split is called within tridia_init.
50!
51!
52! Description:
53! ------------
54! solves the linear system of equations:
55!
56! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
57!   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
58! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
59!
60! by using the Thomas algorithm
61!------------------------------------------------------------------------------!
62
63    USE indices,                                                               &
64        ONLY:  nx, ny, nz
65
66    USE kinds
67
68    USE transpose_indices,                                                     &
69        ONLY:  nxl_z, nyn_z, nxr_z, nys_z
70
71    IMPLICIT NONE
72
73    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddzuw !:
74
75    PRIVATE
76
77    INTERFACE tridia_substi
78       MODULE PROCEDURE tridia_substi
79    END INTERFACE tridia_substi
80
81    INTERFACE tridia_substi_overlap
82       MODULE PROCEDURE tridia_substi_overlap
83    END INTERFACE tridia_substi_overlap
84
85    PUBLIC  tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
86
87 CONTAINS
88
89
90    SUBROUTINE tridia_init
91
92       USE arrays_3d,                                                          &
93           ONLY:  ddzu_pres, ddzw
94
95       USE kinds
96
97       IMPLICIT NONE
98
99       INTEGER(iwp) ::  k !:
100
101       ALLOCATE( ddzuw(0:nz-1,3) )
102
103       DO  k = 0, nz-1
104          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
105          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
106          ddzuw(k,3) = -1.0 * &
107                       ( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
108       ENDDO
109!
110!--    Calculate constant coefficients of the tridiagonal matrix
111#if ! defined ( __check )
112       CALL maketri
113       CALL split
114#endif
115
116    END SUBROUTINE tridia_init
117
118
119    SUBROUTINE maketri
120
121!------------------------------------------------------------------------------!
122! Computes the i- and j-dependent component of the matrix
123!------------------------------------------------------------------------------!
124
125          USE arrays_3d,                                                       &
126              ONLY:  tric
127
128          USE constants,                                                       &
129              ONLY:  pi
130
131          USE control_parameters,                                              &
132              ONLY:  ibc_p_b, ibc_p_t
133
134          USE grid_variables,                                                  &
135              ONLY:  dx, dy
136
137
138          USE kinds
139
140          IMPLICIT NONE
141
142          INTEGER(iwp) ::  i    !:
143          INTEGER(iwp) ::  j    !:
144          INTEGER(iwp) ::  k    !:
145          INTEGER(iwp) ::  nnxh !:
146          INTEGER(iwp) ::  nnyh !:
147
148          REAL(wp)    ::  ll(nxl_z:nxr_z,nys_z:nyn_z) !:
149          !$acc declare create( ll )
150
151
152          nnxh = ( nx + 1 ) / 2
153          nnyh = ( ny + 1 ) / 2
154
155!
156!--       Provide the constant coefficients of the tridiagonal matrix for solution
157!--       of the Poisson equation in Fourier space.
158!--       The coefficients are computed following the method of
159!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
160!--       Siano's original version by discretizing the Poisson equation,
161!--       before it is Fourier-transformed.
162
163          !$acc kernels present( tric )
164          DO  j = nys_z, nyn_z
165             DO  i = nxl_z, nxr_z
166                IF ( j >= 0  .AND.  j <= nnyh )  THEN
167                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
168                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
169                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
170                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
171                                                  REAL( ny+1 ) ) ) / ( dy * dy )
172                   ELSE
173                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
174                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
175                                2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
176                                                  REAL( ny+1 ) ) ) / ( dy * dy )
177                   ENDIF
178                ELSE
179                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
180                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
181                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
182                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
183                                                  REAL( ny+1 ) ) ) / ( dy * dy )
184                   ELSE
185                      ll(i,j) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
186                                                  REAL( nx+1 ) ) ) / ( dx * dx ) + &
187                                2.0 * ( 1.0 - COS( ( 2.0 * pi * ( ny+1-j ) ) / &
188                                                  REAL( ny+1 ) ) ) / ( dy * dy )
189                   ENDIF
190                ENDIF
191             ENDDO
192          ENDDO
193
194          DO  k = 0, nz-1
195             DO  j = nys_z, nyn_z
196                DO  i = nxl_z, nxr_z
197                   tric(i,j,k) = ddzuw(k,3) - ll(i,j)
198                ENDDO
199             ENDDO
200          ENDDO
201          !$acc end kernels
202
203          IF ( ibc_p_b == 1 )  THEN
204             !$acc kernels present( tric )
205             DO  j = nys_z, nyn_z
206                DO  i = nxl_z, nxr_z
207                   tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
208                ENDDO
209             ENDDO
210             !$acc end kernels
211          ENDIF
212          IF ( ibc_p_t == 1 )  THEN
213             !$acc kernels present( tric )
214             DO  j = nys_z, nyn_z
215                DO  i = nxl_z, nxr_z
216                   tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
217                ENDDO
218             ENDDO
219             !$acc end kernels
220          ENDIF
221
222    END SUBROUTINE maketri
223
224
225    SUBROUTINE tridia_substi( ar )
226
227!------------------------------------------------------------------------------!
228! Substitution (Forward and Backward) (Thomas algorithm)
229!------------------------------------------------------------------------------!
230
231          USE arrays_3d,                                                       & 
232              ONLY:  tri
233
234          USE control_parameters,                                              &
235              ONLY:  ibc_p_b, ibc_p_t
236
237          USE kinds
238
239          IMPLICIT NONE
240
241          INTEGER(iwp) ::  i !:
242          INTEGER(iwp) ::  j !:
243          INTEGER(iwp) ::  k !:
244
245          REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !:
246
247          REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1 !:
248          !$acc declare create( ar1 )
249
250!
251!--       Forward substitution
252          DO  k = 0, nz - 1
253             !$acc kernels present( ar, tri )
254             DO  j = nys_z, nyn_z
255                DO  i = nxl_z, nxr_z
256
257                   IF ( k == 0 )  THEN
258                      ar1(i,j,k) = ar(i,j,k+1)
259                   ELSE
260                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
261                   ENDIF
262
263                ENDDO
264             ENDDO
265             !$acc end kernels
266          ENDDO
267
268!
269!--       Backward substitution
270!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
271!--       by zero appearing if the pressure bc is set to neumann at the top of
272!--       the model domain.
273          DO  k = nz-1, 0, -1
274             !$acc kernels present( ar, tri )
275             DO  j = nys_z, nyn_z
276                DO  i = nxl_z, nxr_z
277
278                   IF ( k == nz-1 )  THEN
279                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20 )
280                   ELSE
281                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
282                              / tri(i,j,k,1)
283                   ENDIF
284                ENDDO
285             ENDDO
286             !$acc end kernels
287          ENDDO
288
289!
290!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
291!--       The respective values of ar should be zero at all k-levels if
292!--       acceleration of horizontally averaged vertical velocity is zero.
293          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
294             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
295                !$acc kernels loop present( ar )
296                DO  k = 1, nz
297                   ar(nxl_z,nys_z,k) = 0.0
298                ENDDO
299                !$acc end kernels loop
300             ENDIF
301          ENDIF
302
303    END SUBROUTINE tridia_substi
304
305
306    SUBROUTINE tridia_substi_overlap( ar, jj )
307
308!------------------------------------------------------------------------------!
309! Substitution (Forward and Backward) (Thomas algorithm)
310!------------------------------------------------------------------------------!
311
312          USE arrays_3d,                                                       &
313              ONLY:  tri
314
315          USE control_parameters,                                              &
316              ONLY:  ibc_p_b, ibc_p_t
317
318          USE kinds
319
320          IMPLICIT NONE
321
322          INTEGER(iwp) ::  i  !:
323          INTEGER(iwp) ::  j  !:
324          INTEGER(iwp) ::  jj !:
325          INTEGER(iwp) ::  k  !:
326
327          REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !:
328
329          !$acc declare create( ar1 )
330          REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) ::  ar1 !:
331
332!
333!--       Forward substitution
334          DO  k = 0, nz - 1
335             !$acc kernels present( ar, tri )
336             !$acc loop
337             DO  j = nys_z, nyn_z
338                DO  i = nxl_z, nxr_z
339
340                   IF ( k == 0 )  THEN
341                      ar1(i,j,k) = ar(i,j,k+1)
342                   ELSE
343                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
344                   ENDIF
345
346                ENDDO
347             ENDDO
348             !$acc end kernels
349          ENDDO
350
351!
352!--       Backward substitution
353!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
354!--       by zero appearing if the pressure bc is set to neumann at the top of
355!--       the model domain.
356          DO  k = nz-1, 0, -1
357             !$acc kernels present( ar, tri )
358             !$acc loop
359             DO  j = nys_z, nyn_z
360                DO  i = nxl_z, nxr_z
361
362                   IF ( k == nz-1 )  THEN
363                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20 )
364                   ELSE
365                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
366                              / tri(i,jj,k,1)
367                   ENDIF
368                ENDDO
369             ENDDO
370             !$acc end kernels
371          ENDDO
372
373!
374!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
375!--       The respective values of ar should be zero at all k-levels if
376!--       acceleration of horizontally averaged vertical velocity is zero.
377          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
378             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
379                !$acc kernels loop present( ar )
380                DO  k = 1, nz
381                   ar(nxl_z,nys_z,k) = 0.0
382                ENDDO
383             ENDIF
384          ENDIF
385
386    END SUBROUTINE tridia_substi_overlap
387
388
389    SUBROUTINE split
390
391!------------------------------------------------------------------------------!
392! Splitting of the tridiagonal matrix (Thomas algorithm)
393!------------------------------------------------------------------------------!
394
395          USE arrays_3d,                                                       & 
396              ONLY:  tri, tric
397
398          USE kinds
399
400          IMPLICIT NONE
401
402          INTEGER(iwp) ::  i !:
403          INTEGER(iwp) ::  j !:
404          INTEGER(iwp) ::  k !:
405!
406!--       Splitting
407          !$acc kernels present( tri, tric )
408          !$acc loop
409          DO  j = nys_z, nyn_z
410             !$acc loop vector( 32 )
411             DO  i = nxl_z, nxr_z
412                tri(i,j,0,1) = tric(i,j,0)
413             ENDDO
414          ENDDO
415          !$acc end kernels
416
417          DO  k = 1, nz-1
418             !$acc kernels present( tri, tric )
419             !$acc loop
420             DO  j = nys_z, nyn_z
421                !$acc loop vector( 32 )
422                DO  i = nxl_z, nxr_z
423                   tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
424                   tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
425                ENDDO
426             ENDDO
427             !$acc end kernels
428          ENDDO
429
430    END SUBROUTINE split
431
432
433    SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
434
435!------------------------------------------------------------------------------!
436! Solves the linear system of equations for a 1d-decomposition along x (see
437! tridia)
438!
439! Attention:  when using the intel compilers older than 12.0, array tri must
440!             be passed as an argument to the contained subroutines. Otherwise
441!             addres faults will occur. This feature can be activated with
442!             cpp-switch __intel11
443!             On NEC, tri should not be passed (except for routine substi_1dd)
444!             because this causes very bad performance.
445!------------------------------------------------------------------------------!
446
447       USE arrays_3d,                                                          &
448           ONLY:  ddzu_pres, ddzw
449
450       USE control_parameters,                                                 &
451           ONLY:  ibc_p_b, ibc_p_t
452
453       USE kinds
454
455       IMPLICIT NONE
456
457       INTEGER(iwp) ::  i                  !:
458       INTEGER(iwp) ::  j                  !:
459       INTEGER(iwp) ::  k                  !:
460       INTEGER(iwp) ::  nnyh               !:
461       INTEGER(iwp) ::  nx                 !:
462       INTEGER(iwp) ::  ny                 !:
463       INTEGER(iwp) ::  omp_get_thread_num !:
464       INTEGER(iwp) ::  tn                 !:
465
466       REAL(wp)     ::  ddx2 !:
467       REAL(wp)     ::  ddy2 !:
468
469       REAL(wp), DIMENSION(0:nx,1:nz)     ::  ar         !:
470       REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !:
471
472
473       nnyh = ( ny + 1 ) / 2
474
475!
476!--    Define constant elements of the tridiagonal matrix.
477!--    The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
478!--    the exchanged loops create bank conflicts. The following directive
479!--    prohibits loop exchange and the loops perform much better.
480!       tn = omp_get_thread_num()
481!       WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
482!       CALL local_flush( 120+tn )
483!CDIR NOLOOPCHG
484       DO  k = 0, nz-1
485          DO  i = 0,nx
486             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
487             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
488          ENDDO
489       ENDDO
490!       WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop   thread=', omp_get_thread_num()
491!       CALL local_flush( 120+tn )
492
493       IF ( j <= nnyh )  THEN
494#if defined( __intel11 )
495          CALL maketri_1dd( j, tri_for_1d )
496#else
497          CALL maketri_1dd( j )
498#endif
499       ELSE
500#if defined( __intel11 )
501          CALL maketri_1dd( ny+1-j, tri_for_1d )
502#else
503          CALL maketri_1dd( ny+1-j )
504#endif
505       ENDIF
506#if defined( __intel11 )
507       CALL split_1dd( tri_for_1d )
508#else
509       CALL split_1dd
510#endif
511       CALL substi_1dd( ar, tri_for_1d )
512
513    CONTAINS
514
515#if defined( __intel11 )
516       SUBROUTINE maketri_1dd( j, tri_for_1d )
517#else
518       SUBROUTINE maketri_1dd( j )
519#endif
520
521!------------------------------------------------------------------------------!
522! computes the i- and j-dependent component of the matrix
523!------------------------------------------------------------------------------!
524
525          USE constants,                                                       &
526              ONLY:  pi
527
528          USE kinds
529
530          IMPLICIT NONE
531
532          INTEGER(iwp) ::  i    !:
533          INTEGER(iwp) ::  j    !:
534          INTEGER(iwp) ::  k    !:
535          INTEGER(iwp) ::  nnxh !:
536
537          REAL(wp)     ::  a !:
538          REAL(wp)     ::  c !:
539
540          REAL(wp), DIMENSION(0:nx) ::  l !:
541
542#if defined( __intel11 )
543          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !:
544#endif
545
546
547          nnxh = ( nx + 1 ) / 2
548!
549!--       Provide the tridiagonal matrix for solution of the Poisson equation in
550!--       Fourier space. The coefficients are computed following the method of
551!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
552!--       Siano's original version by discretizing the Poisson equation,
553!--       before it is Fourier-transformed
554          DO  i = 0, nx
555             IF ( i >= 0 .AND. i <= nnxh ) THEN
556                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
557                                         REAL( nx+1 ) ) ) * ddx2 + &
558                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
559                                         REAL( ny+1 ) ) ) * ddy2
560             ELSE
561                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
562                                         REAL( nx+1 ) ) ) * ddx2 + &
563                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
564                                         REAL( ny+1 ) ) ) * ddy2
565             ENDIF
566          ENDDO
567
568          DO  k = 0, nz-1
569             DO  i = 0, nx
570                a = -1.0 * ddzu_pres(k+2) * ddzw(k+1)
571                c = -1.0 * ddzu_pres(k+1) * ddzw(k+1)
572                tri_for_1d(1,i,k) = a + c - l(i)
573             ENDDO
574          ENDDO
575          IF ( ibc_p_b == 1 )  THEN
576             DO  i = 0, nx
577                tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
578             ENDDO
579          ENDIF
580          IF ( ibc_p_t == 1 )  THEN
581             DO  i = 0, nx
582                tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
583             ENDDO
584          ENDIF
585
586       END SUBROUTINE maketri_1dd
587
588
589#if defined( __intel11 )
590       SUBROUTINE split_1dd( tri_for_1d )
591#else
592       SUBROUTINE split_1dd
593#endif
594
595!------------------------------------------------------------------------------!
596! Splitting of the tridiagonal matrix (Thomas algorithm)
597!------------------------------------------------------------------------------!
598
599          IMPLICIT NONE
600
601          INTEGER(iwp) ::  i !:
602          INTEGER(iwp) ::  k !:
603
604#if defined( __intel11 )
605          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !:
606#endif
607
608
609!
610!--       Splitting
611          DO  i = 0, nx
612             tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
613          ENDDO
614          DO  k = 1, nz-1
615             DO  i = 0, nx
616                tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
617                tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k)
618             ENDDO
619          ENDDO
620
621       END SUBROUTINE split_1dd
622
623
624       SUBROUTINE substi_1dd( ar, tri_for_1d )
625
626!------------------------------------------------------------------------------!
627! Substitution (Forward and Backward) (Thomas algorithm)
628!------------------------------------------------------------------------------!
629
630          IMPLICIT NONE
631
632          INTEGER(iwp) ::  i !:
633          INTEGER(iwp) ::  k !:
634
635          REAL(wp), DIMENSION(0:nx,nz)       ::  ar         !:
636          REAL(wp), DIMENSION(0:nx,0:nz-1)   ::  ar1        !:
637          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !:
638
639!
640!--       Forward substitution
641          DO  i = 0, nx
642             ar1(i,0) = ar(i,1)
643          ENDDO
644          DO  k = 1, nz-1
645             DO  i = 0, nx
646                ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
647             ENDDO
648          ENDDO
649
650!
651!--       Backward substitution
652!--       Note, the add of 1.0E-20 in the denominator is due to avoid divisions
653!--       by zero appearing if the pressure bc is set to neumann at the top of
654!--       the model domain.
655          DO  i = 0, nx
656             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20 )
657          ENDDO
658          DO  k = nz-2, 0, -1
659             DO  i = 0, nx
660                ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
661                            / tri_for_1d(4,i,k)
662             ENDDO
663          ENDDO
664
665!
666!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
667!--       The respective values of ar should be zero at all k-levels if
668!--       acceleration of horizontally averaged vertical velocity is zero.
669          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
670             IF ( j == 0 )  THEN
671                DO  k = 1, nz
672                   ar(0,k) = 0.0
673                ENDDO
674             ENDIF
675          ENDIF
676
677       END SUBROUTINE substi_1dd
678
679    END SUBROUTINE tridia_1dd
680
681
682 END MODULE tridia_solver
Note: See TracBrowser for help on using the repository browser.