Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |