Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |