Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |