Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |