implicit none
integer, parameter :: nc=640    !! number of cells in fine grid
integer, parameter :: nn=8      !! number of nodes
real, parameter :: lpss=nc/nn   !! height of each slab
integer, parameter :: hc=nc/2   !! half-grid length
integer, parameter :: np=hc**3  !! number of particles total
real, dimension(6,np) :: xv     !! particle list
integer :: nploc(nn)            !! each nodes particle count
real, dimension(6,np/nn*2,nn) :: xv_out  !! nodal particle list
character(len=40) fn
integer i,k
real z

!! Read in particle positions

fn='xvp.init'
write(*,*) 'Reading ',fn
open(12,file=fn,form='binary')
read(12) xv
close(12)

!! Allocate particles to correct node

nploc=0
do i=1,np
  k=int(xv(3,i)/lpss)+1
  nploc(k)=nploc(k)+1
  xv_out(:,nploc(k),k)=xv(:,i)
  xv_out(3,nploc(k),k)=xv_out(3,nploc(k),k)-(k-1)*lpss
enddo

write(*,*) 'number of particles=',sum(nploc)

!! Write out each nodes initial conditions 

do i=1,nn
  fn=''
  write(fn,'("xvp",i1,".init")') i-1
  open(12,file=fn,form='binary')
  write(12) nploc(i)
  write(12) xv_out(1:6,1:nploc(i),i)
  close(12)
enddo
end