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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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