source: palm/trunk/SOURCE/tridia_solver_mod.f90 @ 1851

Last change on this file since 1851 was 1851, checked in by maronga, 5 years ago

last commit documented

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