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

Last change on this file since 1805 was 1805, checked in by maronga, 5 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 23.0 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 1805 2016-04-05 16:32:34Z maronga $
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!       tn = omp_get_thread_num()
512!       WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
513!       CALL local_flush( 120+tn )
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!       WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop   thread=', omp_get_thread_num()
522!       CALL local_flush( 120+tn )
523
524       IF ( j <= nnyh )  THEN
525#if defined( __intel11 )
526          CALL maketri_1dd( j, tri_for_1d )
527#else
528          CALL maketri_1dd( j )
529#endif
530       ELSE
531#if defined( __intel11 )
532          CALL maketri_1dd( ny+1-j, tri_for_1d )
533#else
534          CALL maketri_1dd( ny+1-j )
535#endif
536       ENDIF
537#if defined( __intel11 )
538       CALL split_1dd( tri_for_1d )
539#else
540       CALL split_1dd
541#endif
542       CALL substi_1dd( ar, tri_for_1d )
543
544    CONTAINS
545
546
547!------------------------------------------------------------------------------!
548! Description:
549! ------------
550!> computes the i- and j-dependent component of the matrix
551!------------------------------------------------------------------------------!
552#if defined( __intel11 )
553       SUBROUTINE maketri_1dd( j, tri_for_1d )
554#else
555       SUBROUTINE maketri_1dd( j )
556#endif
557
558          USE constants,                                                       &
559              ONLY:  pi
560
561          USE kinds
562
563          IMPLICIT NONE
564
565          INTEGER(iwp) ::  i    !<
566          INTEGER(iwp) ::  j    !<
567          INTEGER(iwp) ::  k    !<
568          INTEGER(iwp) ::  nnxh !<
569
570          REAL(wp)     ::  a !<
571          REAL(wp)     ::  c !<
572
573          REAL(wp), DIMENSION(0:nx) ::  l !<
574
575#if defined( __intel11 )
576          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
577#endif
578
579
580          nnxh = ( nx + 1 ) / 2
581!
582!--       Provide the tridiagonal matrix for solution of the Poisson equation in
583!--       Fourier space. The coefficients are computed following the method of
584!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
585!--       Siano's original version by discretizing the Poisson equation,
586!--       before it is Fourier-transformed
587          DO  i = 0, nx
588             IF ( i >= 0 .AND. i <= nnxh ) THEN
589                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * 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             ELSE
594                l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
595                                   REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
596                       2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
597                                   REAL( ny+1, KIND=wp ) ) ) * ddy2
598             ENDIF
599          ENDDO
600
601          DO  k = 0, nz-1
602             DO  i = 0, nx
603                a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
604                c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
605                tri_for_1d(1,i,k) = a + c - l(i)
606             ENDDO
607          ENDDO
608          IF ( ibc_p_b == 1 )  THEN
609             DO  i = 0, nx
610                tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
611             ENDDO
612          ENDIF
613          IF ( ibc_p_t == 1 )  THEN
614             DO  i = 0, nx
615                tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
616             ENDDO
617          ENDIF
618
619       END SUBROUTINE maketri_1dd
620
621
622!------------------------------------------------------------------------------!
623! Description:
624! ------------
625!> Splitting of the tridiagonal matrix (Thomas algorithm)
626!------------------------------------------------------------------------------!
627#if defined( __intel11 )
628       SUBROUTINE split_1dd( tri_for_1d )
629#else
630       SUBROUTINE split_1dd
631#endif
632
633
634          IMPLICIT NONE
635
636          INTEGER(iwp) ::  i !<
637          INTEGER(iwp) ::  k !<
638
639#if defined( __intel11 )
640          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
641#endif
642
643
644!
645!--       Splitting
646          DO  i = 0, nx
647             tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
648          ENDDO
649          DO  k = 1, nz-1
650             DO  i = 0, nx
651                tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
652                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)
653             ENDDO
654          ENDDO
655
656       END SUBROUTINE split_1dd
657
658
659!------------------------------------------------------------------------------!
660! Description:
661! ------------
662!> Substitution (Forward and Backward) (Thomas algorithm)
663!------------------------------------------------------------------------------!
664       SUBROUTINE substi_1dd( ar, tri_for_1d )
665
666
667          IMPLICIT NONE
668
669          INTEGER(iwp) ::  i !<
670          INTEGER(iwp) ::  k !<
671
672          REAL(wp), DIMENSION(0:nx,nz)       ::  ar         !<
673          REAL(wp), DIMENSION(0:nx,0:nz-1)   ::  ar1        !<
674          REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
675
676!
677!--       Forward substitution
678          DO  i = 0, nx
679             ar1(i,0) = ar(i,1)
680          ENDDO
681          DO  k = 1, nz-1
682             DO  i = 0, nx
683                ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
684             ENDDO
685          ENDDO
686
687!
688!--       Backward substitution
689!--       Note, the add of 1.0E-20 in the denominator is due to avoid divisions
690!--       by zero appearing if the pressure bc is set to neumann at the top of
691!--       the model domain.
692          DO  i = 0, nx
693             ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
694          ENDDO
695          DO  k = nz-2, 0, -1
696             DO  i = 0, nx
697                ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
698                            / tri_for_1d(4,i,k)
699             ENDDO
700          ENDDO
701
702!
703!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
704!--       The respective values of ar should be zero at all k-levels if
705!--       acceleration of horizontally averaged vertical velocity is zero.
706          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
707             IF ( j == 0 )  THEN
708                DO  k = 1, nz
709                   ar(0,k) = 0.0_wp
710                ENDDO
711             ENDIF
712          ENDIF
713
714       END SUBROUTINE substi_1dd
715
716    END SUBROUTINE tridia_1dd
717
718
719 END MODULE tridia_solver
Note: See TracBrowser for help on using the repository browser.