Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |