Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |