Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver_mod.f90 @ 4327

Last change on this file since 4327 was 4236, checked in by knoop, 5 years ago
Added missing OpenMP directives within "transpose" and "tridia_solver"
Property svn:keywords set to `Id`
File size: 21.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |