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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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