Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |