Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |