Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |