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

Last change on this file since 3301 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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