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

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