Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 1998

Last change on this file since 1998 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |