Home

Context Navigation

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

Last change on this file since 1127 was 1112, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 45.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |