Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |