Home

Context Navigation

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

Last change on this file since 2062 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |