Home

Context Navigation

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

Last change on this file since 1804 was 1804, checked in by maronga, 8 years ago
removed parameter file check. update of mrungui for compilation with qt5
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	! Removed code for parameter file check (__check)
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: tridia_solver.f90 1804 2016-04-05 16:30:18Z maronga $
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	CALL maketri
133	CALL split
134
135	END SUBROUTINE tridia_init
136
137
138	!------------------------------------------------------------------------------!
139	! Description:
140	! ------------
141	!> Computes the i- and j-dependent component of the matrix
142	!> Provide the constant coefficients of the tridiagonal matrix for solution
143	!> of the Poisson equation in Fourier space.
144	!> The coefficients are computed following the method of
145	!> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
146	!> Siano's original version by discretizing the Poisson equation,
147	!> before it is Fourier-transformed.
148	!------------------------------------------------------------------------------!
149	SUBROUTINE maketri
150
151
152	USE arrays_3d, &
153	ONLY: tric
154
155	USE constants, &
156	ONLY: pi
157
158	USE control_parameters, &
159	ONLY: ibc_p_b, ibc_p_t
160
161	USE grid_variables, &
162	ONLY: dx, dy
163
164
165	USE kinds
166
167	IMPLICIT NONE
168
169	INTEGER(iwp) :: i !<
170	INTEGER(iwp) :: j !<
171	INTEGER(iwp) :: k !<
172	INTEGER(iwp) :: nnxh !<
173	INTEGER(iwp) :: nnyh !<
174
175	REAL(wp) :: ll(nxl_z:nxr_z,nys_z:nyn_z) !<
176	!$acc declare create( ll )
177
178
179	nnxh = ( nx + 1 ) / 2
180	nnyh = ( ny + 1 ) / 2
181
182	!$acc kernels present( tric )
183	DO j = nys_z, nyn_z
184	DO i = nxl_z, nxr_z
185	IF ( j >= 0 .AND. j <= nnyh ) THEN
186	IF ( i >= 0 .AND. i <= nnxh ) THEN
187	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
188	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
189	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
190	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
191	ELSE
192	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
193	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
194	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
195	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
196	ENDIF
197	ELSE
198	IF ( i >= 0 .AND. i <= nnxh ) THEN
199	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
200	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
201	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
202	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
203	ELSE
204	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
205	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
206	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
207	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
208	ENDIF
209	ENDIF
210	ENDDO
211	ENDDO
212
213	DO k = 0, nz-1
214	DO j = nys_z, nyn_z
215	DO i = nxl_z, nxr_z
216	tric(i,j,k) = ddzuw(k,3) - ll(i,j)
217	ENDDO
218	ENDDO
219	ENDDO
220	!$acc end kernels
221
222	IF ( ibc_p_b == 1 ) THEN
223	!$acc kernels present( tric )
224	DO j = nys_z, nyn_z
225	DO i = nxl_z, nxr_z
226	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
227	ENDDO
228	ENDDO
229	!$acc end kernels
230	ENDIF
231	IF ( ibc_p_t == 1 ) THEN
232	!$acc kernels present( tric )
233	DO j = nys_z, nyn_z
234	DO i = nxl_z, nxr_z
235	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
236	ENDDO
237	ENDDO
238	!$acc end kernels
239	ENDIF
240
241	END SUBROUTINE maketri
242
243
244	!------------------------------------------------------------------------------!
245	! Description:
246	! ------------
247	!> Substitution (Forward and Backward) (Thomas algorithm)
248	!------------------------------------------------------------------------------!
249	SUBROUTINE tridia_substi( ar )
250
251
252	USE arrays_3d, &
253	ONLY: tri
254
255	USE control_parameters, &
256	ONLY: ibc_p_b, ibc_p_t
257
258	USE kinds
259
260	IMPLICIT NONE
261
262	INTEGER(iwp) :: i !<
263	INTEGER(iwp) :: j !<
264	INTEGER(iwp) :: k !<
265
266	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
267
268	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
269	!$acc declare create( ar1 )
270
271	!
272	!-- Forward substitution
273	DO k = 0, nz - 1
274	!$acc kernels present( ar, tri )
275	DO j = nys_z, nyn_z
276	DO i = nxl_z, nxr_z
277
278	IF ( k == 0 ) THEN
279	ar1(i,j,k) = ar(i,j,k+1)
280	ELSE
281	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
282	ENDIF
283
284	ENDDO
285	ENDDO
286	!$acc end kernels
287	ENDDO
288
289	!
290	!-- Backward substitution
291	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
292	!-- by zero appearing if the pressure bc is set to neumann at the top of
293	!-- the model domain.
294	DO k = nz-1, 0, -1
295	!$acc kernels present( ar, tri )
296	DO j = nys_z, nyn_z
297	DO i = nxl_z, nxr_z
298
299	IF ( k == nz-1 ) THEN
300	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
301	ELSE
302	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
303	/ tri(i,j,k,1)
304	ENDIF
305	ENDDO
306	ENDDO
307	!$acc end kernels
308	ENDDO
309
310	!
311	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
312	!-- The respective values of ar should be zero at all k-levels if
313	!-- acceleration of horizontally averaged vertical velocity is zero.
314	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
315	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
316	!$acc kernels loop present( ar )
317	DO k = 1, nz
318	ar(nxl_z,nys_z,k) = 0.0_wp
319	ENDDO
320	!$acc end kernels loop
321	ENDIF
322	ENDIF
323
324	END SUBROUTINE tridia_substi
325
326
327	!------------------------------------------------------------------------------!
328	! Description:
329	! ------------
330	!> Substitution (Forward and Backward) (Thomas algorithm)
331	!------------------------------------------------------------------------------!
332	SUBROUTINE tridia_substi_overlap( ar, jj )
333
334
335	USE arrays_3d, &
336	ONLY: tri
337
338	USE control_parameters, &
339	ONLY: ibc_p_b, ibc_p_t
340
341	USE kinds
342
343	IMPLICIT NONE
344
345	INTEGER(iwp) :: i !<
346	INTEGER(iwp) :: j !<
347	INTEGER(iwp) :: jj !<
348	INTEGER(iwp) :: k !<
349
350	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
351
352	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
353	!$acc declare create( ar1 )
354
355	!
356	!-- Forward substitution
357	DO k = 0, nz - 1
358	!$acc kernels present( ar, tri )
359	!$acc loop
360	DO j = nys_z, nyn_z
361	DO i = nxl_z, nxr_z
362
363	IF ( k == 0 ) THEN
364	ar1(i,j,k) = ar(i,j,k+1)
365	ELSE
366	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
367	ENDIF
368
369	ENDDO
370	ENDDO
371	!$acc end kernels
372	ENDDO
373
374	!
375	!-- Backward substitution
376	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
377	!-- by zero appearing if the pressure bc is set to neumann at the top of
378	!-- the model domain.
379	DO k = nz-1, 0, -1
380	!$acc kernels present( ar, tri )
381	!$acc loop
382	DO j = nys_z, nyn_z
383	DO i = nxl_z, nxr_z
384
385	IF ( k == nz-1 ) THEN
386	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
387	ELSE
388	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
389	/ tri(i,jj,k,1)
390	ENDIF
391	ENDDO
392	ENDDO
393	!$acc end kernels
394	ENDDO
395
396	!
397	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
398	!-- The respective values of ar should be zero at all k-levels if
399	!-- acceleration of horizontally averaged vertical velocity is zero.
400	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
401	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
402	!$acc kernels loop present( ar )
403	DO k = 1, nz
404	ar(nxl_z,nys_z,k) = 0.0_wp
405	ENDDO
406	ENDIF
407	ENDIF
408
409	END SUBROUTINE tridia_substi_overlap
410
411
412	!------------------------------------------------------------------------------!
413	! Description:
414	! ------------
415	!> Splitting of the tridiagonal matrix (Thomas algorithm)
416	!------------------------------------------------------------------------------!
417	SUBROUTINE split
418
419
420	USE arrays_3d, &
421	ONLY: tri, tric
422
423	USE kinds
424
425	IMPLICIT NONE
426
427	INTEGER(iwp) :: i !<
428	INTEGER(iwp) :: j !<
429	INTEGER(iwp) :: k !<
430	!
431	!-- Splitting
432	!$acc kernels present( tri, tric )
433	!$acc loop
434	DO j = nys_z, nyn_z
435	!$acc loop vector( 32 )
436	DO i = nxl_z, nxr_z
437	tri(i,j,0,1) = tric(i,j,0)
438	ENDDO
439	ENDDO
440	!$acc end kernels
441
442	DO k = 1, nz-1
443	!$acc kernels present( tri, tric )
444	!$acc loop
445	DO j = nys_z, nyn_z
446	!$acc loop vector( 32 )
447	DO i = nxl_z, nxr_z
448	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
449	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
450	ENDDO
451	ENDDO
452	!$acc end kernels
453	ENDDO
454
455	END SUBROUTINE split
456
457
458	!------------------------------------------------------------------------------!
459	! Description:
460	! ------------
461	!> Solves the linear system of equations for a 1d-decomposition along x (see
462	!> tridia)
463	!>
464	!> @attention when using the intel compilers older than 12.0, array tri must
465	!> be passed as an argument to the contained subroutines. Otherwise
466	!> addres faults will occur. This feature can be activated with
467	!> cpp-switch __intel11
468	!> On NEC, tri should not be passed (except for routine substi_1dd)
469	!> because this causes very bad performance.
470	!------------------------------------------------------------------------------!
471
472	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
473
474
475	USE arrays_3d, &
476	ONLY: ddzu_pres, ddzw
477
478	USE control_parameters, &
479	ONLY: ibc_p_b, ibc_p_t
480
481	USE kinds
482
483	IMPLICIT NONE
484
485	INTEGER(iwp) :: i !<
486	INTEGER(iwp) :: j !<
487	INTEGER(iwp) :: k !<
488	INTEGER(iwp) :: nnyh !<
489	INTEGER(iwp) :: nx !<
490	INTEGER(iwp) :: ny !<
491	INTEGER(iwp) :: omp_get_thread_num !<
492	INTEGER(iwp) :: tn !<
493
494	REAL(wp) :: ddx2 !<
495	REAL(wp) :: ddy2 !<
496
497	REAL(wp), DIMENSION(0:nx,1:nz) :: ar !<
498	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
499
500
501	nnyh = ( ny + 1 ) / 2
502
503	!
504	!-- Define constant elements of the tridiagonal matrix.
505	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
506	!-- the exchanged loops create bank conflicts. The following directive
507	!-- prohibits loop exchange and the loops perform much better.
508	! tn = omp_get_thread_num()
509	! WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
510	! CALL local_flush( 120+tn )
511	!CDIR NOLOOPCHG
512	DO k = 0, nz-1
513	DO i = 0,nx
514	tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
515	tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
516	ENDDO
517	ENDDO
518	! WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop thread=', omp_get_thread_num()
519	! CALL local_flush( 120+tn )
520
521	IF ( j <= nnyh ) THEN
522	#if defined( __intel11 )
523	CALL maketri_1dd( j, tri_for_1d )
524	#else
525	CALL maketri_1dd( j )
526	#endif
527	ELSE
528	#if defined( __intel11 )
529	CALL maketri_1dd( ny+1-j, tri_for_1d )
530	#else
531	CALL maketri_1dd( ny+1-j )
532	#endif
533	ENDIF
534	#if defined( __intel11 )
535	CALL split_1dd( tri_for_1d )
536	#else
537	CALL split_1dd
538	#endif
539	CALL substi_1dd( ar, tri_for_1d )
540
541	CONTAINS
542
543
544	!------------------------------------------------------------------------------!
545	! Description:
546	! ------------
547	!> computes the i- and j-dependent component of the matrix
548	!------------------------------------------------------------------------------!
549	#if defined( __intel11 )
550	SUBROUTINE maketri_1dd( j, tri_for_1d )
551	#else
552	SUBROUTINE maketri_1dd( j )
553	#endif
554
555	USE constants, &
556	ONLY: pi
557
558	USE kinds
559
560	IMPLICIT NONE
561
562	INTEGER(iwp) :: i !<
563	INTEGER(iwp) :: j !<
564	INTEGER(iwp) :: k !<
565	INTEGER(iwp) :: nnxh !<
566
567	REAL(wp) :: a !<
568	REAL(wp) :: c !<
569
570	REAL(wp), DIMENSION(0:nx) :: l !<
571
572	#if defined( __intel11 )
573	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
574	#endif
575
576
577	nnxh = ( nx + 1 ) / 2
578	!
579	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
580	!-- Fourier space. The coefficients are computed following the method of
581	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
582	!-- Siano's original version by discretizing the Poisson equation,
583	!-- before it is Fourier-transformed
584	DO i = 0, nx
585	IF ( i >= 0 .AND. i <= nnxh ) THEN
586	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
587	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
588	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
589	REAL( ny+1, KIND=wp ) ) ) * ddy2
590	ELSE
591	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
592	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
593	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
594	REAL( ny+1, KIND=wp ) ) ) * ddy2
595	ENDIF
596	ENDDO
597
598	DO k = 0, nz-1
599	DO i = 0, nx
600	a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
601	c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
602	tri_for_1d(1,i,k) = a + c - l(i)
603	ENDDO
604	ENDDO
605	IF ( ibc_p_b == 1 ) THEN
606	DO i = 0, nx
607	tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
608	ENDDO
609	ENDIF
610	IF ( ibc_p_t == 1 ) THEN
611	DO i = 0, nx
612	tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
613	ENDDO
614	ENDIF
615
616	END SUBROUTINE maketri_1dd
617
618
619	!------------------------------------------------------------------------------!
620	! Description:
621	! ------------
622	!> Splitting of the tridiagonal matrix (Thomas algorithm)
623	!------------------------------------------------------------------------------!
624	#if defined( __intel11 )
625	SUBROUTINE split_1dd( tri_for_1d )
626	#else
627	SUBROUTINE split_1dd
628	#endif
629
630
631	IMPLICIT NONE
632
633	INTEGER(iwp) :: i !<
634	INTEGER(iwp) :: k !<
635
636	#if defined( __intel11 )
637	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
638	#endif
639
640
641	!
642	!-- Splitting
643	DO i = 0, nx
644	tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
645	ENDDO
646	DO k = 1, nz-1
647	DO i = 0, nx
648	tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
649	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)
650	ENDDO
651	ENDDO
652
653	END SUBROUTINE split_1dd
654
655
656	!------------------------------------------------------------------------------!
657	! Description:
658	! ------------
659	!> Substitution (Forward and Backward) (Thomas algorithm)
660	!------------------------------------------------------------------------------!
661	SUBROUTINE substi_1dd( ar, tri_for_1d )
662
663
664	IMPLICIT NONE
665
666	INTEGER(iwp) :: i !<
667	INTEGER(iwp) :: k !<
668
669	REAL(wp), DIMENSION(0:nx,nz) :: ar !<
670	REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !<
671	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
672
673	!
674	!-- Forward substitution
675	DO i = 0, nx
676	ar1(i,0) = ar(i,1)
677	ENDDO
678	DO k = 1, nz-1
679	DO i = 0, nx
680	ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
681	ENDDO
682	ENDDO
683
684	!
685	!-- Backward substitution
686	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
687	!-- by zero appearing if the pressure bc is set to neumann at the top of
688	!-- the model domain.
689	DO i = 0, nx
690	ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
691	ENDDO
692	DO k = nz-2, 0, -1
693	DO i = 0, nx
694	ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
695	/ tri_for_1d(4,i,k)
696	ENDDO
697	ENDDO
698
699	!
700	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
701	!-- The respective values of ar should be zero at all k-levels if
702	!-- acceleration of horizontally averaged vertical velocity is zero.
703	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
704	IF ( j == 0 ) THEN
705	DO k = 1, nz
706	ar(0,k) = 0.0_wp
707	ENDDO
708	ENDIF
709	ENDIF
710
711	END SUBROUTINE substi_1dd
712
713	END SUBROUTINE tridia_1dd
714
715
716	END MODULE tridia_solver

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |