program main
  use sep
  use base_wavelt_mod
  implicit none

  integer                             :: nx, ny, nz
  real                                :: dx, dy, dz
  real                                :: ox, oy, oz
  real                                :: x1, x2, x3
  real                                :: x, y, z, depth, dplnz
  real                                :: alpha, beta, decay, fm
  real, dimension(:,:,:), allocatable :: vmod, refl, bvel
  integer                             :: ix, iy, iz, nref, iref
  integer, dimension(:) , allocatable :: plnz
 
  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.)

  call from_param("x1", x1, 0.)
  call from_param("x2", x2, 0.)
  call from_param("x3", x3, 1000.)

  call from_param("alpha", alpha, 400.)

  call from_param("nref", nref, 1)

  call from_param("fm", fm, 0.5/dz)

  allocate(vmod(nx,ny,nz))
  allocate(refl(nx,ny,nz))
  allocate(bvel(nx,ny,nz))
  allocate(plnz(nref))

  do iz=1, nz; z = oz+(iz-1)*dz
  do iy=1, ny; y = oy+(iy-1)*dy
  do ix=1, nx; x = ox+(ix-1)*dx
     decay = -9*((x-x1)**2 + (y-x2)**2 + (z-x3)**2)/1000000.
     vmod(ix,iy,iz) = 2000 - alpha*exp(decay) 
  end do
  end do
  end do

  bvel = 2000.

  dplnz = (oz+(nz-1)*dz - oz)/(nref+1)
  do iref=1, nref
     plnz(iref) = 1.5+(iref*dplnz-oz)/dz
  end do

  do iz=1, nz
     z = oz+(iz-1)*dz
     do iref=1, nref
        depth = oz + (plnz(iref)-1)*dz
        refl(:,1,iz) = refl(:,1,iz) + wavelet(z-depth, fm, 1)
     end do
  end do
  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")

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

  call sep_write(vmod)
  call sep_write(refl, "refl")
  call sep_write(bvel, "bvel")
end program
