Home

Context Navigation

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

Last change on this file since 3755 was 3690, checked in by knoop, 6 years ago
Enabled OpenACC usage without using the cudaFFT library. Added respective palmtest build configuration and testcase.
Property svn:keywords set to `Id`
File size: 22.6 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 3690 2019-01-22 22:56:42Z suehring $
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	#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
105
106	!
107	! Description:
108	! ------------
109	!> solves the linear system of equations:
110	!>
111	!> -(4 pi^2(i^2/(dx^2nnx^2)+j^2/(dy^2nny^2))+
112	!> 1/(dzu(k)dzw(k))+1/(dzu(k-1)dzw(k)))*p(i,j,k)+
113	!> 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)
114	!>
115	!> by using the Thomas algorithm
116	!------------------------------------------------------------------------------!
117	MODULE tridia_solver
118
119
120	USE basic_constants_and_equations_mod, &
121	ONLY: pi
122
123	USE indices, &
124	ONLY: nx, ny, nz
125
126	USE kinds
127
128	USE transpose_indices, &
129	ONLY: nxl_z, nyn_z, nxr_z, nys_z
130
131	IMPLICIT NONE
132
133	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddzuw !<
134
135	PRIVATE
136
137	INTERFACE tridia_substi
138	MODULE PROCEDURE tridia_substi
139	END INTERFACE tridia_substi
140
141	INTERFACE tridia_substi_overlap
142	MODULE PROCEDURE tridia_substi_overlap
143	END INTERFACE tridia_substi_overlap
144
145	PUBLIC tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
146
147	CONTAINS
148
149
150	!------------------------------------------------------------------------------!
151	! Description:
152	! ------------
153	!> @todo Missing subroutine description.
154	!------------------------------------------------------------------------------!
155	SUBROUTINE tridia_init
156
157	USE arrays_3d, &
158	ONLY: ddzu_pres, ddzw, rho_air_zw, tri
159
160	IMPLICIT NONE
161
162	INTEGER(iwp) :: k !<
163
164	ALLOCATE( ddzuw(0:nz-1,3) )
165
166	DO k = 0, nz-1
167	ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
168	ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
169	ddzuw(k,3) = -1.0_wp * &
170	( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) + &
171	ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
172	ENDDO
173	!
174	!-- Calculate constant coefficients of the tridiagonal matrix
175	CALL maketri
176	CALL split
177
178	#if __acc_fft_device
179	!$ACC ENTER DATA &
180	!$ACC COPYIN(ddzuw(0:nz-1,1:3)) &
181	!$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2))
182	#endif
183
184	END SUBROUTINE tridia_init
185
186
187	!------------------------------------------------------------------------------!
188	! Description:
189	! ------------
190	!> Computes the i- and j-dependent component of the matrix
191	!> Provide the constant coefficients of the tridiagonal matrix for solution
192	!> of the Poisson equation in Fourier space.
193	!> The coefficients are computed following the method of
194	!> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
195	!> Siano's original version by discretizing the Poisson equation,
196	!> before it is Fourier-transformed.
197	!------------------------------------------------------------------------------!
198	SUBROUTINE maketri
199
200
201	USE arrays_3d, &
202	ONLY: tric, rho_air
203
204	USE control_parameters, &
205	ONLY: ibc_p_b, ibc_p_t
206
207	USE grid_variables, &
208	ONLY: dx, dy
209
210
211	IMPLICIT NONE
212
213	INTEGER(iwp) :: i !<
214	INTEGER(iwp) :: j !<
215	INTEGER(iwp) :: k !<
216	INTEGER(iwp) :: nnxh !<
217	INTEGER(iwp) :: nnyh !<
218
219	REAL(wp) :: ll(nxl_z:nxr_z,nys_z:nyn_z) !<
220
221
222	nnxh = ( nx + 1 ) / 2
223	nnyh = ( ny + 1 ) / 2
224
225	DO j = nys_z, nyn_z
226	DO i = nxl_z, nxr_z
227	IF ( j >= 0 .AND. j <= nnyh ) THEN
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 * 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 * j ) / &
237	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
238	ENDIF
239	ELSE
240	IF ( i >= 0 .AND. i <= nnxh ) THEN
241	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * 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	ELSE
246	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
247	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
248	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
249	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
250	ENDIF
251	ENDIF
252	ENDDO
253	ENDDO
254
255	DO k = 0, nz-1
256	DO j = nys_z, nyn_z
257	DO i = nxl_z, nxr_z
258	tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
259	ENDDO
260	ENDDO
261	ENDDO
262
263	IF ( ibc_p_b == 1 ) THEN
264	DO j = nys_z, nyn_z
265	DO i = nxl_z, nxr_z
266	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
267	ENDDO
268	ENDDO
269	ENDIF
270	IF ( ibc_p_t == 1 ) THEN
271	DO j = nys_z, nyn_z
272	DO i = nxl_z, nxr_z
273	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
274	ENDDO
275	ENDDO
276	ENDIF
277
278	END SUBROUTINE maketri
279
280
281	!------------------------------------------------------------------------------!
282	! Description:
283	! ------------
284	!> Substitution (Forward and Backward) (Thomas algorithm)
285	!------------------------------------------------------------------------------!
286	SUBROUTINE tridia_substi( ar )
287
288
289	USE arrays_3d, &
290	ONLY: tri
291
292	USE control_parameters, &
293	ONLY: ibc_p_b, ibc_p_t
294
295	IMPLICIT NONE
296
297	INTEGER(iwp) :: i !<
298	INTEGER(iwp) :: j !<
299	INTEGER(iwp) :: k !<
300
301	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
302
303	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
304	#if __acc_fft_device
305	!$ACC DECLARE CREATE(ar1)
306	#endif
307
308	!
309	!-- Forward substitution
310	#if __acc_fft_device
311	!$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k)
312	#endif
313	DO k = 0, nz - 1
314	#if __acc_fft_device
315	!$ACC LOOP COLLAPSE(2)
316	#endif
317	DO j = nys_z, nyn_z
318	DO i = nxl_z, nxr_z
319
320	IF ( k == 0 ) THEN
321	ar1(i,j,k) = ar(i,j,k+1)
322	ELSE
323	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
324	ENDIF
325
326	ENDDO
327	ENDDO
328	ENDDO
329	#if __acc_fft_device
330	!$ACC END PARALLEL
331	#endif
332
333	!
334	!-- Backward substitution
335	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
336	!-- by zero appearing if the pressure bc is set to neumann at the top of
337	!-- the model domain.
338	#if __acc_fft_device
339	!$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k)
340	#endif
341	DO k = nz-1, 0, -1
342	#if __acc_fft_device
343	!$ACC LOOP COLLAPSE(2)
344	#endif
345	DO j = nys_z, nyn_z
346	DO i = nxl_z, nxr_z
347
348	IF ( k == nz-1 ) THEN
349	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
350	ELSE
351	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
352	/ tri(i,j,k,1)
353	ENDIF
354	ENDDO
355	ENDDO
356	ENDDO
357	#if __acc_fft_device
358	!$ACC END PARALLEL
359	#endif
360
361	!
362	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
363	!-- The respective values of ar should be zero at all k-levels if
364	!-- acceleration of horizontally averaged vertical velocity is zero.
365	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
366	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
367	#if __acc_fft_device
368	!$ACC PARALLEL LOOP PRESENT(ar)
369	#endif
370	DO k = 1, nz
371	ar(nxl_z,nys_z,k) = 0.0_wp
372	ENDDO
373	ENDIF
374	ENDIF
375
376	END SUBROUTINE tridia_substi
377
378
379	!------------------------------------------------------------------------------!
380	! Description:
381	! ------------
382	!> Substitution (Forward and Backward) (Thomas algorithm)
383	!------------------------------------------------------------------------------!
384	SUBROUTINE tridia_substi_overlap( ar, jj )
385
386
387	USE arrays_3d, &
388	ONLY: tri
389
390	USE control_parameters, &
391	ONLY: ibc_p_b, ibc_p_t
392
393	IMPLICIT NONE
394
395	INTEGER(iwp) :: i !<
396	INTEGER(iwp) :: j !<
397	INTEGER(iwp) :: jj !<
398	INTEGER(iwp) :: k !<
399
400	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
401
402	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
403
404	!
405	!-- Forward substitution
406	DO k = 0, nz - 1
407	DO j = nys_z, nyn_z
408	DO i = nxl_z, nxr_z
409
410	IF ( k == 0 ) THEN
411	ar1(i,j,k) = ar(i,j,k+1)
412	ELSE
413	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
414	ENDIF
415
416	ENDDO
417	ENDDO
418	ENDDO
419
420	!
421	!-- Backward substitution
422	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
423	!-- by zero appearing if the pressure bc is set to neumann at the top of
424	!-- the model domain.
425	DO k = nz-1, 0, -1
426	DO j = nys_z, nyn_z
427	DO i = nxl_z, nxr_z
428
429	IF ( k == nz-1 ) THEN
430	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
431	ELSE
432	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
433	/ tri(i,jj,k,1)
434	ENDIF
435	ENDDO
436	ENDDO
437	ENDDO
438
439	!
440	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
441	!-- The respective values of ar should be zero at all k-levels if
442	!-- acceleration of horizontally averaged vertical velocity is zero.
443	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
444	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
445	DO k = 1, nz
446	ar(nxl_z,nys_z,k) = 0.0_wp
447	ENDDO
448	ENDIF
449	ENDIF
450
451	END SUBROUTINE tridia_substi_overlap
452
453
454	!------------------------------------------------------------------------------!
455	! Description:
456	! ------------
457	!> Splitting of the tridiagonal matrix (Thomas algorithm)
458	!------------------------------------------------------------------------------!
459	SUBROUTINE split
460
461
462	USE arrays_3d, &
463	ONLY: tri, tric
464
465	IMPLICIT NONE
466
467	INTEGER(iwp) :: i !<
468	INTEGER(iwp) :: j !<
469	INTEGER(iwp) :: k !<
470	!
471	!-- Splitting
472	DO j = nys_z, nyn_z
473	DO i = nxl_z, nxr_z
474	tri(i,j,0,1) = tric(i,j,0)
475	ENDDO
476	ENDDO
477
478	DO k = 1, nz-1
479	DO j = nys_z, nyn_z
480	DO i = nxl_z, nxr_z
481	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
482	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
483	ENDDO
484	ENDDO
485	ENDDO
486
487	END SUBROUTINE split
488
489
490	!------------------------------------------------------------------------------!
491	! Description:
492	! ------------
493	!> Solves the linear system of equations for a 1d-decomposition along x (see
494	!> tridia)
495	!>
496	!> @attention when using the intel compilers older than 12.0, array tri must
497	!> be passed as an argument to the contained subroutines. Otherwise
498	!> addres faults will occur. This feature can be activated with
499	!> cpp-switch __intel11
500	!> On NEC, tri should not be passed (except for routine substi_1dd)
501	!> because this causes very bad performance.
502	!------------------------------------------------------------------------------!
503
504	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
505
506
507	USE arrays_3d, &
508	ONLY: ddzu_pres, ddzw, rho_air, rho_air_zw
509
510	USE control_parameters, &
511	ONLY: ibc_p_b, ibc_p_t
512
513	IMPLICIT NONE
514
515	INTEGER(iwp) :: i !<
516	INTEGER(iwp) :: j !<
517	INTEGER(iwp) :: k !<
518	INTEGER(iwp) :: nnyh !<
519	INTEGER(iwp) :: nx !<
520	INTEGER(iwp) :: ny !<
521
522	REAL(wp) :: ddx2 !<
523	REAL(wp) :: ddy2 !<
524
525	REAL(wp), DIMENSION(0:nx,1:nz) :: ar !<
526	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
527
528
529	nnyh = ( ny + 1 ) / 2
530
531	!
532	!-- Define constant elements of the tridiagonal matrix.
533	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
534	!-- the exchanged loops create bank conflicts. The following directive
535	!-- prohibits loop exchange and the loops perform much better.
536	!CDIR NOLOOPCHG
537	DO k = 0, nz-1
538	DO i = 0,nx
539	tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
540	tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
541	ENDDO
542	ENDDO
543
544	IF ( j <= nnyh ) THEN
545	CALL maketri_1dd( j )
546	ELSE
547	CALL maketri_1dd( ny+1-j )
548	ENDIF
549
550	CALL split_1dd
551	CALL substi_1dd( ar, tri_for_1d )
552
553	CONTAINS
554
555
556	!------------------------------------------------------------------------------!
557	! Description:
558	! ------------
559	!> computes the i- and j-dependent component of the matrix
560	!------------------------------------------------------------------------------!
561	SUBROUTINE maketri_1dd( j )
562
563	IMPLICIT NONE
564
565	INTEGER(iwp) :: i !<
566	INTEGER(iwp) :: j !<
567	INTEGER(iwp) :: k !<
568	INTEGER(iwp) :: nnxh !<
569
570	REAL(wp) :: a !<
571	REAL(wp) :: c !<
572
573	REAL(wp), DIMENSION(0:nx) :: l !<
574
575
576	nnxh = ( nx + 1 ) / 2
577	!
578	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
579	!-- Fourier space. The coefficients are computed following the method of
580	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
581	!-- Siano's original version by discretizing the Poisson equation,
582	!-- before it is Fourier-transformed
583	DO i = 0, nx
584	IF ( i >= 0 .AND. i <= nnxh ) THEN
585	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
586	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
587	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
588	REAL( ny+1, KIND=wp ) ) ) * ddy2
589	ELSE
590	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
591	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
592	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
593	REAL( ny+1, KIND=wp ) ) ) * ddy2
594	ENDIF
595	ENDDO
596
597	DO k = 0, nz-1
598	DO i = 0, nx
599	a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
600	c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
601	tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
602	ENDDO
603	ENDDO
604	IF ( ibc_p_b == 1 ) THEN
605	DO i = 0, nx
606	tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
607	ENDDO
608	ENDIF
609	IF ( ibc_p_t == 1 ) THEN
610	DO i = 0, nx
611	tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
612	ENDDO
613	ENDIF
614
615	END SUBROUTINE maketri_1dd
616
617
618	!------------------------------------------------------------------------------!
619	! Description:
620	! ------------
621	!> Splitting of the tridiagonal matrix (Thomas algorithm)
622	!------------------------------------------------------------------------------!
623	SUBROUTINE split_1dd
624
625	IMPLICIT NONE
626
627	INTEGER(iwp) :: i !<
628	INTEGER(iwp) :: k !<
629
630
631	!
632	!-- Splitting
633	DO i = 0, nx
634	tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
635	ENDDO
636	DO k = 1, nz-1
637	DO i = 0, nx
638	tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
639	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)
640	ENDDO
641	ENDDO
642
643	END SUBROUTINE split_1dd
644
645
646	!------------------------------------------------------------------------------!
647	! Description:
648	! ------------
649	!> Substitution (Forward and Backward) (Thomas algorithm)
650	!------------------------------------------------------------------------------!
651	SUBROUTINE substi_1dd( ar, tri_for_1d )
652
653
654	IMPLICIT NONE
655
656	INTEGER(iwp) :: i !<
657	INTEGER(iwp) :: k !<
658
659	REAL(wp), DIMENSION(0:nx,nz) :: ar !<
660	REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !<
661	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
662
663	!
664	!-- Forward substitution
665	DO i = 0, nx
666	ar1(i,0) = ar(i,1)
667	ENDDO
668	DO k = 1, nz-1
669	DO i = 0, nx
670	ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
671	ENDDO
672	ENDDO
673
674	!
675	!-- Backward substitution
676	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
677	!-- by zero appearing if the pressure bc is set to neumann at the top of
678	!-- the model domain.
679	DO i = 0, nx
680	ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
681	ENDDO
682	DO k = nz-2, 0, -1
683	DO i = 0, nx
684	ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
685	/ tri_for_1d(4,i,k)
686	ENDDO
687	ENDDO
688
689	!
690	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
691	!-- The respective values of ar should be zero at all k-levels if
692	!-- acceleration of horizontally averaged vertical velocity is zero.
693	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
694	IF ( j == 0 ) THEN
695	DO k = 1, nz
696	ar(0,k) = 0.0_wp
697	ENDDO
698	ENDIF
699	ENDIF
700
701	END SUBROUTINE substi_1dd
702
703	END SUBROUTINE tridia_1dd
704
705
706	END MODULE tridia_solver

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |