Home

Context Navigation

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

Last change on this file since 1985 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |