!
  !-----------------------------------------------------------------------
  subroutine cfft3(f,nr1,nr2,nr3,isign)
  !-----------------------------------------------------------------------
  !
  !  driver routine for 3d fft using fftw libraries (PG)
  !  FG - stripped from CPE, and modified 
  !
  !-----------------------------------------------------------------------
  !
  use parameters, only : DP
  implicit none
  !
  integer FFTW_FORWARD,FFTW_BACKWARD
  parameter (FFTW_FORWARD=-1,FFTW_BACKWARD=1)
  !   
  integer FFTW_REAL_TO_COMPLEX,FFTW_COMPLEX_TO_REAL
  parameter (FFTW_REAL_TO_COMPLEX=-1,FFTW_COMPLEX_TO_REAL=1)
  !   
  integer FFTW_ESTIMATE,FFTW_MEASURE
  parameter (FFTW_ESTIMATE=0,FFTW_MEASURE=1)
  !   
  integer FFTW_OUT_OF_PLACE,FFTW_IN_PLACE,FFTW_USE_WISDOM
  parameter (FFTW_OUT_OF_PLACE=0)
  parameter (FFTW_IN_PLACE=8,FFTW_USE_WISDOM=16)
  !   
  integer FFTW_THREADSAFE
  parameter (FFTW_THREADSAFE=128)
  !
  integer :: nr1, nr2, nr3, isign
  complex(DP) f(nr1*nr2*nr3)
  real(DP) :: fac
  integer :: ibid, plan(2)
  save plan
  data plan /0,0/
  !
  if       (isign.eq.1) then
     ibid=1
  else if (isign.eq.-1) then
     ibid=2
  else 
     call error('cfft3','isign unexpected',isign)
  end if
  !
  if (plan(ibid).eq.0) call fftw3d_f77_create_plan  &
     &   (plan(ibid),nr1,nr2,nr3,isign,FFTW_ESTIMATE+FFTW_IN_PLACE)
  !
  call fftwnd_f77_one(plan(ibid), f, 0)
  !
  if (isign.eq.-1) then
     fac=1.0/float(nr1*nr2*nr3)
     call DSCAL(2*nr1*nr2*nr3, fac, f, 1)
  end if
  !
  return
  end
  !