program main
  use sep
  use base_wavelt_mod
  implicit none

  integer :: nx, ny, nz, iz
  integer :: index_dfr_x1, index_dfr_x2, index_dfr_x3
  integer :: index_pln_xmin, index_pln_xmax, index_pln_z
  real    :: ox, oy, oz
  real    :: dx, dy, dz
  real    :: dfr_x1, dfr_x2, dfr_x3
  real    :: dfr_z1, dfr_z2, dfr_z3
  real    :: pln_xmin, pln_xmax, pln_z
  real    :: fm, z 
  real, dimension(:,:,:), allocatable :: modl

  call initpar()

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

  call from_param("fm", fm, 0.01)

  allocate(modl(nx,ny,nz)) 

  dfr_x1 = 22000; dfr_z1 = 10200
  dfr_x2 = 30000; dfr_z2 = 8500
  dfr_x3 = 37000; dfr_z3 = 9000
  pln_xmin = 21000; pln_xmax = 40000; pln_z = 11000
  index_dfr_x1   = 1.5 + (dfr_x1-ox)/dx
  index_dfr_x2   = 1.5 + (dfr_x2-ox)/dx
  index_dfr_x3   = 1.5 + (dfr_x3-ox)/dx
  index_pln_xmin = 1.5 + (pln_xmin-ox)/dx
  index_pln_xmax = 1.5 + (pln_xmax-ox)/dx
  index_pln_z    = 1.5 + (pln_z-oz)/dz


  modl = 0.
  do iz=1, nz
     z = oz + (iz-1)*dz
     modl(index_dfr_x1,1,iz) = wavelet(z-dfr_z1,fm,1) * 5.
     modl(index_dfr_x2,1,iz) = modl(index_dfr_x2,1,iz)+wavelet(z-dfr_z2,fm,1) * 10.
     modl(index_dfr_x3,1,iz) = modl(index_dfr_x3,1,iz)+wavelet(z-dfr_z3,fm,1) * 2.
     modl(index_pln_xmin:index_pln_xmax,1,iz) = modl(index_pln_xmin:index_pln_xmax,1,iz)+wavelet(z-pln_z, fm, 1) 
  end do

  call sep_write(modl)
end program
