source: palm/trunk/SOURCE/tridia_solver.f90 @ 1406

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