Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |