Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |