Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |