Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |