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

Last change on this file since 4218 was 4182, checked in by scharf, 5 years ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
File size: 20.9 KB
Line 
1!> @file tridia_solver_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: tridia_solver_mod.f90 4182 2019-08-22 15:20:23Z knoop $
27! Corrected "Former revisions" section
28!
29! 3761 2019-02-25 15:31:42Z raasch
30! OpenACC modification to prevent compiler warning about unused variable
31!
32! 3690 2019-01-22 22:56:42Z knoop
33! OpenACC port for SPEC
34!
35! 1212 2013-08-15 08:46:27Z raasch
36! Initial revision.
37! Routines have been moved to seperate module from former file poisfft to here.
38! The tridiagonal matrix coefficients of array tri are calculated only once at
39! the beginning, i.e. routine split is called within tridia_init.
40!
41!
42! Description:
43! ------------
44!> solves the linear system of equations:
45!>
46!> -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
47!>   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
48!> 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)
49!>
50!> by using the Thomas algorithm
51!------------------------------------------------------------------------------!
52
53#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
54
55 MODULE tridia_solver
56 
57
58    USE basic_constants_and_equations_mod,                                     &
59        ONLY:  pi
60
61    USE indices,                                                               &
62        ONLY:  nx, ny, nz
63
64    USE kinds
65
66    USE transpose_indices,                                                     &
67        ONLY:  nxl_z, nyn_z, nxr_z, nys_z
68
69    IMPLICIT NONE
70
71    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddzuw !<
72
73    PRIVATE
74
75    INTERFACE tridia_substi
76       MODULE PROCEDURE tridia_substi
77    END INTERFACE tridia_substi
78
79    INTERFACE tridia_substi_overlap
80       MODULE PROCEDURE tridia_substi_overlap
81    END INTERFACE tridia_substi_overlap
82
83    PUBLIC  tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
84
85 CONTAINS
86
87
88!------------------------------------------------------------------------------!
89! Description:
90! ------------
91!> @todo Missing subroutine description.
92!------------------------------------------------------------------------------!
93    SUBROUTINE tridia_init
94
95       USE arrays_3d,                                                          &
96           ONLY:  ddzu_pres, ddzw, rho_air_zw
97
98#if defined( _OPENACC )
99       USE arrays_3d,                                                          &
100           ONLY:  tri
101#endif
102
103       IMPLICIT NONE
104
105       INTEGER(iwp) ::  k !<
106
107       ALLOCATE( ddzuw(0:nz-1,3) )
108
109       DO  k = 0, nz-1
110          ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
111          ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
112          ddzuw(k,3) = -1.0_wp * &
113                       ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) +        &
114                         ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
115       ENDDO
116!
117!--    Calculate constant coefficients of the tridiagonal matrix
118       CALL maketri
119       CALL split
120
121#if __acc_fft_device
122       !$ACC ENTER DATA &
123       !$ACC COPYIN(ddzuw(0:nz-1,1:3)) &
124       !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2))
125#endif
126
127    END SUBROUTINE tridia_init
128
129
130!------------------------------------------------------------------------------!
131! Description:
132! ------------
133!> Computes the i- and j-dependent component of the matrix
134!> Provide the constant coefficients of the tridiagonal matrix for solution
135!> of the Poisson equation in Fourier space.
136!> The coefficients are computed following the method of
137!> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
138!> Siano's original version by discretizing the Poisson equation,
139!> before it is Fourier-transformed.
140!------------------------------------------------------------------------------!
141    SUBROUTINE maketri
142
143
144          USE arrays_3d,                                                       &
145              ONLY:  tric, rho_air
146
147          USE control_parameters,                                              &
148              ONLY:  ibc_p_b, ibc_p_t
149
150          USE grid_variables,                                                  &
151              ONLY:  dx, dy
152
153
154          IMPLICIT NONE
155
156          INTEGER(iwp) ::  i    !<
157          INTEGER(iwp) ::  j    !<
158          INTEGER(iwp) ::  k    !<
159          INTEGER(iwp) ::  nnxh !<
160          INTEGER(iwp) ::  nnyh !<
161
162          REAL(wp)    ::  ll(nxl_z:nxr_z,nys_z:nyn_z) !<
163
164
165          nnxh = ( nx + 1 ) / 2
166          nnyh = ( ny + 1 ) / 2
167
168          DO  j = nys_z, nyn_z
169             DO  i = nxl_z, nxr_z
170                IF ( j >= 0  .AND.  j <= nnyh )  THEN
171                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
172                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
173                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
174                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
175                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
176                   ELSE
177                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
178                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
179                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
180                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
181                   ENDIF
182                ELSE
183                   IF ( i >= 0  .AND.  i <= nnxh )  THEN
184                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
185                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
186                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
187                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
188                   ELSE
189                      ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
190                                            REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
191                                2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
192                                            REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
193                   ENDIF
194                ENDIF
195             ENDDO
196          ENDDO
197
198          DO  k = 0, nz-1
199             DO  j = nys_z, nyn_z
200                DO  i = nxl_z, nxr_z
201                   tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
202                ENDDO
203             ENDDO
204          ENDDO
205
206          IF ( ibc_p_b == 1 )  THEN
207             DO  j = nys_z, nyn_z
208                DO  i = nxl_z, nxr_z
209                   tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
210                ENDDO
211             ENDDO
212          ENDIF
213          IF ( ibc_p_t == 1 )  THEN
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          ENDIF
220
221    END SUBROUTINE maketri
222
223
224!------------------------------------------------------------------------------!
225! Description:
226! ------------
227!> Substitution (Forward and Backward) (Thomas algorithm)
228!------------------------------------------------------------------------------!
229    SUBROUTINE tridia_substi( ar )
230
231
232          USE arrays_3d,                                                       & 
233              ONLY:  tri
234
235          USE control_parameters,                                              &
236              ONLY:  ibc_p_b, ibc_p_t
237
238          IMPLICIT NONE
239
240          INTEGER(iwp) ::  i !<
241          INTEGER(iwp) ::  j !<
242          INTEGER(iwp) ::  k !<
243
244          REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
245
246          REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1 !<
247#if __acc_fft_device
248          !$ACC DECLARE CREATE(ar1)
249#endif
250
251!
252!--       Forward substitution
253#if __acc_fft_device
254          !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k)
255#endif
256          DO  k = 0, nz - 1
257#if __acc_fft_device
258             !$ACC LOOP COLLAPSE(2)
259#endif
260             DO  j = nys_z, nyn_z
261                DO  i = nxl_z, nxr_z
262
263                   IF ( k == 0 )  THEN
264                      ar1(i,j,k) = ar(i,j,k+1)
265                   ELSE
266                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
267                   ENDIF
268
269                ENDDO
270             ENDDO
271          ENDDO
272#if __acc_fft_device
273          !$ACC END PARALLEL
274#endif
275
276!
277!--       Backward substitution
278!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
279!--       by zero appearing if the pressure bc is set to neumann at the top of
280!--       the model domain.
281#if __acc_fft_device
282          !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k)
283#endif
284          DO  k = nz-1, 0, -1
285#if __acc_fft_device
286             !$ACC LOOP COLLAPSE(2)
287#endif
288             DO  j = nys_z, nyn_z
289                DO  i = nxl_z, nxr_z
290
291                   IF ( k == nz-1 )  THEN
292                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
293                   ELSE
294                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
295                              / tri(i,j,k,1)
296                   ENDIF
297                ENDDO
298             ENDDO
299          ENDDO
300#if __acc_fft_device
301          !$ACC END PARALLEL
302#endif
303
304!
305!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
306!--       The respective values of ar should be zero at all k-levels if
307!--       acceleration of horizontally averaged vertical velocity is zero.
308          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
309             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
310#if __acc_fft_device
311                !$ACC PARALLEL LOOP PRESENT(ar)
312#endif
313                DO  k = 1, nz
314                   ar(nxl_z,nys_z,k) = 0.0_wp
315                ENDDO
316             ENDIF
317          ENDIF
318
319    END SUBROUTINE tridia_substi
320
321
322!------------------------------------------------------------------------------!
323! Description:
324! ------------
325!> Substitution (Forward and Backward) (Thomas algorithm)
326!------------------------------------------------------------------------------!
327    SUBROUTINE tridia_substi_overlap( ar, jj )
328
329
330          USE arrays_3d,                                                       &
331              ONLY:  tri
332
333          USE control_parameters,                                              &
334              ONLY:  ibc_p_b, ibc_p_t
335
336          IMPLICIT NONE
337
338          INTEGER(iwp) ::  i  !<
339          INTEGER(iwp) ::  j  !<
340          INTEGER(iwp) ::  jj !<
341          INTEGER(iwp) ::  k  !<
342
343          REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
344
345          REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) ::  ar1 !<
346
347!
348!--       Forward substitution
349          DO  k = 0, nz - 1
350             DO  j = nys_z, nyn_z
351                DO  i = nxl_z, nxr_z
352
353                   IF ( k == 0 )  THEN
354                      ar1(i,j,k) = ar(i,j,k+1)
355                   ELSE
356                      ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
357                   ENDIF
358
359                ENDDO
360             ENDDO
361          ENDDO
362
363!
364!--       Backward substitution
365!--       Note, the 1.0E-20 in the denominator is due to avoid divisions
366!--       by zero appearing if the pressure bc is set to neumann at the top of
367!--       the model domain.
368          DO  k = nz-1, 0, -1
369             DO  j = nys_z, nyn_z
370                DO  i = nxl_z, nxr_z
371
372                   IF ( k == nz-1 )  THEN
373                      ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
374                   ELSE
375                      ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
376                              / tri(i,jj,k,1)
377                   ENDIF
378                ENDDO
379             ENDDO
380          ENDDO
381
382!
383!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
384!--       The respective values of ar should be zero at all k-levels if
385!--       acceleration of horizontally averaged vertical velocity is zero.
386          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
387             IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
388                DO  k = 1, nz
389                   ar(nxl_z,nys_z,k) = 0.0_wp
390                ENDDO
391             ENDIF
392          ENDIF
393
394    END SUBROUTINE tridia_substi_overlap
395
396
397!------------------------------------------------------------------------------!
398! Description:
399! ------------
400!> Splitting of the tridiagonal matrix (Thomas algorithm)
401!------------------------------------------------------------------------------!
402    SUBROUTINE split
403
404
405          USE arrays_3d,                                                       & 
406              ONLY:  tri, tric
407
408          IMPLICIT NONE
409
410          INTEGER(iwp) ::  i !<
411          INTEGER(iwp) ::  j !<
412          INTEGER(iwp) ::  k !<
413!
414!--       Splitting
415          DO  j = nys_z, nyn_z
416             DO  i = nxl_z, nxr_z
417                tri(i,j,0,1) = tric(i,j,0)
418             ENDDO
419          ENDDO
420
421          DO  k = 1, nz-1
422             DO  j = nys_z, nyn_z
423                DO  i = nxl_z, nxr_z
424                   tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
425                   tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
426                ENDDO
427             ENDDO
428          ENDDO
429
430    END SUBROUTINE split
431
432
433!------------------------------------------------------------------------------!
434! Description:
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    SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
448
449
450       USE arrays_3d,                                                          &
451           ONLY:  ddzu_pres, ddzw, rho_air, rho_air_zw
452
453       USE control_parameters,                                                 &
454           ONLY:  ibc_p_b, ibc_p_t
455
456       IMPLICIT NONE
457
458       INTEGER(iwp) ::  i                  !<
459       INTEGER(iwp) ::  j                  !<
460       INTEGER(iwp) ::  k                  !<
461       INTEGER(iwp) ::  nnyh               !<
462       INTEGER(iwp) ::  nx                 !<
463       INTEGER(iwp) ::  ny                 !<
464
465       REAL(wp)     ::  ddx2 !<
466       REAL(wp)     ::  ddy2 !<
467
468       REAL(wp), DIMENSION(0:nx,1:nz)     ::  ar         !<
469       REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
470
471
472       nnyh = ( ny + 1 ) / 2
473
474!
475!--    Define constant elements of the tridiagonal matrix.
476!--    The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
477!--    the exchanged loops create bank conflicts. The following directive
478!--    prohibits loop exchange and the loops perform much better.
479!CDIR NOLOOPCHG
480       DO  k = 0, nz-1
481          DO  i = 0,nx
482             tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
483             tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
484          ENDDO
485       ENDDO
486
487       IF ( j <= nnyh )  THEN
488          CALL maketri_1dd( j )
489       ELSE
490          CALL maketri_1dd( ny+1-j )
491       ENDIF
492
493       CALL split_1dd
494       CALL substi_1dd( ar, tri_for_1d )
495
496    CONTAINS
497
498
499!------------------------------------------------------------------------------!
500! Description:
501! ------------
502!> computes the i- and j-dependent component of the matrix
503!------------------------------------------------------------------------------!
504       SUBROUTINE maketri_1dd( j )
505
506          IMPLICIT NONE
507
508          INTEGER(iwp) ::  i    !<
509          INTEGER(iwp) ::  j    !<
510          INTEGER(iwp) ::  k    !<
511          INTEGER(iwp) ::  nnxh !<
512
513          REAL(wp)     ::  a !<
514          REAL(wp)     ::  c !<
515
516          REAL(wp), DIMENSION(0:nx) ::  l !<
517
518
519          nnxh = ( nx + 1 ) / 2
520!
521!--       Provide the tridiagonal matrix for solution of the Poisson equation in
522!--       Fourier space. The coefficients are computed following the method of
523!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
524!--       Siano's original version by discretizing the Poisson equation,
525!--       before it is Fourier-transformed
526          DO  i = 0, nx
527             IF ( i >= 0 .AND. i <= nnxh ) THEN
528                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
529                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
530                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
531                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
532             ELSE
533                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
534                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
535                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
536                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
537             ENDIF
538          ENDDO
539
540          DO  k = 0, nz-1
541             DO  i = 0, nx
542                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
543                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
544                tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
545             ENDDO
546          ENDDO
547          IF ( ibc_p_b == 1 )  THEN
548             DO  i = 0, nx
549                tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
550             ENDDO
551          ENDIF
552          IF ( ibc_p_t == 1 )  THEN
553             DO  i = 0, nx
554                tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
555             ENDDO
556          ENDIF
557
558       END SUBROUTINE maketri_1dd
559
560
561!------------------------------------------------------------------------------!
562! Description:
563! ------------
564!> Splitting of the tridiagonal matrix (Thomas algorithm)
565!------------------------------------------------------------------------------!
566       SUBROUTINE split_1dd
567
568          IMPLICIT NONE
569
570          INTEGER(iwp) ::  i !<
571          INTEGER(iwp) ::  k !<
572
573
574!
575!--       Splitting
576          DO  i = 0, nx
577             tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
578          ENDDO
579          DO  k = 1, nz-1
580             DO  i = 0, nx
581                tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
582                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)
583             ENDDO
584          ENDDO
585
586       END SUBROUTINE split_1dd
587
588
589!------------------------------------------------------------------------------!
590! Description:
591! ------------
592!> Substitution (Forward and Backward) (Thomas algorithm)
593!------------------------------------------------------------------------------!
594       SUBROUTINE substi_1dd( ar, tri_for_1d )
595
596
597          IMPLICIT NONE
598
599          INTEGER(iwp) ::  i !<
600          INTEGER(iwp) ::  k !<
601
602          REAL(wp), DIMENSION(0:nx,nz)       ::  ar         !<
603          REAL(wp), DIMENSION(0:nx,0:nz-1)   ::  ar1        !<
604          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
605
606!
607!--       Forward substitution
608          DO  i = 0, nx
609             ar1(i,0) = ar(i,1)
610          ENDDO
611          DO  k = 1, nz-1
612             DO  i = 0, nx
613                ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
614             ENDDO
615          ENDDO
616
617!
618!--       Backward substitution
619!--       Note, the add of 1.0E-20 in the denominator is due to avoid divisions
620!--       by zero appearing if the pressure bc is set to neumann at the top of
621!--       the model domain.
622          DO  i = 0, nx
623             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
624          ENDDO
625          DO  k = nz-2, 0, -1
626             DO  i = 0, nx
627                ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
628                            / tri_for_1d(4,i,k)
629             ENDDO
630          ENDDO
631
632!
633!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
634!--       The respective values of ar should be zero at all k-levels if
635!--       acceleration of horizontally averaged vertical velocity is zero.
636          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
637             IF ( j == 0 )  THEN
638                DO  k = 1, nz
639                   ar(0,k) = 0.0_wp
640                ENDDO
641             ENDIF
642          ENDIF
643
644       END SUBROUTINE substi_1dd
645
646    END SUBROUTINE tridia_1dd
647
648
649 END MODULE tridia_solver
Note: See TracBrowser for help on using the repository browser.