Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |