Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |