Home

Context Navigation

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

Last change on this file since 1759 was 1750, checked in by raasch, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |