Home

Context Navigation

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

Last change on this file since 332 was 274, checked in by heinze, 16 years ago
Indentation of the message calls corrected
Property svn:keywords set to `Id`
File size: 24.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |