Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |