fourier2d.f90

来自「Ground state of the time-independent Gro」· F90 代码 · 共 88 行

!***
! The following code is for use with the IMSL library
!***

module fft_data
  integer, parameter :: dp=kind(1.d0)
  real(dp) :: fact
  real(dp), dimension(:,:), allocatable :: wsave
end module fft_data


subroutine fft_init(ng)
  !
  ! Initialize the Fourier transform
  !
  use fft_data
  implicit none

  integer, dimension(0:2), intent(in) :: ng

  integer :: m
  
  allocate(wsave(2*ng(0)+15,2))

  do m=1,2
     call dffti(ng(m),wsave(1,m))
  end do

  fact = 1.d0/real(ng(1)*ng(2),dp)

  return
end subroutine fft_init


subroutine fourier2D(ng,fin,fout,direction)
  !
  ! 2D Fourier transform
  !
  use fft_data
  implicit none

  integer, dimension(0:2), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2)), intent(in) :: fin
  real(dp), dimension(ng(1),ng(2)), intent(out) :: fout
  integer, intent(in) :: direction

  integer :: i,k
  real(dp), dimension(ng(0)) :: f

  if(direction >= 0) then
     ! Forward transform
     fout = fin

     ! transform in x
     do k = 1,ng(2)
        call dfftf(ng(1),fout(1,k),wsave(1,1))
     end do

     ! transform in z
     do i = 1,ng(1)
        f(1:ng(2)) = fout(i,:)
        call dfftf(ng(2),f,wsave(1,2))
        fout(i,:) = f(1:ng(2))
     end do

  else
     ! Backward transform
     fout = fin

     ! transform in x
     do k = 1,ng(2)
        call dfftb(ng(1),fout(1,k),wsave(1,1))
     end do

     ! transform in z
     do i = 1,ng(1)
        f(1:ng(2)) = fout(i,:)
        call dfftb(ng(2),f,wsave(1,2))
        fout(i,:) = f(1:ng(2))
     end do
     
     fout = fout*fact

  end if

  return
end subroutine fourier2D