Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |