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

Last change on this file since 3730 was 3690, checked in by knoop, 5 years ago

Enabled OpenACC usage without using the cudaFFT library.
Added respective palmtest build configuration and testcase.

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