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

Last change on this file since 4173 was 3761, checked in by raasch, 5 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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