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

Last change on this file since 3248 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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