Home

Context Navigation

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

Last change on this file since 1219 was 1219, checked in by heinze, 11 years ago
bugfix: use own branch for fftw setting
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	! bugfix: use own branch for fftw
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: fft_xy.f90 1219 2013-08-30 09:33:18Z heinze $
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	#endif
126
127	#if defined( __fftw )
128	INCLUDE 'fftw3.f03'
129	INTEGER(KIND=C_INT) :: nx_c, ny_c
130	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: x_out, y_out
131	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: x_in, y_in
132	TYPE(C_PTR), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
133	#endif
134
135	!
136	!-- Public interfaces
137	INTERFACE fft_init
138	MODULE PROCEDURE fft_init
139	END INTERFACE fft_init
140
141	INTERFACE fft_x
142	MODULE PROCEDURE fft_x
143	END INTERFACE fft_x
144
145	INTERFACE fft_x_1d
146	MODULE PROCEDURE fft_x_1d
147	END INTERFACE fft_x_1d
148
149	INTERFACE fft_y
150	MODULE PROCEDURE fft_y
151	END INTERFACE fft_y
152
153	INTERFACE fft_y_1d
154	MODULE PROCEDURE fft_y_1d
155	END INTERFACE fft_y_1d
156
157	INTERFACE fft_x_m
158	MODULE PROCEDURE fft_x_m
159	END INTERFACE fft_x_m
160
161	INTERFACE fft_y_m
162	MODULE PROCEDURE fft_y_m
163	END INTERFACE fft_y_m
164
165	CONTAINS
166
167
168	SUBROUTINE fft_init
169
170	USE cuda_fft_interfaces
171
172	IMPLICIT NONE
173
174	!
175	!-- The following temporary working arrays have to be on stack or private
176	!-- in OpenMP sense
177	#if defined( __ibm )
178	REAL, DIMENSION(0:nx+2) :: workx
179	REAL, DIMENSION(0:ny+2) :: worky
180	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
181	#elif defined( __nec )
182	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
183	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
184	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
185	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
186	#endif
187
188	!
189	!-- Return, if already called
190	IF ( init_fft ) THEN
191	RETURN
192	ELSE
193	init_fft = .TRUE.
194	ENDIF
195
196	IF ( fft_method == 'system-specific' ) THEN
197
198	dnx = 1.0 / ( nx + 1.0 )
199	dny = 1.0 / ( ny + 1.0 )
200	sqr_dnx = SQRT( dnx )
201	sqr_dny = SQRT( dny )
202	#if defined( __ibm ) && ! defined( __ibmy_special )
203	!
204	!-- Initialize tables for fft along x
205	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
206	aux2, nau2 )
207	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
208	aux4, nau2 )
209	!
210	!-- Initialize tables for fft along y
211	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
212	auy2, nau2 )
213	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
214	auy4, nau2 )
215	#elif defined( __nec )
216	message_string = 'fft method "' // TRIM( fft_method) // &
217	'" currently does not work on NEC'
218	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
219
220	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
221	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
222
223	work_x = 0.0
224	work_y = 0.0
225	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
226	! when using the NEC ffts
227
228	!
229	!-- Initialize tables for fft along x (non-vector and vector case (M))
230	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
231	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
232	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
233	trig_xf, workx, 0 )
234	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
235	trig_xb, workx, 0 )
236	!
237	!-- Initialize tables for fft along y (non-vector and vector case (M))
238	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
239	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
240	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
241	trig_yf, worky, 0 )
242	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
243	trig_yb, worky, 0 )
244	#elif defined( __cuda_fft )
245	total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
246	total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
247	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
248	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
249	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
250	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
251	#else
252	message_string = 'no system-specific fft-call available'
253	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
254	#endif
255	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
256	!
257	!-- Temperton-algorithm
258	!-- Initialize tables for fft along x and y
259	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
260
261	CALL set99( trigs_x, ifax_x, nx+1 )
262	CALL set99( trigs_y, ifax_y, ny+1 )
263
264	ELSEIF ( fft_method == 'fftw' ) THEN
265	!
266	!-- FFTW
267	#if defined( __fftw )
268	nx_c = nx+1
269	ny_c = ny+1
270	ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), &
271	y_out(0:(ny+1)/2) )
272	plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
273	plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
274	plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
275	plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
276	#else
277	message_string = 'preprocessor switch for fftw is missing'
278	CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
279	#endif
280
281	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
282
283	CONTINUE
284
285	ELSE
286
287	message_string = 'fft method "' // TRIM( fft_method) // &
288	'" not available'
289	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
290	ENDIF
291
292	END SUBROUTINE fft_init
293
294
295	SUBROUTINE fft_x( ar, direction, ar_2d )
296
297	!----------------------------------------------------------------------!
298	! fft_x !
299	! !
300	! Fourier-transformation along x-direction !
301	! Version for 2D-decomposition !
302	! !
303	! fft_x uses internal algorithms (Singleton or Temperton) or !
304	! system-specific routines, if they are available !
305	!----------------------------------------------------------------------!
306
307	USE cuda_fft_interfaces
308	#if defined( __cuda_fft )
309	USE ISO_C_BINDING
310	#endif
311
312	IMPLICIT NONE
313
314	CHARACTER (LEN=*) :: direction
315	INTEGER :: i, ishape(1), j, k
316
317	LOGICAL :: forward_fft
318
319	REAL, DIMENSION(0:nx+2) :: work
320	REAL, DIMENSION(nx+2) :: work1
321	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
322	#if defined( __ibm )
323	REAL, DIMENSION(nau2) :: aux2, aux4
324	#elif defined( __nec )
325	REAL, DIMENSION(6*(nx+1)) :: work2
326	#elif defined( __cuda_fft )
327	!$acc declare create( ar_tmp )
328	COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp
329	#endif
330	REAL, DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL :: ar_2d
331	REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar
332
333	IF ( direction == 'forward' ) THEN
334	forward_fft = .TRUE.
335	ELSE
336	forward_fft = .FALSE.
337	ENDIF
338
339	IF ( fft_method == 'singleton-algorithm' ) THEN
340
341	!
342	!-- Performing the fft with singleton's software works on every system,
343	!-- since it is part of the model
344	ALLOCATE( cwork(0:nx) )
345
346	IF ( forward_fft ) then
347
348	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
349	!$OMP DO
350	DO k = nzb_x, nzt_x
351	DO j = nys_x, nyn_x
352
353	DO i = 0, nx
354	cwork(i) = CMPLX( ar(i,j,k) )
355	ENDDO
356
357	ishape = SHAPE( cwork )
358	CALL FFTN( cwork, ishape )
359
360	DO i = 0, (nx+1)/2
361	ar(i,j,k) = REAL( cwork(i) )
362	ENDDO
363	DO i = 1, (nx+1)/2 - 1
364	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
365	ENDDO
366
367	ENDDO
368	ENDDO
369	!$OMP END PARALLEL
370
371	ELSE
372
373	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
374	!$OMP DO
375	DO k = nzb_x, nzt_x
376	DO j = nys_x, nyn_x
377
378	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
379	DO i = 1, (nx+1)/2 - 1
380	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
381	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
382	ENDDO
383	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
384
385	ishape = SHAPE( cwork )
386	CALL FFTN( cwork, ishape, inv = .TRUE. )
387
388	DO i = 0, nx
389	ar(i,j,k) = REAL( cwork(i) )
390	ENDDO
391
392	ENDDO
393	ENDDO
394	!$OMP END PARALLEL
395
396	ENDIF
397
398	DEALLOCATE( cwork )
399
400	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
401
402	!
403	!-- Performing the fft with Temperton's software works on every system,
404	!-- since it is part of the model
405	IF ( forward_fft ) THEN
406
407	!$OMP PARALLEL PRIVATE ( work, i, j, k )
408	!$OMP DO
409	DO k = nzb_x, nzt_x
410	DO j = nys_x, nyn_x
411
412	work(0:nx) = ar(0:nx,j,k)
413	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
414
415	DO i = 0, (nx+1)/2
416	ar(i,j,k) = work(2*i)
417	ENDDO
418	DO i = 1, (nx+1)/2 - 1
419	ar(nx+1-i,j,k) = work(2*i+1)
420	ENDDO
421
422	ENDDO
423	ENDDO
424	!$OMP END PARALLEL
425
426	ELSE
427
428	!$OMP PARALLEL PRIVATE ( work, i, j, k )
429	!$OMP DO
430	DO k = nzb_x, nzt_x
431	DO j = nys_x, nyn_x
432
433	DO i = 0, (nx+1)/2
434	work(2*i) = ar(i,j,k)
435	ENDDO
436	DO i = 1, (nx+1)/2 - 1
437	work(2*i+1) = ar(nx+1-i,j,k)
438	ENDDO
439	work(1) = 0.0
440	work(nx+2) = 0.0
441
442	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
443	ar(0:nx,j,k) = work(0:nx)
444
445	ENDDO
446	ENDDO
447	!$OMP END PARALLEL
448
449	ENDIF
450
451	ELSEIF ( fft_method == 'fftw' ) THEN
452
453	#if defined( __fftw )
454	IF ( forward_fft ) THEN
455
456	!$OMP PARALLEL PRIVATE ( work, i, j, k )
457	!$OMP DO
458	DO k = nzb_x, nzt_x
459	DO j = nys_x, nyn_x
460
461	x_in(0:nx) = ar(0:nx,j,k)
462	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
463
464	IF ( PRESENT( ar_2d ) ) THEN
465
466	DO i = 0, (nx+1)/2
467	ar_2d(i,j) = REAL( x_out(i) ) / ( nx+1 )
468	ENDDO
469	DO i = 1, (nx+1)/2 - 1
470	ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 )
471	ENDDO
472
473	ELSE
474
475	DO i = 0, (nx+1)/2
476	ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
477	ENDDO
478	DO i = 1, (nx+1)/2 - 1
479	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
480	ENDDO
481
482	ENDIF
483
484	ENDDO
485	ENDDO
486	!$OMP END PARALLEL
487
488	ELSE
489	!$OMP PARALLEL PRIVATE ( work, i, j, k )
490	!$OMP DO
491	DO k = nzb_x, nzt_x
492	DO j = nys_x, nyn_x
493
494	IF ( PRESENT( ar_2d ) ) THEN
495
496	x_out(0) = CMPLX( ar_2d(0,j), 0.0 )
497	DO i = 1, (nx+1)/2 - 1
498	x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) )
499	ENDDO
500	x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 )
501
502	ELSE
503
504	x_out(0) = CMPLX( ar(0,j,k), 0.0 )
505	DO i = 1, (nx+1)/2 - 1
506	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
507	ENDDO
508	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
509
510	ENDIF
511
512	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
513	ar(0:nx,j,k) = x_in(0:nx)
514
515	ENDDO
516	ENDDO
517	!$OMP END PARALLEL
518
519	ENDIF
520	#endif
521
522	ELSEIF ( fft_method == 'system-specific' ) THEN
523
524	#if defined( __ibm ) && ! defined( __ibmy_special )
525	IF ( forward_fft ) THEN
526
527	!$OMP PARALLEL PRIVATE ( work, i, j, k )
528	!$OMP DO
529	DO k = nzb_x, nzt_x
530	DO j = nys_x, nyn_x
531
532	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
533	aux2, nau2 )
534
535	DO i = 0, (nx+1)/2
536	ar(i,j,k) = work(2*i)
537	ENDDO
538	DO i = 1, (nx+1)/2 - 1
539	ar(nx+1-i,j,k) = work(2*i+1)
540	ENDDO
541
542	ENDDO
543	ENDDO
544	!$OMP END PARALLEL
545
546	ELSE
547
548	!$OMP PARALLEL PRIVATE ( work, i, j, k )
549	!$OMP DO
550	DO k = nzb_x, nzt_x
551	DO j = nys_x, nyn_x
552
553	DO i = 0, (nx+1)/2
554	work(2*i) = ar(i,j,k)
555	ENDDO
556	DO i = 1, (nx+1)/2 - 1
557	work(2*i+1) = ar(nx+1-i,j,k)
558	ENDDO
559	work(1) = 0.0
560	work(nx+2) = 0.0
561
562	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
563	aux4, nau2 )
564
565	DO i = 0, nx
566	ar(i,j,k) = work(i)
567	ENDDO
568
569	ENDDO
570	ENDDO
571	!$OMP END PARALLEL
572
573	ENDIF
574
575	#elif defined( __nec )
576
577	IF ( forward_fft ) THEN
578
579	!$OMP PARALLEL PRIVATE ( work, i, j, k )
580	!$OMP DO
581	DO k = nzb_x, nzt_x
582	DO j = nys_x, nyn_x
583
584	work(0:nx) = ar(0:nx,j,k)
585
586	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
587
588	DO i = 0, (nx+1)/2
589	ar(i,j,k) = work(2*i)
590	ENDDO
591	DO i = 1, (nx+1)/2 - 1
592	ar(nx+1-i,j,k) = work(2*i+1)
593	ENDDO
594
595	ENDDO
596	ENDDO
597	!$END OMP PARALLEL
598
599	ELSE
600
601	!$OMP PARALLEL PRIVATE ( work, i, j, k )
602	!$OMP DO
603	DO k = nzb_x, nzt_x
604	DO j = nys_x, nyn_x
605
606	DO i = 0, (nx+1)/2
607	work(2*i) = ar(i,j,k)
608	ENDDO
609	DO i = 1, (nx+1)/2 - 1
610	work(2*i+1) = ar(nx+1-i,j,k)
611	ENDDO
612	work(1) = 0.0
613	work(nx+2) = 0.0
614
615	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
616
617	ar(0:nx,j,k) = work(0:nx)
618
619	ENDDO
620	ENDDO
621	!$OMP END PARALLEL
622
623	ENDIF
624
625	#elif defined( __cuda_fft )
626
627	IF ( forward_fft ) THEN
628
629	!$acc data present( ar )
630	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
631
632	!$acc kernels
633	!$acc loop
634	DO k = nzb_x, nzt_x
635	DO j = nys_x, nyn_x
636
637	!$acc loop vector( 32 )
638	DO i = 0, (nx+1)/2
639	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
640	ENDDO
641
642	!$acc loop vector( 32 )
643	DO i = 1, (nx+1)/2 - 1
644	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
645	ENDDO
646
647	ENDDO
648	ENDDO
649	!$acc end kernels
650	!$acc end data
651
652	ELSE
653
654	!$acc data present( ar )
655	!$acc kernels
656	!$acc loop
657	DO k = nzb_x, nzt_x
658	DO j = nys_x, nyn_x
659
660	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
661
662	!$acc loop vector( 32 )
663	DO i = 1, (nx+1)/2 - 1
664	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
665	ENDDO
666	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
667
668	ENDDO
669	ENDDO
670	!$acc end kernels
671
672	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
673	!$acc end data
674
675	ENDIF
676
677	#else
678	message_string = 'no system-specific fft-call available'
679	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
680	#endif
681
682	ELSE
683
684	message_string = 'fft method "' // TRIM( fft_method) // &
685	'" not available'
686	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
687
688	ENDIF
689
690	END SUBROUTINE fft_x
691
692	SUBROUTINE fft_x_1d( ar, direction )
693
694	!----------------------------------------------------------------------!
695	! fft_x_1d !
696	! !
697	! Fourier-transformation along x-direction !
698	! Version for 1D-decomposition !
699	! !
700	! fft_x uses internal algorithms (Singleton or Temperton) or !
701	! system-specific routines, if they are available !
702	!----------------------------------------------------------------------!
703
704	IMPLICIT NONE
705
706	CHARACTER (LEN=*) :: direction
707	INTEGER :: i, ishape(1)
708
709	LOGICAL :: forward_fft
710
711	REAL, DIMENSION(0:nx) :: ar
712	REAL, DIMENSION(0:nx+2) :: work
713	REAL, DIMENSION(nx+2) :: work1
714	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
715	#if defined( __ibm )
716	REAL, DIMENSION(nau2) :: aux2, aux4
717	#elif defined( __nec )
718	REAL, DIMENSION(6*(nx+1)) :: work2
719	#endif
720
721	IF ( direction == 'forward' ) THEN
722	forward_fft = .TRUE.
723	ELSE
724	forward_fft = .FALSE.
725	ENDIF
726
727	IF ( fft_method == 'singleton-algorithm' ) THEN
728
729	!
730	!-- Performing the fft with singleton's software works on every system,
731	!-- since it is part of the model
732	ALLOCATE( cwork(0:nx) )
733
734	IF ( forward_fft ) then
735
736	DO i = 0, nx
737	cwork(i) = CMPLX( ar(i) )
738	ENDDO
739	ishape = SHAPE( cwork )
740	CALL FFTN( cwork, ishape )
741	DO i = 0, (nx+1)/2
742	ar(i) = REAL( cwork(i) )
743	ENDDO
744	DO i = 1, (nx+1)/2 - 1
745	ar(nx+1-i) = -AIMAG( cwork(i) )
746	ENDDO
747
748	ELSE
749
750	cwork(0) = CMPLX( ar(0), 0.0 )
751	DO i = 1, (nx+1)/2 - 1
752	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
753	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
754	ENDDO
755	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
756
757	ishape = SHAPE( cwork )
758	CALL FFTN( cwork, ishape, inv = .TRUE. )
759
760	DO i = 0, nx
761	ar(i) = REAL( cwork(i) )
762	ENDDO
763
764	ENDIF
765
766	DEALLOCATE( cwork )
767
768	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
769
770	!
771	!-- Performing the fft with Temperton's software works on every system,
772	!-- since it is part of the model
773	IF ( forward_fft ) THEN
774
775	work(0:nx) = ar
776	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
777
778	DO i = 0, (nx+1)/2
779	ar(i) = work(2*i)
780	ENDDO
781	DO i = 1, (nx+1)/2 - 1
782	ar(nx+1-i) = work(2*i+1)
783	ENDDO
784
785	ELSE
786
787	DO i = 0, (nx+1)/2
788	work(2*i) = ar(i)
789	ENDDO
790	DO i = 1, (nx+1)/2 - 1
791	work(2*i+1) = ar(nx+1-i)
792	ENDDO
793	work(1) = 0.0
794	work(nx+2) = 0.0
795
796	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
797	ar = work(0:nx)
798
799	ENDIF
800
801	ELSEIF ( fft_method == 'fftw' ) THEN
802
803	#if defined( __fftw )
804	IF ( forward_fft ) THEN
805
806	x_in(0:nx) = ar(0:nx)
807	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
808
809	DO i = 0, (nx+1)/2
810	ar(i) = REAL( x_out(i) ) / ( nx+1 )
811	ENDDO
812	DO i = 1, (nx+1)/2 - 1
813	ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 )
814	ENDDO
815
816	ELSE
817
818	x_out(0) = CMPLX( ar(0), 0.0 )
819	DO i = 1, (nx+1)/2 - 1
820	x_out(i) = CMPLX( ar(i), ar(nx+1-i) )
821	ENDDO
822	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
823
824	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
825	ar(0:nx) = x_in(0:nx)
826
827	ENDIF
828	#endif
829
830	ELSEIF ( fft_method == 'system-specific' ) THEN
831
832	#if defined( __ibm ) && ! defined( __ibmy_special )
833	IF ( forward_fft ) THEN
834
835	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
836	aux2, nau2 )
837
838	DO i = 0, (nx+1)/2
839	ar(i) = work(2*i)
840	ENDDO
841	DO i = 1, (nx+1)/2 - 1
842	ar(nx+1-i) = work(2*i+1)
843	ENDDO
844
845	ELSE
846
847	DO i = 0, (nx+1)/2
848	work(2*i) = ar(i)
849	ENDDO
850	DO i = 1, (nx+1)/2 - 1
851	work(2*i+1) = ar(nx+1-i)
852	ENDDO
853	work(1) = 0.0
854	work(nx+2) = 0.0
855
856	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
857	aux4, nau2 )
858
859	DO i = 0, nx
860	ar(i) = work(i)
861	ENDDO
862
863	ENDIF
864	#elif defined( __nec )
865	IF ( forward_fft ) THEN
866
867	work(0:nx) = ar(0:nx)
868
869	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
870
871	DO i = 0, (nx+1)/2
872	ar(i) = work(2*i)
873	ENDDO
874	DO i = 1, (nx+1)/2 - 1
875	ar(nx+1-i) = work(2*i+1)
876	ENDDO
877
878	ELSE
879
880	DO i = 0, (nx+1)/2
881	work(2*i) = ar(i)
882	ENDDO
883	DO i = 1, (nx+1)/2 - 1
884	work(2*i+1) = ar(nx+1-i)
885	ENDDO
886	work(1) = 0.0
887	work(nx+2) = 0.0
888
889	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
890
891	ar(0:nx) = work(0:nx)
892
893	ENDIF
894	#else
895	message_string = 'no system-specific fft-call available'
896	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
897	#endif
898	ELSE
899	message_string = 'fft method "' // TRIM( fft_method) // &
900	'" not available'
901	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
902
903	ENDIF
904
905	END SUBROUTINE fft_x_1d
906
907	SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
908	nxr_y_l )
909
910	!----------------------------------------------------------------------!
911	! fft_y !
912	! !
913	! Fourier-transformation along y-direction !
914	! Version for 2D-decomposition !
915	! !
916	! fft_y uses internal algorithms (Singleton or Temperton) or !
917	! system-specific routines, if they are available !
918	! !
919	! direction: 'forward' or 'backward' !
920	! ar, ar_tr: 3D data arrays !
921	! forward: ar: before ar_tr: after transformation !
922	! backward: ar_tr: before ar: after transfosition !
923	! !
924	! In case of non-overlapping transposition/transformation: !
925	! nxl_y_bound = nxl_y_l = nxl_y !
926	! nxr_y_bound = nxr_y_l = nxr_y !
927	! !
928	! In case of overlapping transposition/transformation !
929	! - nxl_y_bound and nxr_y_bound have the original values of !
930	! nxl_y, nxr_y. ar_tr is dimensioned using these values. !
931	! - nxl_y_l = nxr_y_r. ar is dimensioned with these values, so that !
932	! transformation is carried out for a 2D-plane only. !
933	!----------------------------------------------------------------------!
934
935	USE cuda_fft_interfaces
936	#if defined( __cuda_fft )
937	USE ISO_C_BINDING
938	#endif
939
940	IMPLICIT NONE
941
942	CHARACTER (LEN=*) :: direction
943	INTEGER :: i, j, jshape(1), k
944	INTEGER :: nxl_y_bound, nxl_y_l, nxr_y_bound, nxr_y_l
945
946	LOGICAL :: forward_fft
947
948	REAL, DIMENSION(0:ny+2) :: work
949	REAL, DIMENSION(ny+2) :: work1
950	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
951	#if defined( __ibm )
952	REAL, DIMENSION(nau2) :: auy2, auy4
953	#elif defined( __nec )
954	REAL, DIMENSION(6*(ny+1)) :: work2
955	#elif defined( __cuda_fft )
956	!$acc declare create( ar_tmp )
957	COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
958	#endif
959	REAL, DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y) :: ar
960	REAL, DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) :: ar_tr
961
962	IF ( direction == 'forward' ) THEN
963	forward_fft = .TRUE.
964	ELSE
965	forward_fft = .FALSE.
966	ENDIF
967
968	IF ( fft_method == 'singleton-algorithm' ) THEN
969
970	!
971	!-- Performing the fft with singleton's software works on every system,
972	!-- since it is part of the model
973	ALLOCATE( cwork(0:ny) )
974
975	IF ( forward_fft ) then
976
977	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
978	!$OMP DO
979	DO k = nzb_y, nzt_y
980	DO i = nxl_y_l, nxr_y_l
981
982	DO j = 0, ny
983	cwork(j) = CMPLX( ar(j,i,k) )
984	ENDDO
985
986	jshape = SHAPE( cwork )
987	CALL FFTN( cwork, jshape )
988
989	DO j = 0, (ny+1)/2
990	ar_tr(j,i,k) = REAL( cwork(j) )
991	ENDDO
992	DO j = 1, (ny+1)/2 - 1
993	ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
994	ENDDO
995
996	ENDDO
997	ENDDO
998	!$OMP END PARALLEL
999
1000	ELSE
1001
1002	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
1003	!$OMP DO
1004	DO k = nzb_y, nzt_y
1005	DO i = nxl_y_l, nxr_y_l
1006
1007	cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )
1008	DO j = 1, (ny+1)/2 - 1
1009	cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) )
1010	cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
1011	ENDDO
1012	cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
1013
1014	jshape = SHAPE( cwork )
1015	CALL FFTN( cwork, jshape, inv = .TRUE. )
1016
1017	DO j = 0, ny
1018	ar(j,i,k) = REAL( cwork(j) )
1019	ENDDO
1020
1021	ENDDO
1022	ENDDO
1023	!$OMP END PARALLEL
1024
1025	ENDIF
1026
1027	DEALLOCATE( cwork )
1028
1029	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1030
1031	!
1032	!-- Performing the fft with Temperton's software works on every system,
1033	!-- since it is part of the model
1034	IF ( forward_fft ) THEN
1035
1036	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1037	!$OMP DO
1038	DO k = nzb_y, nzt_y
1039	DO i = nxl_y_l, nxr_y_l
1040
1041	work(0:ny) = ar(0:ny,i,k)
1042	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1043
1044	DO j = 0, (ny+1)/2
1045	ar_tr(j,i,k) = work(2*j)
1046	ENDDO
1047	DO j = 1, (ny+1)/2 - 1
1048	ar_tr(ny+1-j,i,k) = work(2*j+1)
1049	ENDDO
1050
1051	ENDDO
1052	ENDDO
1053	!$OMP END PARALLEL
1054
1055	ELSE
1056
1057	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1058	!$OMP DO
1059	DO k = nzb_y, nzt_y
1060	DO i = nxl_y_l, nxr_y_l
1061
1062	DO j = 0, (ny+1)/2
1063	work(2*j) = ar_tr(j,i,k)
1064	ENDDO
1065	DO j = 1, (ny+1)/2 - 1
1066	work(2*j+1) = ar_tr(ny+1-j,i,k)
1067	ENDDO
1068	work(1) = 0.0
1069	work(ny+2) = 0.0
1070
1071	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1072	ar(0:ny,i,k) = work(0:ny)
1073
1074	ENDDO
1075	ENDDO
1076	!$OMP END PARALLEL
1077
1078	ENDIF
1079
1080	ELSEIF ( fft_method == 'fftw' ) THEN
1081
1082	#if defined( __fftw )
1083	IF ( forward_fft ) THEN
1084
1085	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1086	!$OMP DO
1087	DO k = nzb_y, nzt_y
1088	DO i = nxl_y_l, nxr_y_l
1089
1090	y_in(0:ny) = ar(0:ny,i,k)
1091	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1092
1093	DO j = 0, (ny+1)/2
1094	ar_tr(j,i,k) = REAL( y_out(j) ) / (ny+1)
1095	ENDDO
1096	DO j = 1, (ny+1)/2 - 1
1097	ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
1098	ENDDO
1099
1100	ENDDO
1101	ENDDO
1102	!$OMP END PARALLEL
1103
1104	ELSE
1105
1106	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1107	!$OMP DO
1108	DO k = nzb_y, nzt_y
1109	DO i = nxl_y_l, nxr_y_l
1110
1111	y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )
1112	DO j = 1, (ny+1)/2 - 1
1113	y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
1114	ENDDO
1115	y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
1116
1117	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1118	ar(0:ny,i,k) = y_in(0:ny)
1119
1120	ENDDO
1121	ENDDO
1122	!$OMP END PARALLEL
1123
1124	ENDIF
1125	#endif
1126
1127	ELSEIF ( fft_method == 'system-specific' ) THEN
1128
1129	#if defined( __ibm ) && ! defined( __ibmy_special )
1130	IF ( forward_fft) THEN
1131
1132	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1133	!$OMP DO
1134	DO k = nzb_y, nzt_y
1135	DO i = nxl_y_l, nxr_y_l
1136
1137	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1138	auy2, nau2 )
1139
1140	DO j = 0, (ny+1)/2
1141	ar_tr(j,i,k) = work(2*j)
1142	ENDDO
1143	DO j = 1, (ny+1)/2 - 1
1144	ar_tr(ny+1-j,i,k) = work(2*j+1)
1145	ENDDO
1146
1147	ENDDO
1148	ENDDO
1149	!$OMP END PARALLEL
1150
1151	ELSE
1152
1153	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1154	!$OMP DO
1155	DO k = nzb_y, nzt_y
1156	DO i = nxl_y_l, nxr_y_l
1157
1158	DO j = 0, (ny+1)/2
1159	work(2*j) = ar_tr(j,i,k)
1160	ENDDO
1161	DO j = 1, (ny+1)/2 - 1
1162	work(2*j+1) = ar_tr(ny+1-j,i,k)
1163	ENDDO
1164	work(1) = 0.0
1165	work(ny+2) = 0.0
1166
1167	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1168	auy4, nau2 )
1169
1170	DO j = 0, ny
1171	ar(j,i,k) = work(j)
1172	ENDDO
1173
1174	ENDDO
1175	ENDDO
1176	!$OMP END PARALLEL
1177
1178	ENDIF
1179	#elif defined( __nec )
1180	IF ( forward_fft ) THEN
1181
1182	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1183	!$OMP DO
1184	DO k = nzb_y, nzt_y
1185	DO i = nxl_y_l, nxr_y_l
1186
1187	work(0:ny) = ar(0:ny,i,k)
1188
1189	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1190
1191	DO j = 0, (ny+1)/2
1192	ar_tr(j,i,k) = work(2*j)
1193	ENDDO
1194	DO j = 1, (ny+1)/2 - 1
1195	ar_tr(ny+1-j,i,k) = work(2*j+1)
1196	ENDDO
1197
1198	ENDDO
1199	ENDDO
1200	!$END OMP PARALLEL
1201
1202	ELSE
1203
1204	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1205	!$OMP DO
1206	DO k = nzb_y, nzt_y
1207	DO i = nxl_y_l, nxr_y_l
1208
1209	DO j = 0, (ny+1)/2
1210	work(2*j) = ar_tr(j,i,k)
1211	ENDDO
1212	DO j = 1, (ny+1)/2 - 1
1213	work(2*j+1) = ar_tr(ny+1-j,i,k)
1214	ENDDO
1215	work(1) = 0.0
1216	work(ny+2) = 0.0
1217
1218	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1219
1220	ar(0:ny,i,k) = work(0:ny)
1221
1222	ENDDO
1223	ENDDO
1224	!$OMP END PARALLEL
1225
1226	ENDIF
1227	#elif defined( __cuda_fft )
1228
1229	IF ( forward_fft ) THEN
1230
1231	!$acc data present( ar )
1232	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1233
1234	!$acc kernels
1235	!$acc loop
1236	DO k = nzb_y, nzt_y
1237	DO i = nxl_y, nxr_y
1238
1239	!$acc loop vector( 32 )
1240	DO j = 0, (ny+1)/2
1241	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
1242	ENDDO
1243
1244	!$acc loop vector( 32 )
1245	DO j = 1, (ny+1)/2 - 1
1246	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1247	ENDDO
1248
1249	ENDDO
1250	ENDDO
1251	!$acc end kernels
1252	!$acc end data
1253
1254	ELSE
1255
1256	!$acc data present( ar )
1257	!$acc kernels
1258	!$acc loop
1259	DO k = nzb_y, nzt_y
1260	DO i = nxl_y, nxr_y
1261
1262	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1263
1264	!$acc loop vector( 32 )
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 |