Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |