Home

Context Navigation

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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago
code has been put under the GNU General Public License (v3)
Property svn:keywords set to `Id`
File size: 25.5 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 1036 2012-10-22 13:43:42Z raasch $
27	!
28	! 274 2009-03-26 15:11:21Z heinze
29	! Output of messages replaced by message handling routine.
30	!
31	! Feb. 2007
32	! RCS Log replace by Id keyword, revision history cleaned up
33	!
34	! Revision 1.4 2006/03/28 12:27:09 raasch
35	! Stop when system-specific fft is selected on NEC. For unknown reasons this
36	! causes a program abort during first allocation in init_grid.
37	!
38	! Revision 1.2 2004/04/30 11:44:27 raasch
39	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
40	! fft_x and fft_y,
41	! function FFT replaced by subroutine FFTN due to problems with 64-bit
42	! mode on ibm,
43	! shape of array cwork is explicitly stored in ishape/jshape and handled
44	! to routine FFTN instead of shape-function (due to compiler error on
45	! decalpha),
46	! non vectorized FFT for nec included
47	!
48	! Revision 1.1 2002/06/11 13:00:49 raasch
49	! Initial revision
50	!
51	!
52	! Description:
53	! ------------
54	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
55	! Original version: Klaus Ketelsen (May 2002)
56	!------------------------------------------------------------------------------!
57
58	USE array_kind
59	USE control_parameters
60	USE indices
61	USE singleton
62	USE temperton_fft
63
64	IMPLICIT NONE
65
66	PRIVATE
67	PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
68
69	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
70
71	LOGICAL, SAVE :: init_fft = .FALSE.
72
73	REAL, SAVE :: sqr_nx, sqr_ny
74	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
75
76	#if defined( __ibm )
77	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
78	!
79	!-- The following working arrays contain tables and have to be "save" and
80	!-- shared in OpenMP sense
81	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
82	#elif defined( __nec )
83	INTEGER, SAVE :: nz1
84	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
85	trig_yf
86	#endif
87
88	!
89	!-- Public interfaces
90	INTERFACE fft_init
91	MODULE PROCEDURE fft_init
92	END INTERFACE fft_init
93
94	INTERFACE fft_x
95	MODULE PROCEDURE fft_x
96	END INTERFACE fft_x
97
98	INTERFACE fft_y
99	MODULE PROCEDURE fft_y
100	END INTERFACE fft_y
101
102	INTERFACE fft_x_m
103	MODULE PROCEDURE fft_x_m
104	END INTERFACE fft_x_m
105
106	INTERFACE fft_y_m
107	MODULE PROCEDURE fft_y_m
108	END INTERFACE fft_y_m
109
110	CONTAINS
111
112
113	SUBROUTINE fft_init
114
115	IMPLICIT NONE
116
117	!
118	!-- The following temporary working arrays have to be on stack or private
119	!-- in OpenMP sense
120	#if defined( __ibm )
121	REAL, DIMENSION(0:nx+2) :: workx
122	REAL, DIMENSION(0:ny+2) :: worky
123	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
124	#elif defined( __nec )
125	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
126	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
127	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
128	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
129	#endif
130
131	!
132	!-- Return, if already called
133	IF ( init_fft ) THEN
134	RETURN
135	ELSE
136	init_fft = .TRUE.
137	ENDIF
138
139	IF ( fft_method == 'system-specific' ) THEN
140
141	sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
142	sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
143	#if defined( __ibm ) && ! defined( __ibmy_special )
144	!
145	!-- Initialize tables for fft along x
146	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
147	aux2, nau2 )
148	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
149	aux4, nau2 )
150	!
151	!-- Initialize tables for fft along y
152	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
153	auy2, nau2 )
154	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
155	auy4, nau2 )
156	#elif defined( __nec )
157	message_string = 'fft method "' // TRIM( fft_method) // &
158	'" currently does not work on NEC'
159	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
160
161	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
162	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
163
164	work_x = 0.0
165	work_y = 0.0
166	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
167	! when using the NEC ffts
168
169	!
170	!-- Initialize tables for fft along x (non-vector and vector case (M))
171	CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
172	CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
173	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
174	trig_xf, workx, 0 )
175	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
176	trig_xb, workx, 0 )
177	!
178	!-- Initialize tables for fft along y (non-vector and vector case (M))
179	CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
180	CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
181	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
182	trig_yf, worky, 0 )
183	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
184	trig_yb, worky, 0 )
185	#else
186	message_string = 'no system-specific fft-call available'
187	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
188	#endif
189	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
190	!
191	!-- Temperton-algorithm
192	!-- Initialize tables for fft along x and y
193	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
194
195	CALL set99( trigs_x, ifax_x, nx+1 )
196	CALL set99( trigs_y, ifax_y, ny+1 )
197
198	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
199
200	CONTINUE
201
202	ELSE
203
204	message_string = 'fft method "' // TRIM( fft_method) // &
205	'" not available'
206	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
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	message_string = 'no system-specific fft-call available'
380	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
381	#endif
382	ELSE
383	message_string = 'fft method "' // TRIM( fft_method) // &
384	'" not available'
385	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
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	message_string = 'no system-specific fft-call available'
560	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
561
562	#endif
563
564	ELSE
565
566	message_string = 'fft method "' // TRIM( fft_method) // &
567	'" not available'
568	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
569
570	ENDIF
571
572	END SUBROUTINE fft_y
573
574	SUBROUTINE fft_x_m( ar, direction )
575
576	!----------------------------------------------------------------------!
577	! fft_x_m !
578	! !
579	! Fourier-transformation along x-direction !
580	! Version for 1d domain decomposition !
581	! using multiple 1D FFT from Math Keisan on NEC !
582	! or Temperton-algorithm !
583	! (no singleton-algorithm on NEC because it does not vectorize) !
584	! !
585	!----------------------------------------------------------------------!
586
587	IMPLICIT NONE
588
589	CHARACTER (LEN=*) :: direction
590	INTEGER :: i, k, siza, sizw
591
592	REAL, DIMENSION(0:nx,nz) :: ar
593	REAL, DIMENSION(0:nx+3,nz+1) :: ai
594	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
595	#if defined( __nec )
596	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
597	#endif
598
599	IF ( fft_method == 'temperton-algorithm' ) THEN
600
601	siza = SIZE( ai, 1 )
602
603	IF ( direction == 'forward') THEN
604
605	ai(0:nx,1:nz) = ar(0:nx,1:nz)
606	ai(nx+1:,:) = 0.0
607
608	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
609
610	DO k = 1, nz
611	DO i = 0, (nx+1)/2
612	ar(i,k) = ai(2*i,k)
613	ENDDO
614	DO i = 1, (nx+1)/2 - 1
615	ar(nx+1-i,k) = ai(2*i+1,k)
616	ENDDO
617	ENDDO
618
619	ELSE
620
621	DO k = 1, nz
622	DO i = 0, (nx+1)/2
623	ai(2*i,k) = ar(i,k)
624	ENDDO
625	DO i = 1, (nx+1)/2 - 1
626	ai(2*i+1,k) = ar(nx+1-i,k)
627	ENDDO
628	ai(1,k) = 0.0
629	ai(nx+2,k) = 0.0
630	ENDDO
631
632	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
633
634	ar(0:nx,1:nz) = ai(0:nx,1:nz)
635
636	ENDIF
637
638	ELSEIF ( fft_method == 'system-specific' ) THEN
639
640	#if defined( __nec )
641	siza = SIZE( ai, 1 )
642	sizw = SIZE( work, 1 )
643
644	IF ( direction == 'forward') THEN
645
646	!
647	!-- Tables are initialized once more. This call should not be
648	!-- necessary, but otherwise program aborts in asymmetric case
649	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
650	trig_xf, work1, 0 )
651
652	ai(0:nx,1:nz) = ar(0:nx,1:nz)
653	IF ( nz1 > nz ) THEN
654	ai(:,nz1) = 0.0
655	ENDIF
656
657	CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
658	trig_xf, work1, 0 )
659
660	DO k = 1, nz
661	DO i = 0, (nx+1)/2
662	ar(i,k) = REAL( work(i+1,k) )
663	ENDDO
664	DO i = 1, (nx+1)/2 - 1
665	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
666	ENDDO
667	ENDDO
668
669	ELSE
670
671	!
672	!-- Tables are initialized once more. This call should not be
673	!-- necessary, but otherwise program aborts in asymmetric case
674	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
675	trig_xb, work1, 0 )
676
677	IF ( nz1 > nz ) THEN
678	work(:,nz1) = 0.0
679	ENDIF
680	DO k = 1, nz
681	work(1,k) = CMPLX( ar(0,k), 0.0 )
682	DO i = 1, (nx+1)/2 - 1
683	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
684	ENDDO
685	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
686	ENDDO
687
688	CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
689	trig_xb, work1, 0 )
690
691	ar(0:nx,1:nz) = ai(0:nx,1:nz)
692
693	ENDIF
694
695	#else
696	message_string = 'no system-specific fft-call available'
697	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
698	#endif
699
700	ELSE
701
702	message_string = 'fft method "' // TRIM( fft_method) // &
703	'" not available'
704	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
705
706	ENDIF
707
708	END SUBROUTINE fft_x_m
709
710	SUBROUTINE fft_y_m( ar, ny1, direction )
711
712	!----------------------------------------------------------------------!
713	! fft_y_m !
714	! !
715	! Fourier-transformation along y-direction !
716	! Version for 1d domain decomposition !
717	! using multiple 1D FFT from Math Keisan on NEC !
718	! or Temperton-algorithm !
719	! (no singleton-algorithm on NEC because it does not vectorize) !
720	! !
721	!----------------------------------------------------------------------!
722
723	IMPLICIT NONE
724
725	CHARACTER (LEN=*) :: direction
726	INTEGER :: j, k, ny1, siza, sizw
727
728	REAL, DIMENSION(0:ny1,nz) :: ar
729	REAL, DIMENSION(0:ny+3,nz+1) :: ai
730	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
731	#if defined( __nec )
732	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
733	#endif
734
735	IF ( fft_method == 'temperton-algorithm' ) THEN
736
737	siza = SIZE( ai, 1 )
738
739	IF ( direction == 'forward') THEN
740
741	ai(0:ny,1:nz) = ar(0:ny,1:nz)
742	ai(ny+1:,:) = 0.0
743
744	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
745
746	DO k = 1, nz
747	DO j = 0, (ny+1)/2
748	ar(j,k) = ai(2*j,k)
749	ENDDO
750	DO j = 1, (ny+1)/2 - 1
751	ar(ny+1-j,k) = ai(2*j+1,k)
752	ENDDO
753	ENDDO
754
755	ELSE
756
757	DO k = 1, nz
758	DO j = 0, (ny+1)/2
759	ai(2*j,k) = ar(j,k)
760	ENDDO
761	DO j = 1, (ny+1)/2 - 1
762	ai(2*j+1,k) = ar(ny+1-j,k)
763	ENDDO
764	ai(1,k) = 0.0
765	ai(ny+2,k) = 0.0
766	ENDDO
767
768	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
769
770	ar(0:ny,1:nz) = ai(0:ny,1:nz)
771
772	ENDIF
773
774	ELSEIF ( fft_method == 'system-specific' ) THEN
775
776	#if defined( __nec )
777	siza = SIZE( ai, 1 )
778	sizw = SIZE( work, 1 )
779
780	IF ( direction == 'forward') THEN
781
782	!
783	!-- Tables are initialized once more. This call should not be
784	!-- necessary, but otherwise program aborts in asymmetric case
785	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
786	trig_yf, work1, 0 )
787
788	ai(0:ny,1:nz) = ar(0:ny,1:nz)
789	IF ( nz1 > nz ) THEN
790	ai(:,nz1) = 0.0
791	ENDIF
792
793	CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
794	trig_yf, work1, 0 )
795
796	DO k = 1, nz
797	DO j = 0, (ny+1)/2
798	ar(j,k) = REAL( work(j+1,k) )
799	ENDDO
800	DO j = 1, (ny+1)/2 - 1
801	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
802	ENDDO
803	ENDDO
804
805	ELSE
806
807	!
808	!-- Tables are initialized once more. This call should not be
809	!-- necessary, but otherwise program aborts in asymmetric case
810	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
811	trig_yb, work1, 0 )
812
813	IF ( nz1 > nz ) THEN
814	work(:,nz1) = 0.0
815	ENDIF
816	DO k = 1, nz
817	work(1,k) = CMPLX( ar(0,k), 0.0 )
818	DO j = 1, (ny+1)/2 - 1
819	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
820	ENDDO
821	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
822	ENDDO
823
824	CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
825	trig_yb, work1, 0 )
826
827	ar(0:ny,1:nz) = ai(0:ny,1:nz)
828
829	ENDIF
830
831	#else
832	message_string = 'no system-specific fft-call available'
833	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
834	#endif
835
836	ELSE
837
838	message_string = 'fft method "' // TRIM( fft_method) // &
839	'" not available'
840	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
841
842	ENDIF
843
844	END SUBROUTINE fft_y_m
845
846	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |