Home

Context Navigation

tridia_solver_mod.f90 @ 4883

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

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

Context Navigation

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

Download in other formats: