Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |