Home

Context Navigation

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

Last change on this file since 3254 was 3241, checked in by raasch, 6 years ago
various changes to avoid compiler warnings (mainly removal of unused variables)
Property svn:keywords set to `Id`
File size: 21.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |