Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |