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

Last change on this file since 1524 was 1407, checked in by raasch, 10 years ago

last commit documented

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