Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |