program main
  use sep
  use base_wavelt_mod
  implicit none

  real, dimension(:,:,:), allocatable :: vmod, refl
  integer                             :: nx, ny, nz
  real                                :: ox, oy, oz
  real                                :: dx, dy, dz
  real                                :: xmin, xmax, ymin, ymax, zmin, zmax
  real                                :: zpln, z, fm, vcons, vanor
  integer                             :: index_xmin, index_xmax
  integer                             :: index_zmin, index_zmax
  integer                             :: iz

  call initpar()

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

  call from_param("vcons", vcons, 2000.)
  call from_param("vanor", vanor, 2200.)

  call from_param("xmin", xmin, -500.)
  call from_param("xmax", xmax, +500.)
  call from_param("zmin", zmin, 1500.)
  call from_param("zmax", zmax, 2500.)

  call from_param("zpln", zpln, 3500.)
  call from_param("fm"  , fm  , 0.01 )

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

  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")


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

  index_xmin = 1.5 + (xmin-ox)/dx
  index_xmax = 1.5 + (xmax-ox)/dx
  index_zmin = 1.5 + (zmin-oz)/dz
  index_zmax = 1.5 + (zmax-oz)/dz

  vmod = vcons

  vmod(index_xmin:index_xmax,:,index_zmin:index_zmax) = vanor

  refl = 0.
  do iz=1, nz
     z = oz + (iz-1)*dz
     refl(:,1,iz) = wavelet(z-zpln, fm, 1)
  end do

  call sep_write(vmod)
  call sep_write(refl, "refl")

end program
