source: LMDZ4/trunk/libf/filtrez/mod_fft_fftw.F90 @ 986

MODULE mod_fft_fftw

#ifdef FFT_FFTW

  REAL,SAVE,ALLOCATABLE    :: Table_forward(:)
  REAL,SAVE,ALLOCATABLE    :: Table_backward(:)
  REAL,SAVE                :: scale_factor
  INTEGER,SAVE             :: vsize
  INTEGER,PARAMETER        :: inc=1
  
  INTEGER,SAVE             :: plan_forward
  INTEGER,SAVE             :: plan_backward
  
CONTAINS
  
  SUBROUTINE Init_fft(iim)
  IMPLICIT NONE
#include <rfftw.h>
    INTEGER :: iim
    REAL    :: rtmp=1.
    COMPLEX*16 :: ctmp
    INTEGER :: itmp=1
    INTEGER :: isign=0
    INTEGER :: ierr
    
    vsize=iim
    scale_factor=1./SQRT(1.*vsize)
    ALLOCATE(Table_forward(2*vsize+64))
    ALLOCATE(Table_backward(2*vsize+64))
    
!    CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr)
    
!    CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr)

    CALL rfftw_f77_create_plan(plan_forward,iim,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
    CALL rfftw_f77_create_plan(plan_backward,iim,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)
    
  END SUBROUTINE Init_fft
  
  
  SUBROUTINE fft_forward(vect,TF_vect,nb_vect)
    IMPLICIT NONE
#include <rfftw.h>
    INTEGER,INTENT(IN)  :: nb_vect
    REAL,INTENT(IN)     :: vect(vsize+inc,nb_vect)
    COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect)
    REAL                :: work(4*vsize*nb_vect)
    INTEGER             :: ierr
    INTEGER, PARAMETER :: isign=-1
  
!    CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr)
     CALL rfftwnd_f77_real_to_complex(plan_forward,nb_vect,vect, 1, vsize+inc , TF_vect, 1, vsize/2+1);  
    
  END SUBROUTINE fft_forward
  
  SUBROUTINE fft_backward(TF_vect,vect,nb_vect)
    IMPLICIT NONE
#include <rfftw.h>
    INTEGER,INTENT(IN)  :: nb_vect
    REAL,INTENT(OUT)    :: vect(vsize+inc,nb_vect)
    COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect)
    REAL                :: work(4*vsize*nb_vect)
    INTEGER             :: ierr
    INTEGER, PARAMETER :: isign=1
  
!    CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr)
    CALL rfftwnd_f77_complex_to_real(plan_forward,nb_vect,TF_vect, 1, vsize/2+1 , vect, 1, vsize+inc);  

  END SUBROUTINE fft_backward

#endif
  
END MODULE mod_fft_fftw