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

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

removed parameter file check. update of mrungui for compilation with qt5

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