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

Last change on this file since 4181 was 4181, checked in by raasch, 2 years ago

bugfix: define directive added, which has been deleted by mistake in r4180

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