program propaga1

  use sep
  use propaga_mod

  implicit none

  real   , allocatable, dimension (:,:)		 :: v ! interval velocity field
  complex, allocatable, dimension (:,:)	 	 :: mod ! wavefield in Fourier domain
  real,    allocatable, dimension (:,:)	 	 :: data 
  integer                                	 :: n1, n2, nx, nt, ntn, stat
  real						 :: d1, o1, d2, o2, dx, ox, dt, ot

  call sep_init()
  call sep_begin_prog()

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! get info from history file
  call from_aux('vel','n1',n1);  call from_aux('vel','d1',d1);  call from_aux('vel','o1',o1);
  call from_aux('vel','n2',n2);  call from_aux('vel','d2',d2);  call from_aux('vel','o2',o2);

  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_aux('data_t','n1',nt);  call from_aux('data_t','d1',dt);  call from_aux('data_t','o1',ot);
     call from_aux('data_t','n2',nx);  call from_aux('data_t','d2',dx);  call from_aux('data_t','o2',ox);
  else
     ntn = int(maxval(sum(v,1)/n1)*d1)
     call from_param("n1",nt,ntn);  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

  allocate(data(nt,nx))

!  !! 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 prop_phsh(mod, v, data) ! split-step Downward continuation step
  call sep_write(data) 
  deallocate(mod,v,data)
  call clean_prop_phsh
  call sep_end_prog()
  end program propaga1
