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

Last change on this file since 2699 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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