Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |