Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |