Home

Context Navigation

tridia_solver.f90 @ 1842

Last change on this file since 1842 was 1818, checked in by maronga, 8 years ago
last commit documented / copyright update
Property svn:keywords set to `Id`
File size: 22.3 KB

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

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

Context Navigation

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

Download in other formats: