program gr_Post3DVz

  use sep, only: sep_init, sep_close, sep_read, sep_write, from_param
!  use ooio, only: axis, readaxis, writeaxis, chk_dim_corresp
  use ooio
!  use gr_post3dvz_mod, only: post3dvz_init, post3dvz_mixdom, clean_post3dvz, rft_fwd_axis1, rft_bck_axis1, post3dvz_fd,post3dvz_ffd
  use gr_post3dvz_mod
!  use Shan_post3dvz_mod

  implicit none

  real   , allocatable, dimension (:,:) :: U ! input wavefield, 2-D, offset=0
  real   , allocatable, dimension (:,:) :: v ! interval velocity field
  complex, allocatable, dimension (:,:) :: C ! wavefield in Fourier domain
  complex, allocatable, dimension (:,:) :: D ! wavefield in 2D-Fourier domain
  integer                               :: ix,iz,i
  integer                               :: op
  type(axis)                            ::  x, z, t, w


  call sep_init()

  call from_param('operator',op,1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Read/write axes
  t = readaxis(1, chk_even=.true.)
  if(t%o/=0) call seperr('For migration t must start at zero (o1=0 in input)')

  z = readaxis(1, from='vel')
  x = readaxis(2, chk_even=.true.)

  call writeaxis(x, 1, axisname='x')
  call writeaxis(z, 2, axisname='z')

  write(0,*) z%n,x%n

  w = post3dvz_init(x, z, t)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Velocity stuff
  call chk_dim_corresp('vel' ,(/z, x/))
  call sep_close()
  allocate( v(z%n,x%n))
  call sep_read(v,'vel')
  if(any(v<=0.)) call seperr('zero or neg value in the velocity model!')
  v = v/2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Input data
  allocate (U(t%n,x%n))
  call sep_read(U)
  allocate (C(w%n,x%n),D(w%n,x%n)); C=0; D=0 ! nt/2 frequencies, no zero freq
  C = rft_fwd_axis1(U)

  deallocate(U)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Migrate
  depth:do iz = 1, z%n
     if (op==1) then
!       call post3dvz_mixdom(C, v(iz,:),D) ! split-step Downward continuation step
       call post3dvz_mixdom(C, v(iz,:)) ! split-step Downward continuation step
!     else if (op==2) then
!       call post3dvz_fd(C, v(iz,:),D) ! implicit finite-difference
!     else if (op==3) then
!       call post3dvz_ffd(C, v(iz,:),D) ! Foutier finite-difference
     end if
     call sep_write(sum(real(C),1)) ! imag cond sum w
    write(0,*) iz*100./z%n, '% migrated'
  end do depth

  deallocate(C,v,D)
  call clean_post3dvz
end program gr_Post3DVz
