Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |