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

Last change on this file since 1809 was 1809, checked in by raasch, 8 years ago

last commit documented

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