# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program propaga1

  use sep, only: sep_init, sep_close, sep_read, sep_write, from_param
  use ooio
  use propag_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)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! get info from history file
  call from_history("n1",n1,"vel"); call from_history("d1",d1,"vel"); call from_history("o1",o1,"vel");
  call from_history("n2",n2,"vel"); call from_history("d2",d2,"vel"); call from_history("o2",o2,"vel");
  allocate(v(n1,n2))
  call sep_read(v,'vel')

  allocate(v(n1,n2),mod(n1,n2))
  call sep_read(mod,'in')
  call sep_read(v,'vel')
  if(any(v<=0.)) call seperr('zero or neg value in the velocity model!')
  v = v/2.

 if (exist_file("data_t")) then
     call from_history("n1",nt,"vel"); call from_history("d1",dt,"vel"); call from_history("o1",ot,"vel");
     call from_history("n2",nx,"vel"); call from_history("d2",dx,"vel"); call from_history("o2",ox,"vel");
  else
     call from_param("n1",nt,int(max(sum(v,1)/n1)*d1)); call from_param("d1",dt,0.004); call from_param("o1",ot,0.);
     call from_param("n2",nx,n2); call from_param("d2",dx,d2); call from_param("o2",ox,o2);
  end if

! !! get info from parameter
! call from_param("n1",n1); call from_param("n2",n2); call from_param("n3",n3)
! call from_param("d1",d1); call from_param("d2",d2); call from_param("d3",d3)
! call from_param("o1",o1); call from_param("o2",o2); call from_param("o3",o3)
! call from_param("n4",n4); call from_param("d4",d4); call from_param("o4",o4)

  !! write info out to history file
  call to_history("n1",nt); call to_history("d1",dt); call to_history("o1",ot)
  call to_history("n2",nx); call to_history("d2",dx); call to_history("o2",ox)
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  call sep_close()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  stat = prop_init(n1, d1, o1, n2, d2, o2, nx, dx, ox, nt, dt, ot)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Modeling
  call post3dvz_mixdom(C, v, D) ! split-step Downward continuation step
  call sep_write(D)
  deallocate(C,v,D)
  call clean_clean_prop_phsh
end program propaga1
