Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |