C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE List(Box,N,Ncp)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMlists2.h'

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
             Do i=1,Np
                Lst(i)=-1
             EndDo
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,j,k)
             Do k=Nm,Nb
             Do j=Nm,Nb
             Do i=Nm,Nb
                Label(i,j,k)=0
             EndDo
             EndDo
             EndDo  
      Do jp=1,Ncp
         i=INT(Xp(jp)/Cell)
         j=INT(Yp(jp)/Cell)
         k=INT(Zp(jp)/Cell)
         i=MIN(MAX(Nm,i),Nb)
         j=MIN(MAX(Nm,j),Nb)
         k=MIN(MAX(Nm,k),Nb)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      Return
c             rearange particles to increase locality
c                       Take only the first N particles
      nn = 0
         Do k=Nm,Nb
         Do j=Nm,Nb
         Do i=Nm,Nb
            jp =Label(i,j,k)
! 10         If(jp.ne.0)Then
            Do while(jp.ne.0)
               If(jp.le.N)Then
                  nn =nn +1
                  Wpp(nn) =iWeight(jp)
	          Xpp(nn) =Xp(jp)
	          Ypp(nn) =Yp(jp)
	          Zpp(nn) =Zp(jp)
	          VXpp(nn) =VXp(jp)
	          VYpp(nn) =VYp(jp)
	          VZpp(nn) =VZp(jp)
               EndIf    
                  jp =Lst(jp)
!               GoTo10
!            EndIf
             End DO
         EndDo
         EndDo
         EndDo
c         write (*,*)  ' new number of poits=',nn,' must be=',N
c                       Take only the last  particles
         Do k=Nm,Nb
         Do j=Nm,Nb
         Do i=Nm,Nb
            jp =Label(i,j,k)
! 20         If(jp.ne.0)Then
!               If(jp.gt.N)Then
!                  nn =nn +1
!                  Wpp(nn) =iWeight(jp)
!	          Xpp(nn) =Xp(jp)
!	          Ypp(nn) =Yp(jp)
!	          Zpp(nn) =Zp(jp)
!	          VXpp(nn) =VXp(jp)
!	          VYpp(nn) =VYp(jp)
!	          VZpp(nn) =VZp(jp)
!               EndIf    
!                  jp =Lst(jp)
!               GoTo20
!            EndIf
             Do while(jp.ne.0)
               If(jp.gt.N)Then
                  nn =nn +1
                  Wpp(nn) =iWeight(jp)
	          Xpp(nn) =Xp(jp)
	          Ypp(nn) =Yp(jp)
	          Zpp(nn) =Zp(jp)
	          VXpp(nn) =VXp(jp)
	          VYpp(nn) =VYp(jp)
	          VZpp(nn) =VZp(jp)
               EndIf    
                  jp =Lst(jp)
             End Do

         EndDo
         EndDo
         EndDo
c         write (*,*)  ' new number of poits=',nn,' must be=',Ncp

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
         Do i=1,Np         ! (1,Ncp) changed
                           ! to Np for parallelization
                  iWeight(i) =Wpp(i)
	          Xp(i) =Xpp(i)
	          Yp(i) =Ypp(i)
	          Zp(i) =Zpp(i)
	          VXp(i) =VXpp(i)
	          VYp(i) =VYpp(i)
	          VZp(i) =VZpp(i)
         EndDo          

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
             Do i=1,Np
                Lst(i)=-1
             EndDo

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,j,k)
             Do k=Nm,Nb
             Do j=Nm,Nb
             Do i=Nm,Nb
                Label(i,j,k)=0
             EndDo
             EndDo
             EndDo  
      Do jp=1,Ncp
         i=INT(Xp(jp)/Cell)
         j=INT(Yp(jp)/Cell)
         k=INT(Zp(jp)/Cell)
         i=MIN(MAX(Nm,i),Nb)
         j=MIN(MAX(Nm,j),Nb)
         k=MIN(MAX(Nm,k),Nb)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      xmin =1.e10
      xmax =-1.e10
      ymin =1.e10
      ymax =-1.e10
      zmin =1.e10
      zmax =-1.e10

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i) REDUCTION(MIN:xmin,ymin,zmin) 
C$OMP+REDUCTION(MAX:xmax,ymax,zmax)
         Do i=1,N
               xmin =min(Xp(i),xmin)
               xmax =max(Xp(i),xmax)
               ymin =min(Yp(i),ymin)
               ymax =max(Yp(i),ymax)
               zmin =min(Zp(i),zmin)
               zmax =max(Zp(i),zmax)
         EndDo          
      write(*,*) '       Particles:',N
      write(*,*) '               x:',xmin,xmax
      write(*,*) '               y:',ymin,ymax
      write(*,*) '               z:',zmin,zmax
      write (*,*)  ' Nmaxpart=',Nmaxpart,Ncp
      Return
      End