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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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