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

Last change on this file since 4281 was 4236, checked in by knoop, 21 months ago

Added missing OpenMP directives within "transpose" and "tridia_solver"

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