module fft ! Fast fourier transform module ! Used to supress high frequency instabilities use prec_const use array implicit none private ! Public Functions PUBLIC :: fft_recursive, initialize_filter, suppress_high_freq contains ! In place Cooley-Tukey FFT, tsource : rosettacode.org/wiki/Fast_Fourier_transform recursive subroutine fft_recursive(x) complex(kind=dp), dimension(:), intent(inout) :: x complex(kind=dp) :: t integer :: N integer :: i complex(kind=dp), dimension(:), allocatable :: even, odd N=size(x) if(N .le. 1) return allocate(odd((N+1)/2)) allocate(even(N/2)) ! divide odd =x(1:N:2) even=x(2:N:2) ! conquer call fft_recursive(odd) call fft_recursive(even) ! combine do i=1,N/2 t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i) x(i) = odd(i) + t x(i+N/2) = odd(i) - t end do deallocate(odd) deallocate(even) end subroutine fft_recursive subroutine initialize_filter use space_grid use array use model implicit none integer :: i integer :: N integer :: center N = (ize-izs+1) center = int( ceiling( (N+1)*0.5 ) ) do i=izs,ize fft_filter(i) = 1._dp- exp( -1._dp*(i-center)*(i-center)/ (2*fft_sigma*fft_sigma) ) ! or just a step function? ! write(*,'("(", F20.15, ",", F20.15, "i )")') fft_filter(i) end do end subroutine initialize_filter subroutine suppress_high_freq(x) use space_grid use array implicit none integer :: i real(dp), dimension(izs:ize), intent(inout) :: x real(dp) :: invN do i=izs,ize fft_input(i) = x(i) end do call fft_recursive(fft_input) ! Fast Fourier Transform do i=izs,ize ! Apply filter fft_input(i) = fft_input(i)*fft_filter(i) end do call fft_recursive(fft_input) ! Invert the FFT invN = 1._dp/(ize-izs+1._dp) do i=izs,ize x(i) = real(fft_input(ize-i+1))*invN end do end subroutine suppress_high_freq end module fft