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

Last change on this file since 4310 was 4236, checked in by knoop, 5 years ago

Added missing OpenMP directives within "transpose" and "tridia_solver"

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