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

Last change on this file since 235 was 4, checked in by raasch, 18 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 24.1 KB
Line 
1 MODULE fft_xy
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: fft_xy.f90 4 2007-02-13 11:33:16Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.4  2006/03/28 12:27:09  raasch
14! Stop when system-specific fft is selected on NEC. For unknown reasons this
15! causes a program abort during first allocation in init_grid.
16!
17! Revision 1.2  2004/04/30 11:44:27  raasch
18! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
19! fft_x and fft_y,
20! function FFT replaced by subroutine FFTN due to problems with 64-bit
21! mode on ibm,
22! shape of array cwork is explicitly stored in ishape/jshape and handled
23! to routine FFTN instead of shape-function (due to compiler error on
24! decalpha),
25! non vectorized FFT for nec included
26!
27! Revision 1.1  2002/06/11 13:00:49  raasch
28! Initial revision
29!
30!
31! Description:
32! ------------
33! Fast Fourier transformation along x and y for 1d domain decomposition along x.
34! Original version: Klaus Ketelsen (May 2002)
35!------------------------------------------------------------------------------!
36
37    USE array_kind
38    USE control_parameters
39    USE indices
40    USE singleton
41    USE temperton_fft
42
43    IMPLICIT NONE
44
45    PRIVATE
46    PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
47
48    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
49
50    LOGICAL, SAVE                            ::  init_fft = .FALSE.
51
52    REAL, SAVE ::  sqr_nx, sqr_ny
53    REAL, DIMENSION(:), ALLOCATABLE, SAVE    ::  trigs_x, trigs_y
54
55#if defined( __ibm )
56    INTEGER, PARAMETER ::  nau1 = 20000, nau2 = 22000
57!
58!-- The following working arrays contain tables and have to be "save" and
59!-- shared in OpenMP sense
60    REAL, DIMENSION(nau1), SAVE ::  aux1, auy1, aux3, auy3
61#elif defined( __nec )
62    INTEGER, SAVE ::  nz1
63    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb, trig_xf, trig_yb, &
64                                              trig_yf
65#endif
66
67!
68!-- Public interfaces
69    INTERFACE fft_init
70       MODULE PROCEDURE fft_init
71    END INTERFACE fft_init
72
73    INTERFACE fft_x
74       MODULE PROCEDURE fft_x
75    END INTERFACE fft_x
76
77    INTERFACE fft_y
78       MODULE PROCEDURE fft_y
79    END INTERFACE fft_y
80
81    INTERFACE fft_x_m
82       MODULE PROCEDURE fft_x_m
83    END INTERFACE fft_x_m
84
85    INTERFACE fft_y_m
86       MODULE PROCEDURE fft_y_m
87    END INTERFACE fft_y_m
88
89 CONTAINS
90
91
92    SUBROUTINE fft_init
93
94       IMPLICIT NONE
95
96!
97!--    The following temporary working arrays have to be on stack or private
98!--    in OpenMP sense
99#if defined( __ibm )
100       REAL, DIMENSION(0:nx+2) :: workx
101       REAL, DIMENSION(0:ny+2) :: worky
102       REAL, DIMENSION(nau2)   :: aux2, auy2, aux4, auy4
103#elif defined( __nec )
104       REAL, DIMENSION(0:nx+3,nz+1)   ::  work_x
105       REAL, DIMENSION(0:ny+3,nz+1)   ::  work_y
106       REAL, DIMENSION(6*(nx+3),nz+1) ::  workx
107       REAL, DIMENSION(6*(ny+3),nz+1) ::  worky
108#endif 
109
110!
111!--    Return, if already called
112       IF ( init_fft )  THEN
113          RETURN
114       ELSE
115          init_fft = .TRUE.
116       ENDIF
117
118       IF ( fft_method == 'system-specific' )  THEN
119
120          sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
121          sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
122#if defined( __ibm )  &&  ! defined( __ibmy_special )
123!
124!--       Initialize tables for fft along x
125          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_nx, aux1, nau1, &
126                      aux2, nau2 )
127          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
128                      aux4, nau2 )
129!
130!--       Initialize tables for fft along y
131          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_ny, auy1, nau1, &
132                      auy2, nau2 )
133          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
134                      auy4, nau2 )
135#elif defined( __nec )
136          PRINT*, '+++ fft_init: fft method "', fft_method, &
137                  '" currently does not work on NEC'
138          CALL local_stop
139
140          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)), &
141                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
142
143          work_x = 0.0
144          work_y = 0.0
145          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
146                                      ! when using the NEC ffts
147
148!
149!--       Initialize tables for fft along x (non-vector and vector case (M))
150          CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
151          CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
152          CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
153                       trig_xf, workx, 0 )
154          CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
155                       trig_xb, workx, 0 )
156!
157!--       Initialize tables for fft along y (non-vector and vector case (M))
158          CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
159          CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
160          CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
161                       trig_yf, worky, 0 )
162          CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
163                       trig_yb, worky, 0 )
164#else
165          PRINT*, '+++ fft_init: no system-specific fft-call available'
166          CALL local_stop
167#endif
168       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
169!
170!--       Temperton-algorithm
171!--       Initialize tables for fft along x and y
172          ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
173
174          CALL set99( trigs_x, ifax_x, nx+1 )
175          CALL set99( trigs_y, ifax_y, ny+1 )
176
177       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
178
179          CONTINUE
180
181       ELSE
182
183          PRINT*, '+++ fft_init: fft method "', fft_method, &
184                  '" not available'
185          CALL local_stop
186
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          PRINT*, '+++ fft_x: no system-specific fft-call available'
360          CALL local_stop
361#endif
362       ELSE
363
364          PRINT*, '+++ fft_x: fft method "', fft_method, '" not available'
365          CALL local_stop
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          PRINT*, '+++ fft_y: no system-specific fft-call available'
540          CALL local_stop
541#endif
542
543       ELSE
544
545          PRINT*, '+++ fft_y: fft method "', fft_method, '" not available'
546          CALL local_stop
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          PRINT*, '+++ fft_x_m: no system-specific fft-call available'
675          STOP
676#endif
677
678       ELSE
679
680          PRINT*, '+++ fft_x_m: fft method "', fft_method, '" not available'
681          CALL local_stop
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          PRINT*, '+++ fft_y_m: no system-specific fft-call available'
810          STOP
811#endif
812
813       ELSE
814
815          PRINT*, '+++ fft_y_m: fft method "', fft_method, '" not available'
816          CALL local_stop
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.