program propaga2a

  use sep
!  use propaga2_mod
  use propaga2a_mod
  use ooio

  implicit none

  integer					 :: op,nv
  real, allocatable, dimension (:,:)		 :: v ! interval velocity field
  real, allocatable, dimension (:,:)	 	 :: mod ! wavefield in Fourier domain
  real, allocatable, dimension (:,:)	 	 :: data 
  type(axis)                         		 :: x, z, t, w, vel

  call sep_init()
  call sep_begin_prog()
op=3 !call from_param('oper',op,2) 
nv=8
call from_param('nvel',nv,3) 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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.)
  t = readaxis(1, from='data_t')
!  vel = readaxis(1,from='velax')

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

  write(0,*) z%n,z%o,z%d
  write(0,*) x%n,x%o,x%d
  write(0,*) t%n,t%o,t%d

  w = prop_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 (mod(z%n,x%n),data(t%n,x%n))
  call sep_read(mod)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Modeling
if (op.eq.1) then
  call prop_phsh(mod, v, data) ! phase-shift modeling
else if (op.eq.2) then
  call prop_splstep(mod, v, data) ! split-step modeling
else if (op.eq.3) then
  call prop_pspi(nv, mod, v, data) ! ext.split-step modeling
end if
  call sep_write(data) 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Finishing
  deallocate(mod,v,data)
  call clean_prop
  call sep_end_prog()
  end program propaga2a
