Home

Context Navigation

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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 56.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |