2D Discrete Fourier Transform using Fortran

Question

I am solving the 2D Wave Equation using Fourier Transform. The Discrete Fourier Transform (and the inverse also) is done inside the kx-loop and ky-loop. But this code runs slow, is there anyway to make it much more efficient? I guess the kx-loop, ky-loop inside the i-loop and j-loop makes it slow. I am also open for external package suggestion.

program fourier_transform_solution
implicit none

integer :: i,j,kx, ky, it
real, parameter :: PI = 3.1415926
integer, parameter :: N=100, N_t=100
integer, parameter :: N_k = N/2
real, dimension(N):: x, y
real, dimension(N,N) :: f, ut_init, u_r
real, dimension(N_t) :: t
complex, dimension(N,N) :: u, u_hat, u_hat_prev, ut_hat, ut_hat_prev
complex, dimension(N) :: wx, wy
complex, parameter :: I_complex = complex(0, 1)
real :: dx, xmin=0, xmax=3*PI, dy, ymin = 0, ymax = 3*PI
real :: d, dt = 0.001

do i=1,N_t,1
    t(i) = 0 + (i-1)*dt
end do

dx = (xmax-xmin)/(N-1); dy = (ymax-ymin)/(N-1); d = dx/(2*PI);

do i=1,N_k,1
    wx(i) = 0 + -(N_k-i)/(d*(N-1)); wx(N_k+i) = 0 + (i)/(d*(N-1))
    wy(i) = 0 + -(N_k-i)/(d*(N-1)); wy(N_k+i) = 0 + (i)/(d*(N-1))
end do

do i=1,N,1
    x(i) = xmin+ (i-1)*dx
    y(i) = ymin+ (i-1)*dy
end do

do i=1,N,1
    do j=1,N,1
        f(i,j) = 0.6*exp( -((x(i)- 1.5*PI)**2) + -((y(j)- 1.5*PI)**2)  )
        u(i,j) = f(i,j)
        u_r(i,j) = f(i,j)
        ut_init(i,j) = 0
    end do
end do

!Fourier Transform
do i=1,(N),1
    do j=1,N,1
        u_hat_prev(i,j) = 0
        ut_hat_prev(i,j) = 0
        do kx=1,N,1
            do ky=1,N,1
                u_hat_prev(i,j) = u_hat_prev(i,j) + u_r(kx, ky)*exp(-(wx(i)*x(kx) + wy(j)*y(ky))*(I_complex))
                ut_hat_prev(i,j) = ut_hat_prev(i, j) + ut_init(kx, ky)*exp(-(wx(i)*x(kx) + wy(j)*y(ky))*(I_complex))
            end do
        end do
    end do
end do


do it=2,N_t,1

    do i=1,(N),1
        do j=1, N, 1
            ut_hat(i, j) = ut_hat_prev(i, j) + dt*(-(wx(i) + wy(j))**2)*(u_hat_prev(i, j))
        end do
    end do

    do i=1,(N),1
        do j=1,N,1
            u_hat(i,j) = u_hat_prev(i,j) + dt*ut_hat(i,j)
        end do
    end do

    !Inverse Fourier
    do i=1,N,1
        do j = 1,N,1
            u(i, j)=0
            do kx=1,N,1
                do ky=1,N,1
                    u(i,j)  = u(i,j) + u_hat(kx, ky)*exp((I_complex)*(wx(kx)*x(i) + wy(ky)*y(j)))
                end do
            end do
        u_r(i,j) = (real(u(i,j)))/(N*N)
        ut_hat_prev(i,j) = ut_hat(i,j)
        u_hat_prev(i,j) = u_hat(i,j)
        end do
    end do
    print *, it
end do


end program fourier_transform_solution
```