Home

Context Navigation

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

Last change on this file since 1482 was 1482, checked in by raasch, 10 years ago
adjustments for using CUDA-aware MPI
Property svn:keywords set to `Id`
File size: 58.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |