source: palm/trunk/SOURCE/tridia_solver_mod.f90 @ 1850

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

added _mod string to several filenames to meet the naming convection for modules

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