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

Last change on this file since 2118 was 2118, checked in by raasch, 5 years ago

all OpenACC directives and related parts removed from the code

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