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

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

Anelastic approximation implemented

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