Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |