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

Last change on this file since 1322 was 1322, checked in by raasch, 8 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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