Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver_mod.f90 @ 2104

Last change on this file since 2104 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |