Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |