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

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

preprocessor directives using machine dependent flags (lc, ibm, etc.) mostly removed from the code

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