Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |