c
      subroutine migravelo(nt,nx,nh,dt,dx,dh,ifrmax,cda,nvel,velo,
     +                     vkzz,nz,out)
c
      parameter(nw=64,nd1=199,nw1=(2*nw)+1,d1=0.005,dmax=1.0-d1)
      parameter(pi=3.14159265,pi2=2.0*pi,piv=1./pi)
c
      real velo(nvel), vkzz(nz,nvel,4)
      complex cxx(nd1,nw1), out(nz,nx,nvel), cda(ifrmax,nh,nx)
      complex cav1, cav2, cav3, cav4
c
      nx1      = (nx/2)+1
      nh1      = (nh/2)+1
      dkx      = pi2/(nx*dx)
      dkh      = pi2/(nh*dh)
      dw       = pi2/(nt*dt)
      vkxmax   = dkx*nx
      vkhmax   = dkh*nh
      dwin     = 1./dw
      ifrmh    = ifrmax-2
c
      do id=1,nd1
         d=d1*id
         temp=pi*d
         assin=sin(temp)
         cav1 = piv*assin*cmplx(cos(temp),-1.*assin)
         do m=1,nw1
            x    =  d+nw - m
            aav  =  (nw-abs(x))/(nw*x)
            cxx(id,m)=aav*cav1 
         enddo
      enddo
c
      do iv=1,nvel
         dz       = velo(iv)*dt*0.5
         dkz      = pi2/(nt*dz)
         do ikz=2,nz
            vkz            = (ikz-1)*dkz
            vkzz(ikz,iv,1) = vkz
            vkzz(ikz,iv,2) = vkz*vkz 
            vkzz(ikz,iv,3) = vkz*vkz*vkz*vkz 
            vkzz(ikz,iv,4) = 1./(vkz*vkz)
         enddo
      enddo
c
      do k=1,nx
         do j=1,nz
            do iv=1,nvel
               out(j,k,iv) = cmplx(0.,0.)
            enddo
         enddo
      enddo
c
c     ******************************************************
c     *            f-k migration-velocity                  *
c     ******************************************************
c
      do ikx=1,nx1
         vkx  = (ikx-1)*dkx
         if(ikx.eq.nx1) vkx=vkx-vkxmax
         vkx2 = vkx*vkx
         ikxd = nx-ikx+2
         do ikh=1,nh1
            vkh  = (ikh-1)*dkh
            if(ikh.eq.nh1) vkh=vkh-vkhmax
            vkh2 = vkh*vkh
            ikhd = nh-ikh+2
c
            vmh  = abs(vkx*vkh)
            vmh2 = vmh*vmh
c
            do iv=1,nvel
               do ikz=2,nz
                  if(vmh.ge.vkzz(ikz,iv,2)) goto 405
                  gm=(1.0+vkx2*vkzz(ikz,iv,4))*(1.0+vkh2*vkzz(ikz,iv,4))   
                  sgm=sqrt(gm)
                  aa11=4.0*(vkzz(ikz,iv,3)-vmh2)
                  bb11=1.0/((vkzz(ikz,iv,2)+vkx2)+(vkzz(ikz,iv,2)+vkh2))
                  dwdkz= 0.125*velo(iv)*sgm*bb11*aa11
C
                  om   = (0.5*velo(iv)*vkzz(ikz,iv,1)*sgm)+0.00001
                  iom=om*dwin
                  if(iom.gt.ifrmax.or.iom.le.2) goto 405
                  wh=(om*dwin)-iom
                  iom1 = iom+1
                  iom2 = iom+2
c
c **** near the boundries the interpolation is not necessary ****************
c
                  cav1 = cmplx(0.,0.)
                  cav2 = cmplx(0.,0.)
                  cav3 = cmplx(0.,0.)
                  cav4 = cmplx(0.,0.)
                  if(wh.le.d1) then 
                     cav1=cda(iom1,ikh,ikx)
                     if(ikh.ne.1.and.ikh.ne.nh1) cav2=cda(iom1,ikhd,ikx)
                     if(ikx.ne.1.and.ikx.ne.nx1) then
                     cav3=cda(iom1,ikh,ikxd)
                    if(ikh.ne.1.and.ikh.ne.nh1) cav4=cda(iom1,ikhd,ikxd)
                     endif
                   endif
c
                   if(wh.ge.dmax) then
                   cav1=cda(iom2,ikh,ikx)
                   if(ikh.ne.1.and.ikh.ne.nh1) cav2=cda(iom2,ikhd,ikx)
                   if(ikx.ne.1.and.ikx.ne.nx1) then
                   cav3=cda(iom2,ikh,ikxd)
                   if(ikh.ne.1.and.ikh.ne.nh1) cav4=cda(iom2,ikhd,ikxd)
                   endif
                   endif
c
c **** far from the boundries the interpolation is necessary ****************
c
                   if(wh.gt.d1.and.wh.lt.dmax) then
                      it=iom+2-nw
                      ie=iom+nw+1
                      itt=max0(1,it)
                      iee=min0(ifrmh,ie)
                      it=it-1
                      id=int((wh/d1)+0.000001)
c
                      do m=itt,iee
                         cav1=cav1+(cda(m,ikh,ikx)*cxx(id,m-it)) 
                         if(ikh.ne.1.and.ikh.ne.nh1) then
                            cav2=cav2+(cda(m,ikhd,ikx)*cxx(id,m-it))
                         endif
                         if(ikx.ne.1.and.ikx.ne.nx1) then
                            cav3=cav3+(cda(m,ikh,ikxd)*cxx(id,m-it))
                            if(ikh.ne.1.and.ikh.ne.nh1) then
                               cav4=cav4+(cda(m,ikhd,ikxd)*cxx(id,m-it))
                            endif 
                          endif 
                      enddo
                   endif 
c
c **** intigral over kh *******************************************
c
                   out(ikz,ikx,iv)=out(ikz,ikx,iv)+dwdkz*(cav1+cav2)
                   if(ikx.ne.1.and.ikx.ne.nx1) then
                   out(ikz,ikxd,iv)=out(ikz,ikxd,iv)+dwdkz*(cav3+cav4)
                   endif
405                continue
                enddo
             enddo
          enddo
      enddo
c
      return
      end