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

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

last commit documented / copyright update

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