subroutine phasemod( image, nz, nx, dz, dx, up, nt, dt, vel)

integer :: nz, nx, nt, iw,nw,ikx,iz
real :: dt,dx,dz, w,w0,dw, kx,kx0,dkx, vel, sig1,sig2,pi, w2, vkx2
complex :: up(nt,nx), image(nz,nx), cc

pi = 3.14159265; sig1 = +1.; sig2 = -1.

call ft2axis( 0, sig2, nz, nx, image)

nw = nt; w0 = -pi/dt; dw = 2.*pi/(nt*dt)
kx0 = -pi/dx; dkx= 2.*pi/(nx*dx)

up = 0., nw*nx*2)
do iw = 2, nw 
   w = w0 + (iw-1) * dw
   do ikx = 2, nx 
      kx = kx0 + (ikx-1) * dkx
      w2 = w * w
      vkx2 = vel*vel * kx*kx / 4.
      if( w2 > vkx2 )
          cc = cexp( cmplx( 0., w * dz * sqrt(1. - vkx2/w2) ))
          do iz = nz, 1, -1
             up(iw,ikx) = up(iw,ikx) * cc + image(iz,ikx)
          end do
   end do  
end do
call ft1axis( 1, sig1, nt, nx, up)
call ft2axis( 1, sig2, nt, nx, up)
return; end

# 1D Fourier transform on a 2D data set along the 1-axis
#
subroutine  ft1axis( conj, sign1, n1,n2, cx)

integer :: conj, n1,n2, i2
complex :: cx(n1,n2)
real sign1

do i2= 1, n2
   call fth( conj, sign1, 1,n1, cx(1,i2))
end do
return; end

# 1D Fourier transform on a 2D data set along the 2-axis
#
subroutine  ft2axis( adj, sign2, n1,n2, cx)
integer i1,          adj,        n1,n2
complex cx(n1,n2)
real    sign2
do i1= 1, n1
   call fth( adj, sign2, n1,n2, cx(i1,1))
end do
return; end

# FT a vector in a matrix, with first omega = - pi
#
subroutine  fth( adj,sign, m1, n12, cx)

integer :: i, adj, m1, n12
real :: sign
complex :: cx(m1,n12)
temporary complex :: temp(n12)

do i= 1, n12
        temp(i) = cx(1,i)
if( adj == 0) then
   do i= 2, n12, 2
      temp(i) =  -temp(i)
      call ftu(  sign, n12, temp)
   end do
else     
   call ftu( -sign, n12, temp)
   do i= 2, n12, 2
      temp(i) =  -temp(i)
   end do
end if
do i= 1, n12
        cx(1,i) = temp(i)
end do
return; end



