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

Last change on this file since 254 was 254, checked in by heinze, 13 years ago

Output of messages replaced by message handling routine.

  • Property svn:keywords set to Id
File size: 24.6 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 254 2009-03-05 15:33:42Z heinze $
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) // '" not available'
364          CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
365
366       ENDIF
367
368    END SUBROUTINE fft_x
369
370    SUBROUTINE fft_y( ar, direction )
371
372!----------------------------------------------------------------------!
373!                                 fft_y                                !
374!                                                                      !
375!               Fourier-transformation along y-direction               !
376!                                                                      !
377!      fft_y uses internal algorithms (Singleton or Temperton) or      !
378!           system-specific routines, if they are available            !
379!----------------------------------------------------------------------!
380
381       IMPLICIT NONE
382
383       CHARACTER (LEN=*) ::  direction
384       INTEGER ::  j, jshape(1)
385
386!kk    REAL, DIMENSION(:)        ::  ar !kk Does NOT work (Bug??)
387       REAL, DIMENSION(0:ny)     ::  ar
388       REAL, DIMENSION(0:ny+2)   ::  work
389       REAL, DIMENSION(ny+2)     ::  work1
390       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
391#if defined( __ibm )
392       REAL, DIMENSION(nau2)     ::  auy2, auy4
393#elif defined( __nec )
394       REAL, DIMENSION(6*(ny+1)) ::  work2
395#endif
396
397       IF ( fft_method == 'singleton-algorithm' )  THEN
398
399!
400!--       Performing the fft with singleton's software works on every system,
401!--       since it is part of the model
402          ALLOCATE( cwork(0:ny) )
403
404          IF ( direction == 'forward')  THEN
405
406             DO  j = 0, ny
407                cwork(j) = CMPLX( ar(j) )
408             ENDDO
409
410             jshape = SHAPE( cwork )
411             CALL FFTN( cwork, jshape )
412
413             DO  j = 0, (ny+1)/2
414                ar(j) = REAL( cwork(j) )
415             ENDDO
416             DO  j = 1, (ny+1)/2 - 1
417                ar(ny+1-j) = -AIMAG( cwork(j) )
418             ENDDO
419
420          ELSE
421
422             cwork(0) = CMPLX( ar(0), 0.0 )
423             DO  j = 1, (ny+1)/2 - 1
424                cwork(j)      = CMPLX( ar(j), -ar(ny+1-j) )
425                cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j) )
426             ENDDO
427             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
428
429             jshape = SHAPE( cwork )
430             CALL FFTN( cwork, jshape, inv = .TRUE. )
431
432             DO  j = 0, ny
433                ar(j) = REAL( cwork(j) )
434             ENDDO
435
436          ENDIF
437
438          DEALLOCATE( cwork )
439
440       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
441
442!
443!--       Performing the fft with Temperton's software works on every system,
444!--       since it is part of the model
445          IF ( direction == 'forward' )  THEN
446
447             work(0:ny) = ar
448             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
449
450             DO  j = 0, (ny+1)/2
451                ar(j) = work(2*j)
452             ENDDO
453             DO  j = 1, (ny+1)/2 - 1
454                ar(ny+1-j) = work(2*j+1)
455             ENDDO
456
457          ELSE
458
459             DO  j = 0, (ny+1)/2
460                work(2*j) = ar(j)
461             ENDDO
462             DO  j = 1, (ny+1)/2 - 1
463                work(2*j+1) = ar(ny+1-j)
464             ENDDO
465             work(1)    = 0.0
466             work(ny+2) = 0.0
467
468             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
469             ar = work(0:ny)
470
471          ENDIF
472
473       ELSEIF ( fft_method == 'system-specific' )  THEN
474
475#if defined( __ibm )  &&  ! defined( __ibmy_special )
476          IF ( direction == 'forward')  THEN
477
478             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
479                         auy2, nau2 )
480
481             DO  j = 0, (ny+1)/2
482                ar(j) = work(2*j)
483             ENDDO
484             DO  j = 1, (ny+1)/2 - 1
485                ar(ny+1-j) = work(2*j+1)
486             ENDDO
487
488          ELSE
489
490             DO  j = 0, (ny+1)/2
491                work(2*j) = ar(j)
492             ENDDO
493             DO  j = 1, (ny+1)/2 - 1
494                work(2*j+1) = ar(ny+1-j)
495             ENDDO
496             work(1)    = 0.0
497             work(ny+2) = 0.0
498
499             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
500                         auy4, nau2 )
501
502             DO  j = 0, ny
503                ar(j) = work(j)
504             ENDDO
505
506          ENDIF
507#elif defined( __nec )
508          IF ( direction == 'forward' )  THEN
509
510             work(0:ny) = ar(0:ny)
511
512             CALL DZFFT( 1, ny+1, sqr_ny, work, work, trig_yf, work2, 0 )
513
514             DO  j = 0, (ny+1)/2
515                ar(j) = work(2*j)
516             ENDDO
517             DO  j = 1, (ny+1)/2 - 1
518                ar(ny+1-j) = work(2*j+1)
519             ENDDO
520
521          ELSE
522
523             DO  j = 0, (ny+1)/2
524                work(2*j) = ar(j)
525             ENDDO
526             DO  j = 1, (ny+1)/2 - 1
527                work(2*j+1) = ar(ny+1-j)
528             ENDDO
529             work(1) = 0.0
530             work(ny+2) = 0.0
531
532             CALL ZDFFT( -1, ny+1, sqr_ny, work, work, trig_yb, work2, 0 )
533
534             ar(0:ny) = work(0:ny)
535
536          ENDIF
537#else
538          message_string = 'no system-specific fft-call available'
539          CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 ) 
540
541#endif
542
543       ELSE
544
545          message_string = 'fft method "' // TRIM( fft_method) // '" not available'
546          CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
547
548       ENDIF
549
550    END SUBROUTINE fft_y
551
552    SUBROUTINE fft_x_m( ar, direction )
553
554!----------------------------------------------------------------------!
555!                               fft_x_m                                !
556!                                                                      !
557!               Fourier-transformation along x-direction               !
558!                 Version for 1d domain decomposition                  !
559!             using multiple 1D FFT from Math Keisan on NEC            !
560!                       or Temperton-algorithm                         !
561!       (no singleton-algorithm on NEC because it does not vectorize)  !
562!                                                                      !
563!----------------------------------------------------------------------!
564
565       IMPLICIT NONE
566
567       CHARACTER (LEN=*)              ::  direction
568       INTEGER                        ::  i, k, siza, sizw
569
570       REAL, DIMENSION(0:nx,nz)       ::  ar
571       REAL, DIMENSION(0:nx+3,nz+1)   ::  ai
572       REAL, DIMENSION(6*(nx+4),nz+1) ::  work1
573#if defined( __nec )
574       COMPLEX, DIMENSION((nx+4)/2+1,nz+1) ::  work
575#endif
576
577       IF ( fft_method == 'temperton-algorithm' )  THEN
578
579          siza = SIZE( ai, 1 )
580
581          IF ( direction == 'forward')  THEN
582
583             ai(0:nx,1:nz) = ar(0:nx,1:nz)
584             ai(nx+1:,:)   = 0.0
585
586             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
587
588             DO  k = 1, nz
589                DO  i = 0, (nx+1)/2
590                   ar(i,k) = ai(2*i,k)
591                ENDDO
592                DO  i = 1, (nx+1)/2 - 1
593                   ar(nx+1-i,k) = ai(2*i+1,k)
594                ENDDO
595             ENDDO
596
597          ELSE
598
599             DO  k = 1, nz
600                DO  i = 0, (nx+1)/2
601                   ai(2*i,k) = ar(i,k)
602                ENDDO
603                DO  i = 1, (nx+1)/2 - 1
604                   ai(2*i+1,k) = ar(nx+1-i,k)
605                ENDDO
606                ai(1,k) = 0.0
607                ai(nx+2,k) = 0.0
608             ENDDO
609
610             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
611
612             ar(0:nx,1:nz) = ai(0:nx,1:nz)
613
614          ENDIF
615
616       ELSEIF ( fft_method == 'system-specific' )  THEN
617
618#if defined( __nec )
619          siza = SIZE( ai, 1 )
620          sizw = SIZE( work, 1 )
621
622          IF ( direction == 'forward')  THEN
623
624!
625!--          Tables are initialized once more. This call should not be
626!--          necessary, but otherwise program aborts in asymmetric case
627             CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
628                          trig_xf, work1, 0 )
629
630             ai(0:nx,1:nz) = ar(0:nx,1:nz)
631             IF ( nz1 > nz )  THEN
632                ai(:,nz1) = 0.0
633             ENDIF
634
635             CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
636                          trig_xf, work1, 0 )
637
638             DO  k = 1, nz
639                DO  i = 0, (nx+1)/2
640                   ar(i,k) = REAL( work(i+1,k) )
641                ENDDO
642                DO  i = 1, (nx+1)/2 - 1
643                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
644                ENDDO
645             ENDDO
646
647          ELSE
648
649!
650!--          Tables are initialized once more. This call should not be
651!--          necessary, but otherwise program aborts in asymmetric case
652             CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
653                          trig_xb, work1, 0 )
654
655             IF ( nz1 > nz )  THEN
656                work(:,nz1) = 0.0
657             ENDIF
658             DO  k = 1, nz
659                work(1,k) = CMPLX( ar(0,k), 0.0 )
660                DO  i = 1, (nx+1)/2 - 1
661                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
662                ENDDO
663                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
664             ENDDO
665
666             CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
667                          trig_xb, work1, 0 )
668
669             ar(0:nx,1:nz) = ai(0:nx,1:nz)
670
671          ENDIF
672
673#else
674          message_string = 'no system-specific fft-call available'
675          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
676#endif
677
678       ELSE
679
680          message_string = 'fft method "' // TRIM( fft_method) // '" not available'
681          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
682
683       ENDIF
684
685    END SUBROUTINE fft_x_m
686
687    SUBROUTINE fft_y_m( ar, ny1, direction )
688
689!----------------------------------------------------------------------!
690!                               fft_y_m                                !
691!                                                                      !
692!               Fourier-transformation along y-direction               !
693!                 Version for 1d domain decomposition                  !
694!             using multiple 1D FFT from Math Keisan on NEC            !
695!                       or Temperton-algorithm                         !
696!       (no singleton-algorithm on NEC because it does not vectorize)  !
697!                                                                      !
698!----------------------------------------------------------------------!
699
700       IMPLICIT NONE
701
702       CHARACTER (LEN=*)              ::  direction
703       INTEGER                        ::  j, k, ny1, siza, sizw
704
705       REAL, DIMENSION(0:ny1,nz)      ::  ar
706       REAL, DIMENSION(0:ny+3,nz+1)   ::  ai
707       REAL, DIMENSION(6*(ny+4),nz+1) ::  work1
708#if defined( __nec )
709       COMPLEX, DIMENSION((ny+4)/2+1,nz+1) ::  work
710#endif
711
712       IF ( fft_method == 'temperton-algorithm' )  THEN
713
714          siza = SIZE( ai, 1 )
715
716          IF ( direction == 'forward')  THEN
717
718             ai(0:ny,1:nz) = ar(0:ny,1:nz)
719             ai(ny+1:,:)   = 0.0
720
721             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
722
723             DO  k = 1, nz
724                DO  j = 0, (ny+1)/2
725                   ar(j,k) = ai(2*j,k)
726                ENDDO
727                DO  j = 1, (ny+1)/2 - 1
728                   ar(ny+1-j,k) = ai(2*j+1,k)
729                ENDDO
730             ENDDO
731
732          ELSE
733
734             DO  k = 1, nz
735                DO  j = 0, (ny+1)/2
736                   ai(2*j,k) = ar(j,k)
737                ENDDO
738                DO  j = 1, (ny+1)/2 - 1
739                   ai(2*j+1,k) = ar(ny+1-j,k)
740                ENDDO
741                ai(1,k) = 0.0
742                ai(ny+2,k) = 0.0
743             ENDDO
744
745             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
746
747             ar(0:ny,1:nz) = ai(0:ny,1:nz)
748
749          ENDIF
750
751       ELSEIF ( fft_method == 'system-specific' )  THEN
752
753#if defined( __nec )
754          siza = SIZE( ai, 1 )
755          sizw = SIZE( work, 1 )
756
757          IF ( direction == 'forward')  THEN
758
759!
760!--          Tables are initialized once more. This call should not be
761!--          necessary, but otherwise program aborts in asymmetric case
762             CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
763                          trig_yf, work1, 0 )
764
765             ai(0:ny,1:nz) = ar(0:ny,1:nz)
766             IF ( nz1 > nz )  THEN
767                ai(:,nz1) = 0.0
768             ENDIF
769
770             CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
771                          trig_yf, work1, 0 )
772
773             DO  k = 1, nz
774                DO  j = 0, (ny+1)/2
775                   ar(j,k) = REAL( work(j+1,k) )
776                ENDDO
777                DO  j = 1, (ny+1)/2 - 1
778                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
779                ENDDO
780             ENDDO
781
782          ELSE
783
784!
785!--          Tables are initialized once more. This call should not be
786!--          necessary, but otherwise program aborts in asymmetric case
787             CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
788                          trig_yb, work1, 0 )
789
790             IF ( nz1 > nz )  THEN
791                work(:,nz1) = 0.0
792             ENDIF
793             DO  k = 1, nz
794                work(1,k) = CMPLX( ar(0,k), 0.0 )
795                DO  j = 1, (ny+1)/2 - 1
796                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
797                ENDDO
798                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
799             ENDDO
800
801             CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
802                          trig_yb, work1, 0 )
803
804             ar(0:ny,1:nz) = ai(0:ny,1:nz)
805
806          ENDIF
807
808#else
809          message_string = 'no system-specific fft-call available'
810          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
811#endif
812
813       ELSE
814         
815          message_string = 'fft method "' // TRIM( fft_method) // '" not available'
816          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
817
818       ENDIF
819
820    END SUBROUTINE fft_y_m
821
822 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.