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

Last change on this file since 1342 was 1342, checked in by kanani, 7 years ago

REAL constants defined as wp-kind

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