source: palm/trunk/SOURCE/tridia_solver_mod.f90

Last change on this file was 4828, checked in by Giersch, 5 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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