Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |