Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |