Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |