diff --git a/src/particles.f90 b/src/particles.f90
index 8fbc53b268dfe277bf253286a7b11c625a922da7..6e53b445f5622d549e6769935d757d0a9721627b 100644
--- a/src/particles.f90
+++ b/src/particles.f90
@@ -6,13 +6,14 @@
 
       if (pnd.ne.0.or.npTot.ne.0) then
         if (myid.eq.0) write(*,'(a)',advance='no') 'particle'
-        if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call partSeed
+        if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
+        call particlesToCells
         call fluidVelocity     
         if (elForceScheme.ne.2) call forcesGauss
         if (elForceScheme.ne.1) call forcesCoulomb
-        call partVelocity
-        call partCollide
-        call partTransport    !particles from neighbour gc added
+        call particlesVelocity
+        call particlesCollide
+        call particlesTransport    !particles from neighbour gc added
         call chargeDensity
         call momentumCoupling !particles from neighbour gc removed
         if (myid.eq.0) write(*,*)
@@ -339,3 +340,29 @@
 
       return
       end
+
+
+!####################################################################
+!> @author Holger Grosshans
+!> @brief locate particles on Eulerian grid
+      subroutine particlesToCells
+      use var
+
+      npic=0
+
+      do 1 n=1,np
+        npic(ip(n),jp(n),lp(n))=npic(ip(n),jp(n),lp(n))+1
+1     enddo
+
+      if (allocated(nic)) deallocate(nic)
+      allocate(nic(ii,jj,ll,maxval(npic)))
+
+      nic=0
+
+      do 2 n=1,np
+        nic(ip(n),jp(n),lp(n),count(nic(ip(n),jp(n),lp(n),:).ne.0)+1)=n
+2     enddo
+
+
+      return
+      end
diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90
index 7584075531da29fb597be98a10ca5fce68aa86f0..36c31b26111ba0c850fe7df0537d32e5e8e2b0d6 100644
--- a/src/particlesTransport.f90
+++ b/src/particlesTransport.f90
@@ -68,7 +68,7 @@
 !####################################################################
 !> @author Holger Grosshans
 !> @brief explicit first order integration of particle velocities
-      subroutine partVelocity
+      subroutine particlesVelocity
       use var      
       real(kind=pr),dimension(3) :: f_fl,f_g
       real(kind=pr) :: fact,dragdt,drag
@@ -105,7 +105,7 @@
 !####################################################################
 !> @author Holger Grosshans
 !> @brief move particles
-      subroutine partTransport
+      subroutine particlesTransport
       use var
       use mpi
       real(kind=pr) :: xpold,ypold,zpold
@@ -347,7 +347,7 @@
 !####################################################################
 !> @author Holger Grosshans
 !> @brief compute particles collisions using ray casting
-      subroutine partCollide
+      subroutine particlesCollide
       use var     
       real(kind=pr),dimension(maxnp) :: xpc,ypc,zpc
       real(kind=pr),dimension(3) :: relvelo,midlink,ncontact
@@ -356,13 +356,19 @@
                        r3n1,r3n2,rxx,projvtox,srel,projvton, &
                        closestrelvelo,fracdt, &
                        vex,vey,vez,vnx,vny,vnz
-      integer :: n1,n2,numcol,ip,jp,lp,m
+      integer :: n1,n2,numcol,ip,jp,lp,m1,m2
       
 
       numcol=0
 
-      do 1 n1=2,np
-      do 2 n2=1,(n1-1)
+      do 3 i=1,ii
+      do 3 j=1,jj
+      do 3 l=1,ll
+      do 1 m1=2,npic(i,j,l)
+      do 2 m2=1,(m1-1)
+
+      n1=nic(i,j,l,m1)
+      n2=nic(i,j,l,m2)
 
 ! condition I to check collision (same or neighbour cell):
       if ((ip(n1).lt.ip(n2)-1.or.ip(n1).gt.ip(n2)+1).or. &
@@ -476,6 +482,7 @@
 
 2     enddo
 1     enddo
+3     enddo
 
       if (myid.eq.0) write(*,'(a,i9)',advance='no') &
         ' |collide =',numcol
@@ -486,7 +493,7 @@
 !####################################################################
 !> @author Holger Grosshans
 !> @brief seed particles
-      subroutine partSeed
+      subroutine particlesSeed
       use var
       real(kind=pr),allocatable,dimension(:) :: randomx,randomy,randomz
       integer :: n,i,j,l,ip,jp,lp,npseed,npseedTot, &
@@ -572,6 +579,7 @@
         write(10,'(3(es10.2e2,x))') radp(n),q_el(n),partn(n)
 
 1     enddo
+
       
       close(10)
 
diff --git a/src/pre.f90 b/src/pre.f90
index 6573ad2f8e1d5082bd940daaa1ea85348dafaace..4fe0c22f08bba9a8832e991055b67d1b0e75c93d 100644
--- a/src/pre.f90
+++ b/src/pre.f90
@@ -191,6 +191,7 @@
       dudt(ii,jj,ll),dvdt(ii,jj,ll),dwdt(ii,jj,ll))
 
 ! particles
+      allocate(npic(ii,jj,ll))
       np=0;npp=0
       call allocateParticleArrays
 
diff --git a/src/var.f90 b/src/var.f90
index b107e4a5240c978495b26bb58ba173a1fc1fe48e..576046cc999aa137ec0eac30815b6199e48c5463 100644
--- a/src/var.f90
+++ b/src/var.f90
@@ -131,6 +131,11 @@
       integer, allocatable, dimension(:) :: &
       wcollnum, &                   !< number of particle-wall collisions
       ppcollnum                     !< number of particle-particle collisions
+      integer, allocatable, dimension(:,:,:) :: &
+      npic                          !< number of particles in a Eulerian cell
+      integer, allocatable, dimension(:,:,:,:) :: &
+      nic                           !< indices of particles in a Eulerian cell
+
 
 ! parallel
       real(kind=pr) :: timebeg,timecom,timenow,timeend