Home

Context Navigation

source: palm/trunk/SOURCE/poisfft.f90 @ 183

Last change on this file since 183 was 164, checked in by raasch, 16 years ago
optimization of transpositions for 2D decompositions, workaround for using -env option with mpiexec, adjustments for lcxt4
Property svn:keywords set to `Id`
File size: 44.4 KB

Line
1	MODULE poisfft_mod
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! Arguments removed from transpose routines
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: poisfft.f90 164 2008-05-15 08:46:15Z letzel $
11	!
12	! 128 2007-10-26 13:11:14Z raasch
13	! Bugfix: wavenumber calculation for even nx in routines maketri
14	!
15	! 85 2007-05-11 09:35:14Z raasch
16	! Bugfix: work_fft*_vec removed from some PRIVATE-declarations
17	!
18	! 76 2007-03-29 00:58:32Z raasch
19	! Tridiagonal coefficients adjusted for Neumann boundary conditions both at
20	! the bottom and the top.
21	!
22	! RCS Log replace by Id keyword, revision history cleaned up
23	!
24	! Revision 1.24 2006/08/04 15:00:24 raasch
25	! Default setting of the thread number tn in case of not using OpenMP
26	!
27	! Revision 1.23 2006/02/23 12:48:38 raasch
28	! Additional compiler directive in routine tridia_1dd for preventing loop
29	! exchange on NEC-SX6
30	!
31	! Revision 1.20 2004/04/30 12:38:09 raasch
32	! Parts of former poisfft_hybrid moved to this subroutine,
33	! former subroutine changed to a module, renaming of FFT-subroutines and
34	! -module, FFTs completely substituted by calls of fft_x and fft_y,
35	! NAG fft used in the non-parallel case completely removed, l in maketri
36	! is now a 1d-array, variables passed by modules instead of using parameter
37	! lists, enlarged transposition arrays introduced
38	!
39	! Revision 1.1 1997/07/24 11:24:14 raasch
40	! Initial revision
41	!
42	!
43	! Description:
44	! ------------
45	! See below.
46	!------------------------------------------------------------------------------!
47
48	!--------------------------------------------------------------------------!
49	! poisfft !
50	! !
51	! Original version: Stephan Siano (pois3d) !
52	! !
53	! Institute of Meteorology and Climatology, University of Hannover !
54	! Germany !
55	! !
56	! Version as of July 23,1996 !
57	! !
58	! !
59	! Version for parallel computers: Siegfried Raasch !
60	! !
61	! Version as of July 03,1997 !
62	! !
63	! Solves the Poisson equation with a 2D spectral method !
64	! d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s !
65	! !
66	! Input: !
67	! real ar contains in the (nnx,nny,nnz) elements, !
68	! starting from the element (1,nys,nxl), the !
69	! values for s !
70	! real work Temporary array !
71	! !
72	! Output: !
73	! real ar contains the solution for p !
74	!--------------------------------------------------------------------------!
75
76	USE fft_xy
77	USE indices
78	USE transpose_indices
79
80	IMPLICIT NONE
81
82	PRIVATE
83	PUBLIC poisfft, poisfft_init
84
85	INTERFACE poisfft
86	MODULE PROCEDURE poisfft
87	END INTERFACE poisfft
88
89	INTERFACE poisfft_init
90	MODULE PROCEDURE poisfft_init
91	END INTERFACE poisfft_init
92
93	CONTAINS
94
95	SUBROUTINE poisfft_init
96
97	CALL fft_init
98
99	END SUBROUTINE poisfft_init
100
101
102	SUBROUTINE poisfft( ar, work )
103
104	USE cpulog
105	USE interfaces
106	USE pegrid
107
108	IMPLICIT NONE
109
110	REAL, DIMENSION(1:nza,nys:nyna,nxl:nxra) :: ar, work
111
112
113	CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
114
115	!
116	!-- Two-dimensional Fourier Transformation in x- and y-direction.
117	#if defined( __parallel )
118	IF ( pdims(2) == 1 ) THEN
119
120	!
121	!-- 1d-domain-decomposition along x:
122	!-- FFT along y and transposition y --> x
123	CALL ffty_tr_yx( ar, work, ar )
124
125	!
126	!-- FFT along x, solving the tridiagonal system and backward FFT
127	CALL fftx_tri_fftx( ar )
128
129	!
130	!-- Transposition x --> y and backward FFT along y
131	CALL tr_xy_ffty( ar, work, ar )
132
133	ELSEIF ( pdims(1) == 1 ) THEN
134
135	!
136	!-- 1d-domain-decomposition along y:
137	!-- FFT along x and transposition x --> y
138	CALL fftx_tr_xy( ar, work, ar )
139
140	!
141	!-- FFT along y, solving the tridiagonal system and backward FFT
142	CALL ffty_tri_ffty( ar )
143
144	!
145	!-- Transposition y --> x and backward FFT along x
146	CALL tr_yx_fftx( ar, work, ar )
147
148	ELSE
149
150	!
151	!-- 2d-domain-decomposition
152	!-- Transposition z --> x
153	CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
154	CALL transpose_zx( ar, work, ar )
155	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
156
157	CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
158	CALL fftxp( ar, 'forward' )
159	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
160
161	!
162	!-- Transposition x --> y
163	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
164	CALL transpose_xy( ar, work, ar )
165	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
166
167	CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
168	CALL fftyp( ar, 'forward' )
169	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
170
171	!
172	!-- Transposition y --> z
173	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
174	CALL transpose_yz( ar, work, ar )
175	CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
176
177	!
178	!-- Solve the Poisson equation in z-direction in cartesian space.
179	CALL cpu_log( log_point_s(6), 'tridia', 'start' )
180	CALL tridia( ar )
181	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
182
183	!
184	!-- Inverse Fourier Transformation
185	!-- Transposition z --> y
186	CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
187	CALL transpose_zy( ar, work, ar )
188	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
189
190	CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
191	CALL fftyp( ar, 'backward' )
192	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
193
194	!
195	!-- Transposition y --> x
196	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
197	CALL transpose_yx( ar, work, ar )
198	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
199
200	CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
201	CALL fftxp( ar, 'backward' )
202	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
203
204	!
205	!-- Transposition x --> z
206	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
207	CALL transpose_xz( ar, work, ar )
208	CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
209
210	ENDIF
211
212	#else
213
214	!
215	!-- Two-dimensional Fourier Transformation along x- and y-direction.
216	CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
217	CALL fftx( ar, 'forward' )
218	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
219	CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
220	CALL ffty( ar, 'forward' )
221	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
222
223	!
224	!-- Solve the Poisson equation in z-direction in cartesian space.
225	CALL cpu_log( log_point_s(6), 'tridia', 'start' )
226	CALL tridia( ar )
227	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
228
229	!
230	!-- Inverse Fourier Transformation.
231	CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
232	CALL ffty( ar, 'backward' )
233	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
234	CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
235	CALL fftx( ar, 'backward' )
236	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
237
238	#endif
239
240	CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
241
242	END SUBROUTINE poisfft
243
244
245
246	SUBROUTINE tridia( ar )
247
248	!------------------------------------------------------------------------------!
249	! solves the linear system of equations:
250	!
251	! -(4 pi^2(i^2/(dx^2nnx^2)+j^2/(dy^2nny^2))+
252	! 1/(dzu(k)dzw(k))+1/(dzu(k-1)dzw(k)))*p(i,j,k)+
253	! 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)
254	!
255	! by using the Thomas algorithm
256	!------------------------------------------------------------------------------!
257
258	USE arrays_3d
259
260	IMPLICIT NONE
261
262	INTEGER :: i, j, k, nnyh
263
264	REAL, DIMENSION(nxl_z:nxr_z,0:nz-1) :: ar1
265	REAL, DIMENSION(5,nxl_z:nxr_z,0:nz-1) :: tri
266
267	#if defined( __parallel )
268	REAL :: ar(nxl_z:nxr_za,nys_z:nyn_za,1:nza)
269	#else
270	REAL :: ar(1:nz,nys_z:nyn_z,nxl_z:nxr_z)
271	#endif
272
273
274	nnyh = (ny+1) / 2
275
276	!
277	!-- Define constant elements of the tridiagonal matrix.
278	DO k = 0, nz-1
279	DO i = nxl_z, nxr_z
280	tri(2,i,k) = ddzu(k+1) * ddzw(k+1)
281	tri(3,i,k) = ddzu(k+2) * ddzw(k+1)
282	ENDDO
283	ENDDO
284
285	#if defined( __parallel )
286	!
287	!-- Repeat for all y-levels.
288	DO j = nys_z, nyn_z
289	IF ( j <= nnyh ) THEN
290	CALL maketri( tri, j )
291	ELSE
292	CALL maketri( tri, ny+1-j )
293	ENDIF
294	CALL split( tri )
295	CALL substi( ar, ar1, tri, j )
296	ENDDO
297	#else
298	!
299	!-- First y-level.
300	CALL maketri( tri, nys_z )
301	CALL split( tri )
302	CALL substi( ar, ar1, tri, 0 )
303
304	!
305	!-- Further y-levels.
306	DO j = 1, nnyh - 1
307	CALL maketri( tri, j )
308	CALL split( tri )
309	CALL substi( ar, ar1, tri, j )
310	CALL substi( ar, ar1, tri, ny+1-j )
311	ENDDO
312	CALL maketri( tri, nnyh )
313	CALL split( tri )
314	CALL substi( ar, ar1, tri, nnyh+nys )
315	#endif
316
317	CONTAINS
318
319	SUBROUTINE maketri( tri, j )
320
321	!------------------------------------------------------------------------------!
322	! Computes the i- and j-dependent component of the matrix
323	!------------------------------------------------------------------------------!
324
325	USE arrays_3d
326	USE constants
327	USE control_parameters
328	USE grid_variables
329
330	IMPLICIT NONE
331
332	INTEGER :: i, j, k, nnxh
333	REAL :: a, c
334	REAL :: ll(nxl_z:nxr_z)
335	REAL :: tri(5,nxl_z:nxr_z,0:nz-1)
336
337
338	nnxh = ( nx + 1 ) / 2
339
340	!
341	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
342	!-- Fourier space. The coefficients are computed following the method of
343	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
344	!-- Siano's original version by discretizing the Poisson equation,
345	!-- before it is Fourier-transformed
346	#if defined( __parallel )
347	DO i = nxl_z, nxr_z
348	IF ( i >= 0 .AND. i <= nnxh ) THEN
349	ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
350	FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
351	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
352	FLOAT( ny+1 ) ) ) / ( dy * dy )
353	ELSE
354	ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
355	FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
356	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
357	FLOAT( ny+1 ) ) ) / ( dy * dy )
358	ENDIF
359	DO k = 0,nz-1
360	a = -1.0 * ddzu(k+2) * ddzw(k+1)
361	c = -1.0 * ddzu(k+1) * ddzw(k+1)
362	tri(1,i,k) = a + c - ll(i)
363	ENDDO
364	ENDDO
365	#else
366	DO i = 0, nnxh
367	ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / FLOAT( nx+1 ) ) ) / &
368	( dx * dx ) + &
369	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / FLOAT( ny+1 ) ) ) / &
370	( dy * dy )
371	DO k = 0, nz-1
372	a = -1.0 * ddzu(k+2) * ddzw(k+1)
373	c = -1.0 * ddzu(k+1) * ddzw(k+1)
374	tri(1,i,k) = a + c - ll(i)
375	IF ( i >= 1 .and. i < nnxh ) THEN
376	tri(1,nx+1-i,k) = tri(1,i,k)
377	ENDIF
378	ENDDO
379	ENDDO
380	#endif
381	IF ( ibc_p_b == 1 .OR. ibc_p_b == 2 ) THEN
382	DO i = nxl_z, nxr_z
383	tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
384	ENDDO
385	ENDIF
386	IF ( ibc_p_t == 1 ) THEN
387	DO i = nxl_z, nxr_z
388	tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
389	ENDDO
390	ENDIF
391
392	END SUBROUTINE maketri
393
394
395	SUBROUTINE substi( ar, ar1, tri, j )
396
397	!------------------------------------------------------------------------------!
398	! Substitution (Forward and Backward) (Thomas algorithm)
399	!------------------------------------------------------------------------------!
400
401	USE control_parameters
402
403	IMPLICIT NONE
404
405	INTEGER :: i, j, k
406	REAL :: ar1(nxl_z:nxr_z,0:nz-1)
407	REAL :: tri(5,nxl_z:nxr_z,0:nz-1)
408	#if defined( __parallel )
409	REAL :: ar(nxl_z:nxr_za,nys_z:nyn_za,1:nza)
410	#else
411	REAL :: ar(1:nz,nys_z:nyn_z,nxl_z:nxr_z)
412	#endif
413
414	!
415	!-- Forward substitution.
416	DO i = nxl_z, nxr_z
417	#if defined( __parallel )
418	ar1(i,0) = ar(i,j,1)
419	#else
420	ar1(i,0) = ar(1,j,i)
421	#endif
422	ENDDO
423	DO k = 1, nz - 1
424	DO i = nxl_z, nxr_z
425	#if defined( __parallel )
426	ar1(i,k) = ar(i,j,k+1) - tri(5,i,k) * ar1(i,k-1)
427	#else
428	ar1(i,k) = ar(k+1,j,i) - tri(5,i,k) * ar1(i,k-1)
429	#endif
430	ENDDO
431	ENDDO
432
433	!
434	!-- Backward substitution.
435	DO i = nxl_z, nxr_z
436	#if defined( __parallel )
437	ar(i,j,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
438	#else
439	ar(nz,j,i) = ar1(i,nz-1) / tri(4,i,nz-1)
440	#endif
441	ENDDO
442	DO k = nz-2, 0, -1
443	DO i = nxl_z, nxr_z
444	#if defined( __parallel )
445	ar(i,j,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,j,k+2) ) &
446	/ tri(4,i,k)
447	#else
448	ar(k+1,j,i) = ( ar1(i,k) - tri(3,i,k) * ar(k+2,j,i) ) &
449	/ tri(4,i,k)
450	#endif
451	ENDDO
452	ENDDO
453
454	!
455	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
456	!-- The respective values of ar should be zero at all k-levels if
457	!-- acceleration of horizontally averaged vertical velocity is zero.
458	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
459	IF ( j == 0 .AND. nxl_z == 0 ) THEN
460	#if defined( __parallel )
461	DO k = 1, nz
462	ar(nxl_z,j,k) = 0.0
463	ENDDO
464	#else
465	DO k = 1, nz
466	ar(k,j,nxl_z) = 0.0
467	ENDDO
468	#endif
469	ENDIF
470	ENDIF
471
472	END SUBROUTINE substi
473
474
475	SUBROUTINE split( tri )
476
477	!------------------------------------------------------------------------------!
478	! Splitting of the tridiagonal matrix (Thomas algorithm)
479	!------------------------------------------------------------------------------!
480
481	IMPLICIT NONE
482
483	INTEGER :: i, k
484	REAL :: tri(5,nxl_z:nxr_z,0:nz-1)
485
486	!
487	!-- Splitting.
488	DO i = nxl_z, nxr_z
489	tri(4,i,0) = tri(1,i,0)
490	ENDDO
491	DO k = 1, nz-1
492	DO i = nxl_z, nxr_z
493	tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
494	tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
495	ENDDO
496	ENDDO
497
498	END SUBROUTINE split
499
500	END SUBROUTINE tridia
501
502
503	#if defined( __parallel )
504	SUBROUTINE fftxp( ar, direction )
505
506	!------------------------------------------------------------------------------!
507	! Fourier-transformation along x-direction Parallelized version
508	!------------------------------------------------------------------------------!
509
510	IMPLICIT NONE
511
512	CHARACTER (LEN=*) :: direction
513	INTEGER :: j, k
514	REAL :: ar(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa)
515
516	!
517	!-- Performing the fft with one of the methods implemented
518	DO k = nzb_x, nzt_x
519	DO j = nys_x, nyn_x
520	CALL fft_x( ar(0:nx,j,k), direction )
521	ENDDO
522	ENDDO
523
524	END SUBROUTINE fftxp
525
526	#else
527	SUBROUTINE fftx( ar, direction )
528
529	!------------------------------------------------------------------------------!
530	! Fourier-transformation along x-direction Non parallel version
531	!------------------------------------------------------------------------------!
532
533	IMPLICIT NONE
534
535	CHARACTER (LEN=*) :: direction
536	INTEGER :: i, j, k
537	REAL :: ar(1:nz,0:ny,0:nx)
538
539	!
540	!-- Performing the fft with one of the methods implemented
541	DO k = 1, nz
542	DO j = 0, ny
543	CALL fft_x( ar(k,j,0:nx), direction )
544	ENDDO
545	ENDDO
546
547	END SUBROUTINE fftx
548	#endif
549
550
551	#if defined( __parallel )
552	SUBROUTINE fftyp( ar, direction )
553
554	!------------------------------------------------------------------------------!
555	! Fourier-transformation along y-direction Parallelized version
556	!------------------------------------------------------------------------------!
557
558	IMPLICIT NONE
559
560	CHARACTER (LEN=*) :: direction
561	INTEGER :: i, k
562	REAL :: ar(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya)
563
564	!
565	!-- Performing the fft with one of the methods implemented
566	DO k = nzb_y, nzt_y
567	DO i = nxl_y, nxr_y
568	CALL fft_y( ar(0:ny,i,k), direction )
569	ENDDO
570	ENDDO
571
572	END SUBROUTINE fftyp
573
574	#else
575	SUBROUTINE ffty( ar, direction )
576
577	!------------------------------------------------------------------------------!
578	! Fourier-transformation along y-direction Non parallel version
579	!------------------------------------------------------------------------------!
580
581	IMPLICIT NONE
582
583	CHARACTER (LEN=*) :: direction
584	INTEGER :: i, k
585	REAL :: ar(1:nz,0:ny,0:nx)
586
587	!
588	!-- Performing the fft with one of the methods implemented
589	DO k = 1, nz
590	DO i = 0, nx
591	CALL fft_y( ar(k,0:ny,i), direction )
592	ENDDO
593	ENDDO
594
595	END SUBROUTINE ffty
596	#endif
597
598	#if defined( __parallel )
599	SUBROUTINE ffty_tr_yx( f_in, work, f_out )
600
601	!------------------------------------------------------------------------------!
602	! Fourier-transformation along y with subsequent transposition y --> x for
603	! a 1d-decomposition along x
604	!
605	! ATTENTION: The performance of this routine is much faster on the NEC-SX6,
606	! if the first index of work_ffty_vec is odd. Otherwise
607	! memory bank conflicts may occur (especially if the index is a
608	! multiple of 128). That's why work_ffty_vec is dimensioned as
609	! 0:ny+1.
610	! Of course, this will not work if users are using an odd number
611	! of gridpoints along y.
612	!------------------------------------------------------------------------------!
613
614	USE control_parameters
615	USE cpulog
616	USE indices
617	USE interfaces
618	USE pegrid
619	USE transpose_indices
620
621	IMPLICIT NONE
622
623	INTEGER :: i, iend, iouter, ir, j, k
624	INTEGER, PARAMETER :: stridex = 4
625
626	REAL, DIMENSION(0:ny,stridex) :: work_ffty
627	#if defined( __nec )
628	REAL, DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec
629	#endif
630	REAL, DIMENSION(1:nza,0:nya,nxl:nxra) :: f_in
631	REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) :: f_out
632	REAL, DIMENSION(nxl:nxra,1:nza,0:nya) :: work
633
634	!
635	!-- Carry out the FFT along y, where all data are present due to the
636	!-- 1d-decomposition along x. Resort the data in a way that x becomes
637	!-- the first index.
638	CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
639
640	IF ( host(1:3) == 'nec' ) THEN
641	#if defined( __nec )
642	!
643	!-- Code optimized for vector processors
644	!$OMP PARALLEL PRIVATE ( i, j, k )
645	!$OMP DO
646	DO i = nxl, nxr
647
648	DO j = 0, ny
649	DO k = 1, nz
650	work_ffty_vec(j,k,i) = f_in(k,j,i)
651	ENDDO
652	ENDDO
653
654	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
655
656	ENDDO
657
658	!$OMP DO
659	DO k = 1, nz
660	DO j = 0, ny
661	DO i = nxl, nxr
662	work(i,k,j) = work_ffty_vec(j,k,i)
663	ENDDO
664	ENDDO
665	ENDDO
666	!$OMP END PARALLEL
667	#endif
668
669	ELSE
670
671	!
672	!-- Cache optimized code.
673	!-- The i-(x-)direction is split into a strided outer loop and an inner
674	!-- loop for better cache performance
675	!$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
676	!$OMP DO
677	DO iouter = nxl, nxr, stridex
678
679	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
680
681	DO k = 1, nz
682
683	DO i = iouter, iend
684
685	ir = i-iouter+1 ! counter within a stride
686	DO j = 0, ny
687	work_ffty(j,ir) = f_in(k,j,i)
688	ENDDO
689	!
690	!-- FFT along y
691	CALL fft_y( work_ffty(:,ir), 'forward' )
692
693	ENDDO
694
695	!
696	!-- Resort
697	DO j = 0, ny
698	DO i = iouter, iend
699	work(i,k,j) = work_ffty(j,i-iouter+1)
700	ENDDO
701	ENDDO
702
703	ENDDO
704
705	ENDDO
706	!$OMP END PARALLEL
707
708	ENDIF
709	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
710
711	!
712	!-- Transpose array
713	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
714	CALL MPI_ALLTOALL( work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
715	f_out(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
716	comm1dx, ierr )
717	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
718
719	END SUBROUTINE ffty_tr_yx
720
721
722	SUBROUTINE tr_xy_ffty( f_in, work, f_out )
723
724	!------------------------------------------------------------------------------!
725	! Transposition x --> y with a subsequent backward Fourier transformation for
726	! a 1d-decomposition along x
727	!------------------------------------------------------------------------------!
728
729	USE control_parameters
730	USE cpulog
731	USE indices
732	USE interfaces
733	USE pegrid
734	USE transpose_indices
735
736	IMPLICIT NONE
737
738	INTEGER :: i, iend, iouter, ir, j, k
739	INTEGER, PARAMETER :: stridex = 4
740
741	REAL, DIMENSION(0:ny,stridex) :: work_ffty
742	#if defined( __nec )
743	REAL, DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec
744	#endif
745	REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) :: f_in
746	REAL, DIMENSION(1:nza,0:nya,nxl:nxra) :: f_out
747	REAL, DIMENSION(nxl:nxra,1:nza,0:nya) :: work
748
749	!
750	!-- Transpose array
751	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
752	CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
753	work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
754	comm1dx, ierr )
755	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
756
757	!
758	!-- Resort the data in a way that y becomes the first index and carry out the
759	!-- backward fft along y.
760	CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
761
762	IF ( host(1:3) == 'nec' ) THEN
763	#if defined( __nec )
764	!
765	!-- Code optimized for vector processors
766	!$OMP PARALLEL PRIVATE ( i, j, k )
767	!$OMP DO
768	DO k = 1, nz
769	DO j = 0, ny
770	DO i = nxl, nxr
771	work_ffty_vec(j,k,i) = work(i,k,j)
772	ENDDO
773	ENDDO
774	ENDDO
775
776	!$OMP DO
777	DO i = nxl, nxr
778
779	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
780
781	DO j = 0, ny
782	DO k = 1, nz
783	f_out(k,j,i) = work_ffty_vec(j,k,i)
784	ENDDO
785	ENDDO
786
787	ENDDO
788	!$OMP END PARALLEL
789	#endif
790
791	ELSE
792
793	!
794	!-- Cache optimized code.
795	!-- The i-(x-)direction is split into a strided outer loop and an inner
796	!-- loop for better cache performance
797	!$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
798	!$OMP DO
799	DO iouter = nxl, nxr, stridex
800
801	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
802
803	DO k = 1, nz
804	!
805	!-- Resort
806	DO j = 0, ny
807	DO i = iouter, iend
808	work_ffty(j,i-iouter+1) = work(i,k,j)
809	ENDDO
810	ENDDO
811
812	DO i = iouter, iend
813
814	!
815	!-- FFT along y
816	ir = i-iouter+1 ! counter within a stride
817	CALL fft_y( work_ffty(:,ir), 'backward' )
818
819	DO j = 0, ny
820	f_out(k,j,i) = work_ffty(j,ir)
821	ENDDO
822	ENDDO
823
824	ENDDO
825
826	ENDDO
827	!$OMP END PARALLEL
828
829	ENDIF
830
831	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
832
833	END SUBROUTINE tr_xy_ffty
834
835
836	SUBROUTINE fftx_tri_fftx( ar )
837
838	!------------------------------------------------------------------------------!
839	! FFT along x, solution of the tridiagonal system and backward FFT for
840	! a 1d-decomposition along x
841	!
842	! WARNING: this subroutine may still not work for hybrid parallelization
843	! with OpenMP (for possible necessary changes see the original
844	! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
845	!------------------------------------------------------------------------------!
846
847	USE control_parameters
848	USE cpulog
849	USE grid_variables
850	USE indices
851	USE interfaces
852	USE pegrid
853	USE transpose_indices
854
855	IMPLICIT NONE
856
857	character(len=3) :: myth_char
858
859	INTEGER :: i, j, k, m, n, omp_get_thread_num, tn
860
861	REAL, DIMENSION(0:nx) :: work_fftx
862	REAL, DIMENSION(0:nx,1:nz) :: work_trix
863	REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) :: ar
864	REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: tri
865
866
867	CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
868
869	ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
870
871	tn = 0 ! Default thread number in case of one thread
872	!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
873	DO j = nys_x, nyn_x
874
875	!$ tn = omp_get_thread_num()
876
877	IF ( host(1:3) == 'nec' ) THEN
878	!
879	!-- Code optimized for vector processors
880	DO k = 1, nz
881
882	m = 0
883	DO n = 1, pdims(1)
884	DO i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
885	work_trix(m,k) = ar(i,k,j,n)
886	m = m + 1
887	ENDDO
888	ENDDO
889
890	ENDDO
891
892	CALL fft_x_m( work_trix, 'forward' )
893
894	ELSE
895	!
896	!-- Cache optimized code
897	DO k = 1, nz
898
899	m = 0
900	DO n = 1, pdims(1)
901	DO i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
902	work_fftx(m) = ar(i,k,j,n)
903	m = m + 1
904	ENDDO
905	ENDDO
906
907	CALL fft_x( work_fftx, 'forward' )
908
909	DO i = 0, nx
910	work_trix(i,k) = work_fftx(i)
911	ENDDO
912
913	ENDDO
914
915	ENDIF
916
917	!
918	!-- Solve the linear equation system
919	CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
920
921	IF ( host(1:3) == 'nec' ) THEN
922	!
923	!-- Code optimized for vector processors
924	CALL fft_x_m( work_trix, 'backward' )
925
926	DO k = 1, nz
927
928	m = 0
929	DO n = 1, pdims(1)
930	DO i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
931	ar(i,k,j,n) = work_trix(m,k)
932	m = m + 1
933	ENDDO
934	ENDDO
935
936	ENDDO
937
938	ELSE
939	!
940	!-- Cache optimized code
941	DO k = 1, nz
942
943	DO i = 0, nx
944	work_fftx(i) = work_trix(i,k)
945	ENDDO
946
947	CALL fft_x( work_fftx, 'backward' )
948
949	m = 0
950	DO n = 1, pdims(1)
951	DO i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
952	ar(i,k,j,n) = work_fftx(m)
953	m = m + 1
954	ENDDO
955	ENDDO
956
957	ENDDO
958
959	ENDIF
960
961	ENDDO
962
963	DEALLOCATE( tri )
964
965	CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
966
967	END SUBROUTINE fftx_tri_fftx
968
969
970	SUBROUTINE fftx_tr_xy( f_in, work, f_out )
971
972	!------------------------------------------------------------------------------!
973	! Fourier-transformation along x with subsequent transposition x --> y for
974	! a 1d-decomposition along y
975	!
976	! ATTENTION: The NEC-branch of this routine may significantly profit from
977	! further optimizations. So far, performance is much worse than
978	! for routine ffty_tr_yx (more than three times slower).
979	!------------------------------------------------------------------------------!
980
981	USE control_parameters
982	USE cpulog
983	USE indices
984	USE interfaces
985	USE pegrid
986	USE transpose_indices
987
988	IMPLICIT NONE
989
990	INTEGER :: i, j, k
991
992	REAL, DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx
993	REAL, DIMENSION(1:nza,nys:nyna,0:nxa) :: f_in
994	REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) :: f_out
995	REAL, DIMENSION(nys:nyna,1:nza,0:nxa) :: work
996
997	!
998	!-- Carry out the FFT along x, where all data are present due to the
999	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1000	!-- the first index.
1001	CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
1002
1003	IF ( host(1:3) == 'nec' ) THEN
1004	!
1005	!-- Code for vector processors
1006	!$OMP PARALLEL PRIVATE ( i, j, k )
1007	!$OMP DO
1008	DO i = 0, nx
1009
1010	DO j = nys, nyn
1011	DO k = 1, nz
1012	work_fftx(i,k,j) = f_in(k,j,i)
1013	ENDDO
1014	ENDDO
1015
1016	ENDDO
1017
1018	!$OMP DO
1019	DO j = nys, nyn
1020
1021	CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1022
1023	DO k = 1, nz
1024	DO i = 0, nx
1025	work(j,k,i) = work_fftx(i,k,j)
1026	ENDDO
1027	ENDDO
1028
1029	ENDDO
1030	!$OMP END PARALLEL
1031
1032	ELSE
1033
1034	!
1035	!-- Cache optimized code (there might be still a potential for better
1036	!-- optimization).
1037	!$OMP PARALLEL PRIVATE (i,j,k,work_fftx)
1038	!$OMP DO
1039	DO i = 0, nx
1040
1041	DO j = nys, nyn
1042	DO k = 1, nz
1043	work_fftx(i,k,j) = f_in(k,j,i)
1044	ENDDO
1045	ENDDO
1046
1047	ENDDO
1048
1049	!$OMP DO
1050	DO j = nys, nyn
1051	DO k = 1, nz
1052
1053	CALL fft_x( work_fftx(0:nx,k,j), 'forward' )
1054
1055	DO i = 0, nx
1056	work(j,k,i) = work_fftx(i,k,j)
1057	ENDDO
1058	ENDDO
1059
1060	ENDDO
1061	!$OMP END PARALLEL
1062
1063	ENDIF
1064	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
1065
1066	!
1067	!-- Transpose array
1068	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1069	CALL MPI_ALLTOALL( work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1070	f_out(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1071	comm1dy, ierr )
1072	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1073
1074	END SUBROUTINE fftx_tr_xy
1075
1076
1077	SUBROUTINE tr_yx_fftx( f_in, work, f_out )
1078
1079	!------------------------------------------------------------------------------!
1080	! Transposition y --> x with a subsequent backward Fourier transformation for
1081	! a 1d-decomposition along x
1082	!------------------------------------------------------------------------------!
1083
1084	USE control_parameters
1085	USE cpulog
1086	USE indices
1087	USE interfaces
1088	USE pegrid
1089	USE transpose_indices
1090
1091	IMPLICIT NONE
1092
1093	INTEGER :: i, j, k
1094
1095	REAL, DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx
1096	REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) :: f_in
1097	REAL, DIMENSION(1:nza,nys:nyna,0:nxa) :: f_out
1098	REAL, DIMENSION(nys:nyna,1:nza,0:nxa) :: work
1099
1100	!
1101	!-- Transpose array
1102	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1103	CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1104	work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1105	comm1dy, ierr )
1106	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1107
1108	!
1109	!-- Carry out the FFT along x, where all data are present due to the
1110	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1111	!-- the first index.
1112	CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
1113
1114	IF ( host(1:3) == 'nec' ) THEN
1115	!
1116	!-- Code optimized for vector processors
1117	!$OMP PARALLEL PRIVATE ( i, j, k )
1118	!$OMP DO
1119	DO j = nys, nyn
1120
1121	DO k = 1, nz
1122	DO i = 0, nx
1123	work_fftx(i,k,j) = work(j,k,i)
1124	ENDDO
1125	ENDDO
1126
1127	CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1128
1129	ENDDO
1130
1131	!$OMP DO
1132	DO i = 0, nx
1133	DO j = nys, nyn
1134	DO k = 1, nz
1135	f_out(k,j,i) = work_fftx(i,k,j)
1136	ENDDO
1137	ENDDO
1138	ENDDO
1139	!$OMP END PARALLEL
1140
1141	ELSE
1142
1143	!
1144	!-- Cache optimized code (there might be still a potential for better
1145	!-- optimization).
1146	!$OMP PARALLEL PRIVATE (i,j,k,work_fftx)
1147	!$OMP DO
1148	DO j = nys, nyn
1149	DO k = 1, nz
1150
1151	DO i = 0, nx
1152	work_fftx(i,k,j) = work(j,k,i)
1153	ENDDO
1154
1155	CALL fft_x( work_fftx(0:nx,k,j), 'backward' )
1156
1157	ENDDO
1158	ENDDO
1159
1160	!$OMP DO
1161	DO i = 0, nx
1162	DO j = nys, nyn
1163	DO k = 1, nz
1164	f_out(k,j,i) = work_fftx(i,k,j)
1165	ENDDO
1166	ENDDO
1167	ENDDO
1168	!$OMP END PARALLEL
1169
1170	ENDIF
1171	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
1172
1173	END SUBROUTINE tr_yx_fftx
1174
1175
1176	SUBROUTINE ffty_tri_ffty( ar )
1177
1178	!------------------------------------------------------------------------------!
1179	! FFT along y, solution of the tridiagonal system and backward FFT for
1180	! a 1d-decomposition along y
1181	!
1182	! WARNING: this subroutine may still not work for hybrid parallelization
1183	! with OpenMP (for possible necessary changes see the original
1184	! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
1185	!------------------------------------------------------------------------------!
1186
1187	USE control_parameters
1188	USE cpulog
1189	USE grid_variables
1190	USE indices
1191	USE interfaces
1192	USE pegrid
1193	USE transpose_indices
1194
1195	IMPLICIT NONE
1196
1197	INTEGER :: i, j, k, m, n, omp_get_thread_num, tn
1198
1199	REAL, DIMENSION(0:ny) :: work_ffty
1200	REAL, DIMENSION(0:ny,1:nz) :: work_triy
1201	REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) :: ar
1202	REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: tri
1203
1204
1205	CALL cpu_log( log_point_s(39), 'fft_y + tridia', 'start' )
1206
1207	ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1208
1209	tn = 0 ! Default thread number in case of one thread
1210	!$OMP PARALLEL PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1211	!$OMP DO
1212	DO i = nxl_y, nxr_y
1213
1214	!$ tn = omp_get_thread_num()
1215
1216	IF ( host(1:3) == 'nec' ) THEN
1217	!
1218	!-- Code optimized for vector processors
1219	DO k = 1, nz
1220
1221	m = 0
1222	DO n = 1, pdims(2)
1223	DO j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1224	work_triy(m,k) = ar(j,k,i,n)
1225	m = m + 1
1226	ENDDO
1227	ENDDO
1228
1229	ENDDO
1230
1231	CALL fft_y_m( work_triy, ny, 'forward' )
1232
1233	ELSE
1234	!
1235	!-- Cache optimized code
1236	DO k = 1, nz
1237
1238	m = 0
1239	DO n = 1, pdims(2)
1240	DO j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1241	work_ffty(m) = ar(j,k,i,n)
1242	m = m + 1
1243	ENDDO
1244	ENDDO
1245
1246	CALL fft_y( work_ffty, 'forward' )
1247
1248	DO j = 0, ny
1249	work_triy(j,k) = work_ffty(j)
1250	ENDDO
1251
1252	ENDDO
1253
1254	ENDIF
1255
1256	!
1257	!-- Solve the linear equation system
1258	CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1259
1260	IF ( host(1:3) == 'nec' ) THEN
1261	!
1262	!-- Code optimized for vector processors
1263	CALL fft_y_m( work_triy, ny, 'backward' )
1264
1265	DO k = 1, nz
1266
1267	m = 0
1268	DO n = 1, pdims(2)
1269	DO j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1270	ar(j,k,i,n) = work_triy(m,k)
1271	m = m + 1
1272	ENDDO
1273	ENDDO
1274
1275	ENDDO
1276
1277	ELSE
1278	!
1279	!-- Cache optimized code
1280	DO k = 1, nz
1281
1282	DO j = 0, ny
1283	work_ffty(j) = work_triy(j,k)
1284	ENDDO
1285
1286	CALL fft_y( work_ffty, 'backward' )
1287
1288	m = 0
1289	DO n = 1, pdims(2)
1290	DO j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1291	ar(j,k,i,n) = work_ffty(m)
1292	m = m + 1
1293	ENDDO
1294	ENDDO
1295
1296	ENDDO
1297
1298	ENDIF
1299
1300	ENDDO
1301	!$OMP END PARALLEL
1302
1303	DEALLOCATE( tri )
1304
1305	CALL cpu_log( log_point_s(39), 'fft_y + tridia', 'stop' )
1306
1307	END SUBROUTINE ffty_tri_ffty
1308
1309
1310	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
1311
1312	!------------------------------------------------------------------------------!
1313	! Solves the linear system of equations for a 1d-decomposition along x (see
1314	! tridia)
1315	!
1316	! Attention: when using the intel compiler, array tri must be passed as an
1317	! argument to the contained subroutines. Otherwise addres faults
1318	! will occur.
1319	! On NEC, tri should not be passed (except for routine substi_1dd)
1320	! because this causes very bad performance.
1321	!------------------------------------------------------------------------------!
1322
1323	USE arrays_3d
1324	USE control_parameters
1325
1326	USE pegrid
1327
1328	IMPLICIT NONE
1329
1330	INTEGER :: i, j, k, nnyh, nx, ny, omp_get_thread_num, tn
1331
1332	REAL :: ddx2, ddy2
1333
1334	REAL, DIMENSION(0:nx,1:nz) :: ar
1335	REAL, DIMENSION(0:nx,0:nz-1) :: ar1
1336	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
1337
1338
1339	nnyh = ( ny + 1 ) / 2
1340
1341	!
1342	!-- Define constant elements of the tridiagonal matrix.
1343	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
1344	!-- the exchanged loops create bank conflicts. The following directive
1345	!-- prohibits loop exchange and the loops perform much better.
1346	! tn = omp_get_thread_num()
1347	! WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
1348	! CALL local_flush( 120+tn )
1349	!CDIR NOLOOPCHG
1350	DO k = 0, nz-1
1351	DO i = 0,nx
1352	tri(2,i,k) = ddzu(k+1) * ddzw(k+1)
1353	tri(3,i,k) = ddzu(k+2) * ddzw(k+1)
1354	ENDDO
1355	ENDDO
1356	! WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop thread=', omp_get_thread_num()
1357	! CALL local_flush( 120+tn )
1358
1359	IF ( j <= nnyh ) THEN
1360	#if defined( __lcmuk )
1361	CALL maketri_1dd( j, tri )
1362	#else
1363	CALL maketri_1dd( j )
1364	#endif
1365	ELSE
1366	#if defined( __lcmuk )
1367	CALL maketri_1dd( ny+1-j, tri )
1368	#else
1369	CALL maketri_1dd( ny+1-j )
1370	#endif
1371	ENDIF
1372	#if defined( __lcmuk )
1373	CALL split_1dd( tri )
1374	#else
1375	CALL split_1dd
1376	#endif
1377	CALL substi_1dd( ar, tri )
1378
1379	CONTAINS
1380
1381	#if defined( __lcmuk )
1382	SUBROUTINE maketri_1dd( j, tri )
1383	#else
1384	SUBROUTINE maketri_1dd( j )
1385	#endif
1386
1387	!------------------------------------------------------------------------------!
1388	! computes the i- and j-dependent component of the matrix
1389	!------------------------------------------------------------------------------!
1390
1391	USE constants
1392
1393	IMPLICIT NONE
1394
1395	INTEGER :: i, j, k, nnxh
1396	REAL :: a, c
1397
1398	REAL, DIMENSION(0:nx) :: l
1399
1400	#if defined( __lcmuk )
1401	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
1402	#endif
1403
1404
1405	nnxh = ( nx + 1 ) / 2
1406	!
1407	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
1408	!-- Fourier space. The coefficients are computed following the method of
1409	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
1410	!-- Siano's original version by discretizing the Poisson equation,
1411	!-- before it is Fourier-transformed
1412	DO i = 0, nx
1413	IF ( i >= 0 .AND. i <= nnxh ) THEN
1414	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
1415	FLOAT( nx+1 ) ) ) * ddx2 + &
1416	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
1417	FLOAT( ny+1 ) ) ) * ddy2
1418	ELSE
1419	l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
1420	FLOAT( nx+1 ) ) ) * ddx2 + &
1421	2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
1422	FLOAT( ny+1 ) ) ) * ddy2
1423	ENDIF
1424	ENDDO
1425
1426	DO k = 0, nz-1
1427	DO i = 0, nx
1428	a = -1.0 * ddzu(k+2) * ddzw(k+1)
1429	c = -1.0 * ddzu(k+1) * ddzw(k+1)
1430	tri(1,i,k) = a + c - l(i)
1431	ENDDO
1432	ENDDO
1433	IF ( ibc_p_b == 1 .OR. ibc_p_b == 2 ) THEN
1434	DO i = 0, nx
1435	tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
1436	ENDDO
1437	ENDIF
1438	IF ( ibc_p_t == 1 ) THEN
1439	DO i = 0, nx
1440	tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
1441	ENDDO
1442	ENDIF
1443
1444	END SUBROUTINE maketri_1dd
1445
1446
1447	#if defined( __lcmuk )
1448	SUBROUTINE split_1dd( tri )
1449	#else
1450	SUBROUTINE split_1dd
1451	#endif
1452
1453	!------------------------------------------------------------------------------!
1454	! Splitting of the tridiagonal matrix (Thomas algorithm)
1455	!------------------------------------------------------------------------------!
1456
1457	IMPLICIT NONE
1458
1459	INTEGER :: i, k
1460
1461	#if defined( __lcmuk )
1462	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
1463	#endif
1464
1465
1466	!
1467	!-- Splitting
1468	DO i = 0, nx
1469	tri(4,i,0) = tri(1,i,0)
1470	ENDDO
1471	DO k = 1, nz-1
1472	DO i = 0, nx
1473	tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
1474	tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
1475	ENDDO
1476	ENDDO
1477
1478	END SUBROUTINE split_1dd
1479
1480
1481	SUBROUTINE substi_1dd( ar, tri )
1482
1483	!------------------------------------------------------------------------------!
1484	! Substitution (Forward and Backward) (Thomas algorithm)
1485	!------------------------------------------------------------------------------!
1486
1487	IMPLICIT NONE
1488
1489	INTEGER :: i, k
1490
1491	REAL, DIMENSION(0:nx,nz) :: ar
1492	REAL, DIMENSION(0:nx,0:nz-1) :: ar1
1493	REAL, DIMENSION(5,0:nx,0:nz-1) :: tri
1494
1495	!
1496	!-- Forward substitution
1497	DO i = 0, nx
1498	ar1(i,0) = ar(i,1)
1499	ENDDO
1500	DO k = 1, nz-1
1501	DO i = 0, nx
1502	ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
1503	ENDDO
1504	ENDDO
1505
1506	!
1507	!-- Backward substitution
1508	DO i = 0, nx
1509	ar(i,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
1510	ENDDO
1511	DO k = nz-2, 0, -1
1512	DO i = 0, nx
1513	ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
1514	/ tri(4,i,k)
1515	ENDDO
1516	ENDDO
1517
1518	!
1519	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
1520	!-- The respective values of ar should be zero at all k-levels if
1521	!-- acceleration of horizontally averaged vertical velocity is zero.
1522	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
1523	IF ( j == 0 ) THEN
1524	DO k = 1, nz
1525	ar(0,k) = 0.0
1526	ENDDO
1527	ENDIF
1528	ENDIF
1529
1530	END SUBROUTINE substi_1dd
1531
1532	END SUBROUTINE tridia_1dd
1533
1534	#endif
1535
1536	END MODULE poisfft_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |