program lloyd
use sep3d_struct_mod
use lloyd_mod
use sep
type(sep3d) :: slow_s
real, allocatable :: slow(:),bound(:),ref_slow(:),slice(:),junk(:,:)
integer :: nref,i,iref,izs,ia,ib,imap
logical :: vel_map

call sep_init(SOURCE)
call from_param("nref",nref,4)
call from_param("vel_map",vel_map,.false.)
call init_lloyd();

call init_sep3d("in",slow_s,"INPUT")
if(sep3d_ndims(slow_s)/=3) &
  call seperr("expecting velocity to be 3-D")
call sreed("in",slow,size(slow)*4)
slow=1./slow

if(.not. vel_map) then
  call to_history("n1",nref)
  call to_history("n2",2)
  call to_history("n3",slow_s%n(3))
end if

allocate(bound(nref),ref_slow(nref), &
  junk(nref,2), slow(slow_s%n(1)*slow_s%n(2)), slice(slow_s%n(1)*slow_s%n(2)))

do izs=1,slow_s%n(3)
   call sreed("in",slow,size(slice)*4)
   slow=1./slow
   slice=slow
   iref=lloyd_go(slice,ref_slow,nref,bound)
   if(vel_map) then
     do ib=1,size(slice)
       imap=1
       do while(slow(ib) > bound(imap))
         imap=imap+1
       end do
       slow(ib)=1./ref_slow(imap)
     end do
     call srite("out",slow,size(slow)*4)
   else
     junk=-1
     junk(1:iref,1)=ref_slow
     junk(1:iref,2)=bound
     call srite("out",junk,size(junk)*4)
   end if
     
end do


















end program
