program main
  use sep
  implicit none

  integer :: nx, ny, nz, nhx, nhy, nhz
  integer :: iz, ix, ihx, jx, jz, jhx
  real, dimension(:,:,:,:), allocatable :: data


  call initpar()

  call from_history("n1", nhx)
  call from_history("n2", nhy)
  call from_history("n3", nx)
  call from_history("n4", ny)
  call from_history("n5", nz)

  call from_param("ix" , ix , nx/2+1)
  call from_param("ihx", ihx, nhx/2+1)
  call from_param("iz" , iz , nz/2+1)
  allocate(data(nhx,nhy,nx,ny))


  do jz=1, nz
     call sep_read(data)
     do jx=1, nx
     do jhx=1, nhx
        if (jx==ix .and. jz==iz .and. jhx==ihx) then
           data(jhx,:,jx,:) = 1
        else
           data(jhx,:,jx,:) = 0.
        end if
     end do 
     end do
     call sep_write(data)
  end do

end program
