!-------------------------------------------------
!
!     
MODULE   setHaloParticles
PUBLIC      

Integer, PARAMETER    :: NROW  =2048,        &
                         NGRID =256,         &
                         nbyteword = 4,      &
                         NPAGE = NROW**2,    & ! # particles in a record
                         NRECL = NPAGE*6,    & 
                         Nrad  = 40,        & ! # of radial bins
                         Nbf   = 10,         & ! width of buffer for linker list
                         Nmaxpart = 1024**3+200e6, & ! max # of particles
                         Nlinker = 200,       &  ! linker-list dimension 
                         FracSearch = 5.,     &  ! 
                         Nhmax   = 150,       & ! max # of selected particles in a halo
                         Srms    = 0.2,       & ! fraction of rms velocity for bound particles 
                         nwant   = 20000,      &    ! pool of preselected particles 
                         Kdynamic= 100             ! OMP chunk                     
Real                  :: INPUT,Mhalo,Ovlim,Ovcnt
Real, DIMENSION(Nrad) :: Rad(Nrad)=0.
Real, DIMENSION(Nmaxpart)    :: Xp,Yp,Zp,VXp,VYp,Vzp,iWeight
Integer, DIMENSION(Nmaxpart) ::Lagrang
Real                  :: Box=0.,Xscale=0.,Vscale=0.,Dscale=0.,Cell=0., weightSmall
Integer               :: Np=0,Ncp=0,Ncenter=0, &
                         Nmx,Nbx,Nmy,Nby,Nmz,Nbz
Real,    ALLOCATABLE, DIMENSION (:)     :: Xpp,Ypp,Zpp,VXpp,VYpp,Vzpp,Wpp
Integer, ALLOCATABLE, DIMENSION (:)     :: LagrangP
COMMON / ROW /	        XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), &
     			VX(NPAGE),VY(NPAGE),VZ(NPAGE)
COMMON / CONTROL/      AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,  & 
                        TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,     &
                        NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, &
                        Ocurv,extras(100)
COMMON / HEADDR/       HEADER
CHARACTER*45           HEADER

Real                ::  RECDAT(NRECL),wspecies(10)
Integer             ::  lspecies(10)
Integer, ALLOCATABLE, DIMENSION (:,:,:) :: Label
Integer, ALLOCATABLE, DIMENSION (:)     :: Lst,nBins,nSelect
Real,    ALLOCATABLE, DIMENSION (:)     :: Xc,Yc,Zc,         &
                                           VXc,VYc,Vzc,      &
                                           Amc,Rmc,Vrms,Vcirc
                                           
Real,    ALLOCATABLE, DIMENSION (:,:)   :: aMr,Vvrms,GravPot
Integer, ALLOCATABLE, DIMENSION (:,:)   :: nPr,nBound
EQUIVALENCE           (wspecies(1),extras(1)), &
                       (lspecies(1),extras(11))



!$OMP THREADPRIVATE(/ROW/)
Contains
!--------------------------------------------------------------
!                    Make linker lists of particles in each cell
      SUBROUTINE List
!--------------------------------------------------------------

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
             Do i=1,Ncp
                Lst(i)=-1
             EndDo

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i,j,k)
       Do k=Nmz,Nbz
       Do j=Nmy,Nby
       Do i=Nmx,Nbx
          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(Nmx,i),Nbx)
         j=MIN(MAX(Nmy,j),Nby)
         k=MIN(MAX(Nmz,k),Nbz)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      write(*,*) ' Made first list'
      ALLOCATE(Xpp(Ncp),Ypp(Ncp),Zpp(Ncp),Stat=iStat)
      If(iStat.ne.0)Then
            write(*,*)'Not enough memory for Xpp in List. No Rearrangment'
            return
      EndIf
      ALLOCATE(VXpp(Ncp),VYpp(Ncp),VZpp(Ncp),Stat=iStat)
      If(iStat.ne.0)Then
            write(*,*)'Not enough memory for Vpp in List. No Rearrangment'
            DEALLOCATE (Xpp,Ypp,Zpp)
            return
      EndIf

      ALLOCATE(LagrangP(Ncp),Wpp(Ncp),Stat=iStat)
      If(iStat.ne.0)Then
            write(*,*)'Not enough memory for LagrangP in List. No Rearrangment'
            DEALLOCATE (Xpp,Ypp,Zpp,VXpp,VYpp,VZpp)
            return
      EndIf
   write(*,*) '  --- working on rearranging '

      nn = 0           ! main particles should be first
         Do k=Nmz,Nbz
            !write (*,*) ' re-arrange k=',k
         Do j=Nmy,Nby
         Do i=Nmx,Nbx
            jp =Label(i,j,k)
            Do while(jp.ne.0)
               If(jp.le.Np)Then
                  nn =nn +1
                  LagrangP(nn) =Lagrang(jp)
                  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
         If(nn.ne.Np)Then
               write(*,*) ' error in rearranging particles in List'
               write(*,*) ' New number of particles is=',nn
               write(*,*) ' It must be equal to old   =',Np
               Stop
         EndIf

     write(*,*) ' --- restore particles '
!                        Restore particles
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
         Do i=1,Np         ! (1,Ncp) changed
                           ! to Np for parallelization
                  Lagrang(i) =LagrangP(i)
                  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          

      DEALLOCATE (Xpp,Ypp,Zpp,VXpp,VYpp,VZpp,LagrangP,Wpp)
      write(*,*) ' re-arranged particles'
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)
             Do i=1,Ncp
                Lst(i)=-1
             EndDo

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i,j,k)
       Do k=Nmz,Nbz
       Do j=Nmy,Nby
       Do i=Nmx,Nbx
          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(Nmx,i),Nbx)
         j=MIN(MAX(Nmy,j),Nby)
         k=MIN(MAX(Nmz,k),Nbz)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      write(*,*) ' Made second list'
      Return
      End SUBROUTINE List
!--------------------------------------------------------------
!                                Find all neibhours for a center
!                                Xc,Yc,Zc - center; 
!                               
      SUBROUTINE SNeib(iHalo)
!--------------------------------------------------------------
Real,    DIMENSION(nwant) :: Rpart,Ek,Et,Ec 
Integer, DIMENSION(nwant) :: Lpart,ind
Integer :: Labpart(1)

         x  = Xc(iHalo) ; y  =  Yc(iHalo); z  =  Zc(iHalo)
         Vx1=VXc(iHalo) ; Vy1= VYc(iHalo); Vz1= VZc(iHalo)
         
         Lpart   = 0
         Labpart = 0
         Ek      = 1.e10
         Et      = 1.e10
         nn =0
         mm =0
         kbin     =Nbins(iHalo)
         RhaloMax =Rad(kbin)
         Do i=1,kbin
            nn = nPr(i,iHalo)
            mm = i
            If(Rad(i).gt.Rmc(iHalo)/2.)exit
         EndDo
         If(nn.eq.0)Then
            Nbound(1,iHalo) =0  ! no particles selected
            return
         EndIf
         Frac  = min(nwant/float(nn),1.)   ! fraction of particles,
                                           ! selected as candidates
         vthresh2 = (max(fSigma(Frac),0.2)*Vrms(iHalo))**2  ! take particles 
                                                   ! with small velocity
         Radmax = Rad(mm)
         dmax =   Radmax**2
         d0   =   Rad(1)
!         write(*,'(3i6,3g12.3)')iHalo,nn,mm,Radmax,Cell
!         write(*,'(10x,2(3x,3f8.2))') x,y,z,Vx1,Vy1,Vz1
         nn   =   0            ! counter for pre-selected particles
!                      limits for Label
         i1=INT((x-Radmax)/Cell)
         j1=INT((y-Radmax)/Cell)
         k1=INT((z-Radmax)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((x+Radmax)/Cell)
         j2=INT((y+Radmax)/Cell)
         k2=INT((z+Radmax)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
!                                        Look for neibhours
         Do k=k1,k2
         Do j=j1,j2
         Do i=i1,i2
            jp =Label(i,j,k)   !  Dark matter
 10         If(jp.ne.0.and.nn<nwant)Then
               dd =(x-Xp(jp))**2+(y-Yp(jp))**2+(z-Zp(jp))**2
!               write(*,'(10x,1p7g12.3)')Xp(jp),Yp(jp),Zp(jp),sqrt(dd),sqrt(dmax)
               If(dd.lt.dmax)Then
                  d1 =sqrt(dd)
                  ir =INT(max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1 ! bin for this particle
                  ir = min(max(1,ir),Nrad)
                  ee =(Vx1-VXp(jp))**2+(Vy1-VYp(jp))**2+(Vz1-VZp(jp))**2
                  if(ee.le.vthresh2)Then    ! take small velocity particles
                     nn = nn +1
!                     Ek(nn) = ee/2.                        ! kinetic energy
                     Et(nn) = ee/2. +GravPot(ir,iHalo) ! total relative energy

!                     Et(nn) = ee/2.  -4.333e-9*aMr(ir,iHalo)/Rad(ir) ! total relative energy
                                                           ! 4.33e-9 = GMsun/Mpc/(km/s)**2
!                     Rpart(nn) = d1
                     Lpart(nn) = Lagrang(jp)
!                     write(*,'(10x,2i8,1p5g12.3,i9)')nn,ir,d1*1000.,Ek(nn),Et(nn), &
!                               ee,Vvrms(ir,iHalo),Lpart(nn) !,nPr(ir,iHalo) 
                  endIf
               EndIf
               jp =Lst(jp)
               GoTo 10
            EndIf
         EndDo
         EndDo
         EndDo

         If(nn.le.Nhmax)Then
                       ! take all nn particles
             nBound(1:nn,iHalo) =Lpart(1:nn)
         Else
!            Ec = Et
             do j=1,Nhmax       ! take most bound
                 Labpart         =MINLOC(Et)
                 kmin            =Labpart(1)
!                 ind(j)          =kmin
                 nBound(j,iHalo) =Lpart(kmin)
                 Et(kmin)        =1.e10
             enddo
         EndIf
       nSelect(iHalo) = min(nn,Nhmax)
!       write(*,'(10x,2i7,/(i4,i9,2x,1p3g12.3))')    &
!                     iHalo,nn,(i,Lpart(ind(i)),sqrt(2.*Ek(ind(i))), &
!                     Ec(ind(i)),Rpart(ind(i))*1000.,i=1,nn)
!       write(*,'(10i10)') nBound(1:nn,iHalo)

!      If(iHalo>10)STOP
      Return
      End SUBROUTINE SNeib
!--------------------------------------------------------------
!                       Add points around the Box to take into account periodical conditions
!  It replicates points periodically and keeps those which are at  distance
!  less than Radius from the Box. 
!                        N     = number of points inside the Box
!                        Ncp = total number of points
      SUBROUTINE Points(Radius)
!--------------------------------------------------------------
      Logical    Inside
         B1  =-Radius
         B2  =Box+Radius
         B3  =Box-Radius
         Ncp =Np
      Do i=1,Np               ! Add dark matter particles
         x0 =xp(i)
         y0 =yp(i)
         z0 =zp(i)
         Inside =            (x0.gt.Radius).and.(x0.lt.B3)
         Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3))
         Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3))
         If(.not.Inside)THEN         
         ux =VXp(i)
         uy =VYp(i)
         uz =VZp(i)
         ww=iWeight(i)
            Do k1 = -1,1
               zr =z0+k1*Box
               If(zr.GT.B1.AND.zr.LT.B2)Then
               kk =k1**2
               Do j1 = -1,1
                  yr =y0+j1*Box
                  If(yr.GT.B1.AND.yr.LT.B2)Then
                     jj = kk +j1**2
                     Do i1 = -1,1
                        xr =x0+i1*Box
                        ii =jj +i1**2
                        If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)Then
                             Ncp =Ncp +1
                             If(Ncp.le.Np)Then
                               xp(Ncp)  =xr
                               yp(Ncp)  =yr
                               zp(Ncp)  =zr
                               VXp(Ncp)  =ux
                               VYp(Ncp)  =uy
                               VZp(Ncp)  =uz
                               iWeight(Ncp) =ww
                               Lagrang(Ncp) =i
                             EndIf
                        EndIf
                     EndDo
                  EndIf
               EndDo
               EndIf
            EndDo
         EndIf
      EndDo
      If(Ncp.Gt.Nmaxpart)Then
         write(*,*)' Too many points:',Ncp,' Increase Nmaxpart =', Nmaxpart
         STOP
      EndIf
      Return
      End SUBROUTINE Points
!-------------------------------- Update Statistics of pairs 
!                   
      SUBROUTINE SPairs
!---------------------------------------------------

!$OMP PARALLEL DO DEFAULT(SHARED) & 
!$OMP PRIVATE (i)                 &
!$OMP SCHEDULE (DYNAMIC,kDynamic)
      Do i=1,Ncenter
         Call Sneib(i)
      EndDo

      Return
      End SUBROUTINE SPairs
!-----------------------------------------------------------------
!               Inverse cumulative chi-2 distribution with 
!               3 degrees of freedom: chi_3(x)=frac. Need it for finding
!               velocity, which has given fraction of particles.
!                    Finds x, 
!                    which has given fraction of particles frac.
!                    fSigma = v/rms_v = [0-infty]
!                    frac   = [0-1]
!                    error of approximation is less than
!                    7.5 percent for v<4.5sigma
!                    15 percent at 6 sigma, 
!                    30 percent at 10sigma
  FUNCTION fSigma(frac)
  real, parameter :: beta =1.,thresh=0.9
  real, parameter :: C =log(1.-thresh)
   If(frac<thresh)Then
     x = frac**0.666667
     fSigma = 0.79788*x*(1.+x*(0.3+1.75*x**2))
   Else if(frac<0.99999)Then
     fSigma =2.*sqrt(1.+beta*(C- log(1.-frac)))
   Else
     fsigma = 10.
   EndIf
  return
  end FUNCTION fSigma

!-------------------------------- Manipulate halo profiles
!                   
      SUBROUTINE SmHalo(iHalo)
!---------------------------------------------------
Real, Dimension(Nrad) :: tmp
      kbin =Nbins(iHalo)
      Vv = Vrms(iHalo)             ! smooth rms profile
      Do i=1,kbin
         if(Vvrms(i,iHalo)<=1.e-5)Vvrms(i,iHalo)=Vv
      EndDo
      Vvrms(1,iHalo) = (Vvrms(1,iHalo)+Vvrms(2,iHalo))/2.
      Do i=2,kbin-1
         tmp(i) =(Vvrms(i-1,iHalo)+Vvrms(i-1,iHalo)+Vvrms(i+1,iHalo))/3.
      EndDo
      Vvrms(2:kbin-1,iHalo) = tmp(2:kbin-1)
                                   ! Gravitational Potential spherical density
                                   ! Units : (km/s)**2
      GravPot(kbin,iHalo) = -4.333e-9*aMr(kbin,iHalo)/Rad(kbin)
      Do k=kbin-1,1,-1
         GravPot(k,iHalo) =GravPot(k+1,iHalo) &
                         -4.333e-9*(aMr(k+1,iHalo)-aMr(k,iHalo)) &
                                /(Rad(k)+Rad(k+1))*2.
      EndDo

      Return
      End SUBROUTINE SmHalo

!-------------------------------- Manipulate Halos 
!                   
      SUBROUTINE SmoothHalos
!---------------------------------------------------

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (i)                 &
!$OMP SCHEDULE (DYNAMIC,kDynamic)
      Do i=1,Ncenter
         Call SmHalo(i)
      EndDo

      Return
      End SUBROUTINE SmoothHalos

!-------------------------------- Update Statistics of pairs 
!                   
      SUBROUTINE WriteData
!---------------------------------------------------

      Do i=1,Ncenter
        write(2,10)i,Xc(i),Yc(i),Zc(i),VXc(i),VYc(i),VZc(i), &
                     Amc(i),Rmc(i)*1000.,Vrms(i),Vcirc(i),     &
                     nSelect(i),(nBound(j,i),j=1,nSelect(i))
      EndDo
10    format(i7,3F9.4,3F8.1,g11.3,f8.2,2f7.1,I6,500i11)
      Return
      End SUBROUTINE WriteData

!---------------------------------------------------------
!                         Read halos catalog
!                       
      SUBROUTINE ReadHalos(iFlag)
!--------------------------------------------------------------
Character  :: Line*80, Txt*8='R(kpc/h)', Txt2*8,txt1*4

   Do i=1,2
      Read(3,'(a)')Line
      If(iFlag==0)write(*,'(a)') Line
   EndDo
      Read(3,'(a4,i4)')txt1,ii
   If(ii.ne.ISTEP)Then
          write(*,*)' Error: Snapshot and Catalog are for different time-steps'
          write(*,'(2(a,i5))') ' Catalog is for ',ii,' Particles are for ',ISTEP
          stop
   EndIf
   Do i=1,11
      Read(3,'(a)')Line
   EndDo
   Read(3,*,err=10) Rad
   Rad =Rad /1000.     ! convert to Mpch
   If(iFlag==0)Then
         write(*,*) '----- Radial bins(Mpch) ----'
         write(*,'(10f8.3)') Rad
   EndIf
   Do
     read(3,'(a)',err=20)Txt2
     If(Txt2==Txt)Exit
   EndDo

   Do
     read (3,45,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ,   &
                   aM,rr,vrmss,vcrc,Npart,Rmax,Conc,Nb
45          Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4)

        If(iStat.ne.0)Exit
        Ncenter = Ncenter +1
        If(iFlag==1)Then
              Xc(Ncenter) = X0
              Yc(Ncenter) = Y0
              Zc(Ncenter) = Z0
              VXc(Ncenter) = VvX
              VYc(Ncenter) = VvY
              VZc(Ncenter) = VvZ
              Amc(Ncenter) = aM
              Rmc(Ncenter) = rr/1000.
              Vcirc(Ncenter) = vcrc
              Vrms(Ncenter) = vrmss
              Nbins(Ncenter) = Nb
        EndIf
!         write(*,'(a,i7,3g12.3,i4)') ' N=',Ncenter,X0, Y0 , Z0,Nbins
          Do j=1,Nb
            read (3,*,err=30,end=30) r2,nn,amass,odens,vv,vc,densd
                 ! write(*,'(10x,i7,3g12.3,i4)')j,r2,amass,vc,nn 
            If(iFlag==1)Then
              nPr(j,Ncenter)    = nn
              aMr(j,Ncenter)    = amass
              Vvrms(j,Ncenter)  = vv
            EndIf                          
          EndDo
   Enddo

    If(iFlag==0)write(*,*) ' Read Halos=',Ncenter
    rewind(3)
      Return
10    write(*,*) ' Could not read radial bins'
      stop
20    write(*,*) ' Could not read text R(kpc/h)'
      stop
30    write(*,*) ' Unexpected End of file'
      stop

      End SUBROUTINE ReadHalos
!---------------------------------------------------------
!                                        Read particles
      SUBROUTINE ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs)
!         N= number of particles
!         x(min/max)s = min/max coordinates of small particles    
!----------------------------------------------------------
      Np =0

      xmin =1.e10 ; xmax =-1.e10
      ymin =1.e10 ; ymax =-1.e10
      zmin =1.e10 ; zmax =-1.e10
      xmins =1.e10 ; xmaxs =-1.e10
      ymins =1.e10 ; ymaxs =-1.e10
      zmins =1.e10 ; zmaxs =-1.e10

      Lsmall =lspecies(1)
      If(Nspecies.eq.0)Then ! old case: 1 complete set
         Npages =NROW          ! Number of pages
         N_in_last=NPAGE       ! Number of particles in last page
      Else                             ! multiple masses
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
      write (*,*) ' N_in_last=',N_in_last
!$OMP PARALLEL DO DEFAULT(SHARED)  &
!$OMP PRIVATE (IROW,In_page,iL,IN) SCHEDULE(DYNAMIC)
      DO  IROW=1, Npages         ! Loop over particle pages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            iL = NPAGE*(IROW-1)
            if(mod(IROW,10).eq.0)write (*,*)' Read page=',IROW,' Np=',iL
!				       Loop over particles
         READ  (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ
         If(iL+In_page > Nmaxpart)Stop ' Too many particles. Increase Nmaxpart'
	   DO   IN=1,In_page
                Xp(iL+IN)  =(XPAR(IN)-1.)*Xscale
                Yp(iL+IN)  =(YPAR(IN)-1.)*Xscale
                Zp(iL+IN)  =(ZPAR(IN)-1.)*Xscale
                VXp(iL+IN)  =VX(IN) *Vscale
                VYp(iL+IN)  =VY(IN) *Vscale
                VZp(iL+IN)  =VZ(IN) *Vscale
           EndDo
      EndDo
      Np      = lspecies(Nspecies)
      write(*,*) ' Done Reading Points. N=',Np
 
!$OMP PARALLEL DO DEFAULT(SHARED)  &
!$OMP PRIVATE (i)                 &
!$OMP REDUCTION(MAX:xmax,ymax,zmax,xmaxs,ymaxs,zmaxs) &
!$OMP REDUCTION(MIN:xmin,ymin,zmin,xmins,ymins,zmins)
      Do i   =1,Np 
          Lagrang(i) = i
          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)
            If(i<=Lsmall)Then
               xmins =min(Xp(i),xmins)
               xmaxs =max(Xp(i),xmaxs)
               ymins =min(Yp(i),ymins)
               ymaxs =max(Yp(i),ymaxs)
               zmins =min(Zp(i),zmins)
               zmaxs =max(Zp(i),zmaxs)
            EndIf  
      ENDDO 


      write(*,*) ' All Particles  :',Np
      write(*,*) '               x:',xmin,xmax
      write(*,*) '               y:',ymin,ymax
      write(*,*) '               z:',zmin,zmax
      write(*,*) ' Small Particles:',Lsmall
      write(*,*) '               x:',xmins,xmaxs
      write(*,*) '               y:',ymins,ymaxs
      write(*,*) '               z:',zmins,zmaxs
      write (*,*)  ' Nmaxpart=',Nmaxpart
      do i=1,1000,200
      write (*,'(i8," Test Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i)
      EndDo 
      Close(21)
      Return
      End Subroutine ReadPnt

!--------------------------------------------------
!                             Set Weights of particles for
!                             fast access
      SUBROUTINE SetWeights
!--------------------------------------------------

      Wpr =NGRID**3/FLOAT(NROW)**3
         N_particles =min(lspecies(Nspecies),Nmaxpart)
         jstart =1
         If(N_particles.le.0)Then
            write (*,*) ' Wrong number of particles:',N_particles
            write (*,*) ' Nspecies=',Nspecies,' N=',lspecies
            STop
         endif 
         Do j=1,Nspecies
            jend =min(lspecies(j),Nmaxpart)
            ww   = wspecies(j)/Wpr
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (k) 
            Do k=jstart ,jend
               iWeight(k) =ww
            EndDo
            jstart =jend
         EndDo
      weightSmall = iWeight(1)
      RETURN
      END SUBROUTINE SetWeights

!     ------------------------
      real function seconds ()
!     ------------------------
!
!     purpose: returns elapsed time in seconds
      Integer, SAVE :: first=0,rate=0
      Real   , SAVE :: t0=0.
      !real*8 dummy
      !real tarray(2)
      If(first==0)Then
         CALL SYSTEM_CLOCK(i,rate)
         first =1
         t0    =float(i)/float(rate)
         seconds = 0.
      Else
!        seconds = etime(tarray)
         CALL SYSTEM_CLOCK(i)
         seconds = float(i)/float(rate)-t0
      EndIf
      return
      end function seconds
!---------------------------------------------------
!                                  Read  current data from disk/tape,
!                                  Open files
!                                  Nrecl is the number of values in a record
!                                  Npage is the number of particles in a page
      SUBROUTINE RDTAPE(file1,file2)
!---------------------------------------------------
!                                     Open file on a tape
Character (LEN=*) :: file1,file2
      OPEN(UNIT=9,FILE=file1,  &
                      FORM='UNFORMATTED', STATUS='UNKNOWN')
!                                 Read control information
!                                 and see whether it has proper
!                                 structure
      READ  (9,err=20,end=20) HEADER,         &
                       AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, &
                       TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,     &
                       NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 &
                        ,Ocurv,extras
     WRITE (*,100) HEADER, &
                       AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, &
                       EKIN,EKIN1,EKIN2,                    &
                       NROWC,NGRID,NRECL,Om0,Oml0,hubble,   &
                       Ocurv
     WRITE (2,100) HEADER,                                 &
                       AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, &
                       EKIN,EKIN1,EKIN2,                    &
                       NROWC,NGRID,NRECL,Om0,Oml0,hubble,   &
                       Ocurv
100  FORMAT(1X,'Header=>',A45,/ &
               1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ &
                 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/        &
                 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/            &
                 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/&
                 1x,' Omega_curvature=',F7.3)
      IF(NROW.NE.NROWC) THEN
         WRITE (*,*) &
                 ' NROW in PARAMETER and in TAPE-FILE are different'
         write (*,*) ' Nrow=',NROW,' NROWC=',NROWC
      ENDIF
      IF(NGRID.NE.NGRIDC) THEN
         WRITE (*,*) &
                ' NGRID in PARAMETER and in TAPE-FILE are different:'
         write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC
      ENDIF
!                                         Open work file on disk
      write (*,*)  ' start openning PMcrs0'
 10   NBYTE = NRECL*4
      NACCES= NBYTE / nbyteword

      OPEN(UNIT=21,FILE=file2,ACCESS='DIRECT',STATUS='UNKNOWN',RECL=NACCES)

      REWIND 9
       write (*,*)  ' done openning PMcrs0'
      RETURN
 20   write (*,*) ' Error reading the header file: Abort'
      stop
      END SUBROUTINE RDTAPE
!--------------------------------------------
!                            Read current PAGE of particles  from disk
!                            NRECL - length of ROW block in words
      SUBROUTINE GETROW(IROW,Ifile)
!--------------------------------------------
         READ  (20+Ifile,REC=IROW) RECDAT
      RETURN
      END SUBROUTINE GETROW


end module SetHaloParticles
!-------------------------------------------------
!
!     
PROGRAM HaloParticles

use     SetHaloParticles 
Character*120 :: file1,file2,catname,listname,txt*8     

!     ------------------------------------------ Read data
      write(*,'(a,$)')' Enter snapshot='
      read(*,'(a)')txt
      file1 ='FILES/PMcrd_'//TRIM(txt)//'.DAT'
      file2 ='FILES/PMcrs0_'//TRIM(txt)//'.DAT'
      catname='CATALOGS/CatalogA.'//TRIM(txt)//'.DAT'
      listname='Halolist.'//TRIM(txt)//'.DAT'
      Open( 2,file=listname,Status='UNKNOWN')! short list of halos
      Open( 3,file=catname,Status='UNKNOWN')! short list of halos

       write(*,'(a)') listname 
       write(*,'(a)') catname 
       tt = seconds()
      CALL RDTAPE(file1,file2)
      CALL SetWeights

      If(extras(100).gt.0.)Then
         Box =extras(100)
      Else
         write(*,'(a,$)')' Enter box size in comoving Mpc/h ='
         read(*,*) Box
      EndIf

      WRITE (*,120)  &
                      AEXPN,ISTEP,NROW,NGRID
      write (*,*) '     Number of mass Species   =',Nspecies
      write (*,*) '     Total number of particles=',lspecies(Nspecies)
      write (*,*) '     Number of small particles=',lspecies(1)
      write (*,*) '     Box size (Mpc/h)         =  ',Box
 120      FORMAT(1X, &
           5x,'Expansion Parameter= ',F8.3,' Step= ',i5,/ &
           5x,'Nrow= ',I4,' Ngrid= ',I4)
           hsmall= hubble            ! Hubble/100.
           Xscale= Box/Ngrid        ! Scale for comoving coordinates
           Vscale= 100.*Xscale/AEXPN ! Scale for velocities
           Dscale= 2.746e+11*Om0*(Box/NROW)**3 *weightSmall ! mass scale
      write (*,'(6x,"Mass of smallest particle",/9x,  &
          "in units M_sun/h is   =",3x,g10.3)')    Dscale
      write(*,'(6x,"first species=",2g11.3)')      wspecies(1),weightSmall
      If(Nspecies.gt.1)Then
      write (*,'(6x,"Mass of largest particle",/9x, &
          "in units M_sun/h is   =",3x,g10.3)')     &
         Dscale*wspecies(Nspecies)/(NGRID**3/FLOAT(NROW)**3)
      EndIf

!     ------------------------------------------ Read data
      CALL ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs)   ! particles
      CALL ReadHalos(0)                                     ! count halos
       Rmax        = Rad(Nrad)
      CALL Points(Rmax)      ! Make more points periodically
           write (*,*)  'Points made'
           !write(*,*) '******* time=',seconds()
                                        ! set limits for linker list    

       Celloptimal = Rad(Nrad)/FracSearch
       Cell        = max(Box/Nlinker,Celloptimal)
       Nmx         = INT(xmins/Cell)-Nbf         ! boundaries for linker list
       Nbx         = INT(xmaxs/Cell)+Nbf
       Nmy         = INT(ymins/Cell)-Nbf
       Nby         = INT(ymaxs/Cell)+Nbf
       Nmz         = INT(zmins/Cell)-Nbf
       Nbz         = INT(zmaxs/Cell)+Nbf
       write(*,'(2(a,f8.4))') ' Cell for linker list(Mpch)=',Cell,' Optimal=',Celloptimal
       write(*,'(a,3(2x,2i4))') ' Dimensions for the linker list=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz
       ww = 4.*((Nbx-Nmx)*(Nby-Nmy)*(Nbz-Nmz)+Ncp)/(1024.**2)
       write (*,'(a,f9.2)') ' Requested memory for linker-list(Mb)=',ww

       ALLOCATE (Lst(Ncp),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz),Stat=iStat)
       If(iStat.ne.0)Then
            write (*,*) ' Not enough memory to allocate linker list tables'
            write (*,'(a,8i4)') ' requested limits:', &
                     Nmx,Nbx,Nmy,Nby,Nmz,Nbz
            ww = -4.*(Nmx-Nbx)*(Nmy-Nby)*(Nmz-Nbz)/1024.**2
            write (*,'(a,f9.2)') ' Requested memory(Mb)=',ww
            Stop
       EndIf
           !write(*,*) '******* time=',seconds()
      Call List             ! make linker List, Label
           !write(*,*) '******* time=',seconds()

      ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter), Stat =iStat0)
      ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter),Nbins(Ncenter), Stat =iStat1)
      ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter), Stat =iStat2)
      If(iStat0**2+iStat1**2+iStat2**2.ne.0)Then
           write(*,*) ' Not enough memory for halo headers.',iStat0,iStat1,iStat2
          stop
      EndIf
      ALLOCATE (aMr(Nrad,Ncenter),Vvrms(Nrad,Ncenter), Stat =iStat0)
      ALLOCATE (nPr(Nrad,Ncenter),GravPot(Nrad,Ncenter), Stat =iStat1)
      If(iStat0**2+iStat1**2.ne.0)Then
           write(*,*) ' Not enough memory for halo profiles.',iStat0,iStat1
          stop
      EndIf
      Ncenter = 0                                           ! initialize halo counter
      CALL ReadHalos(1)                                     ! read halos
      Close(3)
      CALL SmoothHalos

       write (*,20) Np,Ncenter,Box,Cell
       write (2,20) Np,Ncenter,Box,Cell
 20    Format(5x,'Statistics for ',i10,' points. Maxima=',i9, &
          3x,'Box(Mpc)=',f6.1,' Cell for linker list(Mpch)=',f9.5)

       ALLOCATE (nBound(Nhmax,Ncenter),nSelect(Ncenter),Stat=iStat)
       If(iStat.ne.0)Then
           write(*,*) ' Not enough memory for halo profiles.',iStat0,iStat1
          stop
      EndIf
           !write(*,*) '******* time=',seconds()
      CALL SPairs  ! Accumulate  statistics 
           !write(*,*) '******* time=',seconds()
      CALL WriteData
           write(*,*) '******* time(sec)=',seconds()
End PROGRAM HaloParticles