program main
  use sep
  implicit none

  real, dimension(:,:,:), allocatable :: refl, vmod
  integer                             :: iz, iz1, iz2, iz3, iz4
  integer                             :: nx, ny, nz
  real                                :: dx, dy, dz
  real                                :: ox, oy, oz
  character(len=256)                  :: vmod_, refl_

  vmod_ = "vmod"
  refl_ = "refl"

  call initpar()

  call from_param("nx",nx,1); call from_param("dx",dx,1.); call from_param("ox",ox,0.)
  call from_param("ny",ny,1); call from_param("dy",dy,1.); call from_param("oy",oy,0.)
  call from_param("nz",nz,1); call from_param("dz",dz,1.); call from_param("oz",oz,0.)

  allocate(refl(nx,ny,nz))
  allocate(vmod(nx,ny,nz))

  do iz=1, nz
     vmod(:,:,iz) = 2000 !+ (iz-1)*dz*0.3
  end do
  iz1 = 1.5+(500-oz)/dz
  iz2 = 1.5+(1000-oz)/dz
  iz3 = 1.5+(1500-oz)/dz
!  iz4 = 1.5+(4000-oz)/dz

  refl(:,:,iz1) = 1.
  refl(:,:,iz2) = 1.
  refl(:,:,iz3) = 1.
!  refl(:,:,iz4) = 1.

  call to_history("n1",nx,vmod_); call to_history("d1",dx,vmod_); call to_history("o1",ox,vmod_)
  call to_history("n2",ny,vmod_); call to_history("d2",dy,vmod_); call to_history("o2",oy,vmod_)
  call to_history("n3",nz,vmod_); call to_history("d3",dz,vmod_); call to_history("o3",oz,vmod_)

  call to_history("n1",nx,refl_); call to_history("d1",dx,refl_); call to_history("o1",ox,refl_)
  call to_history("n2",ny,refl_); call to_history("d2",dy,refl_); call to_history("o2",oy,refl_)
  call to_history("n3",nz,refl_); call to_history("d3",dz,refl_); call to_history("o3",oz,refl_)

  call sep_write(vmod, vmod_)
  call sep_write(refl, refl_)
end program
