Home

Context Navigation

fft_xy.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: