subroutine phasemod( image, nz, nx, dz, dx, u, nt, dt, v)
integer nt,nw,nx,nz, iw,ikx,iz
real    dt,dx,dz, w,w0,dw, kx,kx0,dkx, kzkz, v, sig1,sig2,pi, signum
complex u(nt,nx), image(nz,nx), c
		pi = 3.14159265;	sig1=+1.;	sig2=-1.
call ft2axis( 0, sig2, nz, nx, image)

w0  = -pi/dt;	dw = 2.*pi/(nt*dt);	nw = nt
kx0 = -pi/dx;	dkx= 2.*pi/(nx*dx)
call zero( nw*nx, u)
do iw  = 2, nw     {  w  = w0  + (iw-1)  * dw     
do ikx = 2, nx     {  kx = kx0 + (ikx-1) * dkx
do iz  = nz, 1, -1 {
	    kzkz = w*w - (v*v) * kx*kx
	    if (kzkz >0) { c = cexp(cmplx(0., signum(w)*dz*sqrt(kzkz))) }
	    else         { c = 0.                                       }
	    u(iw,ikx) = u(iw,ikx) * c
	    u(iw,ikx) = u(iw,ikx)   +  image(iz,ikx)
	    }}}
call ft1axis( 1, sig1, nt, nx, u)
call ft2axis( 1, sig2, nt, nx, u)
return; end