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

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

last commit documented

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