Home

Context Navigation

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

Last change on this file since 2118 was 2118, checked in by raasch, 7 years ago
all OpenACC directives and related parts removed from the code
Property svn:keywords set to `Id`
File size: 21.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	! OpenACC directives removed
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: tridia_solver_mod.f90 2118 2017-01-17 16:38:49Z raasch $
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
198
199	nnxh = ( nx + 1 ) / 2
200	nnyh = ( ny + 1 ) / 2
201
202	DO j = nys_z, nyn_z
203	DO i = nxl_z, nxr_z
204	IF ( j >= 0 .AND. j <= nnyh ) THEN
205	IF ( i >= 0 .AND. i <= nnxh ) THEN
206	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
207	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
208	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
209	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
210	ELSE
211	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
212	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
213	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
214	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
215	ENDIF
216	ELSE
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 * ( ny+1-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 * ( ny+1-j ) ) / &
226	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
227	ENDIF
228	ENDIF
229	ENDDO
230	ENDDO
231
232	DO k = 0, nz-1
233	DO j = nys_z, nyn_z
234	DO i = nxl_z, nxr_z
235	tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
236	ENDDO
237	ENDDO
238	ENDDO
239
240	IF ( ibc_p_b == 1 ) THEN
241	DO j = nys_z, nyn_z
242	DO i = nxl_z, nxr_z
243	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
244	ENDDO
245	ENDDO
246	ENDIF
247	IF ( ibc_p_t == 1 ) THEN
248	DO j = nys_z, nyn_z
249	DO i = nxl_z, nxr_z
250	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
251	ENDDO
252	ENDDO
253	ENDIF
254
255	END SUBROUTINE maketri
256
257
258	!------------------------------------------------------------------------------!
259	! Description:
260	! ------------
261	!> Substitution (Forward and Backward) (Thomas algorithm)
262	!------------------------------------------------------------------------------!
263	SUBROUTINE tridia_substi( ar )
264
265
266	USE arrays_3d, &
267	ONLY: tri
268
269	USE control_parameters, &
270	ONLY: ibc_p_b, ibc_p_t
271
272	USE kinds
273
274	IMPLICIT NONE
275
276	INTEGER(iwp) :: i !<
277	INTEGER(iwp) :: j !<
278	INTEGER(iwp) :: k !<
279
280	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
281
282	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
283
284	!
285	!-- Forward substitution
286	DO k = 0, nz - 1
287	DO j = nys_z, nyn_z
288	DO i = nxl_z, nxr_z
289
290	IF ( k == 0 ) THEN
291	ar1(i,j,k) = ar(i,j,k+1)
292	ELSE
293	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
294	ENDIF
295
296	ENDDO
297	ENDDO
298	ENDDO
299
300	!
301	!-- Backward substitution
302	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
303	!-- by zero appearing if the pressure bc is set to neumann at the top of
304	!-- the model domain.
305	DO k = nz-1, 0, -1
306	DO j = nys_z, nyn_z
307	DO i = nxl_z, nxr_z
308
309	IF ( k == nz-1 ) THEN
310	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
311	ELSE
312	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
313	/ tri(i,j,k,1)
314	ENDIF
315	ENDDO
316	ENDDO
317	ENDDO
318
319	!
320	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
321	!-- The respective values of ar should be zero at all k-levels if
322	!-- acceleration of horizontally averaged vertical velocity is zero.
323	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
324	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
325	DO k = 1, nz
326	ar(nxl_z,nys_z,k) = 0.0_wp
327	ENDDO
328	ENDIF
329	ENDIF
330
331	END SUBROUTINE tridia_substi
332
333
334	!------------------------------------------------------------------------------!
335	! Description:
336	! ------------
337	!> Substitution (Forward and Backward) (Thomas algorithm)
338	!------------------------------------------------------------------------------!
339	SUBROUTINE tridia_substi_overlap( ar, jj )
340
341
342	USE arrays_3d, &
343	ONLY: tri
344
345	USE control_parameters, &
346	ONLY: ibc_p_b, ibc_p_t
347
348	USE kinds
349
350	IMPLICIT NONE
351
352	INTEGER(iwp) :: i !<
353	INTEGER(iwp) :: j !<
354	INTEGER(iwp) :: jj !<
355	INTEGER(iwp) :: k !<
356
357	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
358
359	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
360
361	!
362	!-- Forward substitution
363	DO k = 0, nz - 1
364	DO j = nys_z, nyn_z
365	DO i = nxl_z, nxr_z
366
367	IF ( k == 0 ) THEN
368	ar1(i,j,k) = ar(i,j,k+1)
369	ELSE
370	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
371	ENDIF
372
373	ENDDO
374	ENDDO
375	ENDDO
376
377	!
378	!-- Backward substitution
379	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
380	!-- by zero appearing if the pressure bc is set to neumann at the top of
381	!-- the model domain.
382	DO k = nz-1, 0, -1
383	DO j = nys_z, nyn_z
384	DO i = nxl_z, nxr_z
385
386	IF ( k == nz-1 ) THEN
387	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
388	ELSE
389	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
390	/ tri(i,jj,k,1)
391	ENDIF
392	ENDDO
393	ENDDO
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	DO k = 1, nz
403	ar(nxl_z,nys_z,k) = 0.0_wp
404	ENDDO
405	ENDIF
406	ENDIF
407
408	END SUBROUTINE tridia_substi_overlap
409
410
411	!------------------------------------------------------------------------------!
412	! Description:
413	! ------------
414	!> Splitting of the tridiagonal matrix (Thomas algorithm)
415	!------------------------------------------------------------------------------!
416	SUBROUTINE split
417
418
419	USE arrays_3d, &
420	ONLY: tri, tric
421
422	USE kinds
423
424	IMPLICIT NONE
425
426	INTEGER(iwp) :: i !<
427	INTEGER(iwp) :: j !<
428	INTEGER(iwp) :: k !<
429	!
430	!-- Splitting
431	DO j = nys_z, nyn_z
432	DO i = nxl_z, nxr_z
433	tri(i,j,0,1) = tric(i,j,0)
434	ENDDO
435	ENDDO
436
437	DO k = 1, nz-1
438	DO j = nys_z, nyn_z
439	DO i = nxl_z, nxr_z
440	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
441	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
442	ENDDO
443	ENDDO
444	ENDDO
445
446	END SUBROUTINE split
447
448
449	!------------------------------------------------------------------------------!
450	! Description:
451	! ------------
452	!> Solves the linear system of equations for a 1d-decomposition along x (see
453	!> tridia)
454	!>
455	!> @attention when using the intel compilers older than 12.0, array tri must
456	!> be passed as an argument to the contained subroutines. Otherwise
457	!> addres faults will occur. This feature can be activated with
458	!> cpp-switch __intel11
459	!> On NEC, tri should not be passed (except for routine substi_1dd)
460	!> because this causes very bad performance.
461	!------------------------------------------------------------------------------!
462
463	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
464
465
466	USE arrays_3d, &
467	ONLY: ddzu_pres, ddzw, rho_air, rho_air_zw
468
469	USE control_parameters, &
470	ONLY: ibc_p_b, ibc_p_t
471
472	USE kinds
473
474	IMPLICIT NONE
475
476	INTEGER(iwp) :: i !<
477	INTEGER(iwp) :: j !<
478	INTEGER(iwp) :: k !<
479	INTEGER(iwp) :: nnyh !<
480	INTEGER(iwp) :: nx !<
481	INTEGER(iwp) :: ny !<
482	INTEGER(iwp) :: omp_get_thread_num !<
483	INTEGER(iwp) :: tn !<
484
485	REAL(wp) :: ddx2 !<
486	REAL(wp) :: ddy2 !<
487
488	REAL(wp), DIMENSION(0:nx,1:nz) :: ar !<
489	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
490
491
492	nnyh = ( ny + 1 ) / 2
493
494	!
495	!-- Define constant elements of the tridiagonal matrix.
496	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
497	!-- the exchanged loops create bank conflicts. The following directive
498	!-- prohibits loop exchange and the loops perform much better.
499	!CDIR NOLOOPCHG
500	DO k = 0, nz-1
501	DO i = 0,nx
502	tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
503	tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
504	ENDDO
505	ENDDO
506
507	IF ( j <= nnyh ) THEN
508	CALL maketri_1dd( j )
509	ELSE
510	CALL maketri_1dd( ny+1-j )
511	ENDIF
512
513	CALL split_1dd
514	CALL substi_1dd( ar, tri_for_1d )
515
516	CONTAINS
517
518
519	!------------------------------------------------------------------------------!
520	! Description:
521	! ------------
522	!> computes the i- and j-dependent component of the matrix
523	!------------------------------------------------------------------------------!
524	SUBROUTINE maketri_1dd( j )
525
526	USE constants, &
527	ONLY: pi
528
529	USE kinds
530
531	IMPLICIT NONE
532
533	INTEGER(iwp) :: i !<
534	INTEGER(iwp) :: j !<
535	INTEGER(iwp) :: k !<
536	INTEGER(iwp) :: nnxh !<
537
538	REAL(wp) :: a !<
539	REAL(wp) :: c !<
540
541	REAL(wp), DIMENSION(0:nx) :: l !<
542
543
544	nnxh = ( nx + 1 ) / 2
545	!
546	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
547	!-- Fourier space. The coefficients are computed following the method of
548	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
549	!-- Siano's original version by discretizing the Poisson equation,
550	!-- before it is Fourier-transformed
551	DO i = 0, nx
552	IF ( i >= 0 .AND. i <= nnxh ) THEN
553	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
554	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
555	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
556	REAL( ny+1, KIND=wp ) ) ) * ddy2
557	ELSE
558	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
559	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
560	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
561	REAL( ny+1, KIND=wp ) ) ) * ddy2
562	ENDIF
563	ENDDO
564
565	DO k = 0, nz-1
566	DO i = 0, nx
567	a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
568	c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
569	tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
570	ENDDO
571	ENDDO
572	IF ( ibc_p_b == 1 ) THEN
573	DO i = 0, nx
574	tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
575	ENDDO
576	ENDIF
577	IF ( ibc_p_t == 1 ) THEN
578	DO i = 0, nx
579	tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
580	ENDDO
581	ENDIF
582
583	END SUBROUTINE maketri_1dd
584
585
586	!------------------------------------------------------------------------------!
587	! Description:
588	! ------------
589	!> Splitting of the tridiagonal matrix (Thomas algorithm)
590	!------------------------------------------------------------------------------!
591	SUBROUTINE split_1dd
592
593	IMPLICIT NONE
594
595	INTEGER(iwp) :: i !<
596	INTEGER(iwp) :: k !<
597
598
599	!
600	!-- Splitting
601	DO i = 0, nx
602	tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
603	ENDDO
604	DO k = 1, nz-1
605	DO i = 0, nx
606	tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
607	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)
608	ENDDO
609	ENDDO
610
611	END SUBROUTINE split_1dd
612
613
614	!------------------------------------------------------------------------------!
615	! Description:
616	! ------------
617	!> Substitution (Forward and Backward) (Thomas algorithm)
618	!------------------------------------------------------------------------------!
619	SUBROUTINE substi_1dd( ar, tri_for_1d )
620
621
622	IMPLICIT NONE
623
624	INTEGER(iwp) :: i !<
625	INTEGER(iwp) :: k !<
626
627	REAL(wp), DIMENSION(0:nx,nz) :: ar !<
628	REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !<
629	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
630
631	!
632	!-- Forward substitution
633	DO i = 0, nx
634	ar1(i,0) = ar(i,1)
635	ENDDO
636	DO k = 1, nz-1
637	DO i = 0, nx
638	ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
639	ENDDO
640	ENDDO
641
642	!
643	!-- Backward substitution
644	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
645	!-- by zero appearing if the pressure bc is set to neumann at the top of
646	!-- the model domain.
647	DO i = 0, nx
648	ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
649	ENDDO
650	DO k = nz-2, 0, -1
651	DO i = 0, nx
652	ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
653	/ tri_for_1d(4,i,k)
654	ENDDO
655	ENDDO
656
657	!
658	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
659	!-- The respective values of ar should be zero at all k-levels if
660	!-- acceleration of horizontally averaged vertical velocity is zero.
661	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
662	IF ( j == 0 ) THEN
663	DO k = 1, nz
664	ar(0,k) = 0.0_wp
665	ENDDO
666	ENDIF
667	ENDIF
668
669	END SUBROUTINE substi_1dd
670
671	END SUBROUTINE tridia_1dd
672
673
674	END MODULE tridia_solver

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |