Home

Context Navigation

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

Last change on this file since 1817 was 1816, checked in by raasch, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |