Home

Context Navigation

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

Last change on this file since 1 was 1, checked in by raasch, 18 years ago
Initial repository layout and content
File size: 24.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |