Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |