Home

Context Navigation

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

Last change on this file since 1213 was 1211, checked in by raasch, 11 years ago
last commit documented
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	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: fft_xy.f90 1211 2013-08-14 13:46:36Z 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 )
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,nzb_x:nzt_x) :: ar
324
325	IF ( direction == 'forward' ) THEN
326	forward_fft = .TRUE.
327	ELSE
328	forward_fft = .FALSE.
329	ENDIF
330
331	IF ( fft_method == 'singleton-algorithm' ) THEN
332
333	!
334	!-- Performing the fft with singleton's software works on every system,
335	!-- since it is part of the model
336	ALLOCATE( cwork(0:nx) )
337
338	IF ( forward_fft ) then
339
340	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
341	!$OMP DO
342	DO k = nzb_x, nzt_x
343	DO j = nys_x, nyn_x
344
345	DO i = 0, nx
346	cwork(i) = CMPLX( ar(i,j,k) )
347	ENDDO
348
349	ishape = SHAPE( cwork )
350	CALL FFTN( cwork, ishape )
351
352	DO i = 0, (nx+1)/2
353	ar(i,j,k) = REAL( cwork(i) )
354	ENDDO
355	DO i = 1, (nx+1)/2 - 1
356	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
357	ENDDO
358
359	ENDDO
360	ENDDO
361	!$OMP END PARALLEL
362
363	ELSE
364
365	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
366	!$OMP DO
367	DO k = nzb_x, nzt_x
368	DO j = nys_x, nyn_x
369
370	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
371	DO i = 1, (nx+1)/2 - 1
372	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
373	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
374	ENDDO
375	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
376
377	ishape = SHAPE( cwork )
378	CALL FFTN( cwork, ishape, inv = .TRUE. )
379
380	DO i = 0, nx
381	ar(i,j,k) = REAL( cwork(i) )
382	ENDDO
383
384	ENDDO
385	ENDDO
386	!$OMP END PARALLEL
387
388	ENDIF
389
390	DEALLOCATE( cwork )
391
392	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
393
394	!
395	!-- Performing the fft with Temperton's software works on every system,
396	!-- since it is part of the model
397	IF ( forward_fft ) THEN
398
399	!$OMP PARALLEL PRIVATE ( work, i, j, k )
400	!$OMP DO
401	DO k = nzb_x, nzt_x
402	DO j = nys_x, nyn_x
403
404	work(0:nx) = ar(0:nx,j,k)
405	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
406
407	DO i = 0, (nx+1)/2
408	ar(i,j,k) = work(2*i)
409	ENDDO
410	DO i = 1, (nx+1)/2 - 1
411	ar(nx+1-i,j,k) = work(2*i+1)
412	ENDDO
413
414	ENDDO
415	ENDDO
416	!$OMP END PARALLEL
417
418	ELSE
419
420	!$OMP PARALLEL PRIVATE ( work, i, j, k )
421	!$OMP DO
422	DO k = nzb_x, nzt_x
423	DO j = nys_x, nyn_x
424
425	DO i = 0, (nx+1)/2
426	work(2*i) = ar(i,j,k)
427	ENDDO
428	DO i = 1, (nx+1)/2 - 1
429	work(2*i+1) = ar(nx+1-i,j,k)
430	ENDDO
431	work(1) = 0.0
432	work(nx+2) = 0.0
433
434	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
435	ar(0:nx,j,k) = work(0:nx)
436
437	ENDDO
438	ENDDO
439	!$OMP END PARALLEL
440
441	ENDIF
442
443	ELSEIF ( fft_method == 'fftw' ) THEN
444
445	#if defined( __fftw )
446	IF ( forward_fft ) THEN
447
448	!$OMP PARALLEL PRIVATE ( work, i, j, k )
449	!$OMP DO
450	DO k = nzb_x, nzt_x
451	DO j = nys_x, nyn_x
452
453	x_in(0:nx) = ar(0:nx,j,k)
454	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
455
456	DO i = 0, (nx+1)/2
457	ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
458	ENDDO
459	DO i = 1, (nx+1)/2 - 1
460	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
461	ENDDO
462
463	ENDDO
464	ENDDO
465	!$OMP END PARALLEL
466
467	ELSE
468	!$OMP PARALLEL PRIVATE ( work, i, j, k )
469	!$OMP DO
470	DO k = nzb_x, nzt_x
471	DO j = nys_x, nyn_x
472
473	x_out(0) = CMPLX( ar(0,j,k), 0.0 )
474	DO i = 1, (nx+1)/2 - 1
475	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
476	ENDDO
477	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
478
479	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
480	ar(0:nx,j,k) = x_in(0:nx)
481
482	ENDDO
483	ENDDO
484	!$OMP END PARALLEL
485
486	ENDIF
487	#endif
488
489	ELSEIF ( fft_method == 'system-specific' ) THEN
490
491	#if defined( __ibm ) && ! defined( __ibmy_special )
492	IF ( forward_fft ) THEN
493
494	!$OMP PARALLEL PRIVATE ( work, i, j, k )
495	!$OMP DO
496	DO k = nzb_x, nzt_x
497	DO j = nys_x, nyn_x
498
499	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
500	aux2, nau2 )
501
502	DO i = 0, (nx+1)/2
503	ar(i,j,k) = work(2*i)
504	ENDDO
505	DO i = 1, (nx+1)/2 - 1
506	ar(nx+1-i,j,k) = work(2*i+1)
507	ENDDO
508
509	ENDDO
510	ENDDO
511	!$OMP END PARALLEL
512
513	ELSE
514
515	!$OMP PARALLEL PRIVATE ( work, i, j, k )
516	!$OMP DO
517	DO k = nzb_x, nzt_x
518	DO j = nys_x, nyn_x
519
520	DO i = 0, (nx+1)/2
521	work(2*i) = ar(i,j,k)
522	ENDDO
523	DO i = 1, (nx+1)/2 - 1
524	work(2*i+1) = ar(nx+1-i,j,k)
525	ENDDO
526	work(1) = 0.0
527	work(nx+2) = 0.0
528
529	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
530	aux4, nau2 )
531
532	DO i = 0, nx
533	ar(i,j,k) = work(i)
534	ENDDO
535
536	ENDDO
537	ENDDO
538	!$OMP END PARALLEL
539
540	ENDIF
541
542	#elif defined( __nec )
543
544	IF ( forward_fft ) THEN
545
546	!$OMP PARALLEL PRIVATE ( work, i, j, k )
547	!$OMP DO
548	DO k = nzb_x, nzt_x
549	DO j = nys_x, nyn_x
550
551	work(0:nx) = ar(0:nx,j,k)
552
553	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
554
555	DO i = 0, (nx+1)/2
556	ar(i,j,k) = work(2*i)
557	ENDDO
558	DO i = 1, (nx+1)/2 - 1
559	ar(nx+1-i,j,k) = work(2*i+1)
560	ENDDO
561
562	ENDDO
563	ENDDO
564	!$END OMP PARALLEL
565
566	ELSE
567
568	!$OMP PARALLEL PRIVATE ( work, i, j, k )
569	!$OMP DO
570	DO k = nzb_x, nzt_x
571	DO j = nys_x, nyn_x
572
573	DO i = 0, (nx+1)/2
574	work(2*i) = ar(i,j,k)
575	ENDDO
576	DO i = 1, (nx+1)/2 - 1
577	work(2*i+1) = ar(nx+1-i,j,k)
578	ENDDO
579	work(1) = 0.0
580	work(nx+2) = 0.0
581
582	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
583
584	ar(0:nx,j,k) = work(0:nx)
585
586	ENDDO
587	ENDDO
588	!$OMP END PARALLEL
589
590	ENDIF
591
592	#elif defined( __cuda_fft )
593
594	IF ( forward_fft ) THEN
595
596	!$acc data present( ar )
597	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
598
599	!$acc kernels
600	!$acc loop
601	DO k = nzb_x, nzt_x
602	DO j = nys_x, nyn_x
603
604	!$acc loop vector( 32 )
605	DO i = 0, (nx+1)/2
606	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
607	ENDDO
608
609	!$acc loop vector( 32 )
610	DO i = 1, (nx+1)/2 - 1
611	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
612	ENDDO
613
614	ENDDO
615	ENDDO
616	!$acc end kernels
617	!$acc end data
618
619	ELSE
620
621	!$acc data present( ar )
622	!$acc kernels
623	!$acc loop
624	DO k = nzb_x, nzt_x
625	DO j = nys_x, nyn_x
626
627	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
628
629	!$acc loop vector( 32 )
630	DO i = 1, (nx+1)/2 - 1
631	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
632	ENDDO
633	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
634
635	ENDDO
636	ENDDO
637	!$acc end kernels
638
639	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
640	!$acc end data
641
642	ENDIF
643
644	#else
645	message_string = 'no system-specific fft-call available'
646	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
647	#endif
648
649	ELSE
650
651	message_string = 'fft method "' // TRIM( fft_method) // &
652	'" not available'
653	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
654
655	ENDIF
656
657	END SUBROUTINE fft_x
658
659	SUBROUTINE fft_x_1d( ar, direction )
660
661	!----------------------------------------------------------------------!
662	! fft_x_1d !
663	! !
664	! Fourier-transformation along x-direction !
665	! Version for 1D-decomposition !
666	! !
667	! fft_x uses internal algorithms (Singleton or Temperton) or !
668	! system-specific routines, if they are available !
669	!----------------------------------------------------------------------!
670
671	IMPLICIT NONE
672
673	CHARACTER (LEN=*) :: direction
674	INTEGER :: i, ishape(1)
675
676	LOGICAL :: forward_fft
677
678	REAL, DIMENSION(0:nx) :: ar
679	REAL, DIMENSION(0:nx+2) :: work
680	REAL, DIMENSION(nx+2) :: work1
681	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
682	#if defined( __ibm )
683	REAL, DIMENSION(nau2) :: aux2, aux4
684	#elif defined( __nec )
685	REAL, DIMENSION(6*(nx+1)) :: work2
686	#endif
687
688	IF ( direction == 'forward' ) THEN
689	forward_fft = .TRUE.
690	ELSE
691	forward_fft = .FALSE.
692	ENDIF
693
694	IF ( fft_method == 'singleton-algorithm' ) THEN
695
696	!
697	!-- Performing the fft with singleton's software works on every system,
698	!-- since it is part of the model
699	ALLOCATE( cwork(0:nx) )
700
701	IF ( forward_fft ) then
702
703	DO i = 0, nx
704	cwork(i) = CMPLX( ar(i) )
705	ENDDO
706	ishape = SHAPE( cwork )
707	CALL FFTN( cwork, ishape )
708	DO i = 0, (nx+1)/2
709	ar(i) = REAL( cwork(i) )
710	ENDDO
711	DO i = 1, (nx+1)/2 - 1
712	ar(nx+1-i) = -AIMAG( cwork(i) )
713	ENDDO
714
715	ELSE
716
717	cwork(0) = CMPLX( ar(0), 0.0 )
718	DO i = 1, (nx+1)/2 - 1
719	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
720	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
721	ENDDO
722	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
723
724	ishape = SHAPE( cwork )
725	CALL FFTN( cwork, ishape, inv = .TRUE. )
726
727	DO i = 0, nx
728	ar(i) = REAL( cwork(i) )
729	ENDDO
730
731	ENDIF
732
733	DEALLOCATE( cwork )
734
735	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
736
737	!
738	!-- Performing the fft with Temperton's software works on every system,
739	!-- since it is part of the model
740	IF ( forward_fft ) THEN
741
742	work(0:nx) = ar
743	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
744
745	DO i = 0, (nx+1)/2
746	ar(i) = work(2*i)
747	ENDDO
748	DO i = 1, (nx+1)/2 - 1
749	ar(nx+1-i) = work(2*i+1)
750	ENDDO
751
752	ELSE
753
754	DO i = 0, (nx+1)/2
755	work(2*i) = ar(i)
756	ENDDO
757	DO i = 1, (nx+1)/2 - 1
758	work(2*i+1) = ar(nx+1-i)
759	ENDDO
760	work(1) = 0.0
761	work(nx+2) = 0.0
762
763	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
764	ar = work(0:nx)
765
766	ENDIF
767
768	ELSEIF ( fft_method == 'system-specific' ) THEN
769
770	#if defined( __ibm ) && ! defined( __ibmy_special )
771	IF ( forward_fft ) THEN
772
773	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
774	aux2, nau2 )
775
776	DO i = 0, (nx+1)/2
777	ar(i) = work(2*i)
778	ENDDO
779	DO i = 1, (nx+1)/2 - 1
780	ar(nx+1-i) = work(2*i+1)
781	ENDDO
782
783	ELSE
784
785	DO i = 0, (nx+1)/2
786	work(2*i) = ar(i)
787	ENDDO
788	DO i = 1, (nx+1)/2 - 1
789	work(2*i+1) = ar(nx+1-i)
790	ENDDO
791	work(1) = 0.0
792	work(nx+2) = 0.0
793
794	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
795	aux4, nau2 )
796
797	DO i = 0, nx
798	ar(i) = work(i)
799	ENDDO
800
801	ENDIF
802	#elif defined( __nec )
803	IF ( forward_fft ) THEN
804
805	work(0:nx) = ar(0:nx)
806
807	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
808
809	DO i = 0, (nx+1)/2
810	ar(i) = work(2*i)
811	ENDDO
812	DO i = 1, (nx+1)/2 - 1
813	ar(nx+1-i) = work(2*i+1)
814	ENDDO
815
816	ELSE
817
818	DO i = 0, (nx+1)/2
819	work(2*i) = ar(i)
820	ENDDO
821	DO i = 1, (nx+1)/2 - 1
822	work(2*i+1) = ar(nx+1-i)
823	ENDDO
824	work(1) = 0.0
825	work(nx+2) = 0.0
826
827	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
828
829	ar(0:nx) = work(0:nx)
830
831	ENDIF
832	#else
833	message_string = 'no system-specific fft-call available'
834	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
835	#endif
836	ELSE
837	message_string = 'fft method "' // TRIM( fft_method) // &
838	'" not available'
839	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
840
841	ENDIF
842
843	END SUBROUTINE fft_x_1d
844
845	SUBROUTINE fft_y( ar, direction )
846
847	!----------------------------------------------------------------------!
848	! fft_y !
849	! !
850	! Fourier-transformation along y-direction !
851	! Version for 2D-decomposition !
852	! !
853	! fft_y uses internal algorithms (Singleton or Temperton) or !
854	! system-specific routines, if they are available !
855	!----------------------------------------------------------------------!
856
857	USE cuda_fft_interfaces
858	#if defined( __cuda_fft )
859	USE ISO_C_BINDING
860	#endif
861
862	IMPLICIT NONE
863
864	CHARACTER (LEN=*) :: direction
865	INTEGER :: i, j, jshape(1), k
866
867	LOGICAL :: forward_fft
868
869	REAL, DIMENSION(0:ny+2) :: work
870	REAL, DIMENSION(ny+2) :: work1
871	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
872	#if defined( __ibm )
873	REAL, DIMENSION(nau2) :: auy2, auy4
874	#elif defined( __nec )
875	REAL, DIMENSION(6*(ny+1)) :: work2
876	#elif defined( __cuda_fft )
877	!$acc declare create( ar_tmp )
878	COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
879	#endif
880	REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) :: ar
881
882	IF ( direction == 'forward' ) THEN
883	forward_fft = .TRUE.
884	ELSE
885	forward_fft = .FALSE.
886	ENDIF
887
888	IF ( fft_method == 'singleton-algorithm' ) THEN
889
890	!
891	!-- Performing the fft with singleton's software works on every system,
892	!-- since it is part of the model
893	ALLOCATE( cwork(0:ny) )
894
895	IF ( forward_fft ) then
896
897	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
898	!$OMP DO
899	DO k = nzb_y, nzt_y
900	DO i = nxl_y, nxr_y
901
902	DO j = 0, ny
903	cwork(j) = CMPLX( ar(j,i,k) )
904	ENDDO
905
906	jshape = SHAPE( cwork )
907	CALL FFTN( cwork, jshape )
908
909	DO j = 0, (ny+1)/2
910	ar(j,i,k) = REAL( cwork(j) )
911	ENDDO
912	DO j = 1, (ny+1)/2 - 1
913	ar(ny+1-j,i,k) = -AIMAG( cwork(j) )
914	ENDDO
915
916	ENDDO
917	ENDDO
918	!$OMP END PARALLEL
919
920	ELSE
921
922	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
923	!$OMP DO
924	DO k = nzb_y, nzt_y
925	DO i = nxl_y, nxr_y
926
927	cwork(0) = CMPLX( ar(0,i,k), 0.0 )
928	DO j = 1, (ny+1)/2 - 1
929	cwork(j) = CMPLX( ar(j,i,k), -ar(ny+1-j,i,k) )
930	cwork(ny+1-j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
931	ENDDO
932	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
933
934	jshape = SHAPE( cwork )
935	CALL FFTN( cwork, jshape, inv = .TRUE. )
936
937	DO j = 0, ny
938	ar(j,i,k) = REAL( cwork(j) )
939	ENDDO
940
941	ENDDO
942	ENDDO
943	!$OMP END PARALLEL
944
945	ENDIF
946
947	DEALLOCATE( cwork )
948
949	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
950
951	!
952	!-- Performing the fft with Temperton's software works on every system,
953	!-- since it is part of the model
954	IF ( forward_fft ) THEN
955
956	!$OMP PARALLEL PRIVATE ( work, i, j, k )
957	!$OMP DO
958	DO k = nzb_y, nzt_y
959	DO i = nxl_y, nxr_y
960
961	work(0:ny) = ar(0:ny,i,k)
962	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
963
964	DO j = 0, (ny+1)/2
965	ar(j,i,k) = work(2*j)
966	ENDDO
967	DO j = 1, (ny+1)/2 - 1
968	ar(ny+1-j,i,k) = work(2*j+1)
969	ENDDO
970
971	ENDDO
972	ENDDO
973	!$OMP END PARALLEL
974
975	ELSE
976
977	!$OMP PARALLEL PRIVATE ( work, i, j, k )
978	!$OMP DO
979	DO k = nzb_y, nzt_y
980	DO i = nxl_y, nxr_y
981
982	DO j = 0, (ny+1)/2
983	work(2*j) = ar(j,i,k)
984	ENDDO
985	DO j = 1, (ny+1)/2 - 1
986	work(2*j+1) = ar(ny+1-j,i,k)
987	ENDDO
988	work(1) = 0.0
989	work(ny+2) = 0.0
990
991	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
992	ar(0:ny,i,k) = work(0:ny)
993
994	ENDDO
995	ENDDO
996	!$OMP END PARALLEL
997
998	ENDIF
999
1000	ELSEIF ( fft_method == 'fftw' ) THEN
1001
1002	#if defined( __fftw )
1003	IF ( forward_fft ) THEN
1004
1005	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1006	!$OMP DO
1007	DO k = nzb_y, nzt_y
1008	DO i = nxl_y, nxr_y
1009
1010	y_in(0:ny) = ar(0:ny,i,k)
1011	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1012
1013	DO j = 0, (ny+1)/2
1014	ar(j,i,k) = REAL( y_out(j) ) /(ny+1)
1015	ENDDO
1016	DO j = 1, (ny+1)/2 - 1
1017	ar(ny+1-j,i,k) = AIMAG( y_out(j) ) /(ny+1)
1018	ENDDO
1019
1020	ENDDO
1021	ENDDO
1022	!$OMP END PARALLEL
1023
1024	ELSE
1025
1026	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1027	!$OMP DO
1028	DO k = nzb_y, nzt_y
1029	DO i = nxl_y, nxr_y
1030
1031	y_out(0) = CMPLX( ar(0,i,k), 0.0 )
1032	DO j = 1, (ny+1)/2 - 1
1033	y_out(j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1034	ENDDO
1035	y_out((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1036
1037	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1038	ar(0:ny,i,k) = y_in(0:ny)
1039
1040	ENDDO
1041	ENDDO
1042	!$OMP END PARALLEL
1043
1044	ENDIF
1045	#endif
1046
1047	ELSEIF ( fft_method == 'system-specific' ) THEN
1048
1049	#if defined( __ibm ) && ! defined( __ibmy_special )
1050	IF ( forward_fft) THEN
1051
1052	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1053	!$OMP DO
1054	DO k = nzb_y, nzt_y
1055	DO i = nxl_y, nxr_y
1056
1057	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1058	auy2, nau2 )
1059
1060	DO j = 0, (ny+1)/2
1061	ar(j,i,k) = work(2*j)
1062	ENDDO
1063	DO j = 1, (ny+1)/2 - 1
1064	ar(ny+1-j,i,k) = work(2*j+1)
1065	ENDDO
1066
1067	ENDDO
1068	ENDDO
1069	!$OMP END PARALLEL
1070
1071	ELSE
1072
1073	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1074	!$OMP DO
1075	DO k = nzb_y, nzt_y
1076	DO i = nxl_y, nxr_y
1077
1078	DO j = 0, (ny+1)/2
1079	work(2*j) = ar(j,i,k)
1080	ENDDO
1081	DO j = 1, (ny+1)/2 - 1
1082	work(2*j+1) = ar(ny+1-j,i,k)
1083	ENDDO
1084	work(1) = 0.0
1085	work(ny+2) = 0.0
1086
1087	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1088	auy4, nau2 )
1089
1090	DO j = 0, ny
1091	ar(j,i,k) = work(j)
1092	ENDDO
1093
1094	ENDDO
1095	ENDDO
1096	!$OMP END PARALLEL
1097
1098	ENDIF
1099	#elif defined( __nec )
1100	IF ( forward_fft ) THEN
1101
1102	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1103	!$OMP DO
1104	DO k = nzb_y, nzt_y
1105	DO i = nxl_y, nxr_y
1106
1107	work(0:ny) = ar(0:ny,i,k)
1108
1109	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1110
1111	DO j = 0, (ny+1)/2
1112	ar(j,i,k) = work(2*j)
1113	ENDDO
1114	DO j = 1, (ny+1)/2 - 1
1115	ar(ny+1-j,i,k) = work(2*j+1)
1116	ENDDO
1117
1118	ENDDO
1119	ENDDO
1120	!$END OMP PARALLEL
1121
1122	ELSE
1123
1124	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1125	!$OMP DO
1126	DO k = nzb_y, nzt_y
1127	DO i = nxl_y, nxr_y
1128
1129	DO j = 0, (ny+1)/2
1130	work(2*j) = ar(j,i,k)
1131	ENDDO
1132	DO j = 1, (ny+1)/2 - 1
1133	work(2*j+1) = ar(ny+1-j,i,k)
1134	ENDDO
1135	work(1) = 0.0
1136	work(ny+2) = 0.0
1137
1138	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1139
1140	ar(0:ny,i,k) = work(0:ny)
1141
1142	ENDDO
1143	ENDDO
1144	!$OMP END PARALLEL
1145
1146	ENDIF
1147	#elif defined( __cuda_fft )
1148
1149	IF ( forward_fft ) THEN
1150
1151	!$acc data present( ar )
1152	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1153
1154	!$acc kernels
1155	!$acc loop
1156	DO k = nzb_y, nzt_y
1157	DO i = nxl_y, nxr_y
1158
1159	!$acc loop vector( 32 )
1160	DO j = 0, (ny+1)/2
1161	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
1162	ENDDO
1163
1164	!$acc loop vector( 32 )
1165	DO j = 1, (ny+1)/2 - 1
1166	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1167	ENDDO
1168
1169	ENDDO
1170	ENDDO
1171	!$acc end kernels
1172	!$acc end data
1173
1174	ELSE
1175
1176	!$acc data present( ar )
1177	!$acc kernels
1178	!$acc loop
1179	DO k = nzb_y, nzt_y
1180	DO i = nxl_y, nxr_y
1181
1182	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1183
1184	!$acc loop vector( 32 )
1185	DO j = 1, (ny+1)/2 - 1
1186	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1187	ENDDO
1188	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1189
1190	ENDDO
1191	ENDDO
1192	!$acc end kernels
1193
1194	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1195	!$acc end data
1196
1197	ENDIF
1198
1199	#else
1200	message_string = 'no system-specific fft-call available'
1201	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
1202	#endif
1203
1204	ELSE
1205
1206	message_string = 'fft method "' // TRIM( fft_method) // &
1207	'" not available'
1208	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
1209
1210	ENDIF
1211
1212	END SUBROUTINE fft_y
1213
1214	SUBROUTINE fft_y_1d( ar, direction )
1215
1216	!----------------------------------------------------------------------!
1217	! fft_y_1d !
1218	! !
1219	! Fourier-transformation along y-direction !
1220	! Version for 1D-decomposition !
1221	! !
1222	! fft_y uses internal algorithms (Singleton or Temperton) or !
1223	! system-specific routines, if they are available !
1224	!----------------------------------------------------------------------!
1225
1226	IMPLICIT NONE
1227
1228	CHARACTER (LEN=*) :: direction
1229	INTEGER :: j, jshape(1)
1230
1231	LOGICAL :: forward_fft
1232
1233	REAL, DIMENSION(0:ny) :: ar
1234	REAL, DIMENSION(0:ny+2) :: work
1235	REAL, DIMENSION(ny+2) :: work1
1236	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
1237	#if defined( __ibm )
1238	REAL, DIMENSION(nau2) :: auy2, auy4
1239	#elif defined( __nec )
1240	REAL, DIMENSION(6*(ny+1)) :: work2
1241	#endif
1242
1243	IF ( direction == 'forward' ) THEN
1244	forward_fft = .TRUE.
1245	ELSE
1246	forward_fft = .FALSE.
1247	ENDIF
1248
1249	IF ( fft_method == 'singleton-algorithm' ) THEN
1250
1251	!
1252	!-- Performing the fft with singleton's software works on every system,
1253	!-- since it is part of the model
1254	ALLOCATE( cwork(0:ny) )
1255
1256	IF ( forward_fft ) THEN
1257
1258	DO j = 0, ny
1259	cwork(j) = CMPLX( ar(j) )
1260	ENDDO
1261
1262	jshape = SHAPE( cwork )
1263	CALL FFTN( cwork, jshape )
1264
1265	DO j = 0, (ny+1)/2
1266	ar(j) = REAL( cwork(j) )
1267	ENDDO
1268	DO j = 1, (ny+1)/2 - 1
1269	ar(ny+1-j) = -AIMAG( cwork(j) )
1270	ENDDO
1271
1272	ELSE
1273
1274	cwork(0) = CMPLX( ar(0), 0.0 )
1275	DO j = 1, (ny+1)/2 - 1
1276	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
1277	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
1278	ENDDO
1279	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
1280
1281	jshape = SHAPE( cwork )
1282	CALL FFTN( cwork, jshape, inv = .TRUE. )
1283
1284	DO j = 0, ny
1285	ar(j) = REAL( cwork(j) )
1286	ENDDO
1287
1288	ENDIF
1289
1290	DEALLOCATE( cwork )
1291
1292	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1293
1294	!
1295	!-- Performing the fft with Temperton's software works on every system,
1296	!-- since it is part of the model
1297	IF ( forward_fft ) THEN
1298
1299	work(0:ny) = ar
1300	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1301
1302	DO j = 0, (ny+1)/2
1303	ar(j) = work(2*j)
1304	ENDDO
1305	DO j = 1, (ny+1)/2 - 1
1306	ar(ny+1-j) = work(2*j+1)
1307	ENDDO
1308
1309	ELSE
1310
1311	DO j = 0, (ny+1)/2
1312	work(2*j) = ar(j)
1313	ENDDO
1314	DO j = 1, (ny+1)/2 - 1
1315	work(2*j+1) = ar(ny+1-j)
1316	ENDDO
1317	work(1) = 0.0
1318	work(ny+2) = 0.0
1319
1320	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1321	ar = work(0:ny)
1322
1323	ENDIF
1324
1325	ELSEIF ( fft_method == 'system-specific' ) THEN
1326
1327	#if defined( __ibm ) && ! defined( __ibmy_special )
1328	IF ( forward_fft ) THEN
1329
1330	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1331	auy2, nau2 )
1332
1333	DO j = 0, (ny+1)/2
1334	ar(j) = work(2*j)
1335	ENDDO
1336	DO j = 1, (ny+1)/2 - 1
1337	ar(ny+1-j) = work(2*j+1)
1338	ENDDO
1339
1340	ELSE
1341
1342	DO j = 0, (ny+1)/2
1343	work(2*j) = ar(j)
1344	ENDDO
1345	DO j = 1, (ny+1)/2 - 1
1346	work(2*j+1) = ar(ny+1-j)
1347	ENDDO
1348	work(1) = 0.0
1349	work(ny+2) = 0.0
1350
1351	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1352	auy4, nau2 )
1353
1354	DO j = 0, ny
1355	ar(j) = work(j)
1356	ENDDO
1357
1358	ENDIF
1359	#elif defined( __nec )
1360	IF ( forward_fft ) THEN
1361
1362	work(0:ny) = ar(0:ny)
1363
1364	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1365
1366	DO j = 0, (ny+1)/2
1367	ar(j) = work(2*j)
1368	ENDDO
1369	DO j = 1, (ny+1)/2 - 1
1370	ar(ny+1-j) = work(2*j+1)
1371	ENDDO
1372
1373	ELSE
1374
1375	DO j = 0, (ny+1)/2
1376	work(2*j) = ar(j)
1377	ENDDO
1378	DO j = 1, (ny+1)/2 - 1
1379	work(2*j+1) = ar(ny+1-j)
1380	ENDDO
1381	work(1) = 0.0
1382	work(ny+2) = 0.0
1383
1384	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1385
1386	ar(0:ny) = work(0:ny)
1387
1388	ENDIF
1389	#else
1390	message_string = 'no system-specific fft-call available'
1391	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
1392
1393	#endif
1394
1395	ELSE
1396
1397	message_string = 'fft method "' // TRIM( fft_method) // &
1398	'" not available'
1399	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1400
1401	ENDIF
1402
1403	END SUBROUTINE fft_y_1d
1404
1405	SUBROUTINE fft_x_m( ar, direction )
1406
1407	!----------------------------------------------------------------------!
1408	! fft_x_m !
1409	! !
1410	! Fourier-transformation along x-direction !
1411	! Version for 1d domain decomposition !
1412	! using multiple 1D FFT from Math Keisan on NEC !
1413	! or Temperton-algorithm !
1414	! (no singleton-algorithm on NEC because it does not vectorize) !
1415	! !
1416	!----------------------------------------------------------------------!
1417
1418	IMPLICIT NONE
1419
1420	CHARACTER (LEN=*) :: direction
1421	INTEGER :: i, k, siza
1422
1423	REAL, DIMENSION(0:nx,nz) :: ar
1424	REAL, DIMENSION(0:nx+3,nz+1) :: ai
1425	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
1426	#if defined( __nec )
1427	INTEGER :: sizw
1428	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
1429	#endif
1430
1431	IF ( fft_method == 'temperton-algorithm' ) THEN
1432
1433	siza = SIZE( ai, 1 )
1434
1435	IF ( direction == 'forward') THEN
1436
1437	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1438	ai(nx+1:,:) = 0.0
1439
1440	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1441
1442	DO k = 1, nz
1443	DO i = 0, (nx+1)/2
1444	ar(i,k) = ai(2*i,k)
1445	ENDDO
1446	DO i = 1, (nx+1)/2 - 1
1447	ar(nx+1-i,k) = ai(2*i+1,k)
1448	ENDDO
1449	ENDDO
1450
1451	ELSE
1452
1453	DO k = 1, nz
1454	DO i = 0, (nx+1)/2
1455	ai(2*i,k) = ar(i,k)
1456	ENDDO
1457	DO i = 1, (nx+1)/2 - 1
1458	ai(2*i+1,k) = ar(nx+1-i,k)
1459	ENDDO
1460	ai(1,k) = 0.0
1461	ai(nx+2,k) = 0.0
1462	ENDDO
1463
1464	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1465
1466	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1467
1468	ENDIF
1469
1470	ELSEIF ( fft_method == 'system-specific' ) THEN
1471
1472	#if defined( __nec )
1473	siza = SIZE( ai, 1 )
1474	sizw = SIZE( work, 1 )
1475
1476	IF ( direction == 'forward') THEN
1477
1478	!
1479	!-- Tables are initialized once more. This call should not be
1480	!-- necessary, but otherwise program aborts in asymmetric case
1481	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1482	trig_xf, work1, 0 )
1483
1484	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1485	IF ( nz1 > nz ) THEN
1486	ai(:,nz1) = 0.0
1487	ENDIF
1488
1489	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
1490	trig_xf, work1, 0 )
1491
1492	DO k = 1, nz
1493	DO i = 0, (nx+1)/2
1494	ar(i,k) = REAL( work(i+1,k) )
1495	ENDDO
1496	DO i = 1, (nx+1)/2 - 1
1497	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1498	ENDDO
1499	ENDDO
1500
1501	ELSE
1502
1503	!
1504	!-- Tables are initialized once more. This call should not be
1505	!-- necessary, but otherwise program aborts in asymmetric case
1506	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1507	trig_xb, work1, 0 )
1508
1509	IF ( nz1 > nz ) THEN
1510	work(:,nz1) = 0.0
1511	ENDIF
1512	DO k = 1, nz
1513	work(1,k) = CMPLX( ar(0,k), 0.0 )
1514	DO i = 1, (nx+1)/2 - 1
1515	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
1516	ENDDO
1517	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
1518	ENDDO
1519
1520	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1521	trig_xb, work1, 0 )
1522
1523	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1524
1525	ENDIF
1526
1527	#else
1528	message_string = 'no system-specific fft-call available'
1529	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
1530	#endif
1531
1532	ELSE
1533
1534	message_string = 'fft method "' // TRIM( fft_method) // &
1535	'" not available'
1536	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1537
1538	ENDIF
1539
1540	END SUBROUTINE fft_x_m
1541
1542	SUBROUTINE fft_y_m( ar, ny1, direction )
1543
1544	!----------------------------------------------------------------------!
1545	! fft_y_m !
1546	! !
1547	! Fourier-transformation along y-direction !
1548	! Version for 1d domain decomposition !
1549	! using multiple 1D FFT from Math Keisan on NEC !
1550	! or Temperton-algorithm !
1551	! (no singleton-algorithm on NEC because it does not vectorize) !
1552	! !
1553	!----------------------------------------------------------------------!
1554
1555	IMPLICIT NONE
1556
1557	CHARACTER (LEN=*) :: direction
1558	INTEGER :: j, k, ny1, siza
1559
1560	REAL, DIMENSION(0:ny1,nz) :: ar
1561	REAL, DIMENSION(0:ny+3,nz+1) :: ai
1562	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
1563	#if defined( __nec )
1564	INTEGER :: sizw
1565	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
1566	#endif
1567
1568	IF ( fft_method == 'temperton-algorithm' ) THEN
1569
1570	siza = SIZE( ai, 1 )
1571
1572	IF ( direction == 'forward') THEN
1573
1574	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1575	ai(ny+1:,:) = 0.0
1576
1577	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1578
1579	DO k = 1, nz
1580	DO j = 0, (ny+1)/2
1581	ar(j,k) = ai(2*j,k)
1582	ENDDO
1583	DO j = 1, (ny+1)/2 - 1
1584	ar(ny+1-j,k) = ai(2*j+1,k)
1585	ENDDO
1586	ENDDO
1587
1588	ELSE
1589
1590	DO k = 1, nz
1591	DO j = 0, (ny+1)/2
1592	ai(2*j,k) = ar(j,k)
1593	ENDDO
1594	DO j = 1, (ny+1)/2 - 1
1595	ai(2*j+1,k) = ar(ny+1-j,k)
1596	ENDDO
1597	ai(1,k) = 0.0
1598	ai(ny+2,k) = 0.0
1599	ENDDO
1600
1601	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1602
1603	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1604
1605	ENDIF
1606
1607	ELSEIF ( fft_method == 'system-specific' ) THEN
1608
1609	#if defined( __nec )
1610	siza = SIZE( ai, 1 )
1611	sizw = SIZE( work, 1 )
1612
1613	IF ( direction == 'forward') THEN
1614
1615	!
1616	!-- Tables are initialized once more. This call should not be
1617	!-- necessary, but otherwise program aborts in asymmetric case
1618	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1619	trig_yf, work1, 0 )
1620
1621	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1622	IF ( nz1 > nz ) THEN
1623	ai(:,nz1) = 0.0
1624	ENDIF
1625
1626	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
1627	trig_yf, work1, 0 )
1628
1629	DO k = 1, nz
1630	DO j = 0, (ny+1)/2
1631	ar(j,k) = REAL( work(j+1,k) )
1632	ENDDO
1633	DO j = 1, (ny+1)/2 - 1
1634	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1635	ENDDO
1636	ENDDO
1637
1638	ELSE
1639
1640	!
1641	!-- Tables are initialized once more. This call should not be
1642	!-- necessary, but otherwise program aborts in asymmetric case
1643	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1644	trig_yb, work1, 0 )
1645
1646	IF ( nz1 > nz ) THEN
1647	work(:,nz1) = 0.0
1648	ENDIF
1649	DO k = 1, nz
1650	work(1,k) = CMPLX( ar(0,k), 0.0 )
1651	DO j = 1, (ny+1)/2 - 1
1652	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
1653	ENDDO
1654	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
1655	ENDDO
1656
1657	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1658	trig_yb, work1, 0 )
1659
1660	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1661
1662	ENDIF
1663
1664	#else
1665	message_string = 'no system-specific fft-call available'
1666	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
1667	#endif
1668
1669	ELSE
1670
1671	message_string = 'fft method "' // TRIM( fft_method) // &
1672	'" not available'
1673	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1674
1675	ENDIF
1676
1677	END SUBROUTINE fft_y_m
1678
1679
1680	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |