program greenf_Vz1

  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 grf_vz1_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
  real, allocatable, dimension (:,:,:)	 	 :: D ! Green_Fct.
  integer                               	 :: ix,iz,i,it
  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')
  t = readaxis(1, from='in')
  x = readaxis(2, chk_even=.true.)
  x = readaxis(2, from='vel')

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

!  write(0,*) z%n,x%n,t%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 (C(z%n,x%n))
  call sep_read(C)
  allocate (D(t%n,x%n)); D=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Modeling
  call post3dvz_mixdom(C, v, D) ! split-step Downward continuation step
  call sep_write(D) 
  deallocate(C,v,D)
  call clean_post3dvz
end program greenf_Vz1
