Home

Context Navigation

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

Last change on this file since 1153 was 1153, checked in by raasch, 11 years ago
code adjustment of accelerator version for PGI 12.3 / CUDA 5.0
Property svn:keywords set to `Id`
File size: 45.9 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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! code adjustment of data types for CUDA fft required by PGI 12.3 / CUDA 5.0
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: fft_xy.f90 1153 2013-05-10 14:33:08Z raasch $
27	!
28	! 1111 2013-03-08 23:54:10Z raasch
29	! further openACC statements added, CUDA branch completely runs on GPU
30	! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
31	! used total domain)
32	!
33	! 1106 2013-03-04 05:31:38Z raasch
34	! CUDA fft added
35	! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
36	! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
37	!
38	! 1092 2013-02-02 11:24:22Z raasch
39	! variable sizw declared for NEC case only
40	!
41	! 1036 2012-10-22 13:43:42Z raasch
42	! code put under GPL (PALM 3.9)
43	!
44	! 274 2009-03-26 15:11:21Z heinze
45	! Output of messages replaced by message handling routine.
46	!
47	! Feb. 2007
48	! RCS Log replace by Id keyword, revision history cleaned up
49	!
50	! Revision 1.4 2006/03/28 12:27:09 raasch
51	! Stop when system-specific fft is selected on NEC. For unknown reasons this
52	! causes a program abort during first allocation in init_grid.
53	!
54	! Revision 1.2 2004/04/30 11:44:27 raasch
55	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
56	! fft_x and fft_y,
57	! function FFT replaced by subroutine FFTN due to problems with 64-bit
58	! mode on ibm,
59	! shape of array cwork is explicitly stored in ishape/jshape and handled
60	! to routine FFTN instead of shape-function (due to compiler error on
61	! decalpha),
62	! non vectorized FFT for nec included
63	!
64	! Revision 1.1 2002/06/11 13:00:49 raasch
65	! Initial revision
66	!
67	!
68	! Description:
69	! ------------
70	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
71	! Original version: Klaus Ketelsen (May 2002)
72	!------------------------------------------------------------------------------!
73
74	USE control_parameters
75	USE indices
76	#if defined( __cuda_fft )
77	USE ISO_C_BINDING
78	#endif
79	USE precision_kind
80	USE singleton
81	USE temperton_fft
82	USE transpose_indices
83
84	IMPLICIT NONE
85
86	PRIVATE
87	PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
88
89	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
90
91	LOGICAL, SAVE :: init_fft = .FALSE.
92
93	REAL, SAVE :: dnx, dny, sqr_dnx, sqr_dny
94	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
95
96	#if defined( __ibm )
97	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
98	!
99	!-- The following working arrays contain tables and have to be "save" and
100	!-- shared in OpenMP sense
101	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
102	#elif defined( __nec )
103	INTEGER, SAVE :: nz1
104	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
105	trig_yf
106	#elif defined( __cuda_fft )
107	INTEGER(C_INT), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
108	INTEGER, SAVE :: total_points_x_transpo, total_points_y_transpo
109	#endif
110
111	!
112	!-- Public interfaces
113	INTERFACE fft_init
114	MODULE PROCEDURE fft_init
115	END INTERFACE fft_init
116
117	INTERFACE fft_x
118	MODULE PROCEDURE fft_x
119	END INTERFACE fft_x
120
121	INTERFACE fft_x_1d
122	MODULE PROCEDURE fft_x_1d
123	END INTERFACE fft_x_1d
124
125	INTERFACE fft_y
126	MODULE PROCEDURE fft_y
127	END INTERFACE fft_y
128
129	INTERFACE fft_y_1d
130	MODULE PROCEDURE fft_y_1d
131	END INTERFACE fft_y_1d
132
133	INTERFACE fft_x_m
134	MODULE PROCEDURE fft_x_m
135	END INTERFACE fft_x_m
136
137	INTERFACE fft_y_m
138	MODULE PROCEDURE fft_y_m
139	END INTERFACE fft_y_m
140
141	CONTAINS
142
143
144	SUBROUTINE fft_init
145
146	USE cuda_fft_interfaces
147
148	IMPLICIT NONE
149
150	!
151	!-- The following temporary working arrays have to be on stack or private
152	!-- in OpenMP sense
153	#if defined( __ibm )
154	REAL, DIMENSION(0:nx+2) :: workx
155	REAL, DIMENSION(0:ny+2) :: worky
156	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
157	#elif defined( __nec )
158	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
159	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
160	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
161	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
162	#endif
163
164	!
165	!-- Return, if already called
166	IF ( init_fft ) THEN
167	RETURN
168	ELSE
169	init_fft = .TRUE.
170	ENDIF
171
172	IF ( fft_method == 'system-specific' ) THEN
173
174	dnx = 1.0 / ( nx + 1.0 )
175	dny = 1.0 / ( ny + 1.0 )
176	sqr_dnx = SQRT( dnx )
177	sqr_dny = SQRT( dny )
178	#if defined( __ibm ) && ! defined( __ibmy_special )
179	!
180	!-- Initialize tables for fft along x
181	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
182	aux2, nau2 )
183	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
184	aux4, nau2 )
185	!
186	!-- Initialize tables for fft along y
187	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
188	auy2, nau2 )
189	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
190	auy4, nau2 )
191	#elif defined( __nec )
192	message_string = 'fft method "' // TRIM( fft_method) // &
193	'" currently does not work on NEC'
194	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
195
196	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
197	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
198
199	work_x = 0.0
200	work_y = 0.0
201	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
202	! when using the NEC ffts
203
204	!
205	!-- Initialize tables for fft along x (non-vector and vector case (M))
206	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
207	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
208	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
209	trig_xf, workx, 0 )
210	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
211	trig_xb, workx, 0 )
212	!
213	!-- Initialize tables for fft along y (non-vector and vector case (M))
214	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
215	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
216	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
217	trig_yf, worky, 0 )
218	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
219	trig_yb, worky, 0 )
220	#elif defined( __cuda_fft )
221	total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
222	total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
223	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
224	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
225	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
226	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
227	#else
228	message_string = 'no system-specific fft-call available'
229	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
230	#endif
231	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
232	!
233	!-- Temperton-algorithm
234	!-- Initialize tables for fft along x and y
235	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
236
237	CALL set99( trigs_x, ifax_x, nx+1 )
238	CALL set99( trigs_y, ifax_y, ny+1 )
239
240	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
241
242	CONTINUE
243
244	ELSE
245
246	message_string = 'fft method "' // TRIM( fft_method) // &
247	'" not available'
248	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
249	ENDIF
250
251	END SUBROUTINE fft_init
252
253
254	SUBROUTINE fft_x( ar, direction )
255
256	!----------------------------------------------------------------------!
257	! fft_x !
258	! !
259	! Fourier-transformation along x-direction !
260	! Version for 2D-decomposition !
261	! !
262	! fft_x uses internal algorithms (Singleton or Temperton) or !
263	! system-specific routines, if they are available !
264	!----------------------------------------------------------------------!
265
266	USE cuda_fft_interfaces
267	#if defined( __cuda_fft )
268	USE ISO_C_BINDING
269	#endif
270
271	IMPLICIT NONE
272
273	CHARACTER (LEN=*) :: direction
274	INTEGER :: i, ishape(1), j, k
275
276	LOGICAL :: forward_fft
277
278	REAL, DIMENSION(0:nx+2) :: work
279	REAL, DIMENSION(nx+2) :: work1
280	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
281	#if defined( __ibm )
282	REAL, DIMENSION(nau2) :: aux2, aux4
283	#elif defined( __nec )
284	REAL, DIMENSION(6*(nx+1)) :: work2
285	#elif defined( __cuda_fft )
286	!$acc declare create( ar_tmp )
287	COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp
288	#endif
289	REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar
290
291	IF ( direction == 'forward' ) THEN
292	forward_fft = .TRUE.
293	ELSE
294	forward_fft = .FALSE.
295	ENDIF
296
297	IF ( fft_method == 'singleton-algorithm' ) THEN
298
299	!
300	!-- Performing the fft with singleton's software works on every system,
301	!-- since it is part of the model
302	ALLOCATE( cwork(0:nx) )
303
304	IF ( forward_fft ) then
305
306	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
307	!$OMP DO
308	DO k = nzb_x, nzt_x
309	DO j = nys_x, nyn_x
310
311	DO i = 0, nx
312	cwork(i) = CMPLX( ar(i,j,k) )
313	ENDDO
314
315	ishape = SHAPE( cwork )
316	CALL FFTN( cwork, ishape )
317
318	DO i = 0, (nx+1)/2
319	ar(i,j,k) = REAL( cwork(i) )
320	ENDDO
321	DO i = 1, (nx+1)/2 - 1
322	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
323	ENDDO
324
325	ENDDO
326	ENDDO
327	!$OMP END PARALLEL
328
329	ELSE
330
331	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
332	!$OMP DO
333	DO k = nzb_x, nzt_x
334	DO j = nys_x, nyn_x
335
336	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
337	DO i = 1, (nx+1)/2 - 1
338	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
339	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
340	ENDDO
341	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
342
343	ishape = SHAPE( cwork )
344	CALL FFTN( cwork, ishape, inv = .TRUE. )
345
346	DO i = 0, nx
347	ar(i,j,k) = REAL( cwork(i) )
348	ENDDO
349
350	ENDDO
351	ENDDO
352	!$OMP END PARALLEL
353
354	ENDIF
355
356	DEALLOCATE( cwork )
357
358	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
359
360	!
361	!-- Performing the fft with Temperton's software works on every system,
362	!-- since it is part of the model
363	IF ( forward_fft ) THEN
364
365	!$OMP PARALLEL PRIVATE ( work, i, j, k )
366	!$OMP DO
367	DO k = nzb_x, nzt_x
368	DO j = nys_x, nyn_x
369
370	work(0:nx) = ar(0:nx,j,k)
371	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
372
373	DO i = 0, (nx+1)/2
374	ar(i,j,k) = work(2*i)
375	ENDDO
376	DO i = 1, (nx+1)/2 - 1
377	ar(nx+1-i,j,k) = work(2*i+1)
378	ENDDO
379
380	ENDDO
381	ENDDO
382	!$OMP END PARALLEL
383
384	ELSE
385
386	!$OMP PARALLEL PRIVATE ( work, i, j, k )
387	!$OMP DO
388	DO k = nzb_x, nzt_x
389	DO j = nys_x, nyn_x
390
391	DO i = 0, (nx+1)/2
392	work(2*i) = ar(i,j,k)
393	ENDDO
394	DO i = 1, (nx+1)/2 - 1
395	work(2*i+1) = ar(nx+1-i,j,k)
396	ENDDO
397	work(1) = 0.0
398	work(nx+2) = 0.0
399
400	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
401	ar(0:nx,j,k) = work(0:nx)
402
403	ENDDO
404	ENDDO
405	!$OMP END PARALLEL
406
407	ENDIF
408
409	ELSEIF ( fft_method == 'system-specific' ) THEN
410
411	#if defined( __ibm ) && ! defined( __ibmy_special )
412	IF ( forward_fft ) THEN
413
414	!$OMP PARALLEL PRIVATE ( work, i, j, k )
415	!$OMP DO
416	DO k = nzb_x, nzt_x
417	DO j = nys_x, nyn_x
418
419	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
420	aux2, nau2 )
421
422	DO i = 0, (nx+1)/2
423	ar(i,j,k) = work(2*i)
424	ENDDO
425	DO i = 1, (nx+1)/2 - 1
426	ar(nx+1-i,j,k) = work(2*i+1)
427	ENDDO
428
429	ENDDO
430	ENDDO
431	!$OMP END PARALLEL
432
433	ELSE
434
435	!$OMP PARALLEL PRIVATE ( work, i, j, k )
436	!$OMP DO
437	DO k = nzb_x, nzt_x
438	DO j = nys_x, nyn_x
439
440	DO i = 0, (nx+1)/2
441	work(2*i) = ar(i,j,k)
442	ENDDO
443	DO i = 1, (nx+1)/2 - 1
444	work(2*i+1) = ar(nx+1-i,j,k)
445	ENDDO
446	work(1) = 0.0
447	work(nx+2) = 0.0
448
449	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
450	aux4, nau2 )
451
452	DO i = 0, nx
453	ar(i,j,k) = work(i)
454	ENDDO
455
456	ENDDO
457	ENDDO
458	!$OMP END PARALLEL
459
460	ENDIF
461
462	#elif defined( __nec )
463
464	IF ( forward_fft ) THEN
465
466	!$OMP PARALLEL PRIVATE ( work, i, j, k )
467	!$OMP DO
468	DO k = nzb_x, nzt_x
469	DO j = nys_x, nyn_x
470
471	work(0:nx) = ar(0:nx,j,k)
472
473	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
474
475	DO i = 0, (nx+1)/2
476	ar(i,j,k) = work(2*i)
477	ENDDO
478	DO i = 1, (nx+1)/2 - 1
479	ar(nx+1-i,j,k) = work(2*i+1)
480	ENDDO
481
482	ENDDO
483	ENDDO
484	!$END OMP PARALLEL
485
486	ELSE
487
488	!$OMP PARALLEL PRIVATE ( work, i, j, k )
489	!$OMP DO
490	DO k = nzb_x, nzt_x
491	DO j = nys_x, nyn_x
492
493	DO i = 0, (nx+1)/2
494	work(2*i) = ar(i,j,k)
495	ENDDO
496	DO i = 1, (nx+1)/2 - 1
497	work(2*i+1) = ar(nx+1-i,j,k)
498	ENDDO
499	work(1) = 0.0
500	work(nx+2) = 0.0
501
502	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
503
504	ar(0:nx,j,k) = work(0:nx)
505
506	ENDDO
507	ENDDO
508	!$OMP END PARALLEL
509
510	ENDIF
511
512	#elif defined( __cuda_fft )
513
514	IF ( forward_fft ) THEN
515
516	!$acc data present( ar )
517	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
518
519	!$acc kernels
520	!$acc loop
521	DO k = nzb_x, nzt_x
522	DO j = nys_x, nyn_x
523
524	!$acc loop vector( 32 )
525	DO i = 0, (nx+1)/2
526	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
527	ENDDO
528
529	!$acc loop vector( 32 )
530	DO i = 1, (nx+1)/2 - 1
531	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
532	ENDDO
533
534	ENDDO
535	ENDDO
536	!$acc end kernels
537	!$acc end data
538
539	ELSE
540
541	!$acc data present( ar )
542	!$acc kernels
543	!$acc loop
544	DO k = nzb_x, nzt_x
545	DO j = nys_x, nyn_x
546
547	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
548
549	!$acc loop vector( 32 )
550	DO i = 1, (nx+1)/2 - 1
551	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
552	ENDDO
553	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
554
555	ENDDO
556	ENDDO
557	!$acc end kernels
558
559	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
560	!$acc end data
561
562	ENDIF
563
564	#else
565	message_string = 'no system-specific fft-call available'
566	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
567	#endif
568
569	ELSE
570
571	message_string = 'fft method "' // TRIM( fft_method) // &
572	'" not available'
573	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
574
575	ENDIF
576
577	END SUBROUTINE fft_x
578
579	SUBROUTINE fft_x_1d( ar, direction )
580
581	!----------------------------------------------------------------------!
582	! fft_x_1d !
583	! !
584	! Fourier-transformation along x-direction !
585	! Version for 1D-decomposition !
586	! !
587	! fft_x uses internal algorithms (Singleton or Temperton) or !
588	! system-specific routines, if they are available !
589	!----------------------------------------------------------------------!
590
591	IMPLICIT NONE
592
593	CHARACTER (LEN=*) :: direction
594	INTEGER :: i, ishape(1)
595
596	LOGICAL :: forward_fft
597
598	REAL, DIMENSION(0:nx) :: ar
599	REAL, DIMENSION(0:nx+2) :: work
600	REAL, DIMENSION(nx+2) :: work1
601	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
602	#if defined( __ibm )
603	REAL, DIMENSION(nau2) :: aux2, aux4
604	#elif defined( __nec )
605	REAL, DIMENSION(6*(nx+1)) :: work2
606	#endif
607
608	IF ( direction == 'forward' ) THEN
609	forward_fft = .TRUE.
610	ELSE
611	forward_fft = .FALSE.
612	ENDIF
613
614	IF ( fft_method == 'singleton-algorithm' ) THEN
615
616	!
617	!-- Performing the fft with singleton's software works on every system,
618	!-- since it is part of the model
619	ALLOCATE( cwork(0:nx) )
620
621	IF ( forward_fft ) then
622
623	DO i = 0, nx
624	cwork(i) = CMPLX( ar(i) )
625	ENDDO
626	ishape = SHAPE( cwork )
627	CALL FFTN( cwork, ishape )
628	DO i = 0, (nx+1)/2
629	ar(i) = REAL( cwork(i) )
630	ENDDO
631	DO i = 1, (nx+1)/2 - 1
632	ar(nx+1-i) = -AIMAG( cwork(i) )
633	ENDDO
634
635	ELSE
636
637	cwork(0) = CMPLX( ar(0), 0.0 )
638	DO i = 1, (nx+1)/2 - 1
639	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
640	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
641	ENDDO
642	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
643
644	ishape = SHAPE( cwork )
645	CALL FFTN( cwork, ishape, inv = .TRUE. )
646
647	DO i = 0, nx
648	ar(i) = REAL( cwork(i) )
649	ENDDO
650
651	ENDIF
652
653	DEALLOCATE( cwork )
654
655	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
656
657	!
658	!-- Performing the fft with Temperton's software works on every system,
659	!-- since it is part of the model
660	IF ( forward_fft ) THEN
661
662	work(0:nx) = ar
663	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
664
665	DO i = 0, (nx+1)/2
666	ar(i) = work(2*i)
667	ENDDO
668	DO i = 1, (nx+1)/2 - 1
669	ar(nx+1-i) = work(2*i+1)
670	ENDDO
671
672	ELSE
673
674	DO i = 0, (nx+1)/2
675	work(2*i) = ar(i)
676	ENDDO
677	DO i = 1, (nx+1)/2 - 1
678	work(2*i+1) = ar(nx+1-i)
679	ENDDO
680	work(1) = 0.0
681	work(nx+2) = 0.0
682
683	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
684	ar = work(0:nx)
685
686	ENDIF
687
688	ELSEIF ( fft_method == 'system-specific' ) THEN
689
690	#if defined( __ibm ) && ! defined( __ibmy_special )
691	IF ( forward_fft ) THEN
692
693	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
694	aux2, nau2 )
695
696	DO i = 0, (nx+1)/2
697	ar(i) = work(2*i)
698	ENDDO
699	DO i = 1, (nx+1)/2 - 1
700	ar(nx+1-i) = work(2*i+1)
701	ENDDO
702
703	ELSE
704
705	DO i = 0, (nx+1)/2
706	work(2*i) = ar(i)
707	ENDDO
708	DO i = 1, (nx+1)/2 - 1
709	work(2*i+1) = ar(nx+1-i)
710	ENDDO
711	work(1) = 0.0
712	work(nx+2) = 0.0
713
714	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
715	aux4, nau2 )
716
717	DO i = 0, nx
718	ar(i) = work(i)
719	ENDDO
720
721	ENDIF
722	#elif defined( __nec )
723	IF ( forward_fft ) THEN
724
725	work(0:nx) = ar(0:nx)
726
727	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
728
729	DO i = 0, (nx+1)/2
730	ar(i) = work(2*i)
731	ENDDO
732	DO i = 1, (nx+1)/2 - 1
733	ar(nx+1-i) = work(2*i+1)
734	ENDDO
735
736	ELSE
737
738	DO i = 0, (nx+1)/2
739	work(2*i) = ar(i)
740	ENDDO
741	DO i = 1, (nx+1)/2 - 1
742	work(2*i+1) = ar(nx+1-i)
743	ENDDO
744	work(1) = 0.0
745	work(nx+2) = 0.0
746
747	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
748
749	ar(0:nx) = work(0:nx)
750
751	ENDIF
752	#else
753	message_string = 'no system-specific fft-call available'
754	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
755	#endif
756	ELSE
757	message_string = 'fft method "' // TRIM( fft_method) // &
758	'" not available'
759	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
760
761	ENDIF
762
763	END SUBROUTINE fft_x_1d
764
765	SUBROUTINE fft_y( ar, direction )
766
767	!----------------------------------------------------------------------!
768	! fft_y !
769	! !
770	! Fourier-transformation along y-direction !
771	! Version for 2D-decomposition !
772	! !
773	! fft_y uses internal algorithms (Singleton or Temperton) or !
774	! system-specific routines, if they are available !
775	!----------------------------------------------------------------------!
776
777	USE cuda_fft_interfaces
778	#if defined( __cuda_fft )
779	USE ISO_C_BINDING
780	#endif
781
782	IMPLICIT NONE
783
784	CHARACTER (LEN=*) :: direction
785	INTEGER :: i, j, jshape(1), k
786
787	LOGICAL :: forward_fft
788
789	REAL, DIMENSION(0:ny+2) :: work
790	REAL, DIMENSION(ny+2) :: work1
791	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
792	#if defined( __ibm )
793	REAL, DIMENSION(nau2) :: auy2, auy4
794	#elif defined( __nec )
795	REAL, DIMENSION(6*(ny+1)) :: work2
796	#elif defined( __cuda_fft )
797	!$acc declare create( ar_tmp )
798	COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
799	#endif
800	REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) :: ar
801
802	IF ( direction == 'forward' ) THEN
803	forward_fft = .TRUE.
804	ELSE
805	forward_fft = .FALSE.
806	ENDIF
807
808	IF ( fft_method == 'singleton-algorithm' ) THEN
809
810	!
811	!-- Performing the fft with singleton's software works on every system,
812	!-- since it is part of the model
813	ALLOCATE( cwork(0:ny) )
814
815	IF ( forward_fft ) then
816
817	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
818	!$OMP DO
819	DO k = nzb_y, nzt_y
820	DO i = nxl_y, nxr_y
821
822	DO j = 0, ny
823	cwork(j) = CMPLX( ar(j,i,k) )
824	ENDDO
825
826	jshape = SHAPE( cwork )
827	CALL FFTN( cwork, jshape )
828
829	DO j = 0, (ny+1)/2
830	ar(j,i,k) = REAL( cwork(j) )
831	ENDDO
832	DO j = 1, (ny+1)/2 - 1
833	ar(ny+1-j,i,k) = -AIMAG( cwork(j) )
834	ENDDO
835
836	ENDDO
837	ENDDO
838	!$OMP END PARALLEL
839
840	ELSE
841
842	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
843	!$OMP DO
844	DO k = nzb_y, nzt_y
845	DO i = nxl_y, nxr_y
846
847	cwork(0) = CMPLX( ar(0,i,k), 0.0 )
848	DO j = 1, (ny+1)/2 - 1
849	cwork(j) = CMPLX( ar(j,i,k), -ar(ny+1-j,i,k) )
850	cwork(ny+1-j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
851	ENDDO
852	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
853
854	jshape = SHAPE( cwork )
855	CALL FFTN( cwork, jshape, inv = .TRUE. )
856
857	DO j = 0, ny
858	ar(j,i,k) = REAL( cwork(j) )
859	ENDDO
860
861	ENDDO
862	ENDDO
863	!$OMP END PARALLEL
864
865	ENDIF
866
867	DEALLOCATE( cwork )
868
869	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
870
871	!
872	!-- Performing the fft with Temperton's software works on every system,
873	!-- since it is part of the model
874	IF ( forward_fft ) THEN
875
876	!$OMP PARALLEL PRIVATE ( work, i, j, k )
877	!$OMP DO
878	DO k = nzb_y, nzt_y
879	DO i = nxl_y, nxr_y
880
881	work(0:ny) = ar(0:ny,i,k)
882	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
883
884	DO j = 0, (ny+1)/2
885	ar(j,i,k) = work(2*j)
886	ENDDO
887	DO j = 1, (ny+1)/2 - 1
888	ar(ny+1-j,i,k) = work(2*j+1)
889	ENDDO
890
891	ENDDO
892	ENDDO
893	!$OMP END PARALLEL
894
895	ELSE
896
897	!$OMP PARALLEL PRIVATE ( work, i, j, k )
898	!$OMP DO
899	DO k = nzb_y, nzt_y
900	DO i = nxl_y, nxr_y
901
902	DO j = 0, (ny+1)/2
903	work(2*j) = ar(j,i,k)
904	ENDDO
905	DO j = 1, (ny+1)/2 - 1
906	work(2*j+1) = ar(ny+1-j,i,k)
907	ENDDO
908	work(1) = 0.0
909	work(ny+2) = 0.0
910
911	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
912	ar(0:ny,i,k) = work(0:ny)
913
914	ENDDO
915	ENDDO
916	!$OMP END PARALLEL
917
918	ENDIF
919
920	ELSEIF ( fft_method == 'system-specific' ) THEN
921
922	#if defined( __ibm ) && ! defined( __ibmy_special )
923	IF ( forward_fft) THEN
924
925	!$OMP PARALLEL PRIVATE ( work, i, j, k )
926	!$OMP DO
927	DO k = nzb_y, nzt_y
928	DO i = nxl_y, nxr_y
929
930	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
931	auy2, nau2 )
932
933	DO j = 0, (ny+1)/2
934	ar(j,i,k) = work(2*j)
935	ENDDO
936	DO j = 1, (ny+1)/2 - 1
937	ar(ny+1-j,i,k) = work(2*j+1)
938	ENDDO
939
940	ENDDO
941	ENDDO
942	!$OMP END PARALLEL
943
944	ELSE
945
946	!$OMP PARALLEL PRIVATE ( work, i, j, k )
947	!$OMP DO
948	DO k = nzb_y, nzt_y
949	DO i = nxl_y, nxr_y
950
951	DO j = 0, (ny+1)/2
952	work(2*j) = ar(j,i,k)
953	ENDDO
954	DO j = 1, (ny+1)/2 - 1
955	work(2*j+1) = ar(ny+1-j,i,k)
956	ENDDO
957	work(1) = 0.0
958	work(ny+2) = 0.0
959
960	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
961	auy4, nau2 )
962
963	DO j = 0, ny
964	ar(j,i,k) = work(j)
965	ENDDO
966
967	ENDDO
968	ENDDO
969	!$OMP END PARALLEL
970
971	ENDIF
972	#elif defined( __nec )
973	IF ( forward_fft ) THEN
974
975	!$OMP PARALLEL PRIVATE ( work, i, j, k )
976	!$OMP DO
977	DO k = nzb_y, nzt_y
978	DO i = nxl_y, nxr_y
979
980	work(0:ny) = ar(0:ny,i,k)
981
982	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
983
984	DO j = 0, (ny+1)/2
985	ar(j,i,k) = work(2*j)
986	ENDDO
987	DO j = 1, (ny+1)/2 - 1
988	ar(ny+1-j,i,k) = work(2*j+1)
989	ENDDO
990
991	ENDDO
992	ENDDO
993	!$END OMP PARALLEL
994
995	ELSE
996
997	!$OMP PARALLEL PRIVATE ( work, i, j, k )
998	!$OMP DO
999	DO k = nzb_y, nzt_y
1000	DO i = nxl_y, nxr_y
1001
1002	DO j = 0, (ny+1)/2
1003	work(2*j) = ar(j,i,k)
1004	ENDDO
1005	DO j = 1, (ny+1)/2 - 1
1006	work(2*j+1) = ar(ny+1-j,i,k)
1007	ENDDO
1008	work(1) = 0.0
1009	work(ny+2) = 0.0
1010
1011	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1012
1013	ar(0:ny,i,k) = work(0:ny)
1014
1015	ENDDO
1016	ENDDO
1017	!$OMP END PARALLEL
1018
1019	ENDIF
1020	#elif defined( __cuda_fft )
1021
1022	IF ( forward_fft ) THEN
1023
1024	!$acc data present( ar )
1025	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1026
1027	!$acc kernels
1028	!$acc loop
1029	DO k = nzb_y, nzt_y
1030	DO i = nxl_y, nxr_y
1031
1032	!$acc loop vector( 32 )
1033	DO j = 0, (ny+1)/2
1034	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
1035	ENDDO
1036
1037	!$acc loop vector( 32 )
1038	DO j = 1, (ny+1)/2 - 1
1039	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1040	ENDDO
1041
1042	ENDDO
1043	ENDDO
1044	!$acc end kernels
1045	!$acc end data
1046
1047	ELSE
1048
1049	!$acc data present( ar )
1050	!$acc kernels
1051	!$acc loop
1052	DO k = nzb_y, nzt_y
1053	DO i = nxl_y, nxr_y
1054
1055	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1056
1057	!$acc loop vector( 32 )
1058	DO j = 1, (ny+1)/2 - 1
1059	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1060	ENDDO
1061	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1062
1063	ENDDO
1064	ENDDO
1065	!$acc end kernels
1066
1067	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1068	!$acc end data
1069
1070	ENDIF
1071
1072	#else
1073	message_string = 'no system-specific fft-call available'
1074	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
1075	#endif
1076
1077	ELSE
1078
1079	message_string = 'fft method "' // TRIM( fft_method) // &
1080	'" not available'
1081	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
1082
1083	ENDIF
1084
1085	END SUBROUTINE fft_y
1086
1087	SUBROUTINE fft_y_1d( ar, direction )
1088
1089	!----------------------------------------------------------------------!
1090	! fft_y_1d !
1091	! !
1092	! Fourier-transformation along y-direction !
1093	! Version for 1D-decomposition !
1094	! !
1095	! fft_y uses internal algorithms (Singleton or Temperton) or !
1096	! system-specific routines, if they are available !
1097	!----------------------------------------------------------------------!
1098
1099	IMPLICIT NONE
1100
1101	CHARACTER (LEN=*) :: direction
1102	INTEGER :: j, jshape(1)
1103
1104	LOGICAL :: forward_fft
1105
1106	REAL, DIMENSION(0:ny) :: ar
1107	REAL, DIMENSION(0:ny+2) :: work
1108	REAL, DIMENSION(ny+2) :: work1
1109	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
1110	#if defined( __ibm )
1111	REAL, DIMENSION(nau2) :: auy2, auy4
1112	#elif defined( __nec )
1113	REAL, DIMENSION(6*(ny+1)) :: work2
1114	#endif
1115
1116	IF ( direction == 'forward' ) THEN
1117	forward_fft = .TRUE.
1118	ELSE
1119	forward_fft = .FALSE.
1120	ENDIF
1121
1122	IF ( fft_method == 'singleton-algorithm' ) THEN
1123
1124	!
1125	!-- Performing the fft with singleton's software works on every system,
1126	!-- since it is part of the model
1127	ALLOCATE( cwork(0:ny) )
1128
1129	IF ( forward_fft ) THEN
1130
1131	DO j = 0, ny
1132	cwork(j) = CMPLX( ar(j) )
1133	ENDDO
1134
1135	jshape = SHAPE( cwork )
1136	CALL FFTN( cwork, jshape )
1137
1138	DO j = 0, (ny+1)/2
1139	ar(j) = REAL( cwork(j) )
1140	ENDDO
1141	DO j = 1, (ny+1)/2 - 1
1142	ar(ny+1-j) = -AIMAG( cwork(j) )
1143	ENDDO
1144
1145	ELSE
1146
1147	cwork(0) = CMPLX( ar(0), 0.0 )
1148	DO j = 1, (ny+1)/2 - 1
1149	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
1150	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
1151	ENDDO
1152	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
1153
1154	jshape = SHAPE( cwork )
1155	CALL FFTN( cwork, jshape, inv = .TRUE. )
1156
1157	DO j = 0, ny
1158	ar(j) = REAL( cwork(j) )
1159	ENDDO
1160
1161	ENDIF
1162
1163	DEALLOCATE( cwork )
1164
1165	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1166
1167	!
1168	!-- Performing the fft with Temperton's software works on every system,
1169	!-- since it is part of the model
1170	IF ( forward_fft ) THEN
1171
1172	work(0:ny) = ar
1173	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1174
1175	DO j = 0, (ny+1)/2
1176	ar(j) = work(2*j)
1177	ENDDO
1178	DO j = 1, (ny+1)/2 - 1
1179	ar(ny+1-j) = work(2*j+1)
1180	ENDDO
1181
1182	ELSE
1183
1184	DO j = 0, (ny+1)/2
1185	work(2*j) = ar(j)
1186	ENDDO
1187	DO j = 1, (ny+1)/2 - 1
1188	work(2*j+1) = ar(ny+1-j)
1189	ENDDO
1190	work(1) = 0.0
1191	work(ny+2) = 0.0
1192
1193	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1194	ar = work(0:ny)
1195
1196	ENDIF
1197
1198	ELSEIF ( fft_method == 'system-specific' ) THEN
1199
1200	#if defined( __ibm ) && ! defined( __ibmy_special )
1201	IF ( forward_fft ) THEN
1202
1203	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1204	auy2, nau2 )
1205
1206	DO j = 0, (ny+1)/2
1207	ar(j) = work(2*j)
1208	ENDDO
1209	DO j = 1, (ny+1)/2 - 1
1210	ar(ny+1-j) = work(2*j+1)
1211	ENDDO
1212
1213	ELSE
1214
1215	DO j = 0, (ny+1)/2
1216	work(2*j) = ar(j)
1217	ENDDO
1218	DO j = 1, (ny+1)/2 - 1
1219	work(2*j+1) = ar(ny+1-j)
1220	ENDDO
1221	work(1) = 0.0
1222	work(ny+2) = 0.0
1223
1224	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1225	auy4, nau2 )
1226
1227	DO j = 0, ny
1228	ar(j) = work(j)
1229	ENDDO
1230
1231	ENDIF
1232	#elif defined( __nec )
1233	IF ( forward_fft ) THEN
1234
1235	work(0:ny) = ar(0:ny)
1236
1237	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1238
1239	DO j = 0, (ny+1)/2
1240	ar(j) = work(2*j)
1241	ENDDO
1242	DO j = 1, (ny+1)/2 - 1
1243	ar(ny+1-j) = work(2*j+1)
1244	ENDDO
1245
1246	ELSE
1247
1248	DO j = 0, (ny+1)/2
1249	work(2*j) = ar(j)
1250	ENDDO
1251	DO j = 1, (ny+1)/2 - 1
1252	work(2*j+1) = ar(ny+1-j)
1253	ENDDO
1254	work(1) = 0.0
1255	work(ny+2) = 0.0
1256
1257	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1258
1259	ar(0:ny) = work(0:ny)
1260
1261	ENDIF
1262	#else
1263	message_string = 'no system-specific fft-call available'
1264	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
1265
1266	#endif
1267
1268	ELSE
1269
1270	message_string = 'fft method "' // TRIM( fft_method) // &
1271	'" not available'
1272	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1273
1274	ENDIF
1275
1276	END SUBROUTINE fft_y_1d
1277
1278	SUBROUTINE fft_x_m( ar, direction )
1279
1280	!----------------------------------------------------------------------!
1281	! fft_x_m !
1282	! !
1283	! Fourier-transformation along x-direction !
1284	! Version for 1d domain decomposition !
1285	! using multiple 1D FFT from Math Keisan on NEC !
1286	! or Temperton-algorithm !
1287	! (no singleton-algorithm on NEC because it does not vectorize) !
1288	! !
1289	!----------------------------------------------------------------------!
1290
1291	IMPLICIT NONE
1292
1293	CHARACTER (LEN=*) :: direction
1294	INTEGER :: i, k, siza
1295
1296	REAL, DIMENSION(0:nx,nz) :: ar
1297	REAL, DIMENSION(0:nx+3,nz+1) :: ai
1298	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
1299	#if defined( __nec )
1300	INTEGER :: sizw
1301	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
1302	#endif
1303
1304	IF ( fft_method == 'temperton-algorithm' ) THEN
1305
1306	siza = SIZE( ai, 1 )
1307
1308	IF ( direction == 'forward') THEN
1309
1310	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1311	ai(nx+1:,:) = 0.0
1312
1313	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1314
1315	DO k = 1, nz
1316	DO i = 0, (nx+1)/2
1317	ar(i,k) = ai(2*i,k)
1318	ENDDO
1319	DO i = 1, (nx+1)/2 - 1
1320	ar(nx+1-i,k) = ai(2*i+1,k)
1321	ENDDO
1322	ENDDO
1323
1324	ELSE
1325
1326	DO k = 1, nz
1327	DO i = 0, (nx+1)/2
1328	ai(2*i,k) = ar(i,k)
1329	ENDDO
1330	DO i = 1, (nx+1)/2 - 1
1331	ai(2*i+1,k) = ar(nx+1-i,k)
1332	ENDDO
1333	ai(1,k) = 0.0
1334	ai(nx+2,k) = 0.0
1335	ENDDO
1336
1337	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1338
1339	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1340
1341	ENDIF
1342
1343	ELSEIF ( fft_method == 'system-specific' ) THEN
1344
1345	#if defined( __nec )
1346	siza = SIZE( ai, 1 )
1347	sizw = SIZE( work, 1 )
1348
1349	IF ( direction == 'forward') THEN
1350
1351	!
1352	!-- Tables are initialized once more. This call should not be
1353	!-- necessary, but otherwise program aborts in asymmetric case
1354	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1355	trig_xf, work1, 0 )
1356
1357	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1358	IF ( nz1 > nz ) THEN
1359	ai(:,nz1) = 0.0
1360	ENDIF
1361
1362	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
1363	trig_xf, work1, 0 )
1364
1365	DO k = 1, nz
1366	DO i = 0, (nx+1)/2
1367	ar(i,k) = REAL( work(i+1,k) )
1368	ENDDO
1369	DO i = 1, (nx+1)/2 - 1
1370	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1371	ENDDO
1372	ENDDO
1373
1374	ELSE
1375
1376	!
1377	!-- Tables are initialized once more. This call should not be
1378	!-- necessary, but otherwise program aborts in asymmetric case
1379	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1380	trig_xb, work1, 0 )
1381
1382	IF ( nz1 > nz ) THEN
1383	work(:,nz1) = 0.0
1384	ENDIF
1385	DO k = 1, nz
1386	work(1,k) = CMPLX( ar(0,k), 0.0 )
1387	DO i = 1, (nx+1)/2 - 1
1388	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
1389	ENDDO
1390	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
1391	ENDDO
1392
1393	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1394	trig_xb, work1, 0 )
1395
1396	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1397
1398	ENDIF
1399
1400	#else
1401	message_string = 'no system-specific fft-call available'
1402	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
1403	#endif
1404
1405	ELSE
1406
1407	message_string = 'fft method "' // TRIM( fft_method) // &
1408	'" not available'
1409	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1410
1411	ENDIF
1412
1413	END SUBROUTINE fft_x_m
1414
1415	SUBROUTINE fft_y_m( ar, ny1, direction )
1416
1417	!----------------------------------------------------------------------!
1418	! fft_y_m !
1419	! !
1420	! Fourier-transformation along y-direction !
1421	! Version for 1d domain decomposition !
1422	! using multiple 1D FFT from Math Keisan on NEC !
1423	! or Temperton-algorithm !
1424	! (no singleton-algorithm on NEC because it does not vectorize) !
1425	! !
1426	!----------------------------------------------------------------------!
1427
1428	IMPLICIT NONE
1429
1430	CHARACTER (LEN=*) :: direction
1431	INTEGER :: j, k, ny1, siza
1432
1433	REAL, DIMENSION(0:ny1,nz) :: ar
1434	REAL, DIMENSION(0:ny+3,nz+1) :: ai
1435	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
1436	#if defined( __nec )
1437	INTEGER :: sizw
1438	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
1439	#endif
1440
1441	IF ( fft_method == 'temperton-algorithm' ) THEN
1442
1443	siza = SIZE( ai, 1 )
1444
1445	IF ( direction == 'forward') THEN
1446
1447	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1448	ai(ny+1:,:) = 0.0
1449
1450	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1451
1452	DO k = 1, nz
1453	DO j = 0, (ny+1)/2
1454	ar(j,k) = ai(2*j,k)
1455	ENDDO
1456	DO j = 1, (ny+1)/2 - 1
1457	ar(ny+1-j,k) = ai(2*j+1,k)
1458	ENDDO
1459	ENDDO
1460
1461	ELSE
1462
1463	DO k = 1, nz
1464	DO j = 0, (ny+1)/2
1465	ai(2*j,k) = ar(j,k)
1466	ENDDO
1467	DO j = 1, (ny+1)/2 - 1
1468	ai(2*j+1,k) = ar(ny+1-j,k)
1469	ENDDO
1470	ai(1,k) = 0.0
1471	ai(ny+2,k) = 0.0
1472	ENDDO
1473
1474	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1475
1476	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1477
1478	ENDIF
1479
1480	ELSEIF ( fft_method == 'system-specific' ) THEN
1481
1482	#if defined( __nec )
1483	siza = SIZE( ai, 1 )
1484	sizw = SIZE( work, 1 )
1485
1486	IF ( direction == 'forward') THEN
1487
1488	!
1489	!-- Tables are initialized once more. This call should not be
1490	!-- necessary, but otherwise program aborts in asymmetric case
1491	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1492	trig_yf, work1, 0 )
1493
1494	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1495	IF ( nz1 > nz ) THEN
1496	ai(:,nz1) = 0.0
1497	ENDIF
1498
1499	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
1500	trig_yf, work1, 0 )
1501
1502	DO k = 1, nz
1503	DO j = 0, (ny+1)/2
1504	ar(j,k) = REAL( work(j+1,k) )
1505	ENDDO
1506	DO j = 1, (ny+1)/2 - 1
1507	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1508	ENDDO
1509	ENDDO
1510
1511	ELSE
1512
1513	!
1514	!-- Tables are initialized once more. This call should not be
1515	!-- necessary, but otherwise program aborts in asymmetric case
1516	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1517	trig_yb, work1, 0 )
1518
1519	IF ( nz1 > nz ) THEN
1520	work(:,nz1) = 0.0
1521	ENDIF
1522	DO k = 1, nz
1523	work(1,k) = CMPLX( ar(0,k), 0.0 )
1524	DO j = 1, (ny+1)/2 - 1
1525	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
1526	ENDDO
1527	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
1528	ENDDO
1529
1530	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1531	trig_yb, work1, 0 )
1532
1533	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1534
1535	ENDIF
1536
1537	#else
1538	message_string = 'no system-specific fft-call available'
1539	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
1540	#endif
1541
1542	ELSE
1543
1544	message_string = 'fft method "' // TRIM( fft_method) // &
1545	'" not available'
1546	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1547
1548	ENDIF
1549
1550	END SUBROUTINE fft_y_m
1551
1552
1553	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |