Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |