Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |