Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy.f90 @ 1210

Last change on this file since 1210 was 1210, checked in by raasch, 11 years ago
fftw support added; object file list in Makefile replaced by a short one line statement
Property svn:keywords set to `Id`
File size: 49.7 KB

Line
1	MODULE fft_xy
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later 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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! fftw added
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: fft_xy.f90 1210 2013-08-14 10:58:20Z raasch $
27	!
28	! 1166 2013-05-24 13:55:44Z raasch
29	! C_DOUBLE/COMPLEX reset to dpk
30	!
31	! 1153 2013-05-10 14:33:08Z raasch
32	! code adjustment of data types for CUDA fft required by PGI 12.3 / CUDA 5.0
33	!
34	! 1111 2013-03-08 23:54:10Z raasch
35	! further openACC statements added, CUDA branch completely runs on GPU
36	! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
37	! used total domain)
38	!
39	! 1106 2013-03-04 05:31:38Z raasch
40	! CUDA fft added
41	! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
42	! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
43	!
44	! 1092 2013-02-02 11:24:22Z raasch
45	! variable sizw declared for NEC case only
46	!
47	! 1036 2012-10-22 13:43:42Z raasch
48	! code put under GPL (PALM 3.9)
49	!
50	! 274 2009-03-26 15:11:21Z heinze
51	! Output of messages replaced by message handling routine.
52	!
53	! Feb. 2007
54	! RCS Log replace by Id keyword, revision history cleaned up
55	!
56	! Revision 1.4 2006/03/28 12:27:09 raasch
57	! Stop when system-specific fft is selected on NEC. For unknown reasons this
58	! causes a program abort during first allocation in init_grid.
59	!
60	! Revision 1.2 2004/04/30 11:44:27 raasch
61	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
62	! fft_x and fft_y,
63	! function FFT replaced by subroutine FFTN due to problems with 64-bit
64	! mode on ibm,
65	! shape of array cwork is explicitly stored in ishape/jshape and handled
66	! to routine FFTN instead of shape-function (due to compiler error on
67	! decalpha),
68	! non vectorized FFT for nec included
69	!
70	! Revision 1.1 2002/06/11 13:00:49 raasch
71	! Initial revision
72	!
73	!
74	! Description:
75	! ------------
76	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
77	! Original version: Klaus Ketelsen (May 2002)
78	!------------------------------------------------------------------------------!
79
80	USE control_parameters
81	USE indices
82	#if defined( __cuda_fft )
83	USE ISO_C_BINDING
84	#elif defined( __fftw )
85	USE, INTRINSIC :: ISO_C_BINDING
86	#endif
87	USE precision_kind
88	USE singleton
89	USE temperton_fft
90	USE transpose_indices
91
92	IMPLICIT NONE
93
94	PRIVATE
95	PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
96
97	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
98
99	LOGICAL, SAVE :: init_fft = .FALSE.
100
101	REAL, SAVE :: dnx, dny, sqr_dnx, sqr_dny
102	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
103
104	#if defined( __ibm )
105	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
106	!
107	!-- The following working arrays contain tables and have to be "save" and
108	!-- shared in OpenMP sense
109	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
110	#elif defined( __nec )
111	INTEGER, SAVE :: nz1
112	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
113	trig_yf
114	#elif defined( __cuda_fft )
115	INTEGER(C_INT), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
116	INTEGER, SAVE :: total_points_x_transpo, total_points_y_transpo
117	#elif defined( __fftw )
118	INCLUDE 'fftw3.f03'
119	INTEGER(KIND=C_INT) :: nx_c, ny_c
120	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: x_out, y_out
121	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: x_in, y_in
122	TYPE(C_PTR), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
123	#endif
124
125	!
126	!-- Public interfaces
127	INTERFACE fft_init
128	MODULE PROCEDURE fft_init
129	END INTERFACE fft_init
130
131	INTERFACE fft_x
132	MODULE PROCEDURE fft_x
133	END INTERFACE fft_x
134
135	INTERFACE fft_x_1d
136	MODULE PROCEDURE fft_x_1d
137	END INTERFACE fft_x_1d
138
139	INTERFACE fft_y
140	MODULE PROCEDURE fft_y
141	END INTERFACE fft_y
142
143	INTERFACE fft_y_1d
144	MODULE PROCEDURE fft_y_1d
145	END INTERFACE fft_y_1d
146
147	INTERFACE fft_x_m
148	MODULE PROCEDURE fft_x_m
149	END INTERFACE fft_x_m
150
151	INTERFACE fft_y_m
152	MODULE PROCEDURE fft_y_m
153	END INTERFACE fft_y_m
154
155	CONTAINS
156
157
158	SUBROUTINE fft_init
159
160	USE cuda_fft_interfaces
161
162	IMPLICIT NONE
163
164	!
165	!-- The following temporary working arrays have to be on stack or private
166	!-- in OpenMP sense
167	#if defined( __ibm )
168	REAL, DIMENSION(0:nx+2) :: workx
169	REAL, DIMENSION(0:ny+2) :: worky
170	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
171	#elif defined( __nec )
172	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
173	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
174	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
175	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
176	#endif
177
178	!
179	!-- Return, if already called
180	IF ( init_fft ) THEN
181	RETURN
182	ELSE
183	init_fft = .TRUE.
184	ENDIF
185
186	IF ( fft_method == 'system-specific' ) THEN
187
188	dnx = 1.0 / ( nx + 1.0 )
189	dny = 1.0 / ( ny + 1.0 )
190	sqr_dnx = SQRT( dnx )
191	sqr_dny = SQRT( dny )
192	#if defined( __ibm ) && ! defined( __ibmy_special )
193	!
194	!-- Initialize tables for fft along x
195	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
196	aux2, nau2 )
197	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
198	aux4, nau2 )
199	!
200	!-- Initialize tables for fft along y
201	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
202	auy2, nau2 )
203	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
204	auy4, nau2 )
205	#elif defined( __nec )
206	message_string = 'fft method "' // TRIM( fft_method) // &
207	'" currently does not work on NEC'
208	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
209
210	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
211	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
212
213	work_x = 0.0
214	work_y = 0.0
215	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
216	! when using the NEC ffts
217
218	!
219	!-- Initialize tables for fft along x (non-vector and vector case (M))
220	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
221	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
222	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
223	trig_xf, workx, 0 )
224	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
225	trig_xb, workx, 0 )
226	!
227	!-- Initialize tables for fft along y (non-vector and vector case (M))
228	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
229	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
230	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
231	trig_yf, worky, 0 )
232	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
233	trig_yb, worky, 0 )
234	#elif defined( __cuda_fft )
235	total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
236	total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
237	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
238	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
239	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
240	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
241	#else
242	message_string = 'no system-specific fft-call available'
243	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
244	#endif
245	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
246	!
247	!-- Temperton-algorithm
248	!-- Initialize tables for fft along x and y
249	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
250
251	CALL set99( trigs_x, ifax_x, nx+1 )
252	CALL set99( trigs_y, ifax_y, ny+1 )
253
254	ELSEIF ( fft_method == 'fftw' ) THEN
255	!
256	!-- FFTW
257	#if defined( __fftw )
258	nx_c = nx+1
259	ny_c = ny+1
260	ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), &
261	y_out(0:(ny+1)/2) )
262	plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
263	plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
264	plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
265	plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
266	#else
267	message_string = 'preprocessor switch for fftw is missing'
268	CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
269	#endif
270
271	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
272
273	CONTINUE
274
275	ELSE
276
277	message_string = 'fft method "' // TRIM( fft_method) // &
278	'" not available'
279	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
280	ENDIF
281
282	END SUBROUTINE fft_init
283
284
285	SUBROUTINE fft_x( ar, direction )
286
287	!----------------------------------------------------------------------!
288	! fft_x !
289	! !
290	! Fourier-transformation along x-direction !
291	! Version for 2D-decomposition !
292	! !
293	! fft_x uses internal algorithms (Singleton or Temperton) or !
294	! system-specific routines, if they are available !
295	!----------------------------------------------------------------------!
296
297	USE cuda_fft_interfaces
298	#if defined( __cuda_fft )
299	USE ISO_C_BINDING
300	#endif
301
302	IMPLICIT NONE
303
304	CHARACTER (LEN=*) :: direction
305	INTEGER :: i, ishape(1), j, k
306
307	LOGICAL :: forward_fft
308
309	REAL, DIMENSION(0:nx+2) :: work
310	REAL, DIMENSION(nx+2) :: work1
311	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
312	#if defined( __ibm )
313	REAL, DIMENSION(nau2) :: aux2, aux4
314	#elif defined( __nec )
315	REAL, DIMENSION(6*(nx+1)) :: work2
316	#elif defined( __cuda_fft )
317	!$acc declare create( ar_tmp )
318	COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp
319	#endif
320	REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar
321
322	IF ( direction == 'forward' ) THEN
323	forward_fft = .TRUE.
324	ELSE
325	forward_fft = .FALSE.
326	ENDIF
327
328	IF ( fft_method == 'singleton-algorithm' ) THEN
329
330	!
331	!-- Performing the fft with singleton's software works on every system,
332	!-- since it is part of the model
333	ALLOCATE( cwork(0:nx) )
334
335	IF ( forward_fft ) then
336
337	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
338	!$OMP DO
339	DO k = nzb_x, nzt_x
340	DO j = nys_x, nyn_x
341
342	DO i = 0, nx
343	cwork(i) = CMPLX( ar(i,j,k) )
344	ENDDO
345
346	ishape = SHAPE( cwork )
347	CALL FFTN( cwork, ishape )
348
349	DO i = 0, (nx+1)/2
350	ar(i,j,k) = REAL( cwork(i) )
351	ENDDO
352	DO i = 1, (nx+1)/2 - 1
353	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
354	ENDDO
355
356	ENDDO
357	ENDDO
358	!$OMP END PARALLEL
359
360	ELSE
361
362	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
363	!$OMP DO
364	DO k = nzb_x, nzt_x
365	DO j = nys_x, nyn_x
366
367	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
368	DO i = 1, (nx+1)/2 - 1
369	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
370	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
371	ENDDO
372	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
373
374	ishape = SHAPE( cwork )
375	CALL FFTN( cwork, ishape, inv = .TRUE. )
376
377	DO i = 0, nx
378	ar(i,j,k) = REAL( cwork(i) )
379	ENDDO
380
381	ENDDO
382	ENDDO
383	!$OMP END PARALLEL
384
385	ENDIF
386
387	DEALLOCATE( cwork )
388
389	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
390
391	!
392	!-- Performing the fft with Temperton's software works on every system,
393	!-- since it is part of the model
394	IF ( forward_fft ) THEN
395
396	!$OMP PARALLEL PRIVATE ( work, i, j, k )
397	!$OMP DO
398	DO k = nzb_x, nzt_x
399	DO j = nys_x, nyn_x
400
401	work(0:nx) = ar(0:nx,j,k)
402	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
403
404	DO i = 0, (nx+1)/2
405	ar(i,j,k) = work(2*i)
406	ENDDO
407	DO i = 1, (nx+1)/2 - 1
408	ar(nx+1-i,j,k) = work(2*i+1)
409	ENDDO
410
411	ENDDO
412	ENDDO
413	!$OMP END PARALLEL
414
415	ELSE
416
417	!$OMP PARALLEL PRIVATE ( work, i, j, k )
418	!$OMP DO
419	DO k = nzb_x, nzt_x
420	DO j = nys_x, nyn_x
421
422	DO i = 0, (nx+1)/2
423	work(2*i) = ar(i,j,k)
424	ENDDO
425	DO i = 1, (nx+1)/2 - 1
426	work(2*i+1) = ar(nx+1-i,j,k)
427	ENDDO
428	work(1) = 0.0
429	work(nx+2) = 0.0
430
431	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
432	ar(0:nx,j,k) = work(0:nx)
433
434	ENDDO
435	ENDDO
436	!$OMP END PARALLEL
437
438	ENDIF
439
440	ELSEIF ( fft_method == 'fftw' ) THEN
441
442	#if defined( __fftw )
443	IF ( forward_fft ) THEN
444
445	!$OMP PARALLEL PRIVATE ( work, i, j, k )
446	!$OMP DO
447	DO k = nzb_x, nzt_x
448	DO j = nys_x, nyn_x
449
450	x_in(0:nx) = ar(0:nx,j,k)
451	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
452
453	DO i = 0, (nx+1)/2
454	ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
455	ENDDO
456	DO i = 1, (nx+1)/2 - 1
457	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
458	ENDDO
459
460	ENDDO
461	ENDDO
462	!$OMP END PARALLEL
463
464	ELSE
465	!$OMP PARALLEL PRIVATE ( work, i, j, k )
466	!$OMP DO
467	DO k = nzb_x, nzt_x
468	DO j = nys_x, nyn_x
469
470	x_out(0) = CMPLX( ar(0,j,k), 0.0 )
471	DO i = 1, (nx+1)/2 - 1
472	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
473	ENDDO
474	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
475
476	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
477	ar(0:nx,j,k) = x_in(0:nx)
478
479	ENDDO
480	ENDDO
481	!$OMP END PARALLEL
482
483	ENDIF
484	#endif
485
486	ELSEIF ( fft_method == 'system-specific' ) THEN
487
488	#if defined( __ibm ) && ! defined( __ibmy_special )
489	IF ( forward_fft ) THEN
490
491	!$OMP PARALLEL PRIVATE ( work, i, j, k )
492	!$OMP DO
493	DO k = nzb_x, nzt_x
494	DO j = nys_x, nyn_x
495
496	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
497	aux2, nau2 )
498
499	DO i = 0, (nx+1)/2
500	ar(i,j,k) = work(2*i)
501	ENDDO
502	DO i = 1, (nx+1)/2 - 1
503	ar(nx+1-i,j,k) = work(2*i+1)
504	ENDDO
505
506	ENDDO
507	ENDDO
508	!$OMP END PARALLEL
509
510	ELSE
511
512	!$OMP PARALLEL PRIVATE ( work, i, j, k )
513	!$OMP DO
514	DO k = nzb_x, nzt_x
515	DO j = nys_x, nyn_x
516
517	DO i = 0, (nx+1)/2
518	work(2*i) = ar(i,j,k)
519	ENDDO
520	DO i = 1, (nx+1)/2 - 1
521	work(2*i+1) = ar(nx+1-i,j,k)
522	ENDDO
523	work(1) = 0.0
524	work(nx+2) = 0.0
525
526	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
527	aux4, nau2 )
528
529	DO i = 0, nx
530	ar(i,j,k) = work(i)
531	ENDDO
532
533	ENDDO
534	ENDDO
535	!$OMP END PARALLEL
536
537	ENDIF
538
539	#elif defined( __nec )
540
541	IF ( forward_fft ) THEN
542
543	!$OMP PARALLEL PRIVATE ( work, i, j, k )
544	!$OMP DO
545	DO k = nzb_x, nzt_x
546	DO j = nys_x, nyn_x
547
548	work(0:nx) = ar(0:nx,j,k)
549
550	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
551
552	DO i = 0, (nx+1)/2
553	ar(i,j,k) = work(2*i)
554	ENDDO
555	DO i = 1, (nx+1)/2 - 1
556	ar(nx+1-i,j,k) = work(2*i+1)
557	ENDDO
558
559	ENDDO
560	ENDDO
561	!$END OMP PARALLEL
562
563	ELSE
564
565	!$OMP PARALLEL PRIVATE ( work, i, j, k )
566	!$OMP DO
567	DO k = nzb_x, nzt_x
568	DO j = nys_x, nyn_x
569
570	DO i = 0, (nx+1)/2
571	work(2*i) = ar(i,j,k)
572	ENDDO
573	DO i = 1, (nx+1)/2 - 1
574	work(2*i+1) = ar(nx+1-i,j,k)
575	ENDDO
576	work(1) = 0.0
577	work(nx+2) = 0.0
578
579	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
580
581	ar(0:nx,j,k) = work(0:nx)
582
583	ENDDO
584	ENDDO
585	!$OMP END PARALLEL
586
587	ENDIF
588
589	#elif defined( __cuda_fft )
590
591	IF ( forward_fft ) THEN
592
593	!$acc data present( ar )
594	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
595
596	!$acc kernels
597	!$acc loop
598	DO k = nzb_x, nzt_x
599	DO j = nys_x, nyn_x
600
601	!$acc loop vector( 32 )
602	DO i = 0, (nx+1)/2
603	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
604	ENDDO
605
606	!$acc loop vector( 32 )
607	DO i = 1, (nx+1)/2 - 1
608	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
609	ENDDO
610
611	ENDDO
612	ENDDO
613	!$acc end kernels
614	!$acc end data
615
616	ELSE
617
618	!$acc data present( ar )
619	!$acc kernels
620	!$acc loop
621	DO k = nzb_x, nzt_x
622	DO j = nys_x, nyn_x
623
624	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
625
626	!$acc loop vector( 32 )
627	DO i = 1, (nx+1)/2 - 1
628	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
629	ENDDO
630	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
631
632	ENDDO
633	ENDDO
634	!$acc end kernels
635
636	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
637	!$acc end data
638
639	ENDIF
640
641	#else
642	message_string = 'no system-specific fft-call available'
643	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
644	#endif
645
646	ELSE
647
648	message_string = 'fft method "' // TRIM( fft_method) // &
649	'" not available'
650	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
651
652	ENDIF
653
654	END SUBROUTINE fft_x
655
656	SUBROUTINE fft_x_1d( ar, direction )
657
658	!----------------------------------------------------------------------!
659	! fft_x_1d !
660	! !
661	! Fourier-transformation along x-direction !
662	! Version for 1D-decomposition !
663	! !
664	! fft_x uses internal algorithms (Singleton or Temperton) or !
665	! system-specific routines, if they are available !
666	!----------------------------------------------------------------------!
667
668	IMPLICIT NONE
669
670	CHARACTER (LEN=*) :: direction
671	INTEGER :: i, ishape(1)
672
673	LOGICAL :: forward_fft
674
675	REAL, DIMENSION(0:nx) :: ar
676	REAL, DIMENSION(0:nx+2) :: work
677	REAL, DIMENSION(nx+2) :: work1
678	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
679	#if defined( __ibm )
680	REAL, DIMENSION(nau2) :: aux2, aux4
681	#elif defined( __nec )
682	REAL, DIMENSION(6*(nx+1)) :: work2
683	#endif
684
685	IF ( direction == 'forward' ) THEN
686	forward_fft = .TRUE.
687	ELSE
688	forward_fft = .FALSE.
689	ENDIF
690
691	IF ( fft_method == 'singleton-algorithm' ) THEN
692
693	!
694	!-- Performing the fft with singleton's software works on every system,
695	!-- since it is part of the model
696	ALLOCATE( cwork(0:nx) )
697
698	IF ( forward_fft ) then
699
700	DO i = 0, nx
701	cwork(i) = CMPLX( ar(i) )
702	ENDDO
703	ishape = SHAPE( cwork )
704	CALL FFTN( cwork, ishape )
705	DO i = 0, (nx+1)/2
706	ar(i) = REAL( cwork(i) )
707	ENDDO
708	DO i = 1, (nx+1)/2 - 1
709	ar(nx+1-i) = -AIMAG( cwork(i) )
710	ENDDO
711
712	ELSE
713
714	cwork(0) = CMPLX( ar(0), 0.0 )
715	DO i = 1, (nx+1)/2 - 1
716	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
717	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
718	ENDDO
719	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
720
721	ishape = SHAPE( cwork )
722	CALL FFTN( cwork, ishape, inv = .TRUE. )
723
724	DO i = 0, nx
725	ar(i) = REAL( cwork(i) )
726	ENDDO
727
728	ENDIF
729
730	DEALLOCATE( cwork )
731
732	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
733
734	!
735	!-- Performing the fft with Temperton's software works on every system,
736	!-- since it is part of the model
737	IF ( forward_fft ) THEN
738
739	work(0:nx) = ar
740	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
741
742	DO i = 0, (nx+1)/2
743	ar(i) = work(2*i)
744	ENDDO
745	DO i = 1, (nx+1)/2 - 1
746	ar(nx+1-i) = work(2*i+1)
747	ENDDO
748
749	ELSE
750
751	DO i = 0, (nx+1)/2
752	work(2*i) = ar(i)
753	ENDDO
754	DO i = 1, (nx+1)/2 - 1
755	work(2*i+1) = ar(nx+1-i)
756	ENDDO
757	work(1) = 0.0
758	work(nx+2) = 0.0
759
760	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
761	ar = work(0:nx)
762
763	ENDIF
764
765	ELSEIF ( fft_method == 'system-specific' ) THEN
766
767	#if defined( __ibm ) && ! defined( __ibmy_special )
768	IF ( forward_fft ) THEN
769
770	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
771	aux2, nau2 )
772
773	DO i = 0, (nx+1)/2
774	ar(i) = work(2*i)
775	ENDDO
776	DO i = 1, (nx+1)/2 - 1
777	ar(nx+1-i) = work(2*i+1)
778	ENDDO
779
780	ELSE
781
782	DO i = 0, (nx+1)/2
783	work(2*i) = ar(i)
784	ENDDO
785	DO i = 1, (nx+1)/2 - 1
786	work(2*i+1) = ar(nx+1-i)
787	ENDDO
788	work(1) = 0.0
789	work(nx+2) = 0.0
790
791	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
792	aux4, nau2 )
793
794	DO i = 0, nx
795	ar(i) = work(i)
796	ENDDO
797
798	ENDIF
799	#elif defined( __nec )
800	IF ( forward_fft ) THEN
801
802	work(0:nx) = ar(0:nx)
803
804	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
805
806	DO i = 0, (nx+1)/2
807	ar(i) = work(2*i)
808	ENDDO
809	DO i = 1, (nx+1)/2 - 1
810	ar(nx+1-i) = work(2*i+1)
811	ENDDO
812
813	ELSE
814
815	DO i = 0, (nx+1)/2
816	work(2*i) = ar(i)
817	ENDDO
818	DO i = 1, (nx+1)/2 - 1
819	work(2*i+1) = ar(nx+1-i)
820	ENDDO
821	work(1) = 0.0
822	work(nx+2) = 0.0
823
824	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
825
826	ar(0:nx) = work(0:nx)
827
828	ENDIF
829	#else
830	message_string = 'no system-specific fft-call available'
831	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
832	#endif
833	ELSE
834	message_string = 'fft method "' // TRIM( fft_method) // &
835	'" not available'
836	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
837
838	ENDIF
839
840	END SUBROUTINE fft_x_1d
841
842	SUBROUTINE fft_y( ar, direction )
843
844	!----------------------------------------------------------------------!
845	! fft_y !
846	! !
847	! Fourier-transformation along y-direction !
848	! Version for 2D-decomposition !
849	! !
850	! fft_y uses internal algorithms (Singleton or Temperton) or !
851	! system-specific routines, if they are available !
852	!----------------------------------------------------------------------!
853
854	USE cuda_fft_interfaces
855	#if defined( __cuda_fft )
856	USE ISO_C_BINDING
857	#endif
858
859	IMPLICIT NONE
860
861	CHARACTER (LEN=*) :: direction
862	INTEGER :: i, j, jshape(1), k
863
864	LOGICAL :: forward_fft
865
866	REAL, DIMENSION(0:ny+2) :: work
867	REAL, DIMENSION(ny+2) :: work1
868	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
869	#if defined( __ibm )
870	REAL, DIMENSION(nau2) :: auy2, auy4
871	#elif defined( __nec )
872	REAL, DIMENSION(6*(ny+1)) :: work2
873	#elif defined( __cuda_fft )
874	!$acc declare create( ar_tmp )
875	COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
876	#endif
877	REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) :: ar
878
879	IF ( direction == 'forward' ) THEN
880	forward_fft = .TRUE.
881	ELSE
882	forward_fft = .FALSE.
883	ENDIF
884
885	IF ( fft_method == 'singleton-algorithm' ) THEN
886
887	!
888	!-- Performing the fft with singleton's software works on every system,
889	!-- since it is part of the model
890	ALLOCATE( cwork(0:ny) )
891
892	IF ( forward_fft ) then
893
894	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
895	!$OMP DO
896	DO k = nzb_y, nzt_y
897	DO i = nxl_y, nxr_y
898
899	DO j = 0, ny
900	cwork(j) = CMPLX( ar(j,i,k) )
901	ENDDO
902
903	jshape = SHAPE( cwork )
904	CALL FFTN( cwork, jshape )
905
906	DO j = 0, (ny+1)/2
907	ar(j,i,k) = REAL( cwork(j) )
908	ENDDO
909	DO j = 1, (ny+1)/2 - 1
910	ar(ny+1-j,i,k) = -AIMAG( cwork(j) )
911	ENDDO
912
913	ENDDO
914	ENDDO
915	!$OMP END PARALLEL
916
917	ELSE
918
919	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
920	!$OMP DO
921	DO k = nzb_y, nzt_y
922	DO i = nxl_y, nxr_y
923
924	cwork(0) = CMPLX( ar(0,i,k), 0.0 )
925	DO j = 1, (ny+1)/2 - 1
926	cwork(j) = CMPLX( ar(j,i,k), -ar(ny+1-j,i,k) )
927	cwork(ny+1-j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
928	ENDDO
929	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
930
931	jshape = SHAPE( cwork )
932	CALL FFTN( cwork, jshape, inv = .TRUE. )
933
934	DO j = 0, ny
935	ar(j,i,k) = REAL( cwork(j) )
936	ENDDO
937
938	ENDDO
939	ENDDO
940	!$OMP END PARALLEL
941
942	ENDIF
943
944	DEALLOCATE( cwork )
945
946	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
947
948	!
949	!-- Performing the fft with Temperton's software works on every system,
950	!-- since it is part of the model
951	IF ( forward_fft ) THEN
952
953	!$OMP PARALLEL PRIVATE ( work, i, j, k )
954	!$OMP DO
955	DO k = nzb_y, nzt_y
956	DO i = nxl_y, nxr_y
957
958	work(0:ny) = ar(0:ny,i,k)
959	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
960
961	DO j = 0, (ny+1)/2
962	ar(j,i,k) = work(2*j)
963	ENDDO
964	DO j = 1, (ny+1)/2 - 1
965	ar(ny+1-j,i,k) = work(2*j+1)
966	ENDDO
967
968	ENDDO
969	ENDDO
970	!$OMP END PARALLEL
971
972	ELSE
973
974	!$OMP PARALLEL PRIVATE ( work, i, j, k )
975	!$OMP DO
976	DO k = nzb_y, nzt_y
977	DO i = nxl_y, nxr_y
978
979	DO j = 0, (ny+1)/2
980	work(2*j) = ar(j,i,k)
981	ENDDO
982	DO j = 1, (ny+1)/2 - 1
983	work(2*j+1) = ar(ny+1-j,i,k)
984	ENDDO
985	work(1) = 0.0
986	work(ny+2) = 0.0
987
988	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
989	ar(0:ny,i,k) = work(0:ny)
990
991	ENDDO
992	ENDDO
993	!$OMP END PARALLEL
994
995	ENDIF
996
997	ELSEIF ( fft_method == 'fftw' ) THEN
998
999	#if defined( __fftw )
1000	IF ( forward_fft ) THEN
1001
1002	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1003	!$OMP DO
1004	DO k = nzb_y, nzt_y
1005	DO i = nxl_y, nxr_y
1006
1007	y_in(0:ny) = ar(0:ny,i,k)
1008	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1009
1010	DO j = 0, (ny+1)/2
1011	ar(j,i,k) = REAL( y_out(j) ) /(ny+1)
1012	ENDDO
1013	DO j = 1, (ny+1)/2 - 1
1014	ar(ny+1-j,i,k) = AIMAG( y_out(j) ) /(ny+1)
1015	ENDDO
1016
1017	ENDDO
1018	ENDDO
1019	!$OMP END PARALLEL
1020
1021	ELSE
1022
1023	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1024	!$OMP DO
1025	DO k = nzb_y, nzt_y
1026	DO i = nxl_y, nxr_y
1027
1028	y_out(0) = CMPLX( ar(0,i,k), 0.0 )
1029	DO j = 1, (ny+1)/2 - 1
1030	y_out(j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1031	ENDDO
1032	y_out((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1033
1034	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1035	ar(0:ny,i,k) = y_in(0:ny)
1036
1037	ENDDO
1038	ENDDO
1039	!$OMP END PARALLEL
1040
1041	ENDIF
1042	#endif
1043
1044	ELSEIF ( fft_method == 'system-specific' ) THEN
1045
1046	#if defined( __ibm ) && ! defined( __ibmy_special )
1047	IF ( forward_fft) THEN
1048
1049	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1050	!$OMP DO
1051	DO k = nzb_y, nzt_y
1052	DO i = nxl_y, nxr_y
1053
1054	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1055	auy2, nau2 )
1056
1057	DO j = 0, (ny+1)/2
1058	ar(j,i,k) = work(2*j)
1059	ENDDO
1060	DO j = 1, (ny+1)/2 - 1
1061	ar(ny+1-j,i,k) = work(2*j+1)
1062	ENDDO
1063
1064	ENDDO
1065	ENDDO
1066	!$OMP END PARALLEL
1067
1068	ELSE
1069
1070	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1071	!$OMP DO
1072	DO k = nzb_y, nzt_y
1073	DO i = nxl_y, nxr_y
1074
1075	DO j = 0, (ny+1)/2
1076	work(2*j) = ar(j,i,k)
1077	ENDDO
1078	DO j = 1, (ny+1)/2 - 1
1079	work(2*j+1) = ar(ny+1-j,i,k)
1080	ENDDO
1081	work(1) = 0.0
1082	work(ny+2) = 0.0
1083
1084	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1085	auy4, nau2 )
1086
1087	DO j = 0, ny
1088	ar(j,i,k) = work(j)
1089	ENDDO
1090
1091	ENDDO
1092	ENDDO
1093	!$OMP END PARALLEL
1094
1095	ENDIF
1096	#elif defined( __nec )
1097	IF ( forward_fft ) THEN
1098
1099	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1100	!$OMP DO
1101	DO k = nzb_y, nzt_y
1102	DO i = nxl_y, nxr_y
1103
1104	work(0:ny) = ar(0:ny,i,k)
1105
1106	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1107
1108	DO j = 0, (ny+1)/2
1109	ar(j,i,k) = work(2*j)
1110	ENDDO
1111	DO j = 1, (ny+1)/2 - 1
1112	ar(ny+1-j,i,k) = work(2*j+1)
1113	ENDDO
1114
1115	ENDDO
1116	ENDDO
1117	!$END OMP PARALLEL
1118
1119	ELSE
1120
1121	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1122	!$OMP DO
1123	DO k = nzb_y, nzt_y
1124	DO i = nxl_y, nxr_y
1125
1126	DO j = 0, (ny+1)/2
1127	work(2*j) = ar(j,i,k)
1128	ENDDO
1129	DO j = 1, (ny+1)/2 - 1
1130	work(2*j+1) = ar(ny+1-j,i,k)
1131	ENDDO
1132	work(1) = 0.0
1133	work(ny+2) = 0.0
1134
1135	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1136
1137	ar(0:ny,i,k) = work(0:ny)
1138
1139	ENDDO
1140	ENDDO
1141	!$OMP END PARALLEL
1142
1143	ENDIF
1144	#elif defined( __cuda_fft )
1145
1146	IF ( forward_fft ) THEN
1147
1148	!$acc data present( ar )
1149	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1150
1151	!$acc kernels
1152	!$acc loop
1153	DO k = nzb_y, nzt_y
1154	DO i = nxl_y, nxr_y
1155
1156	!$acc loop vector( 32 )
1157	DO j = 0, (ny+1)/2
1158	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
1159	ENDDO
1160
1161	!$acc loop vector( 32 )
1162	DO j = 1, (ny+1)/2 - 1
1163	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1164	ENDDO
1165
1166	ENDDO
1167	ENDDO
1168	!$acc end kernels
1169	!$acc end data
1170
1171	ELSE
1172
1173	!$acc data present( ar )
1174	!$acc kernels
1175	!$acc loop
1176	DO k = nzb_y, nzt_y
1177	DO i = nxl_y, nxr_y
1178
1179	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1180
1181	!$acc loop vector( 32 )
1182	DO j = 1, (ny+1)/2 - 1
1183	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1184	ENDDO
1185	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1186
1187	ENDDO
1188	ENDDO
1189	!$acc end kernels
1190
1191	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1192	!$acc end data
1193
1194	ENDIF
1195
1196	#else
1197	message_string = 'no system-specific fft-call available'
1198	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
1199	#endif
1200
1201	ELSE
1202
1203	message_string = 'fft method "' // TRIM( fft_method) // &
1204	'" not available'
1205	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
1206
1207	ENDIF
1208
1209	END SUBROUTINE fft_y
1210
1211	SUBROUTINE fft_y_1d( ar, direction )
1212
1213	!----------------------------------------------------------------------!
1214	! fft_y_1d !
1215	! !
1216	! Fourier-transformation along y-direction !
1217	! Version for 1D-decomposition !
1218	! !
1219	! fft_y uses internal algorithms (Singleton or Temperton) or !
1220	! system-specific routines, if they are available !
1221	!----------------------------------------------------------------------!
1222
1223	IMPLICIT NONE
1224
1225	CHARACTER (LEN=*) :: direction
1226	INTEGER :: j, jshape(1)
1227
1228	LOGICAL :: forward_fft
1229
1230	REAL, DIMENSION(0:ny) :: ar
1231	REAL, DIMENSION(0:ny+2) :: work
1232	REAL, DIMENSION(ny+2) :: work1
1233	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
1234	#if defined( __ibm )
1235	REAL, DIMENSION(nau2) :: auy2, auy4
1236	#elif defined( __nec )
1237	REAL, DIMENSION(6*(ny+1)) :: work2
1238	#endif
1239
1240	IF ( direction == 'forward' ) THEN
1241	forward_fft = .TRUE.
1242	ELSE
1243	forward_fft = .FALSE.
1244	ENDIF
1245
1246	IF ( fft_method == 'singleton-algorithm' ) THEN
1247
1248	!
1249	!-- Performing the fft with singleton's software works on every system,
1250	!-- since it is part of the model
1251	ALLOCATE( cwork(0:ny) )
1252
1253	IF ( forward_fft ) THEN
1254
1255	DO j = 0, ny
1256	cwork(j) = CMPLX( ar(j) )
1257	ENDDO
1258
1259	jshape = SHAPE( cwork )
1260	CALL FFTN( cwork, jshape )
1261
1262	DO j = 0, (ny+1)/2
1263	ar(j) = REAL( cwork(j) )
1264	ENDDO
1265	DO j = 1, (ny+1)/2 - 1
1266	ar(ny+1-j) = -AIMAG( cwork(j) )
1267	ENDDO
1268
1269	ELSE
1270
1271	cwork(0) = CMPLX( ar(0), 0.0 )
1272	DO j = 1, (ny+1)/2 - 1
1273	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
1274	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
1275	ENDDO
1276	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
1277
1278	jshape = SHAPE( cwork )
1279	CALL FFTN( cwork, jshape, inv = .TRUE. )
1280
1281	DO j = 0, ny
1282	ar(j) = REAL( cwork(j) )
1283	ENDDO
1284
1285	ENDIF
1286
1287	DEALLOCATE( cwork )
1288
1289	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1290
1291	!
1292	!-- Performing the fft with Temperton's software works on every system,
1293	!-- since it is part of the model
1294	IF ( forward_fft ) THEN
1295
1296	work(0:ny) = ar
1297	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1298
1299	DO j = 0, (ny+1)/2
1300	ar(j) = work(2*j)
1301	ENDDO
1302	DO j = 1, (ny+1)/2 - 1
1303	ar(ny+1-j) = work(2*j+1)
1304	ENDDO
1305
1306	ELSE
1307
1308	DO j = 0, (ny+1)/2
1309	work(2*j) = ar(j)
1310	ENDDO
1311	DO j = 1, (ny+1)/2 - 1
1312	work(2*j+1) = ar(ny+1-j)
1313	ENDDO
1314	work(1) = 0.0
1315	work(ny+2) = 0.0
1316
1317	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1318	ar = work(0:ny)
1319
1320	ENDIF
1321
1322	ELSEIF ( fft_method == 'system-specific' ) THEN
1323
1324	#if defined( __ibm ) && ! defined( __ibmy_special )
1325	IF ( forward_fft ) THEN
1326
1327	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1328	auy2, nau2 )
1329
1330	DO j = 0, (ny+1)/2
1331	ar(j) = work(2*j)
1332	ENDDO
1333	DO j = 1, (ny+1)/2 - 1
1334	ar(ny+1-j) = work(2*j+1)
1335	ENDDO
1336
1337	ELSE
1338
1339	DO j = 0, (ny+1)/2
1340	work(2*j) = ar(j)
1341	ENDDO
1342	DO j = 1, (ny+1)/2 - 1
1343	work(2*j+1) = ar(ny+1-j)
1344	ENDDO
1345	work(1) = 0.0
1346	work(ny+2) = 0.0
1347
1348	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1349	auy4, nau2 )
1350
1351	DO j = 0, ny
1352	ar(j) = work(j)
1353	ENDDO
1354
1355	ENDIF
1356	#elif defined( __nec )
1357	IF ( forward_fft ) THEN
1358
1359	work(0:ny) = ar(0:ny)
1360
1361	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1362
1363	DO j = 0, (ny+1)/2
1364	ar(j) = work(2*j)
1365	ENDDO
1366	DO j = 1, (ny+1)/2 - 1
1367	ar(ny+1-j) = work(2*j+1)
1368	ENDDO
1369
1370	ELSE
1371
1372	DO j = 0, (ny+1)/2
1373	work(2*j) = ar(j)
1374	ENDDO
1375	DO j = 1, (ny+1)/2 - 1
1376	work(2*j+1) = ar(ny+1-j)
1377	ENDDO
1378	work(1) = 0.0
1379	work(ny+2) = 0.0
1380
1381	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1382
1383	ar(0:ny) = work(0:ny)
1384
1385	ENDIF
1386	#else
1387	message_string = 'no system-specific fft-call available'
1388	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
1389
1390	#endif
1391
1392	ELSE
1393
1394	message_string = 'fft method "' // TRIM( fft_method) // &
1395	'" not available'
1396	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1397
1398	ENDIF
1399
1400	END SUBROUTINE fft_y_1d
1401
1402	SUBROUTINE fft_x_m( ar, direction )
1403
1404	!----------------------------------------------------------------------!
1405	! fft_x_m !
1406	! !
1407	! Fourier-transformation along x-direction !
1408	! Version for 1d domain decomposition !
1409	! using multiple 1D FFT from Math Keisan on NEC !
1410	! or Temperton-algorithm !
1411	! (no singleton-algorithm on NEC because it does not vectorize) !
1412	! !
1413	!----------------------------------------------------------------------!
1414
1415	IMPLICIT NONE
1416
1417	CHARACTER (LEN=*) :: direction
1418	INTEGER :: i, k, siza
1419
1420	REAL, DIMENSION(0:nx,nz) :: ar
1421	REAL, DIMENSION(0:nx+3,nz+1) :: ai
1422	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
1423	#if defined( __nec )
1424	INTEGER :: sizw
1425	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
1426	#endif
1427
1428	IF ( fft_method == 'temperton-algorithm' ) THEN
1429
1430	siza = SIZE( ai, 1 )
1431
1432	IF ( direction == 'forward') THEN
1433
1434	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1435	ai(nx+1:,:) = 0.0
1436
1437	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1438
1439	DO k = 1, nz
1440	DO i = 0, (nx+1)/2
1441	ar(i,k) = ai(2*i,k)
1442	ENDDO
1443	DO i = 1, (nx+1)/2 - 1
1444	ar(nx+1-i,k) = ai(2*i+1,k)
1445	ENDDO
1446	ENDDO
1447
1448	ELSE
1449
1450	DO k = 1, nz
1451	DO i = 0, (nx+1)/2
1452	ai(2*i,k) = ar(i,k)
1453	ENDDO
1454	DO i = 1, (nx+1)/2 - 1
1455	ai(2*i+1,k) = ar(nx+1-i,k)
1456	ENDDO
1457	ai(1,k) = 0.0
1458	ai(nx+2,k) = 0.0
1459	ENDDO
1460
1461	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1462
1463	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1464
1465	ENDIF
1466
1467	ELSEIF ( fft_method == 'system-specific' ) THEN
1468
1469	#if defined( __nec )
1470	siza = SIZE( ai, 1 )
1471	sizw = SIZE( work, 1 )
1472
1473	IF ( direction == 'forward') THEN
1474
1475	!
1476	!-- Tables are initialized once more. This call should not be
1477	!-- necessary, but otherwise program aborts in asymmetric case
1478	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1479	trig_xf, work1, 0 )
1480
1481	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1482	IF ( nz1 > nz ) THEN
1483	ai(:,nz1) = 0.0
1484	ENDIF
1485
1486	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
1487	trig_xf, work1, 0 )
1488
1489	DO k = 1, nz
1490	DO i = 0, (nx+1)/2
1491	ar(i,k) = REAL( work(i+1,k) )
1492	ENDDO
1493	DO i = 1, (nx+1)/2 - 1
1494	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1495	ENDDO
1496	ENDDO
1497
1498	ELSE
1499
1500	!
1501	!-- Tables are initialized once more. This call should not be
1502	!-- necessary, but otherwise program aborts in asymmetric case
1503	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1504	trig_xb, work1, 0 )
1505
1506	IF ( nz1 > nz ) THEN
1507	work(:,nz1) = 0.0
1508	ENDIF
1509	DO k = 1, nz
1510	work(1,k) = CMPLX( ar(0,k), 0.0 )
1511	DO i = 1, (nx+1)/2 - 1
1512	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
1513	ENDDO
1514	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
1515	ENDDO
1516
1517	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1518	trig_xb, work1, 0 )
1519
1520	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1521
1522	ENDIF
1523
1524	#else
1525	message_string = 'no system-specific fft-call available'
1526	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
1527	#endif
1528
1529	ELSE
1530
1531	message_string = 'fft method "' // TRIM( fft_method) // &
1532	'" not available'
1533	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1534
1535	ENDIF
1536
1537	END SUBROUTINE fft_x_m
1538
1539	SUBROUTINE fft_y_m( ar, ny1, direction )
1540
1541	!----------------------------------------------------------------------!
1542	! fft_y_m !
1543	! !
1544	! Fourier-transformation along y-direction !
1545	! Version for 1d domain decomposition !
1546	! using multiple 1D FFT from Math Keisan on NEC !
1547	! or Temperton-algorithm !
1548	! (no singleton-algorithm on NEC because it does not vectorize) !
1549	! !
1550	!----------------------------------------------------------------------!
1551
1552	IMPLICIT NONE
1553
1554	CHARACTER (LEN=*) :: direction
1555	INTEGER :: j, k, ny1, siza
1556
1557	REAL, DIMENSION(0:ny1,nz) :: ar
1558	REAL, DIMENSION(0:ny+3,nz+1) :: ai
1559	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
1560	#if defined( __nec )
1561	INTEGER :: sizw
1562	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
1563	#endif
1564
1565	IF ( fft_method == 'temperton-algorithm' ) THEN
1566
1567	siza = SIZE( ai, 1 )
1568
1569	IF ( direction == 'forward') THEN
1570
1571	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1572	ai(ny+1:,:) = 0.0
1573
1574	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1575
1576	DO k = 1, nz
1577	DO j = 0, (ny+1)/2
1578	ar(j,k) = ai(2*j,k)
1579	ENDDO
1580	DO j = 1, (ny+1)/2 - 1
1581	ar(ny+1-j,k) = ai(2*j+1,k)
1582	ENDDO
1583	ENDDO
1584
1585	ELSE
1586
1587	DO k = 1, nz
1588	DO j = 0, (ny+1)/2
1589	ai(2*j,k) = ar(j,k)
1590	ENDDO
1591	DO j = 1, (ny+1)/2 - 1
1592	ai(2*j+1,k) = ar(ny+1-j,k)
1593	ENDDO
1594	ai(1,k) = 0.0
1595	ai(ny+2,k) = 0.0
1596	ENDDO
1597
1598	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1599
1600	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1601
1602	ENDIF
1603
1604	ELSEIF ( fft_method == 'system-specific' ) THEN
1605
1606	#if defined( __nec )
1607	siza = SIZE( ai, 1 )
1608	sizw = SIZE( work, 1 )
1609
1610	IF ( direction == 'forward') THEN
1611
1612	!
1613	!-- Tables are initialized once more. This call should not be
1614	!-- necessary, but otherwise program aborts in asymmetric case
1615	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1616	trig_yf, work1, 0 )
1617
1618	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1619	IF ( nz1 > nz ) THEN
1620	ai(:,nz1) = 0.0
1621	ENDIF
1622
1623	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
1624	trig_yf, work1, 0 )
1625
1626	DO k = 1, nz
1627	DO j = 0, (ny+1)/2
1628	ar(j,k) = REAL( work(j+1,k) )
1629	ENDDO
1630	DO j = 1, (ny+1)/2 - 1
1631	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1632	ENDDO
1633	ENDDO
1634
1635	ELSE
1636
1637	!
1638	!-- Tables are initialized once more. This call should not be
1639	!-- necessary, but otherwise program aborts in asymmetric case
1640	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1641	trig_yb, work1, 0 )
1642
1643	IF ( nz1 > nz ) THEN
1644	work(:,nz1) = 0.0
1645	ENDIF
1646	DO k = 1, nz
1647	work(1,k) = CMPLX( ar(0,k), 0.0 )
1648	DO j = 1, (ny+1)/2 - 1
1649	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
1650	ENDDO
1651	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
1652	ENDDO
1653
1654	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1655	trig_yb, work1, 0 )
1656
1657	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1658
1659	ENDIF
1660
1661	#else
1662	message_string = 'no system-specific fft-call available'
1663	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
1664	#endif
1665
1666	ELSE
1667
1668	message_string = 'fft method "' // TRIM( fft_method) // &
1669	'" not available'
1670	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1671
1672	ENDIF
1673
1674	END SUBROUTINE fft_y_m
1675
1676
1677	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |