Home

Context Navigation

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

Last change on this file since 1099 was 1093, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 25.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 1093 2013-02-02 12:58:49Z raasch $
27	!
28	! 1092 2013-02-02 11:24:22Z raasch
29	! variable sizw declared for NEC case only
30	!
31	! 1036 2012-10-22 13:43:42Z raasch
32	! code put under GPL (PALM 3.9)
33	!
34	! 274 2009-03-26 15:11:21Z heinze
35	! Output of messages replaced by message handling routine.
36	!
37	! Feb. 2007
38	! RCS Log replace by Id keyword, revision history cleaned up
39	!
40	! Revision 1.4 2006/03/28 12:27:09 raasch
41	! Stop when system-specific fft is selected on NEC. For unknown reasons this
42	! causes a program abort during first allocation in init_grid.
43	!
44	! Revision 1.2 2004/04/30 11:44:27 raasch
45	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
46	! fft_x and fft_y,
47	! function FFT replaced by subroutine FFTN due to problems with 64-bit
48	! mode on ibm,
49	! shape of array cwork is explicitly stored in ishape/jshape and handled
50	! to routine FFTN instead of shape-function (due to compiler error on
51	! decalpha),
52	! non vectorized FFT for nec included
53	!
54	! Revision 1.1 2002/06/11 13:00:49 raasch
55	! Initial revision
56	!
57	!
58	! Description:
59	! ------------
60	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
61	! Original version: Klaus Ketelsen (May 2002)
62	!------------------------------------------------------------------------------!
63
64	USE array_kind
65	USE control_parameters
66	USE indices
67	USE singleton
68	USE temperton_fft
69
70	IMPLICIT NONE
71
72	PRIVATE
73	PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
74
75	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
76
77	LOGICAL, SAVE :: init_fft = .FALSE.
78
79	REAL, SAVE :: sqr_nx, sqr_ny
80	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
81
82	#if defined( __ibm )
83	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
84	!
85	!-- The following working arrays contain tables and have to be "save" and
86	!-- shared in OpenMP sense
87	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
88	#elif defined( __nec )
89	INTEGER, SAVE :: nz1
90	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
91	trig_yf
92	#endif
93
94	!
95	!-- Public interfaces
96	INTERFACE fft_init
97	MODULE PROCEDURE fft_init
98	END INTERFACE fft_init
99
100	INTERFACE fft_x
101	MODULE PROCEDURE fft_x
102	END INTERFACE fft_x
103
104	INTERFACE fft_y
105	MODULE PROCEDURE fft_y
106	END INTERFACE fft_y
107
108	INTERFACE fft_x_m
109	MODULE PROCEDURE fft_x_m
110	END INTERFACE fft_x_m
111
112	INTERFACE fft_y_m
113	MODULE PROCEDURE fft_y_m
114	END INTERFACE fft_y_m
115
116	CONTAINS
117
118
119	SUBROUTINE fft_init
120
121	IMPLICIT NONE
122
123	!
124	!-- The following temporary working arrays have to be on stack or private
125	!-- in OpenMP sense
126	#if defined( __ibm )
127	REAL, DIMENSION(0:nx+2) :: workx
128	REAL, DIMENSION(0:ny+2) :: worky
129	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
130	#elif defined( __nec )
131	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
132	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
133	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
134	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
135	#endif
136
137	!
138	!-- Return, if already called
139	IF ( init_fft ) THEN
140	RETURN
141	ELSE
142	init_fft = .TRUE.
143	ENDIF
144
145	IF ( fft_method == 'system-specific' ) THEN
146
147	sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
148	sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
149	#if defined( __ibm ) && ! defined( __ibmy_special )
150	!
151	!-- Initialize tables for fft along x
152	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
153	aux2, nau2 )
154	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
155	aux4, nau2 )
156	!
157	!-- Initialize tables for fft along y
158	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
159	auy2, nau2 )
160	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
161	auy4, nau2 )
162	#elif defined( __nec )
163	message_string = 'fft method "' // TRIM( fft_method) // &
164	'" currently does not work on NEC'
165	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
166
167	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
168	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
169
170	work_x = 0.0
171	work_y = 0.0
172	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
173	! when using the NEC ffts
174
175	!
176	!-- Initialize tables for fft along x (non-vector and vector case (M))
177	CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
178	CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
179	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
180	trig_xf, workx, 0 )
181	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
182	trig_xb, workx, 0 )
183	!
184	!-- Initialize tables for fft along y (non-vector and vector case (M))
185	CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
186	CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
187	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
188	trig_yf, worky, 0 )
189	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
190	trig_yb, worky, 0 )
191	#else
192	message_string = 'no system-specific fft-call available'
193	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
194	#endif
195	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
196	!
197	!-- Temperton-algorithm
198	!-- Initialize tables for fft along x and y
199	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
200
201	CALL set99( trigs_x, ifax_x, nx+1 )
202	CALL set99( trigs_y, ifax_y, ny+1 )
203
204	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
205
206	CONTINUE
207
208	ELSE
209
210	message_string = 'fft method "' // TRIM( fft_method) // &
211	'" not available'
212	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
213	ENDIF
214
215	END SUBROUTINE fft_init
216
217
218	SUBROUTINE fft_x( ar, direction )
219
220	!----------------------------------------------------------------------!
221	! fft_x !
222	! !
223	! Fourier-transformation along x-direction !
224	! !
225	! fft_x uses internal algorithms (Singleton or Temperton) or !
226	! system-specific routines, if they are available !
227	!----------------------------------------------------------------------!
228
229	IMPLICIT NONE
230
231	CHARACTER (LEN=*) :: direction
232	INTEGER :: i, ishape(1)
233
234	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
235	REAL, DIMENSION(0:nx) :: ar
236	REAL, DIMENSION(0:nx+2) :: work
237	REAL, DIMENSION(nx+2) :: work1
238	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
239	#if defined( __ibm )
240	REAL, DIMENSION(nau2) :: aux2, aux4
241	#elif defined( __nec )
242	REAL, DIMENSION(6*(nx+1)) :: work2
243	#endif
244
245	IF ( fft_method == 'singleton-algorithm' ) THEN
246
247	!
248	!-- Performing the fft with singleton's software works on every system,
249	!-- since it is part of the model
250	ALLOCATE( cwork(0:nx) )
251
252	IF ( direction == 'forward') then
253
254	DO i = 0, nx
255	cwork(i) = CMPLX( ar(i) )
256	ENDDO
257	ishape = SHAPE( cwork )
258	CALL FFTN( cwork, ishape )
259
260	DO i = 0, (nx+1)/2
261	ar(i) = REAL( cwork(i) )
262	ENDDO
263	DO i = 1, (nx+1)/2 - 1
264	ar(nx+1-i) = -AIMAG( cwork(i) )
265	ENDDO
266
267	ELSE
268
269	cwork(0) = CMPLX( ar(0), 0.0 )
270	DO i = 1, (nx+1)/2 - 1
271	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
272	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
273	ENDDO
274	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
275
276	ishape = SHAPE( cwork )
277	CALL FFTN( cwork, ishape, inv = .TRUE. )
278
279	DO i = 0, nx
280	ar(i) = REAL( cwork(i) )
281	ENDDO
282
283	ENDIF
284
285	DEALLOCATE( cwork )
286
287	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
288
289	!
290	!-- Performing the fft with Temperton's software works on every system,
291	!-- since it is part of the model
292	IF ( direction == 'forward' ) THEN
293
294	work(0:nx) = ar
295	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
296
297	DO i = 0, (nx+1)/2
298	ar(i) = work(2*i)
299	ENDDO
300	DO i = 1, (nx+1)/2 - 1
301	ar(nx+1-i) = work(2*i+1)
302	ENDDO
303
304	ELSE
305
306	DO i = 0, (nx+1)/2
307	work(2*i) = ar(i)
308	ENDDO
309	DO i = 1, (nx+1)/2 - 1
310	work(2*i+1) = ar(nx+1-i)
311	ENDDO
312	work(1) = 0.0
313	work(nx+2) = 0.0
314
315	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
316	ar = work(0:nx)
317
318	ENDIF
319
320	ELSEIF ( fft_method == 'system-specific' ) THEN
321
322	#if defined( __ibm ) && ! defined( __ibmy_special )
323	IF ( direction == 'forward' ) THEN
324
325	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
326	aux2, nau2 )
327
328	DO i = 0, (nx+1)/2
329	ar(i) = work(2*i)
330	ENDDO
331	DO i = 1, (nx+1)/2 - 1
332	ar(nx+1-i) = work(2*i+1)
333	ENDDO
334
335	ELSE
336
337	DO i = 0, (nx+1)/2
338	work(2*i) = ar(i)
339	ENDDO
340	DO i = 1, (nx+1)/2 - 1
341	work(2*i+1) = ar(nx+1-i)
342	ENDDO
343	work(1) = 0.0
344	work(nx+2) = 0.0
345
346	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
347	aux4, nau2 )
348
349	DO i = 0, nx
350	ar(i) = work(i)
351	ENDDO
352
353	ENDIF
354	#elif defined( __nec )
355	IF ( direction == 'forward' ) THEN
356
357	work(0:nx) = ar(0:nx)
358
359	CALL DZFFT( 1, nx+1, sqr_nx, work, work, trig_xf, work2, 0 )
360
361	DO i = 0, (nx+1)/2
362	ar(i) = work(2*i)
363	ENDDO
364	DO i = 1, (nx+1)/2 - 1
365	ar(nx+1-i) = work(2*i+1)
366	ENDDO
367
368	ELSE
369
370	DO i = 0, (nx+1)/2
371	work(2*i) = ar(i)
372	ENDDO
373	DO i = 1, (nx+1)/2 - 1
374	work(2*i+1) = ar(nx+1-i)
375	ENDDO
376	work(1) = 0.0
377	work(nx+2) = 0.0
378
379	CALL ZDFFT( -1, nx+1, sqr_nx, work, work, trig_xb, work2, 0 )
380
381	ar(0:nx) = work(0:nx)
382
383	ENDIF
384	#else
385	message_string = 'no system-specific fft-call available'
386	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
387	#endif
388	ELSE
389	message_string = 'fft method "' // TRIM( fft_method) // &
390	'" not available'
391	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
392
393	ENDIF
394
395	END SUBROUTINE fft_x
396
397	SUBROUTINE fft_y( ar, direction )
398
399	!----------------------------------------------------------------------!
400	! fft_y !
401	! !
402	! Fourier-transformation along y-direction !
403	! !
404	! fft_y uses internal algorithms (Singleton or Temperton) or !
405	! system-specific routines, if they are available !
406	!----------------------------------------------------------------------!
407
408	IMPLICIT NONE
409
410	CHARACTER (LEN=*) :: direction
411	INTEGER :: j, jshape(1)
412
413	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
414	REAL, DIMENSION(0:ny) :: ar
415	REAL, DIMENSION(0:ny+2) :: work
416	REAL, DIMENSION(ny+2) :: work1
417	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
418	#if defined( __ibm )
419	REAL, DIMENSION(nau2) :: auy2, auy4
420	#elif defined( __nec )
421	REAL, DIMENSION(6*(ny+1)) :: work2
422	#endif
423
424	IF ( fft_method == 'singleton-algorithm' ) THEN
425
426	!
427	!-- Performing the fft with singleton's software works on every system,
428	!-- since it is part of the model
429	ALLOCATE( cwork(0:ny) )
430
431	IF ( direction == 'forward') THEN
432
433	DO j = 0, ny
434	cwork(j) = CMPLX( ar(j) )
435	ENDDO
436
437	jshape = SHAPE( cwork )
438	CALL FFTN( cwork, jshape )
439
440	DO j = 0, (ny+1)/2
441	ar(j) = REAL( cwork(j) )
442	ENDDO
443	DO j = 1, (ny+1)/2 - 1
444	ar(ny+1-j) = -AIMAG( cwork(j) )
445	ENDDO
446
447	ELSE
448
449	cwork(0) = CMPLX( ar(0), 0.0 )
450	DO j = 1, (ny+1)/2 - 1
451	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
452	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
453	ENDDO
454	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
455
456	jshape = SHAPE( cwork )
457	CALL FFTN( cwork, jshape, inv = .TRUE. )
458
459	DO j = 0, ny
460	ar(j) = REAL( cwork(j) )
461	ENDDO
462
463	ENDIF
464
465	DEALLOCATE( cwork )
466
467	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
468
469	!
470	!-- Performing the fft with Temperton's software works on every system,
471	!-- since it is part of the model
472	IF ( direction == 'forward' ) THEN
473
474	work(0:ny) = ar
475	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
476
477	DO j = 0, (ny+1)/2
478	ar(j) = work(2*j)
479	ENDDO
480	DO j = 1, (ny+1)/2 - 1
481	ar(ny+1-j) = work(2*j+1)
482	ENDDO
483
484	ELSE
485
486	DO j = 0, (ny+1)/2
487	work(2*j) = ar(j)
488	ENDDO
489	DO j = 1, (ny+1)/2 - 1
490	work(2*j+1) = ar(ny+1-j)
491	ENDDO
492	work(1) = 0.0
493	work(ny+2) = 0.0
494
495	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
496	ar = work(0:ny)
497
498	ENDIF
499
500	ELSEIF ( fft_method == 'system-specific' ) THEN
501
502	#if defined( __ibm ) && ! defined( __ibmy_special )
503	IF ( direction == 'forward') THEN
504
505	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
506	auy2, nau2 )
507
508	DO j = 0, (ny+1)/2
509	ar(j) = work(2*j)
510	ENDDO
511	DO j = 1, (ny+1)/2 - 1
512	ar(ny+1-j) = work(2*j+1)
513	ENDDO
514
515	ELSE
516
517	DO j = 0, (ny+1)/2
518	work(2*j) = ar(j)
519	ENDDO
520	DO j = 1, (ny+1)/2 - 1
521	work(2*j+1) = ar(ny+1-j)
522	ENDDO
523	work(1) = 0.0
524	work(ny+2) = 0.0
525
526	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
527	auy4, nau2 )
528
529	DO j = 0, ny
530	ar(j) = work(j)
531	ENDDO
532
533	ENDIF
534	#elif defined( __nec )
535	IF ( direction == 'forward' ) THEN
536
537	work(0:ny) = ar(0:ny)
538
539	CALL DZFFT( 1, ny+1, sqr_ny, work, work, trig_yf, work2, 0 )
540
541	DO j = 0, (ny+1)/2
542	ar(j) = work(2*j)
543	ENDDO
544	DO j = 1, (ny+1)/2 - 1
545	ar(ny+1-j) = work(2*j+1)
546	ENDDO
547
548	ELSE
549
550	DO j = 0, (ny+1)/2
551	work(2*j) = ar(j)
552	ENDDO
553	DO j = 1, (ny+1)/2 - 1
554	work(2*j+1) = ar(ny+1-j)
555	ENDDO
556	work(1) = 0.0
557	work(ny+2) = 0.0
558
559	CALL ZDFFT( -1, ny+1, sqr_ny, work, work, trig_yb, work2, 0 )
560
561	ar(0:ny) = work(0:ny)
562
563	ENDIF
564	#else
565	message_string = 'no system-specific fft-call available'
566	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
567
568	#endif
569
570	ELSE
571
572	message_string = 'fft method "' // TRIM( fft_method) // &
573	'" not available'
574	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
575
576	ENDIF
577
578	END SUBROUTINE fft_y
579
580	SUBROUTINE fft_x_m( ar, direction )
581
582	!----------------------------------------------------------------------!
583	! fft_x_m !
584	! !
585	! Fourier-transformation along x-direction !
586	! Version for 1d domain decomposition !
587	! using multiple 1D FFT from Math Keisan on NEC !
588	! or Temperton-algorithm !
589	! (no singleton-algorithm on NEC because it does not vectorize) !
590	! !
591	!----------------------------------------------------------------------!
592
593	IMPLICIT NONE
594
595	CHARACTER (LEN=*) :: direction
596	INTEGER :: i, k, siza
597
598	REAL, DIMENSION(0:nx,nz) :: ar
599	REAL, DIMENSION(0:nx+3,nz+1) :: ai
600	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
601	#if defined( __nec )
602	INTEGER :: sizw
603	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
604	#endif
605
606	IF ( fft_method == 'temperton-algorithm' ) THEN
607
608	siza = SIZE( ai, 1 )
609
610	IF ( direction == 'forward') THEN
611
612	ai(0:nx,1:nz) = ar(0:nx,1:nz)
613	ai(nx+1:,:) = 0.0
614
615	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
616
617	DO k = 1, nz
618	DO i = 0, (nx+1)/2
619	ar(i,k) = ai(2*i,k)
620	ENDDO
621	DO i = 1, (nx+1)/2 - 1
622	ar(nx+1-i,k) = ai(2*i+1,k)
623	ENDDO
624	ENDDO
625
626	ELSE
627
628	DO k = 1, nz
629	DO i = 0, (nx+1)/2
630	ai(2*i,k) = ar(i,k)
631	ENDDO
632	DO i = 1, (nx+1)/2 - 1
633	ai(2*i+1,k) = ar(nx+1-i,k)
634	ENDDO
635	ai(1,k) = 0.0
636	ai(nx+2,k) = 0.0
637	ENDDO
638
639	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
640
641	ar(0:nx,1:nz) = ai(0:nx,1:nz)
642
643	ENDIF
644
645	ELSEIF ( fft_method == 'system-specific' ) THEN
646
647	#if defined( __nec )
648	siza = SIZE( ai, 1 )
649	sizw = SIZE( work, 1 )
650
651	IF ( direction == 'forward') THEN
652
653	!
654	!-- Tables are initialized once more. This call should not be
655	!-- necessary, but otherwise program aborts in asymmetric case
656	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
657	trig_xf, work1, 0 )
658
659	ai(0:nx,1:nz) = ar(0:nx,1:nz)
660	IF ( nz1 > nz ) THEN
661	ai(:,nz1) = 0.0
662	ENDIF
663
664	CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
665	trig_xf, work1, 0 )
666
667	DO k = 1, nz
668	DO i = 0, (nx+1)/2
669	ar(i,k) = REAL( work(i+1,k) )
670	ENDDO
671	DO i = 1, (nx+1)/2 - 1
672	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
673	ENDDO
674	ENDDO
675
676	ELSE
677
678	!
679	!-- Tables are initialized once more. This call should not be
680	!-- necessary, but otherwise program aborts in asymmetric case
681	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
682	trig_xb, work1, 0 )
683
684	IF ( nz1 > nz ) THEN
685	work(:,nz1) = 0.0
686	ENDIF
687	DO k = 1, nz
688	work(1,k) = CMPLX( ar(0,k), 0.0 )
689	DO i = 1, (nx+1)/2 - 1
690	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
691	ENDDO
692	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
693	ENDDO
694
695	CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
696	trig_xb, work1, 0 )
697
698	ar(0:nx,1:nz) = ai(0:nx,1:nz)
699
700	ENDIF
701
702	#else
703	message_string = 'no system-specific fft-call available'
704	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
705	#endif
706
707	ELSE
708
709	message_string = 'fft method "' // TRIM( fft_method) // &
710	'" not available'
711	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
712
713	ENDIF
714
715	END SUBROUTINE fft_x_m
716
717	SUBROUTINE fft_y_m( ar, ny1, direction )
718
719	!----------------------------------------------------------------------!
720	! fft_y_m !
721	! !
722	! Fourier-transformation along y-direction !
723	! Version for 1d domain decomposition !
724	! using multiple 1D FFT from Math Keisan on NEC !
725	! or Temperton-algorithm !
726	! (no singleton-algorithm on NEC because it does not vectorize) !
727	! !
728	!----------------------------------------------------------------------!
729
730	IMPLICIT NONE
731
732	CHARACTER (LEN=*) :: direction
733	INTEGER :: j, k, ny1, siza
734
735	REAL, DIMENSION(0:ny1,nz) :: ar
736	REAL, DIMENSION(0:ny+3,nz+1) :: ai
737	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
738	#if defined( __nec )
739	INTEGER :: sizw
740	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
741	#endif
742
743	IF ( fft_method == 'temperton-algorithm' ) THEN
744
745	siza = SIZE( ai, 1 )
746
747	IF ( direction == 'forward') THEN
748
749	ai(0:ny,1:nz) = ar(0:ny,1:nz)
750	ai(ny+1:,:) = 0.0
751
752	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
753
754	DO k = 1, nz
755	DO j = 0, (ny+1)/2
756	ar(j,k) = ai(2*j,k)
757	ENDDO
758	DO j = 1, (ny+1)/2 - 1
759	ar(ny+1-j,k) = ai(2*j+1,k)
760	ENDDO
761	ENDDO
762
763	ELSE
764
765	DO k = 1, nz
766	DO j = 0, (ny+1)/2
767	ai(2*j,k) = ar(j,k)
768	ENDDO
769	DO j = 1, (ny+1)/2 - 1
770	ai(2*j+1,k) = ar(ny+1-j,k)
771	ENDDO
772	ai(1,k) = 0.0
773	ai(ny+2,k) = 0.0
774	ENDDO
775
776	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
777
778	ar(0:ny,1:nz) = ai(0:ny,1:nz)
779
780	ENDIF
781
782	ELSEIF ( fft_method == 'system-specific' ) THEN
783
784	#if defined( __nec )
785	siza = SIZE( ai, 1 )
786	sizw = SIZE( work, 1 )
787
788	IF ( direction == 'forward') THEN
789
790	!
791	!-- Tables are initialized once more. This call should not be
792	!-- necessary, but otherwise program aborts in asymmetric case
793	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
794	trig_yf, work1, 0 )
795
796	ai(0:ny,1:nz) = ar(0:ny,1:nz)
797	IF ( nz1 > nz ) THEN
798	ai(:,nz1) = 0.0
799	ENDIF
800
801	CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
802	trig_yf, work1, 0 )
803
804	DO k = 1, nz
805	DO j = 0, (ny+1)/2
806	ar(j,k) = REAL( work(j+1,k) )
807	ENDDO
808	DO j = 1, (ny+1)/2 - 1
809	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
810	ENDDO
811	ENDDO
812
813	ELSE
814
815	!
816	!-- Tables are initialized once more. This call should not be
817	!-- necessary, but otherwise program aborts in asymmetric case
818	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
819	trig_yb, work1, 0 )
820
821	IF ( nz1 > nz ) THEN
822	work(:,nz1) = 0.0
823	ENDIF
824	DO k = 1, nz
825	work(1,k) = CMPLX( ar(0,k), 0.0 )
826	DO j = 1, (ny+1)/2 - 1
827	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
828	ENDDO
829	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
830	ENDDO
831
832	CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
833	trig_yb, work1, 0 )
834
835	ar(0:ny,1:nz) = ai(0:ny,1:nz)
836
837	ENDIF
838
839	#else
840	message_string = 'no system-specific fft-call available'
841	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
842	#endif
843
844	ELSE
845
846	message_string = 'fft method "' // TRIM( fft_method) // &
847	'" not available'
848	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
849
850	ENDIF
851
852	END SUBROUTINE fft_y_m
853
854	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |