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

Last change on this file since 3643 was 3634, checked in by knoop, 6 years ago

OpenACC port for SPEC

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