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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

Code annotations made doxygen readable

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