Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |