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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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