program main
  use sep
  use base_ioutil_mod
  use base_window_mod
  implicit none

  integer                             :: n1, n2, n3
  real                                :: o1, o2, o3
  real                                :: d1, d2, d3
  real, dimension(:,:,:), allocatable :: data
  real, dimension(:,:,:), allocatable :: mask
  real, dimension(:)    , allocatable :: wind
  real, dimension(:)    , allocatable :: bond
  integer                             :: nbnd, nwnd
  integer                             :: i1, i2, i3
  integer                             :: ind1, ind2, ind3, ind4, ind5, ind6
  real                                :: vcut
  call initpar()

  call from_history("n1", n1); call from_history("d1", d1); call from_history("o1", o1)
  call from_history("n2", n2); call from_history("d2", d2); call from_history("o2", o2)
  call from_history("n3", n3); call from_history("d3", d3); call from_history("o3", o3)

  call from_param("nbnd", nbnd, 0)
  call from_param("nwnd", nwnd, 20)
  call from_param("vcut", vcut, 12000.)
  allocate(data(n1,n2,n3))
  allocate(mask(n1,n2,n3))
  allocate(wind(nwnd))
  allocate(bond(nbnd))

  call hammiwin1d(bond, nbnd)
  call hammiwin1d(wind, nwnd)

  call sep_read(data)

  mask = 1.
  do i3=1, n3
  do i2=1, n2
  do i1=1, n1
     if (data(i1,i2,i3)>vcut) then
        ind1 = max(1, i1-nwnd); ind2 = min(i1+nwnd, n1)
        ind3 = max(1, i2-nwnd); ind4 = min(i2+nwnd, n2) 
        ind5 = max(1, i3-nwnd); ind6 = min(i3+nwnd, n3)
        mask(ind1:ind2,ind3:ind4,ind5:ind6) = 0.
     end if
  end do
  end do
  end do

   


  do i3=1, n3
  do i2=1, n2
     mask(      1:nbnd,i2,i3) = mask(      1:nbnd,i2,i3)*bond(1:nbnd)
     mask(n1-nbnd+1:n1,i2,i3) = mask(n1-nbnd+1:n1,i2,i3)*bond(nbnd:1:-1)
  end do
  end do
  do i2=1, n2
  do i1=1, n1
     mask(i1,i2,1:nbnd      ) = mask(i1,i2,1:nbnd      )*bond(1:nbnd)
     mask(i1,i2,n3-nbnd+1:n3) = mask(i1,i2,n3-nbnd+1:n3)*bond(nbnd:1:-1)
  end do
  end do
  call sep_write(mask)
end program
