Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |