Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |