Home

Context Navigation

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

Last change on this file since 2036 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |