Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |