Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |