Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver.f90 @ 1362

Last change on this file since 1362 was 1343, checked in by kanani, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |