Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |