C-------------------------------------------------
C        A.Klypin 25 Nov.1997 (NMSU)
C                            Bound-Density-Maxima code:
C                            Find  maxima of mass in spheres of radius Rsearch1 
C                            for particles in a periodic cube     of size     Box
C                           Coordinates: 0-NGRID,  Box =NGRID
C       All variables and constants are scaled to the Hubble constant of your
C       model. Only FINAL values are rescled to H=100km/s/Mpc
C       More detailed description of the algorithm and the code can be
C       found in file BDM_describe
C--------------------------------------------------   
C       Arrays Label and Lst provide quick access to neibhor particles
C       Label(i,j,k) = number of the first particle in
c                      the cell (i,j,k), if Label .ne. 0
C                    = 0 - empty cell
C       Cell         = size of a cell for Label
C       Lst(i)       = number of the next particle in the cell,
C                      where i-th was found.
C                      If Lst=0 - last point in the cell was found.
C     
      Module Structures

C       Nc - number of centers, Np - number of dm particles
C       Nrad - number of shells for halos
C
C                       Configuration parameters
      INCLUDE 'PMparameters.h'
      PARAMETER (Nrad    =    60)      ! Number of shells for halo profile
      PARAMETER (FracSearch = 7. )      ! Factor of Rmax for small halos
      PARAMETER (Np      = Nmaxpart)    ! Max Nparticles
      PARAMETER (Npp      = Np      )    ! =1 or Np  If Npp=Np rearrange particles
                                        !    to increase locality. Uses lots of RAM
      PARAMETER (Nrad_lim =Nrad)
      PARAMETER (iDebug  =     0)       ! debugging parameteter
      PARAMETER (iHelp   =     0)       ! help

C                       Parameters to tune halo finder
      PARAMETER (dClose  =    15. )         ! Only one seed within distance
                                           ! dRdubl/dClose will be kept; dRdubl is distance
                                           ! to remove duplicates.
	                                   ! This is used in CATALOG.
                                           ! It is maximum of the radius of the first bin
                                           ! Rad(1) and  dRdubl/dClose

      PARAMETER (nFirst  =    5 )          ! minimum number of particles inside 
                                           ! radius for chosing first bin for overdensity
                                           ! and mass profiles
      PARAMETER (mFirst  =    30)          ! minimum number of particles in side 
                                           ! radius for chosing first bin for radial velocity
					   ! profile Vrad. This is for <Vradial>=0.5Vcirc test,
                                           ! which used to truncate outer radius of halos

      PARAMETER (dRmaxrot =    2.)         ! defines minimum concentration for a halo when
                                           ! deciding on escape velocity.
                                           ! Radius of maximum circular velocity is
                                           ! limited to be smaller than Rhalo/dRmaxrot.
                                           ! Concentration is larger than dRmaxrot*2.15

      PARAMETER (nmin     = 30 )           ! minimum number of particles, which defines
                                           ! radius of a sphere used for bulk velocity of a halo
 
      PARAMETER (iCentr   = 2 )           ! minimum bin number for  bulk velocity of a halo:
                                           ! Radial bin for the velocity is
                                           ! the largest of iCentr and radius containing nmin
 
c                       Optimization parameters
      PARAMETER (Nm      =    -5)          ! define size of linker mesh
      PARAMETER (Nlinker =    400)         ! 
      PARAMETER (Nb      =    Nlinker-Nm)  ! define size of linker mesh
      PARAMETER (mDynamic = 500 )     ! dynamical scheduling for OpenMP
      PARAMETER (kDynamic = 100 )     ! these parameters define length of
                                     ! chunk of halos to assign for 
                                     ! one OMP thread
      Real, DIMENSION(Np) :: Xp,Yp,Zp,VXp,VYp,Vzp
      Integer, DIMENSION(Np) :: Lst

      Real, ALLOCATABLE, DIMENSION(:) :: Xpp,Ypp,Zpp,
     &                          VXpp,VYpp,VZpp,Wpp
      Real, SAVE, ALLOCATABLE, DIMENSION(:,:) :: Wxc,Wyc,Wzc,
     &                                        Vrmc,Rrc,Vrada
      Integer, SAVE, ALLOCATABLE,DIMENSION(:,:) :: Nbc,NbcM,
     &                                            NbcG,NbcS
      Real, ALLOCATABLE, DIMENSION(:) :: Xm,Ym,Zm,Amc,Rmc,Vrm,RmaxV
      Integer, ALLOCATABLE, DIMENSION(:) :: iRadius,Lab


      COMMON /Lims/ Nmx,Nbx,Nmy,Nby,Nmz,Nbz
      Integer,  ALLOCATABLE :: Label(:,:,:)

      COMMON /PNT/ RadSr1,RadSr2,aNvir,Cell,
     +             Rad(Nrad),Rad2(Nrad),weightSmall
!$OMP THREADPRIVATE(/ROW/)

      End Module Structures

!-------------------------------------------------------------------
      Program BDM

      use Structures

      Real             INPUT,Mhalo
      Real             Overdens(Nrad),Mass(Nrad),MnMass(Nrad)
      COMMON /OVER/ Ovlim,Ovcnt
      DIMENSION      Radsc(Nrad)


C     ------------------------------------------ Read data
      If(iHelp.eq.1)Call Help
      CALL RDTAPE
      If(iDebug.eq.1)write (*,*) ' RDTAPE is done'
      CALL SetWeights

      If(extras(100).gt.0.)Then
         Box =extras(100)
      Else
         Box  =INPUT(' Enter box size in comoving Mpc/h =')
      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.
           Box   = Box /hsmall       ! Scale to real Mpc          
           Xscale= Box/Ngrid        ! Scale for comoving coordinates
           Vscale= 100.*hsmall*Xscale/AEXPN ! Scale for velocities
           Dscale= 2.746e+11*hsmall**2*(Box/NROW)**3 *weightSmall ! mass scale
           aMassOne= Om0*Dscale*hsmall
      write (*,'(6x,"Mass of smallest particle",/9x,
     +    "in units M_sun/h is   =",3x,g10.3)') 
     +      aMassOne
      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)') 
     +   Om0*Dscale*wspecies(Nspecies)*hsmall/(NGRID**3/FLOAT(NROW)**3)
      EndIf

      If(iHelp.eq.1)write(*,140)
 140  format(10x,'Give average overdensity inside halo radius',/
     +10x,'= 178 for CDM, Omg_0=1; =340 for LCDM Omg_0=0.3')

      Ovlim   = OverdenVir()        ! set virial overdensity
      Ovcnt   = Ovlim   ! Min. Center Overdensity for Halos 
                                   ! scale factor to get number of virial particles
      aNvir  = 4.1888*Ovlim*(NROW/float(NGRID))**3/weightSmall
      Amassl  = AmassOne*5.  ! Minimum halo mass in Msun/h 
      If(iHelp.eq.1)write(*,150)
 150  format(10x,'Code searches for maxima of mass inside spheres',/
     +10x,'with radii ranging from small radius (second question) to',/
     +10x,'large radius (first question). The radii may be equal.',/
     +10x,'The search radii should be larger than the resolution ',/
     +10x,'of the run. Code will fail if there is too few particles ',/
     +10x,'inside the spheres.',/
     +10x,'The smaller radius also is used for the search of ',/
     +10x,'more accurate position of density peak',/
     +10x,' Never give the value of the small radius which is ',/
     +10x,'smaller than 1.5-2 highest spatial resolution')

      Rsearch1=INPUT('Enter comoving search radius(Mpc/h)       => ')
      Rsearch2=INPUT('Enter smaller radius(Mpc/h) of final halos=> ')
      Rminhalo=0.0001   ! min.radius for halos(Mpc/h) 
      Fract   = 1.            ! fraction of DM particles (1/4,1/2,1)
      Fract   =Min(Max(Fract,0.),1.0)
                             
      If(iHelp.eq.1)write(*,160)
 160  format(10x,'Unbounding particles.',/
     +   10x,' =1 for removing unbound particles; =0 or >10',
     +    ' for no unbounding')
      Toohot  =INPUT('Enter rejection velocity limit (V/Vescape)=> ')
      If(iHelp.eq.1)write(*,170)
 170  format(10x,'Halos with close masses, close velocities and with',/
     +11x,'distances less than some value (Mpc/h) are considered ',/
     +10x,'as the same halo. Only the largest is kept.')

      dRdubl  =INPUT('Distance to check for Velocity duplicates => ')
      fVdubl  =0.15    ! Define duplicates if  (v1-v2)/Vrms   <   
      if(iHelp.eq.1)Then
         write (25,*) Ovlim,'  Overdensity Threshold for Halos'
         write (25,*) Rsearch1,' comoving search radius(Mpc/h)'
         write (25,*) Rsearch2,' smaller radius(Mpc/h) of final halos'
         write (25,*) Toohot,' rejection velocity limit (V/Vescape)'
         write (25,*) dRdubl,
     +             '  Distance to check for Velocity duplicates'     
      EndIf
C     ------------------------------------------ Open files
      Open( 2,file='Catshort.DAT',Status='UNKNOWN')! short list of halos
      Open(20,file='Catalog.DAT',Status='UNKNOWN') ! catalog of halos
C     ------------------------------------------ Read data
      CALL ReadPntP(N,Fract,xshift,yshift,zshift)                     ! particles
      If(iDebug.eq.1)write (*,*)' done Reading Particles'
C       CALL ReadMax(N,Fract)                     ! particles
c      Rsearch1=Rsearch1/AEXPN              ! scale to proper distances
c      Rsearch2=Rsearch2/AEXPN
c      dRdubl  =dRdubl /AEXPN
      Boxx =Box    !  store box size in Mpc

         WRITE (2,100) HEADER,
     +               AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +               EKIN,NROW,NGRID,Om0,Oml0,hubble
         WRITE (20,100) HEADER,
     +               AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +               EKIN,NROW,NGRID,Om0,Oml0,hubble
100      FORMAT(1X,'Header=>',A45,/
     +         ' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +         ' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,
     +         ' Nrow=',I4,' Ngrid=',I4,/' Omega_0=',f6.3,
     +         ' Omega_L=',f6.3,' h=',f6.3)
C     ------------------------------------------ Constants ---------
           Omcold= Om0 
           Omhot = 0.
           Amassl= Amassl/hsmall     ! Scale to real masses
           Radius= Rsearch1/hsmall  ! Comoving Search radius in real Mpc
           Rminhalo=Rminhalo/hsmall
           dRdubl= dRdubl/hsmall
C                         Dscale is a factor used to convert particles to mass 
C                                    =  mass of a particle for Omega=1 model
C                         Xscale is a factor to convert PM coordinates into Mpc
C                                    = Box / Ngrid = length of a PM cell in Mpc
C                        Vscale is a factor to convert PM pomenta to V (km/sec)
C                                    = the Hubble velocity at the distance 
C                                         of 1 PM cell
C                the factor AEXPN is needed to go from momenta to pec.velocity
      If(iDebug.eq.1)write (*,*) ' Vscale=',Vscale, ' X=',Xscale,
     +        ' Box(Mpc)=',Box
           Radius=Radius/Xscale   ! comoving radius in cell units
           RadSr1 =Rsearch1/hsmall/Xscale ! max Min search radii
           RadSr2 =Rsearch2/hsmall/Xscale           
           Box   =NGRID           ! all internal variables are in Grid units
           Riter =0.0001          ! Error for iteration of radius in cell units
           Rmax  =Radius          ! Treat maxima Identical  if dR<Rmax
           Cell   = Box/Nlinker
c           Cell   = 0.25
               If(iDebug.eq.1) write (*,*) ' Cell size=',Cell,
     +               ' Search Radius in cell units=',Radius
         DtoNscale = Om0*Dscale/Fract ! scale dark matter particles
cc         ANlimit =!!!Amassl/(DtoNscale)         ! and this
         ANlimit =2.*aNvir*weightSmall ! factor to calculate
c         ANlimit =8.*aNvir*weightSmall ! factor to calculate
                                      ! central overdensity
                                      ! It will be 8 times larger than
                                      ! virial overdensity
           R0    =0.004
           Rnrad =R0*((Nrad-1)/4.+1.)**2 
           write (*,'(/" The first and the last bins have radii:",
     &                   2f8.4," Mpc/h, comoving")')  R0,Rnrad
           write (*,*) '   If you want to change shells, enter a factor'
           write (*,'(A,$)') '    or 1 to keep current radii : '
           read  (*,*) factor
         if(iHelp.eq.1)write(25,*)factor,'  scale factor for bins'
           If(abs(factor-1.).ge.1.e-3.and.factor.gt.0.)Then
              R0    =0.004 * factor
              Rnrad =R0*((Nrad-1)/4.+1.)**2           
              write (*,*) '   The first and the last radii for halos'
              write (*,*) '   will be',R0,Rnrad,' Mpc/h, comoving'
           EndIf
           Do ir=1,Nrad                            ! Radii for halo shells
              Rad(ir) =R0/hsmall*((ir-1)/4.+1.)**2/Xscale ! comoving
              Rad2(ir)=Rad(ir)**2
           EndDo
           write (*,'(3(a,g12.3))') 
     +                ' Final Radii(Mpc/h): 1st=',Rad(1)*Xscale*hsmall,
     +                ' last=',Rad(Nrad)*Xscale*hsmall,' Cell=',Cell

c------------------------------------------ Prepare Data
           RadOld =Radius
           Radius  =Rad(Nrad)            ! Make radius temporarily bigger
      CALL Points(Box,N,Ncp,Radius)      ! Make more points periodically
                     If(iDebug.eq.1)
     +               write (*,*) ' Nparticles in periodic set=',Ncp,
     +                           '( Maximum =',Np,')'
           Radius =RadOld                ! Restore radius
           write (*,*)  'Points made'

      Call LimList(Box,N,Ncp)            ! find limits for linker-list
      write(*,'(a,6i6)')' Limits=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz
      write(*,'(a,g12.4)')' Memory required (Mb)=',
     &                    4.*(Nbz-Nmz)*(Nby-Nmy)*(Nbx-Nmx)/1024**2
      ALLOCATE (Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz),Stat =iStatus)
         If(iStatus.ne.0)Stop ' not enough memory for Label'
         write (*,*) ' Done allocation'
      Call List(Box,N,Ncp)               ! make linker List, Label
         write (*,*)  'List made'
      CALL ReadCent(N,Ncentr,Boxx,0)     ! count initial centers
      Nc = Ncentr
      ALLOCATE (Xm(Nc),Ym(Nc),Zm(Nc),Amc(Nc),Rmc(Nc),Stat=iStatus)
          If(iStatus.ne.0)Stop ' not enough memory for Centers'
      CALL ReadCent(N,Ncentr,Boxx,1)     ! set initial centers
          write (*,*) ' Nparticles=',N,' Ncenters=',Ncentr
C----------------------------- Find maxima
      ALLOCATE (Lab(Ncentr))
      CALL Pairs(Box,N,Ncentr,Riter)       ! find maxima 
      DEALLOCATE (Lab)
      write (*,16)  Amassl*hsmall,ANlimit,Rminhalo*hsmall,
     &          dRdubl/dClose*hsmall*1.e3,Rsearch1,Rsearch2,Toohot,Cell
      write (2,16)  Amassl*hsmall,ANlimit,Rminhalo*hsmall,
     &          dRdubl/dClose*hsmall*1.e3,Rsearch1,Rsearch2,Toohot,Cell
      write (20,16) Amassl*hsmall,ANlimit,Rminhalo*hsmall,
     &          dRdubl/dClose*hsmall*1.e3,Rsearch1,Rsearch2,Toohot,Cell
             If(iDebug.eq.1)write (*,*) ' Om0=',Om0,' Dscale=',Dscale
 16   format(5x,'Small halos with Mass(M_sun/h) less than',
     .       g10.3,'(N_eff<',g11.3,') were removed',/
     .       5x,'Minimum Radius of Halo(Mpc/h)=',f7.4,
     &          ' Minimum distance beteen halos(kpc/h)=',f7.4,/
     .       5x,'Comov.Search radius(Mpc/h)=',f7.4,
     .         ' Second radius(Mpc/h)=',f7.4,/
     .       5x,'Particles with V>',f5.1,'V_escape were removed',
     .       ' Cell (in grids) for linker list=',f8.3)
           write (*,*)  'Pairs made'
      CALL SmallRem(ANlimit,Ncentr)  ! Remove small halos
           write (*,*)  'Small Rem made. Ncentr=',Ncentr
C                           Final touch: if Rsearch1 and Rsearch2 are different
C                           then  improve position of halos by shrinking
C                           the search Radius to Rsearch2 in Nit steps

      Rsearch2 = Rsearch2/2.          ! reduce the small radius 
c      RadSr2 =Rsearch2/hsmall/Xscale  ! for more accurate search of center
      ALLOCATE (Lab(Ncentr))
      IF(ABS(Rsearch1-Rsearch2).gt.0.001*Rsearch1)Then
         Nit =3
         hNit =(Rsearch1/Rsearch2)**(1./Nit)
         If(iDebug.eq.1)write (*,*) ' Nit =',Nit,' Factor:',hNit
        Do ii =1,Nit
            n1 = 0
            r1  =0.
            Do i=1,Ncentr
      	     If(Rmc(i).gt.RadSr2)Then
                    n1     = n1 +1
                    Rmc(i) = max(Rmc(i)/hNit,RadSr2)
                    r1     = r1 +Rmc(i)
                 Endif
             Enddo 
          write (*,*)'Iterate:',
     &               ' Mean Comoving Radius for halos(Mpc/h)=',
     &                   r1/(max(n1,1))*Xscale*hsmall,' n=',n1
          CALL Pairs(Box,N,Ncentr,Riter) 
        EndDo
        Radius=Rsearch1
      EndIf

       write(*,*) '   ---- do catalog --'
      dRdubl =dRdubl/Xscale                ! scale to cell units
      CALL CATALOG(Box,Ncentr,dRdubl)      ! Remove duplicates
      DEALLOCATE (Lab)

       write(*,*) '   ---- end catalog --'
      DEALLOCATE (Label)
      Celloptimal = (Rad(Nrad)/FracSearch)
      CellOld     = Cell
      Cell = max(Celloptimal,float(NGRID)/float(Nlinker))
         ff   = CellOld/Cell      ! scaling factor to new limits
         Nmx  = Nmx*ff
         Nmy  = Nmy*ff
         Nmz  = Nmz*ff
         Nbx  = Nbx*ff
         Nby  = Nby*ff
         Nbz  = Nbz*ff
      
      ALLOCATE (Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz), Stat =iStat)
      If(iStat.ne.0)Then
           write(*,*) '--- Failed to allocate Label --'
           write(*,*) '--- X:',Nmx,Nbx 
           write(*,*) '--- Y:',Nmy,Nby 
           write(*,*) '--- Z:',Nmz,Nbz 
           stop '--- Stop'
      EndIf
       write(*,*) '   ---- do List --'
      Call List(Box,N,Ncp)             ! make linker List, Label

       write (*,20) N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell
       write (2,20) N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell
       write (20,20)N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell
20         Format(5x,'Statistics for ',i9,' points. Maxima=',i8,/
     .         5x,'Box(Mpc)=',f6.1,
     .            5x,'Overdensity Limit=',f6.1,
     .             ' Cell for linker list(grids)=',f8.3)
          t0 =seconds()
          write (*,*) ' time  before Spairs(sec)=',t0
      ALLOCATE (Wxc(Nrad,Ncentr),Wyc(Nrad,Ncentr),Wzc(Nrad,Ncentr),
     &            Stat=iStat)
      If(iStat.ne.0)Stop ' not enough memory for Wxxx'
      ALLOCATE (Vrmc(Nrad,Ncentr),Rrc(Nrad,Ncentr),Vrada(Nrad,Ncentr),
     &            Stat=iStat)
      If(iStat.ne.0)Stop ' not enough memory for Vrmc'
      ALLOCATE (Nbc(Nrad,Ncentr),NbcM(Nrad,Ncentr),NbcG(Nrad,Ncentr),
     &            NbcS(Nrad,Ncentr),Stat=iStat)
      If(iStat.ne.0)Stop ' not enough memory for Wxxx'
      CALL SPairs(Box,N,Ncentr)  ! Accumulate  statistics 

          t1 =seconds()
          write (*,*) ' time for Spairs (sec)=',t1-t0
      Amassl  = AmassOne*9.  ! Minimum halo mass in Msun/h 
C-----------------------------      Auxiliary Settings
      nfalse =0
      Do ir=1,Nrad                  ! all scaled to hubble, not h^(-1)
         Radsc(ir) =Rad(ir) *Xscale ! scale radii to comoving Mpc
         MnMass(ir)=Om0*4.188*2.746e+11*hubble**2*Radsc(ir)**3 
c                                   mass at the mean density
      EndDo
      write (*,12) (Radsc(i)*hsmall*1000.,i=1,Nrad)
      write (2,12) (Radsc(i)*hsmall*1000.,i=1,Nrad)
      write (20,12) (Radsc(i)*hsmall*1000.,i=1,Nrad)
 12   Format(5x,'Radii of shells in comoving kpc/h:',/20(4x,10f8.2/))
C---------------------------- Process Halos : Unbind, Dublicates 
      ALLOCATE (Vrm(Ncentr),RmaxV(Ncentr),iRadius(Ncentr),Stat=iStat)
          If(iStaus.ne.0)Stop ' not enough memory for Centers'
      CALL CLEAN(Ncentr,Xscale,Vscale,Dscale, 
     +            Fract,Radsc,MnMass,Toohot,Rminhalo)
      write (*,*)  '  CLEAN: Ncenters=',Ncentr
          t2 =seconds()
          write (*,*) ' time for CLEAN(sec)=',t2-t1

      IF(fVdubl .gt.0. .and. dRdubl.gt.0.)Then    !  Remove VelDublicates
      ALLOCATE (Lab(Ncentr)) 
         CALL RemVelDublA(Box,Ncentr,dRdubl,fVdubl,Xscale,Vscale)
      DEALLOCATE (Lab) 
         write (*,25) Ncentr,dRdubl*Xscale*hsmall ,fVdubl
         write (2,25) Ncentr,dRdubl*Xscale*hsmall ,fVdubl
         write (20,25)Ncentr,dRdubl*Xscale*hsmall ,fVdubl
 25       Format(5x,'Number of centers=',i6,
     &                   2x,'Velocity Duplicates: dR(Mpc/h)<',f6.3,
     &                   2x,'dV/Vrms<',f6.3)
       ENDIF
C---------------------------------- Find largest halo
c       NpMax =0
c       inpMax =0
c      Do i=1,Ncentr        
c         write (*,*) ' halo=',i,NbcS(1,i),NbcS(2,i),NbcS(3,i)
c         If(NbcS(2,i).gt. NpMax)Then
c            inpMax =i
c            NpMax =NbcS(2,i)
c         EndIf
c      EndDo
c      If(inpMax.eq.0.or.NpMax.eq.0)Write (*,*) ' Errrror'
c      Do i=1,Ncentr         
cc         If(i.ne.inpMax) iRadius(i) =0
c         If(Vrm(i).lt. 100.) iRadius(i) =0
c      EndDo
c      iRadius(inpMax) =Nrad
C------------------------------- Print Results
      iCount = 0
      write ( 2,13)
      write (20,13)
 13   Format(6x,'Coordinates Mpc',12x,'Velocity km/s',7x,
     & 'Mass     Radius Vrms(3D) Vcirc Npart(<Rad) Rmax(kpch)',
     & ' Conctr Bins')
       write (20,14)
 14   Format('R(kpc/h)',5x,' N(<R)    M/Msun   Overdens   Vrms  Vcirc',
     &                  ' Dens(Msun/Mpc^3) Nsmall Nlrg Nbig Vrad(km/s)')
      Do i=1,Ncentr          ! Scale all to physical units
         Xc    =(Xm(i)-xshift)  *Xscale
         Yc    =(Ym(i)-yshift)  *Xscale
         Zc    =(Zm(i)-zshift)  *Xscale
         irad  =iRadius(i)
         Mhalo =Amc(i)
         Rhalo =Rmc(i)
         RmaxH =RmaxV(i)
         Vcirc =Vrm(i)

         iS      =NbcS(irad,i)   ! small-size particles
         iM      =NbcM(irad,i)   ! Medium-size particles
         iG      =NbcG(irad,i)   ! Giant  particles
         Nnp   =INT(Mhalo/(Dscale*Om0))
         Do ir =1,Nrad
            Summass     = Om0*Nbc(ir,i)/Fract
            Mass(ir)    = Dscale*Summass
            Overdens(ir)= Mass(ir)/MnMass(ir)
         EndDo
         im =Nbc(Nrad,i)
             Do ir =Nrad-1,1,-1
                if(Nbc(ir,i).ne.im)Then
                   iLast =ir
                   goto 220
                EndIf 
             EndDo
 220         iLast = min(Nrad,iLast+1)
       If(irad.eq.0.or.Rhalo.le.Rminhalo
     .      .or.Mhalo.le.Amassl)Then
              nfalse=nfalse+1
            irr =Max(irad,2)         ! don't take V from the first bin
            VXc =Wxc(irr,i)   *Vscale
            VYc =Wyc(irr,i)   *Vscale
            VZc =Wzc(irr,i)   *Vscale
            Nnp =INT(Mhalo/(Dscale*Om0))
       EndIf
         If((irad.ne.0.and.Rhalo.gt.Rminhalo)
     .                .and.Mhalo.gt.Amassl
     .                )Then
            irr = Max(irad,2)
            irr2 = 2  !Max(irad,2)         ! don't take V from the first bin
 10         if(Nbc(irr2,i).lt.nmin.and.irr2.lt.iCentr)Then
               irr2 =irr2+1
               goto 10
            EndIf 

            VXc =Wxc(irr2,i)   *Vscale
            VYc =Wyc(irr2,i)   *Vscale
            VZc =Wzc(irr2,i)   *Vscale
            Nnp =INT(Mhalo/(Dscale*Om0))

             write (20,45) Xc*hsmall,Yc*hsmall,Zc*hsmall,
     &                      VXc,VYc,VZc,Mhalo*hsmall,
     &                      Rhalo*hsmall*1000.,Vrmc(irr,i),
     &                      Vcirc,Nnp,
     &                      RmaxH*hsmall*1000.,Rhalo/RmaxH*2.15,iLast

             write ( 2,46) Xc*hsmall,Yc*hsmall,Zc*hsmall,
     &                      VXc,VYc,VZc,Mhalo*hsmall,
     &                      Rhalo*hsmall*1000.,Vrmc(irr,i),
     &                      Vcirc,Nnp,
     &                      RmaxH*hsmall*1000.,Rhalo/RmaxH*2.15

            Do ir=1,iLast
               If(ir.eq.1)Then
                  Vvc =Vrmc(ir,i)
                  Volume =4.188*Radsc(ir)**3
                  DenDm=Mass(ir)/Volume
                  If(Nbc(ir,i).le.2)Then ! too few particles
                      Radmean =Radsc(ir)/2.**0.3333
                  Else         ! take real mean radius
                      Radmean  =Rrc(ir,i)*Xscale/Nbc(ir,i)
                  EndIf
                  Vcr    =sqrt(Mass(1)/Radsc(1)/AEXPN)*6.58e-5
               Else              ! remove contribution from inner rad
                 Vvc =sqrt(MAX((Vrmc(ir,i)**2*Nbc(ir,i) -
     .              Vrmc(ir-1,i)**2*Nbc(ir-1,i))/
     .              max(Nbc(ir,i) -Nbc(ir-1,i),1),1.e-10))
c                  Vvc =(Vrmc(ir,i)*Nbc(ir,i) -  !radial velocity
c     .              Vrmc(ir-1,i)*Nbc(ir-1,i))/
c     .              max(Nbc(ir,i) -Nbc(ir-1,i),1)
                  Volume =4.188*(Radsc(ir)**3- Radsc(ir-1)**3) 
                  DenDm  =(Mass(ir)-Mass(ir-1)) /Volume
                  If(Nbc(ir,i) -Nbc(ir-1,i).le.2)Then ! too few particles
                     Radmean=((Radsc(ir)**3+Radsc(ir-1)**3)/2.)**0.333
                  Else         ! take real mean radius
                     Radmean=(Rrc(ir,i)-Rrc(ir-1,i))*Xscale/
     .                        (Nbc(ir,i)-Nbc(ir-1,i))
                  EndIf
                  Vcr    =sqrt(Mass(ir)/Radsc(ir)/AEXPN)*6.58e-5
                EndIf
               write (20,40) Radmean*hsmall*1000.,
     +                       Nbc(ir,i),Mass(ir)*hsmall,Overdens(ir),
     +                       Vvc,Vcr,DenDm/hsmall**2,NbcS(ir,i),
     +                       NbcM(ir,i),NbcG(ir,i),
     &                       Vrada(ir,i)*Vscale
            EndDo  ! end loop ir
45          Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4)
46          Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2F7.1,i9,f8.1,f8.2,i3)
47         Format(F8.3,2F8.3,g11.3,f8.2,i9,i8,i7,i5)
40          Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3)
41          Format(f7.4,I7,I6,g11.3,g10.3,2f7.1,f6.3,g10.3,4G10.3)
         Endif
      EndDo    !  end loop Ncentr
      write (*,50) Ncentr-nfalse,nfalse
c      write (2,50) Ncentr,nfalse
 50   Format('---------- Final Catalog has ',i7,' halos ------------',
     .      /'                    ',i7,' were too puffy or too small') 
      Stop
      End Program BDM
C--------------------------------------------------------------
C                        virial overdensity for cosmological model
C                        at different expansion parameter AEXPN
      Function OverdenVir()
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      

      xx =-(1.-Om0)*AEXPN**3/(Om0+(1.-Om0)*AEXPN**3)
      Ovlim =(178.+82.*xx-39.*xx**2)/(1.+xx)
      write (*,*)  'Overdensity Delta =',Ovlim
                              ! FOR DARK ENERGY
      ISUGRA = 0   ! 1- for RP
      If(ISUGRA.eq.1)Then
         Qlambda  = 3    ! log(lambda/Gev) for RP
         a = -2.119e-2*Qlambda-0.259
         b = -1.833e-2*Qlambda+0.975
         c = -6.975e-3*Qlambda+9.771e-2
         z = 1./AEXPN-1.
         Omz = 1.-(1.-Om0)*AEXPN**(a+b*abs(z)**c)
         write (*,*) ' a,b,c=',a,b,c
         write (*,*) ' power =',a+b*abs(z)**c
         a = -1.45e-2*Qlambda+0.186
         b =  -1.1e-2*Qlambda +0.22
         Ovlim =178.*Omz**(a+b*Omz-1)
         write (*,*) ' OverdensLimit=',Ovlim
         write (*,*) ' Omega_matter(z) =',Omz
         write (*,*) ' Omega_matter(0) =',Om0
      EndIf

      OverdenVir = Ovlim
c      Ovlim   =INPUT('Enter Overdensity Threshold for Halos     => ')

 
      RETURN
      END
C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE LimList(Box,N,Ncp)
C--------------------------------------------------------------
      use Structures


      first =iWeight(1)*1.1   ! add small factor to avoid rounding errors
      xmin = 1.e10
      ymin = 1.e10
      zmin = 1.e10
      xmax = -1.e10
      ymax = -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
            If(iWeight(i).lt.first)Then
               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)
            EndIf
         EndDo          

         Nmx=INT(xmin/Cell)-5
         Nmy=INT(ymin/Cell)-5
         Nmz=INT(zmin/Cell)-5
         Nbx=INT(xmax/Cell)+5
         Nby=INT(ymax/Cell)+5
         Nbz=INT(zmax/Cell)+5

      write(*,*) '       Particles:',N
      write(*,*) '               x:',xmin,xmax
      write(*,*) '               y:',ymin,ymax
      write(*,*) '               z:',zmin,zmax
      write(*,*) ' Label limits  x:',Nmx,Nbx
      write(*,*) '               y:',Nmy,Nby
      write(*,*) '               z:',Nmz,Nbz

      Return
      End
C--------------------------------------------------------------
C                          Make linker lists of particles in each cell
      SUBROUTINE List(Box,N,Ncp)
C--------------------------------------------------------------
      use Structures

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
             Do i=1,Ncp
                Lst(i)=-1
             EndDo
C$OMP PARALLEL DO DEFAULT(SHARED)
C$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

      If(Npp.ne.Np)Return  ! Xpp/Ypp/Zpp arrays are used
                           ! for re-arranging particles 
                           ! to increase locality only when Npp=Np
      ALLOCATE (Xpp(Ncp),Ypp(Ncp),Zpp(Ncp),
     &          VXpp(Ncp),VYpp(Ncp),VZpp(Ncp),Wpp(Ncp), Stat =iStat)
      if(iStat.ne.0)return
                          ! Take only the first N particles
      nn = 0
         Do k=Nmz,Nbz
         Do j=Nmy,Nby
         Do i=Nmx,Nbx
            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
         write (*,*)  ' new number of poits=',nn,' must be=',N
                ! Take only the last  particles
         Do k=Nmz,Nbz
         Do j=Nmy,Nby
         Do i=Nmx,Nbx
            jp =Label(i,j,k)
             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
         write(*,*)' new number of points=',nn,' must be=',Ncp
!                        Restore particles
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
      Do i=1,Ncp         ! (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          
      write(*,*) ' done moving particles back'
      DEALLOCATE (Xpp,Ypp,Zpp,VXpp,VYpp,VZpp,Wpp)
      write(*,*) ' Make new Linker List'
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
             Do i=1,Ncp
                Lst(i)=-1
             EndDo

C$OMP PARALLEL DO DEFAULT(SHARED)
C$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

      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
C--------------------------------------------------------------
C                                 Remove dublicates 
      SUBROUTINE RemVelDublA(Box,Ncentr,dRdubl,fVdubl,Xscale,Vscale)
C--------------------------------------------------------------
      use Structures

      Real      Mhalo,Mhalo2
      double precision q1
      Logical  inside

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
             Do i=1,Ncentr
                Lst(i)=-1
             EndDo
C$OMP PARALLEL DO DEFAULT(SHARED)
C$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,Ncentr
         i=INT(Xm(jp)/Cell)
         j=INT(Ym(jp)/Cell)
         k=INT(Zm(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

      dR2 = dRdubl**2
      fV2 = fVdubl**2
      Vol_2 =Rad(2)**3-Rad(1)**3
      nrem= 0


C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,Mhalo,irad) REDUCTION(+:nrem)
      Do i=1,Ncentr
        Mhalo  =Amc(i)
        irad   =Max(iRadius(i),2)         
        If(Mhalo.le.1.e-7 .or. irad.eq.0)Then
          Lab(i) =-1
          nrem = nrem +1
        Else
          Lab(i) =1
        EndIf
c        if(Lab(i).eq.1)
      EndDo
      write(*,*) ' RemVelDubl:',Ncentr,Xscale,Vscale,nrem
      nrem =0
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,jp,Mhalo,Mhalo2,irad,jrad,Rh_1,Rh_2,Rhh,Vrh2,dd,  
C$OMP+dv,VV,irr,q1,Xc,Yc,Zc,i1,j1,k1,i2,j2,k2,ii,jj,kk,inside)
C$OMP+REDUCTION(+:nrem) SCHEDULE(DYNAMIC,mDynamic)
      Do i=1,Ncentr
c        If(i/1000*1000.eq.i)write (*,*)  ' RmDubl i=',i
        If(Lab(i).ne.-1)Then         ! do not take empty centers
        Mhalo  =Amc(i)
        irad   =Max(iRadius(i),2)
        Rh_1   =Rad(irad)
        Vrh2   =(Vrmc(irad,i)/Vscale)**2

        Xc =Xm(i)
        Yc =Ym(i)
        Zc =Zm(i)
c                      limits for Label
         i1=INT((Xc-dRdubl)/Cell)
         j1=INT((Yc-dRdubl)/Cell)
         k1=INT((Zc-dRdubl)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((Xc+dRdubl)/Cell)
         j2=INT((Yc+dRdubl)/Cell)
         k2=INT((Zc+dRdubl)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
C                                        Look for neibhours in nearest cells
         Do kk =k1,k2
         Do jj =j1,j2
         Do ii =i1,i2
            jp =Label(ii,jj,kk)   !  Dark matter Halo 
 10         If(jp.ne.0)Then

          If(i.ne.jp.and.Lab(jp).ne.-1)Then ! do not compare with empty centers
          dd =(Xm(jp)-Xc)**2 + (Ym(jp)-Yc)**2+ (Zm(jp)-Zc)**2     
          Mhalo2 =Amc(jp)
          inside = (dd.lt.dR2).and.
     &               (Mhalo2.gt.Mhalo/1.5.or.Mhalo2.lt.Mhalo*1.5)
          inside = inside.or.(dd.lt.dR2/10.)
          If(inside)Then    ! close distance pair
             jrad=Max(iRadius(jp),2)
             Rh_2 =Rad(jrad)
             Rhh =MAX(Rh_1,Rh_2,1.e-8)
             dv =(Wxc(jrad,jp)-Wxc(irad,i))**2 +
     .           (Wyc(jrad,jp)-Wyc(irad,i))**2 +      
     .           (Wzc(jrad,jp)-Wzc(irad,i))**2
             VV =MAX(Vrh2,(Vrmc(jrad,jp)/Vscale)**2) ! max Vrms of the two
c             If(dv.lt.(6.-5.*dd/dR2)*fV2*VV)Then   ! close in velocity
c             If(dv.lt.(2.-dd/(dRdubl*Rhh))*fV2*VV)Then ! close in velocity
             If(dv.lt.fV2*VV.or.dd.lt.dR2/10.)Then ! close in velocity
                 irr = min(10,Nrad)             ! costruct weight to find
                q1 = (Mhalo/Mhalo2)*     ! which halo is larger
     &               (Vrmc(irr,i)/Vrmc(irr,jp))**2*
     &              (Nbc(3,i)/real(max(Nbc(3,jp),1)))*
     &              Nbc(irr,i)/real(max(Nbc(irr,jp),1))
                if(abs(q1-1.).lt.1.e-5)q1=jp/float(i)
                nrem = nrem +1

31           format(f10.4,f9.3,f11.2,2i14,3f10.3,2g11.3,i4,2i5,2i3)               
 30             format(2i6,f10.4,f9.3,f11.2,4i6,g12.4,6f10.3,2g11.3,i4)                
 35             format(' MM ',2i8,f10.4,4i6,g11.3,6f10.3,2g11.3,i4)                
                If(q1.lt.1.)Then
                   Lab(i)=0     ! kill dublicate
                EndIf 
             EndIf                 ! end dv     check
          EndIf                    ! end dd<dR2 + Mhalo2>Mhalo/1.5 check
          EndIf                    ! end i=jp     check
          jp = Lst(jp)
          GoTo 10
          EndIf                    ! end jp neq 0
        EndDo                      ! end j loop
        EndDo                      ! end j loop
        EndDo                      ! end j loop
      EndIf                        ! end Lab(i)=-1 check
      EndDo                        ! end i loop

      m  =0
      meanmass =0
      write (*,*)  '  packing Halos in remVelDubl '

      Do i=1,Ncentr
         If(Lab(i).eq.1)Then
            m = m +1
            Xm(m)      =Xm(i)
            Ym(m)      =Ym(i)
            Zm(m)      =Zm(i)
            Amc(m)     =Amc(i)   
            Rmc(m)     =Rmc(i)  
            RmaxV(m)   =RmaxV(i)
            Vrm(m)     =Vrm(i) 
            iRadius(m) =iRadius(i)

            Do ir=1,Nrad  !  Store Results 
               Wxc(ir,m)  =Wxc(ir,i) 
               Wyc(ir,m)  =Wyc(ir,i) 
               Wzc(ir,m)  =Wzc(ir,i) 
               Nbc(ir,m)  =Nbc(ir,i) 
               NbcS(ir,m)  =NbcS(ir,i) 
               NbcM(ir,m)  =NbcM(ir,i) 
               NbcG(ir,m)  =NbcG(ir,i) 
               Rrc(ir,m)  =Rrc(ir,i)
               Vrmc(ir,m) =Vrmc(ir,i)
               Vrada(ir,m) =Vrada(ir,i)
            EndDo
            meanmass =meanmass +Nbc(Nrad,m)
         EndIf
      EndDo
      xmeanmass = FLOAT(meanmass)/FLOAT(m)      
      write(*,15) m,Ncentr,xmeanmass
      write(20,15) m,Ncentr,xmeanmass
      write(*,20) dRdubl*Xscale*hubble,fVdubl
 15   format(5x,'Vel.Dublicates: New number of objects=',i7,' Old=',i7,
     .       ' mean N_particles=',g11.3)
 20   format('     Distance Difference=',f8.4,/
     .       '     Vrelative/Vrms <    ',f8.4)
      Ncentr =m

      Return
      End
C--------------------------------------------------------------
C                                Find all neibhours for a center
C                                Xc,Yc,Zc - center; a0 -its weight
C
      SUBROUTINE Neib(Xnew,Ynew,Znew,Amnew,Xc,Yc,Zc,Radius)
C--------------------------------------------------------------
      use Structures

      Real*8  Xw,Yw,Zw,aMw,ww
C                                       Initiate counters
        Xnew =0.
        Ynew =0.
        Znew =0.
  	Amnew=0.
        Xw =0.
        Yw =0.
        Zw =0.
  	aMw=0.
        Radius2=Radius**2
        weig1 = 1.01*weightSmall
c                      limits for Label
         i1=INT((Xc-Radius)/Cell)
         j1=INT((Yc-Radius)/Cell)
         k1=INT((Zc-Radius)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((Xc+Radius)/Cell)
         j2=INT((Yc+Radius)/Cell)
         k2=INT((Zc+Radius)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
C                                        Look for neibhours
         nn =0
         Do k=k1,k2
         Do j=j1,j2
         Do i=i1,i2
            jp =Label(i,j,k)
 10         If(jp.ne.0)Then
               dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
               If(dd.lt.Radius2.and.iWeight(jp).le.weig1)Then
                  nn =nn +1
                  ww =iWeight(jp)
	             Xw =Xw +Xp(jp)*ww
	             Yw =Yw +Yp(jp)*ww
	             Zw =Zw +Zp(jp)*ww
	             aMw=aMw+ww
               EndIf  
                  jp =Lst(jp)
               GoTo10
            EndIf
         EndDo
         EndDo
         EndDo
         	If(aMw.GT.0.)Then
	            Xnew =Xw /aMw
	            Ynew =Yw /aMw
	            Znew =Zw /aMw
            EndIf
            Amnew = aMw
      Return
      End
C--------------------------------------------------------------
C                                Find all neibhours for a center
C                                Xc,Yc,Zc - center; 
C                               
      SUBROUTINE SNeib(Vx1,Vy1,Vz1,Nob1,Xc,Yc,Zc)
C--------------------------------------------------------------
      use Structures

      DIMENSION  Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Nob1(Nrad)
      Real*8     Wx1(Nrad),Wy1(Nrad),Wz1(Nrad)


      iFlag = 0
         Radmax =Rad(Nrad)/FracSearch 
         d0      = Rad(1)
 50      dmax = Radmax**2
         Ninside = aNvir*Radmax**3


         Do i=1,Nrad
              Nob1(i) =0
              Wx1(i)  =0.
              Wy1(i)  =0.
              Wz1(i)  =0.
         EndDo
c                      limits for Label
         i1=INT((Xc-Radmax)/Cell)
         j1=INT((Yc-Radmax)/Cell)
         k1=INT((Zc-Radmax)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((Xc+Radmax)/Cell)
         j2=INT((Yc+Radmax)/Cell)
         k2=INT((Zc+Radmax)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
            nps =0
C                                        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)Then
               dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
               If(dd.lt.dmax)Then
                  d1 =sqrt(dd)
                    ir =INT(max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1
                    ir = min(max(1,ir),Nrad)
                    wu =iWeight(jp)
                    iww = INT(wu/weightSmall+0.05)
                    Nob1(ir)= Nob1(ir)+ iww
c             write(30,'(20x,i9,i5,i4,6g11.3)')jp,iww,ir,wu,dd,Rad2(ir)
c                   If(iww.eq.iw1)NobS(ir)= NobS(ir)+ 1
c                    If(iww.eq.iw8)NobM(ir)= NobM(ir)+ 1
c                    If(iww.eq.iw64)NobG(ir)= NobG(ir)+ 1
c                    Rob1(ir)= Rob1(ir)+ d1*iww
c                    Vxx(ir) = Vxx(ir) +
c     +                            iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2)
                    Wx1(ir) = Wx1(ir) +iww*VXp(jp)  
                    Wy1(ir) = Wy1(ir) +iww*VYp(jp)
                    Wz1(ir) = Wz1(ir) +iww*VZp(jp)
               EndIf
               jp =Lst(jp)
               GoTo10
            EndIf
         EndDo
         EndDo
         EndDo
         Do ir =2,Nrad
            Nob1(ir)= Nob1(ir)   + Nob1(ir-1)
            Wx1(ir) = Wx1(ir)      + Wx1(ir-1)
            Wy1(ir) = Wy1(ir)      + Wy1(ir-1) 
            Wz1(ir) = Wz1(ir)      + Wz1(ir-1) 
         EndDo
 
         Do ir=1,Nrad
           Vx1(ir)  =0.
           Vy1(ir)  =0.
           Vz1(ir)  =0.
           If(Nob1(ir).GT.0)Then
               Nobj    = Nob1(ir)
               Vx1(ir) = Wx1(ir)/Nobj
               Vy1(ir) = Wy1(ir)/Nobj
               Vz1(ir) = Wz1(ir)/Nobj
           EndIf
         EndDo
         If(Nob1(Nrad).gt.Ninside.and.iFlag.eq.0)Then
             Radmax =Rad(Nrad)
             Iflag = 1
             goto50
          EndIf 
      Return
      End
C--------------------------------------------------------------
C                       Add points around the Box to take into account periodical conditions
C  It replicates points periodically and keeps those which are at  distance
C  less than Radius from the Box. 
C                        N     = number of points inside the Box
C                        Ncp = total number of points
      SUBROUTINE Points(Box,N,Ncp,Radius)
C--------------------------------------------------------------
      Use Structures

      Logical    Inside
         B1  =-Radius
         B2  =Box+Radius
         B3  =Box-Radius
         Ncp =N
      Do i=1,N               ! 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
                             EndIf
                        EndIf
                     EndDo
                  EndIf
               EndDo
               EndIf
            EndDo
         EndIf
      EndDo
      If(Ncp.Gt.Np)Then
         write(*,*)' Too many points:',Ncp,' Maximum is =', Np
         STOP
      EndIf
      Return
      End
C-------------------------------------------------------------
C                     Find positions of 'Ncentr' spheres, which maximize number of particles
C                    inside 'radius Radius'. 
C                    Iterations are repeated until 1) mass stops increasing   OR
C                                                       2) the sphere stops moving: dR is less than Riter 
C                    Looks for all particles inside Radius for periodical set 
      SUBROUTINE Pairs(Box,N,Ncentr,Riter)
C--------------------------------------------------------------
      Use Structures

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i)
      Do i=1,Ncentr
      	Lab(i) =1
         Amc(i) =0
         If(Rmc(i).lt.1.e-8)write (*,*) ' Error in Rmc: zero',Rmc(i),i
         If(Rmc(i).le.0.5*RadSr2)
     &   write (*,'(" Error Rmc too small=",g11.3," min=",g11.3,i8)')
     &              Rmc(i),RadSr2,i
         If(Rmc(i).gt.RadSr1)write (*,*) ' Error in Rmc: big ',Rmc(i),i
c         Rmc(i) =Radius
      EndDo
      Niter =0
      iDynamic = mDynamic
10    Iter  =0
      Niter =Niter +1
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,Xc,Yc,Zc,a0,in,Radius)
C$OMP+PRIVATE (Xnew,Ynew,Znew,Amnew,Dr)
C$OMP+SCHEDULE(DYNAMIC,iDynamic) REDUCTION(+:Iter)
      Do i=1,Ncentr
         Xc =Xm(i)
         Yc =Ym(i)
         Zc =Zm(i)
         a0 =Amc(i)
         in =Lab(i)
         Radius =Rmc(i)
         If(Radius.lt.0.999*RadSR2)in =-Niter
         If(in.gt.0)Then         ! If not converged, keep iterating
            Call Neib(Xnew,Ynew,Znew,Amnew,Xc,Yc,Zc,Radius) ! new center of mass
            Dr =MAX(ABS(Xnew-Xc),ABS(Ynew-Yc),ABS(Znew-Zc))
            IF(Dr.LT.Riter.or.Amnew.lt.Amc(i))
     .                  Lab(i)=-Niter     ! iterations converged
            Iter =Iter +1
	         IF(Xnew.lt.0.)Xnew =Xnew+Box
	         IF(Ynew.lt.0.)Ynew =Ynew+Box
	         IF(Znew.lt.0.)Znew =Znew+Box
	         IF(Xnew.gt.Box)Xnew =Xnew-Box
	         IF(Ynew.gt.Box)Ynew =Ynew-Box
	         IF(Znew.gt.Box)Znew =Znew-Box
            Xm(i)  = Xnew
            Ym(i)  = Ynew
            Zm(i)  = Znew
            Amc(i) = Amnew
c            write (*,'(i6,6g12.5)') i,
c     &           Xm(i)*20./128.,Ym(i)*20./128.,Zm(i)*20./128.,Amc(i) 
         EndIf
      EndDo

      Write (*,*) ' Iteration=',Niter,' Centers left=',Iter,iDynamic
      iDynamic =max(min(INT(Iter/float(Ncentr)*mDynamic),
     &                                     mDynamic),kDynamic)
      If(Iter.GE.1 .and. Niter.LT.100)Goto10

      Return
      End

C-------------------------------- Update Statistics of pairs 
C                   
      SUBROUTINE SPairs(Box,N,Ncentr)
C---------------------------------------------------
      Use Structures

      DIMENSION  Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad),
     .             Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad)

C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,ir)
      Do i=1,Ncentr
            Do ir=1,Nrad
               Wxc(ir,i) = 0.
               Wyc(ir,i) = 0.
               Wzc(ir,i) = 0.
               Nbc(ir,i) = 0
               NbcS(ir,i) = 0
               NbcM(ir,i) = 0
               NbcG(ir,i) = 0
               Rrc(ir,i) = 0.
               Vrmc(ir,i)= 0.
            EndDo
      EndDo
C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,ir,Xc,Yc,Zc,Vx1,Vy1,Vz1,Nob1)
C$OMP+SCHEDULE (DYNAMIC,kDynamic)
       Do i=1,Ncentr
!          if(mod(i,10000).eq.0)write (*,*)  ' Sneib halo=',i
         Xc =Xm(i)
         Yc =Ym(i)
         Zc =Zm(i)
         Call SNeib(Vx1,Vy1,Vz1,Nob1,Xc,Yc,Zc)
         Do ir=1,Nrad  !  Store Results 
             Wxc(ir,i) =Vx1(ir)
             Wyc(ir,i) =Vy1(ir)
             Wzc(ir,i) =Vz1(ir)
             Nbc(ir,i) =Nob1(ir)
         EndDo
      EndDo

      Return
      End

C-------------------------------- Remove small halos
C           
C           
      SUBROUTINE SmallRem(ANlimit,Ncentr)
C--------------------------------------------------------------
      Use Structures

      New =0
      Do i=1,Ncentr
         If(Amc(i).GT.ANlimit*Rmc(i)**3)Then
          New =New +1
          Xm(New) =Xm(i)
          Ym(New) =Ym(i)
          Zm(New) =Zm(i)
          Amc(New)=Amc(i)
          Rmc(New)=Rmc(i)
         EndIf
      EndDo
      Ncentr =New
      write (*,*) ' small centers removed. New Ncentr=',Ncentr
      Return
      End
C-------------------------------- Get Catalog of Halos
C                        Remove close centers, keep only the largest
C               
      SUBROUTINE CATALOG(Box,Ncentr,dRdubl)
C--------------------------------------------------------------
      Use Structures

      Real       Mhalo
      double precision q1
      Logical    inside

      dR2 =(max(Rad(1),dRdubl/dClose))**2
      Do i=1,Ncentr
        Lst(i)=-1
      EndDo
C$OMP PARALLEL DO DEFAULT(SHARED)
C$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,Ncentr
         i=INT(Xm(jp)/Cell)
         j=INT(Ym(jp)/Cell)
         k=INT(Zm(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 (*,*) ' CATALOG:',Ncentr
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,Mhalo,irad)
      Do i=1,Ncentr
        If(Amc(i).le.1.e-7)Then
          Lab(i) =-1
        Else
          Lab(i) =1
        EndIf
      EndDo
      write (*,*) ' CATALOG:',Ncentr
      nrem =0
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP+PRIVATE (i,jp,Mhalo,dd,  
C$OMP+q1,Xc,Yc,Zc,i1,j1,k1,i2,j2,k2,ii,jj,kk,inside)
C$OMP+REDUCTION(+:nrem) SCHEDULE(DYNAMIC,500)
      Do i=1,Ncentr
c        If(i/10000*10000.eq.i)write (*,*)  ' RmClose i=',i
        If(Lab(i).ne.-1)Then         ! do not take empty centers
        Mhalo  =Amc(i)

        Xc =Xm(i)
        Yc =Ym(i)
        Zc =Zm(i)
c                      limits for Label
         i1=INT((Xc-dRdubl)/Cell)
         j1=INT((Yc-dRdubl)/Cell)
         k1=INT((Zc-dRdubl)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((Xc+dRdubl)/Cell)
         j2=INT((Yc+dRdubl)/Cell)
         k2=INT((Zc+dRdubl)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
C                                        Look for neibhours in nearest cells
         Do kk =k1,k2
         Do jj =j1,j2
         Do ii =i1,i2
            jp =Label(ii,jj,kk)   !  Dark matter Halo 
 10         If(jp.ne.0)Then

          If(i.ne.jp.and.Lab(jp).ne.-1)Then ! do not compare with empty centers
          dd =(Xm(jp)-Xc)**2 + (Ym(jp)-Yc)**2+ (Zm(jp)-Zc)**2     
          inside = (dd.lt.dR2)
           If(inside)Then    ! close distance pair
             q1   =Mhalo/Amc(jp)
             if(abs(q1-1.d0).lt.1.d-5.and.i.lt.jp)q1=0.9
             if(abs(q1-1.d0).lt.1.d-5.and.i.gt.jp)q1=1.1
 31           format(f10.4,f9.3,f11.2,2i14,3f10.3,2g11.3,i4,2i5,2i3)               
 30             format(2i9,2f10.4,8f9.3)
 35             format('      MM ',2i8,f10.4,g11.3,6f10.3,2g11.3,i4)                

             If(q1.lt.1.)Then
                Lab(i)=0     ! kill dublicate
                nrem  = nrem +1
             EndIf 
          EndIf                    ! end dd<dR2 
          EndIf                    ! end i=jp     check
          jp = Lst(jp)
          GoTo 10
          EndIf                    ! end jp neq 0
        EndDo                      ! end j loop
        EndDo                      ! end j loop
        EndDo                      ! end j loop
      EndIf                        ! end Lab(i)=-1 check
      EndDo                        ! end i loop

      write (*,*)' CATALOG: done removed halos=',nrem 
C                                 Remove maxima with Label =0
      New =0
      Do i=1,Ncentr
        If(Lab(i).eq.1)Then
          New =New +1
          Xm(New) =Xm(i)
          Ym(New) =Ym(i)
          Zm(New) =Zm(i)
          Amc(New)=Amc(i)
          Rmc(New)=Rmc(i)
      !      write (*,'(i6,3f9.3,g12.3,f8.2)')
      !&     New,Xm(i),Ym(i),Zm(i),Amc(i),Rmc(i)  

         EndIf 
      EndDo
      Ncentr =New
      write (*,*) ' dublicates removed. New Ncentr=',Ncentr

      Return
      End
C---------------------------------------------------------
C                              Select particles
c                              as initial seeds of centers
c                              Need to run it twice
c                      iFlag = 0 => count centers
c                            = 1 => get centers
      SUBROUTINE ReadCent(N,Ncentr,Boxx,iFlag)
C--------------------------------------------------------------
      Use Structures

      Real INPUT,Iweigh_first,Iweigh_last,Iweigh_min
      Logical      Ldiff
      SAVE Nline,jneib
c      Xccc =INPUT(' Center x=')
c         write (*,*) Xccc
c      Yccc =INPUT(' Center y=')
c         write (*,*) Yccc
c      Zccc =INPUT(' Center z=')
c         write (*,*) Zccc
c      Rccc =INPUT(' Center r=')
c         write (*,*) Rccc
         hsmall= hubble            ! Hubble/100.
c       Xccc =Xccc/Boxx*NGRID/hsmall
c      Yccc =Yccc/Boxx*NGRID/hsmall
c      Zccc =Zccc/Boxx*NGRID/hsmall
c      Rccc =Rccc/Boxx*NGRID/hsmall
c      write (*,'(a,4f8.2)')' Center for seeds: coords,Rad(Cells)=',
c     &       Xccc,Yccc,Zccc,Rccc 
      write (*,'(a,f8.2,i5)')' Boxx,Ngrid=',Boxx,NGRID
      Ncentr =0
      Iweigh_first =weightSmall
      Iweigh_last  =iWeight(N)
      Iweigh_min  = Iweigh_first
      If(iDebug.eq.1)
     & write (*,*) ' Weights: first species=',Iweigh_first,
     &  ' last=',Iweigh_last,' Weight minimum center=',Iweigh_min
        n1        = 0  
        nhalf    = 0
        n2        = 0
      If(RadSr1.lt.1.1*RadSr2.or.ABS(RadSr2).lt.1.e-8)Then
                                   ! there is only one radius
         Ldiff =.false.
         Radius = RadSr1
         R1        = RadSr1
         R2        = RadSr1
         Rhalf    = RadSr1
      Else
        Ldiff =.true.
        R1       = 1.1*RadSr2
        R2       = 0.9*RadSr1
        Rhalf   =(RadSr1+RadSr2)/2.
        ialpha =2
        RRs  = (RadSr1-RadSr2)
        If(iDebug.eq.1)Then
	  write(*,*) ' Min Radius=',RadSr2,' Max=',RadSr1
	  write(*,*) ' Power index =',ialpha
        EndIf
      Endif 
      If(iHelp.eq.1)Write(*,130)
 130     format(/
     +10x,'These parameters define how and how many spheres are ',/
     +10x,'used to search for density maxima. Tests can be done ',/
     +10x,'with few spheres (answers: 100 0). For final search ',/
     +10x,'a large fraction of particles should be used. ',/
     +10x,'Combinations (np/10 10) or (np/3 5) are recommended',/
     +10x,'Only particles of the smallest mass are used for ',/
     +10x,'initial centers of shperes')

      If(iFlag.eq.0)Then
      Nline=INPUT(' Enter Number of particles for initial seeds=> ')
      jneib=INPUT(' Enter Number of neibhours for a seed       => ')
      EndIf

      If(Nline.eq.0)Goto 50
         fr    = N/FLOAT(Nline)
         xrand =1.00001
         Nstep =MAX(INT(fr+1.),1)
         fr    =1./fr
         Nseed = 121071
         write (*,*) ' Every ',Nstep,' particle is used',
     &                ' as initial center',' Fraction=',fr
         Nplast =0
         Do i=1,N        ! take centers if weight is less than Iweigh_min
           If( iWeight(i).le.Iweigh_min)Then
              Nplast =Nplast +1
           Endif 
         Enddo 
         write (*,*) ' Total pool for particles centers=',Nplast
C         Do i=1,N,Nstep
          Do i=1,Nplast
             If(fr.lt.0.999)xrand =RANDd(Nseed) 
               IF(xrand.lt.fr)Then 
c                 IF(ABS(Xp(i)-Xccc).lt.Rccc.and. 
c     .                ABS(Yp(i)-Yccc).lt.Rccc.and.
c     .                ABS(Zp(i)-Zccc).lt.Rccc)THEN
                  Ncentr           =Ncentr+1
                  IF(Ldiff)Then           ! different radii
                        Rrrr   =RadSr2+
     &                            RRs*(MAX(RANDd(Nseed),1.e-8))**ialpha
                        If(iFlag.eq.1) Rmc(Ncentr) = Rrrr

                        If(Rrrr.lt.R1)n1 = n1+1
                        If(Rrrr.lt.R2.and.Rrrr.ge.R1)nhalf = nhalf+1
                        If(Rrrr.gt.R2)n2 = n2+1
                  Else
                        If(iFlag.eq.1) Rmc(Ncentr) =Radius
                  Endif 
                        If(iFlag.eq.1)Then 
                           Amc(Ncentr) =1.
                           Xm(Ncentr)   =Xp(i)
                           Ym(Ncentr)   =Yp(i)
                           Zm(Ncentr)   =Zp(i)
                        EndIf
                  j=Ncentr

c                 EndIf      ! endif for small box test
               EndIf
         EndDo
         write (*,*) ' Ncenters (random fraction)=',Ncentr
      If(jneib.eq.0)goto 150               ! No more points
 50   Ncnt_max = 0         ! select a fraction of particles from
                           ! each cell
      Ncnt_lrg = 0

      Do k=Nmz,Nbz
      Do j=Nmy,Nby
      Do i=Nmx,Nbx
         jp =Label(i,j,k)
         If(jp.ne.0)Then
           Ncnt  =0
 5         If(jp.ne.0)Then  ! count number of particles in the cell
              Ncnt = Ncnt+1
              jp=Lst(jp) 
              goto 5                                                                                             
           Endif
             !    step for selecting centers  
           jstep = jneib*min(10,max(INT(Ncnt/1e3),1))
             !
           jp        = Label(i,j,k)
           jcount    = 0
           Ncount    = 0
 10        jcount    = jcount+1
              If(jp.ne.0.and.jcount.lt.jstep)Then
                     jp=Lst(jp) ! find if there are enough neib
                     goto 10
              endif 
               jcount    =0
           If(jp.ne.0)Then 
           If(iWeight(jp).le.Iweigh_min)Then !  take this particle 
              Ncount =Ncount +1
c                IF(ABS(Xp(jp)-Xccc).lt.Rccc.and. 
c     .                 ABS(Yp(jp)-Yccc).lt.Rccc.and.
c     .                 ABS(Zp(jp)-Zccc).lt.Rccc)THEN
                Ncentr=Ncentr+1
                  IF(Ldiff)Then           ! different radii
                        Rrrr   =RadSr2+
     &                            RRs*(MAX(RANDd(Nseed),1.e-8))**ialpha
                        If(iFlag.eq.1) Rmc(Ncentr) = Rrrr

                        If(Rrrr.lt.R1)n1 = n1+1
                        If(Rrrr.lt.R2.and.Rrrr.ge.R1)nhalf = nhalf+1
                        If(Rrrr.gt.R2)n2 = n2+1
                  Else
                        If(iFlag.eq.1) Rmc(Ncentr) =Radius
                  Endif 
                        If(iFlag.eq.1) Then 
                           Amc(Ncentr) =1.
                           Xm(Ncentr)  =Xp(jp)
                           Ym(Ncentr)  =Yp(jp)
                           Zm(Ncentr)  =Zp(jp)
                        EndIf
c                If(abs(Xm(Ncentr)-152.).lt.0.5.and.
c     &           abs(Ym(Ncentr)-112.3).lt.0.5)
c                write (30,30)Ncentr,i,j,k,
c     &              Xm(Ncentr), Ym(Ncentr), Zm(Ncentr)
c 30             format(i6,5x,3i4,8f9.2)

c                  EndIf           ! endif for small-box test
c                 If(Ncount.lt.300) Goto 10
                 Goto 10
           EndIf    ! end (iWeight< min) condition
           EndIf    ! end (jp>0)         condition
           EndIf    ! end (jp>0)         condition
           If(Ncount.gt.0)Then
c             write (*,'(a,i8,a,i5)')' Np_in_cell=',Ncnt,
c     &                            ' select each ',jstep  
c             write (*,'(5x,a,i8,3i4,i8)')'Count=',Ncount,i,j,k,Ncentr
             Ncnt_max = max(Ncnt_max,Ncnt)
             If(Ncnt.gt.1e4)Ncnt_lrg = Ncnt_lrg+1
           EndIf 
      EndDo
      EndDo
      EndDo
150   write(*,*) ' Ncentr (total)=',Ncentr  !,' Max allowed=',Nc
      write(*,*) ' Max number of particles in a cell   =',Ncnt_max
      write(*,*) ' N cells with more than 1e4 particles=',Ncnt_lrg
      !If(Ncentr.gt.Nc)Then
      !   write (*,*) ' Too many centers'
      !   STOP ' Too many centers'
      !Endif 
         If(Ncentr.gt.0)Then
            write (*,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2
            write (2,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2
            write (20,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2

 200   format(' Number of centers with radii between ',3x,2f7.4,
     &                ' was ',i8,
     &              /18x,' with radii between ',3x,2f7.4,' was ',i8,
     &              /18x,' with radii between ',3x,2f7.4,' was ',i8)
          Endif 
      Return
      End
C---------------------------------------------------------
C                                        Read list of maxima
c                                     Xnew =X-1: coord. now are 0-NGRID       
      SUBROUTINE ReadMax(N,Xscale)
C--------------------------------------------------------------
      Use Structures

         read(3) 
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,
     +                  NROWC,NGRIDC,NREC0,Om0,Oml0,hubble,
     +                  Ocurv
         write(*,200)
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,
     +                  NROWC,NGRIDC,NREC0,Om0,Oml0,hubble,
     +                  Ocurv
200   FORMAT(
     +           1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.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)
       N =0
      xxm =1000.
      xxn =-1000.
 10   read (3,err=5,end=5) X0, Y0, Z0, VvX,VvY,VvZ
            N=N+1
            If(N.le.Np)Then
            Xp(N)  =X0-1.
            Yp(N)  =Y0-1.
            Zp(N)  =Z0-1.
            VXp(N)=VvX
            VYp(N)=VvY
            VZp(N)=VvZ
            xxm    =min(xxm,Xp(N))
            xxn     =max(xxn,Xp(N))
            EndIf
            GoTo 10
 5    write(*,*) ' Read Particles: X{Min,Max}=',xxm,xxn,' N=',N
      If(N.GT.Np)Then
         write (*,*) ' Too many particles: Max=',Np
         Stop
      EndIf
      Return
      End
C---------------------------------------------------------
C                                        Read particles
      SUBROUTINE ReadPntP(N,Fract,xshift,yshift,zshift)    
C--------------------------------------------------------------
      Use Structures

      N =0
      Iz=1
      Iy=1
      Ff=1.
c      xshift = 100.
c      yshift = 100.
c      zshift = 0.
      xshift = 0.
      yshift = 0.
      zshift = 0.

      xmin =1.e10
      xmax =-1.e10
      ymin =1.e10
      ymax =-1.e10
      zmin =1.e10
      zmax =-1.e10
      Nsmall =0
      Lsmall =lspecies(1)
      If(Fract.le.0.666)Then
         Iz=2
         Ff=0.5
         If(Fract.le.0.333)Then
            Iy=2
            Ff=0.25
         Endif
      EndIf
      Fract =Ff
      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
c      write (*,*) ' Pages=',Npages,' Species=',Nspecies
c      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
C				       Loop over particles
         READ  (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ
	   DO   IN=1,In_page
                Xp(iL+IN)  =XPAR(IN)
                Yp(iL+IN)  =YPAR(IN)
                Zp(iL+IN)  =ZPAR(IN)
                VXp(iL+IN)  =VX(IN)
                VYp(iL+IN)  =VY(IN)
                VZp(iL+IN)  =VZ(IN)
           EndDo
      EndDo
      N      = lspecies(Nspecies)
      write(*,*) ' Done Reading Points. N=',N
 
!$OMP PARALLEL DO DEFAULT(SHARED) 
!$OMP+ PRIVATE (i)                 
!$OMP+ REDUCTION(MAX:xmax,ymax,zmax) 
!$OMP+ REDUCTION(MIN:xmin,ymin,zmin) 
      Do i   =1,N 
            Xp(i)  =Xp(i)-1.
            Yp(i)  =Yp(i)-1.
            Zp(i)  =Zp(i)-1.
            If(Xp(i).gt.NGRID)Xp(i)  =Xp(i)-NGRID
            If(Yp(i).gt.NGRID)Yp(i)  =Yp(i)-NGRID
            If(Zp(i).gt.NGRID)Zp(i)  =Zp(i)-NGRID
            If(Xp(i).le.0.)Xp(i)  =Xp(i)+NGRID
            If(Yp(i).le.0.)Yp(i)  =Yp(i)+NGRID
            If(Zp(i).le.0.)Zp(i)  =Zp(i)+NGRID
            If(Xp(i).ge.NGRID)Xp(i)  =Xp(i)-1.e-4
            If(Yp(i).ge.NGRID)Yp(i)  =Yp(i)-1.e-4
            If(Zp(i).ge.NGRID)Zp(i)  =Zp(i)-1.e-4
       if(Xp(i).le.0.)write (*,'(" X:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
       if(Yp(i).le.0.)write (*,'(" Y:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
       if(Zp(i).le.0.)write (*,'(" Z:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
      if(Xp(i).ge.NGRID)write(*,'(" X:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
      if(Yp(i).ge.NGRID)write(*,'(" Y:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
      if(Zp(i).ge.NGRID)write(*,'(" Z:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i
            If(i.le.Lsmall)Then
               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)
            EndIf  
      ENDDO 


      write(*,*) ' Small Particles:',Lsmall
      write(*,*) '               x:',xmin,xmax
      write(*,*) '               y:',ymin,ymax
      write(*,*) '               z:',zmin,zmax
      write (*,*)  ' Nmaxpart=',Nmaxpart,Np
      do i=1,1000,200
      write (*,'(i8," Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i)
      EndDo 

      Return
      End

C---------------------------------------------------------
C                                        Read particles
      SUBROUTINE ReadPnt(N,Fract,xshift,yshift,zshift)    
C--------------------------------------------------------------
      Use Structures

      N =0
      Iz=1
      Iy=1
      Ff=1.
c      xshift = 100.
c      yshift = 100.
c      zshift = 0.
      xshift = 0.
      yshift = 0.
      zshift = 0.

      xmin =1.e10
      xmax =-1.e10
      ymin =1.e10
      ymax =-1.e10
      zmin =1.e10
      zmax =-1.e10
      Nsmall =0
      Lsmall =lspecies(1)
      If(Fract.le.0.666)Then
         Iz=2
         Ff=0.5
         If(Fract.le.0.333)Then
            Iy=2
            Ff=0.25
         Endif
      EndIf
      Fract =Ff
      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
c      write (*,*) ' Pages=',Npages,' Species=',Nspecies
c      write (*,*) ' N_in_last=',N_in_last
      DO  IROW=1, Npages         ! Loop over particle pages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            iL = NPAGE*(IROW-1)
c            write (*,*)  ' page=',IROW,' Np=',N
C				       Loop over particles
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,In_page
            N   = N +1
            Xp(N)  =XPAR(IN)-1.+xshift
            Yp(N)  =YPAR(IN)-1.+yshift
            Zp(N)  =ZPAR(IN)-1.+zshift
            If(Xp(N).gt.NGRID)Xp(N)  =Xp(N)-NGRID
            If(Yp(N).gt.NGRID)Yp(N)  =Yp(N)-NGRID
            If(Zp(N).gt.NGRID)Zp(N)  =Zp(N)-NGRID
            If(Xp(N).le.0.)Xp(N)  =Xp(N)+NGRID
            If(Yp(N).le.0.)Yp(N)  =Yp(N)+NGRID
            If(Zp(N).le.0.)Zp(N)  =Zp(N)+NGRID
            If(Xp(N).ge.NGRID)Xp(N)  =Xp(N)-1.e-4
            If(Yp(N).ge.NGRID)Yp(N)  =Yp(N)-1.e-4
            If(Zp(N).ge.NGRID)Zp(N)  =Zp(N)-1.e-4
       if(Xp(N).le.0.)write (*,'(" X:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
       if(Yp(N).le.0.)write (*,'(" Y:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
       if(Zp(N).le.0.)write (*,'(" Z:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
      if(Xp(N).ge.NGRID)write(*,'(" X:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
      if(Yp(N).ge.NGRID)write(*,'(" Y:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
      if(Zp(N).ge.NGRID)write(*,'(" Z:",3g11.4,i10)')Xp(N),Yp(N),Zp(N),N
            If(N.le.Lsmall)Then
               Nsmall =Nsmall + 1
               xmin =min(Xp(N),xmin)
               xmax =max(Xp(N),xmax)
               ymin =min(Yp(N),ymin)
               ymax =max(Yp(N),ymax)
               zmin =min(Zp(N),zmin)
               zmax =max(Zp(N),zmax)
            EndIf
            VXp(N)=VX(IN)
            VYp(N)=VY(IN)
            VZp(N)=VZ(IN)

                Ipart =min(IN+iL,Nmaxpart)                     ! current particle number
             iWeight(N) =iWeight(Ipart)   ! particles weight 
         ENDDO
      ENDDO 

      write(*,*) ' Small Particles:',Nsmall
      write(*,*) '               x:',xmin,xmax
      write(*,*) '               y:',ymin,ymax
      write(*,*) '               z:',zmin,zmax
      write (*,*)  ' Nmaxpart=',Nmaxpart,Np
      do i=1,1000,200
      write (*,'(i8," Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i)
      EndDo 

      Return
      End

C---------------------------------------------------
C                                Remove Unbound particles, find radius and mass
C                                No removal of unbound particles If:
C                                       Toohot =< 0 or Toohot =>5
C                             Vscale =100.*hsmall*Xscale/AEXPN 
C                              Xscale=Box/Ngrid 
C                              Dscale=2.75e+11*hsmall**2*(Box/NROW)**3 
      SUBROUTINE CLEAN(Ncentr,Xscale,Vscale,Dscale,
     +     Fract,Radsc,MnMass,Toohot,Rminhalo)  
C--------------------------------------------------------------
      Use Structures

      dimension Mass(Nrad),MnMass(Nrad),Overdens(Nrad),Radsc(Nrad),
     &          Vrad(Nrad)
      Real          Mass,MnMass,Mhalo

C.... Open_MP
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,ir,iter,Xc,Yc,Zc,Summass,Mass,Overdens)
C$OMP+PRIVATE (irad,Mhalo,Rhalo,Vmaxrot,Rmaxrot,Vcirc,Vcc)
C$OMP+PRIVATE (NobG,Rob1,Vxx,ifirst,iflag,Vrad )
C$OMP+SCHEDULE (DYNAMIC,kDynamic)
      Do i=1,Ncentr          
c         if(mod(i,10000).eq.0)write(*,*)' Clean halo=',i,Ncentr
         Xc =Xm(i)  *Xscale !          Scale all to physical units: comoving
         Yc =Ym(i)  *Xscale
         Zc =Zm(i)  *Xscale
         iflag =0
         ifirst=0
         Do ir =1,Nrad
            Summass     = Om0*Nbc(ir,i)/Fract
            If(Nbc(ir,i).gt.nFirst.and.iflag.eq.0)Then
               ifirst = ir            ! find first significant bin
               iflag  = 1
            EndIf 
            Mass(ir)    = Dscale*Summass
            Overdens(ir)= Mass(ir)/MnMass(ir)
            Vrad(ir) =0.
         EndDo  

C                                     find size, mass, rad ... for a halo
         Call Decide(Mass,Radsc,Overdens,Nrad,irad,Mhalo,Rhalo,
     +               Vrad,Vmaxrot,Rmaxrot,ifirst,dRmaxrot )
         Vcirc=sqrt(Vmaxrot/AEXPN)*6.58e-5 ! =10^-5*sqrt(G 2e33/3.08e24) 

         If(Toohot.gt.0. .and. Toohot.lt.5.)Then
         Do iter =0,5      !  Remove High Energy particles from halos
c         If(irad .lt.1 .or. Rhalo.le.Rminhalo)Then
c           iic = iic +1
c            write (30,'(i7,i4,2g11.3,3f9.3,3x,10i6)')
c     &         i,irad,Rhalo,Mhalo,Xc,Yc,Zc,(Nbc(ir,i),ir=1,10) 
c         else 
c            write (32,'(i7,2i4,2g11.3,3f9.3,3x,10i6)')
c     &         i,irad,iter,Rhalo,Mhalo,Xc,Yc,Zc,(Nbc(ir,i),ir=1,10) 
c         EndIf 
         If(irad .ge.1 .and. Rhalo.gt.Rminhalo)Then
           Vcc =Vcirc/Vscale          ! back to dimensionless units
           Vmaxrot=Vcc*Toohot*1.03**(4-iter) ! Increase the Rotatioanl Vel.
c           Vmaxrot=Vcc*Toohot ! Increase the Rotatioanl Vel.

           Call RemEng(i,irad,Vmaxrot,Rmaxrot/Xscale)  ! Remove Unbound Particles
           iflag =0
           ifirst=0
           Do ir =1,Nrad
            If(Nbc(ir,i).gt.nFirst.and.iflag.eq.0)Then
               ifirst = ir            ! find first significant bin
               iflag  = 1
            EndIf 
            If(Nbc(ir,i).gt.mFirst.and.ir.gt.2)then
               Vrad(ir)    = abs(Vrada(ir,i))*Vscale/
     &           (sqrt(Mass(ir)/Rad(ir)/Xscale/AEXPN)*6.58e-5)
            else
               Vrad(ir) =0.
            EndIf
              Vrmc(ir,i)=Vrmc(ir,i) *Vscale
             Summass     =Om0*Nbc(ir,i)/Fract
             Mass(ir)    = Dscale*Summass
             Overdens(ir)=Mass(ir)/MnMass(ir)
           EndDo  

          Call Decide(Mass,Radsc,Overdens,Nrad,irad,Mhalo,Rhalo,
     +                Vrad, Vmaxrot,Rmaxrot,ifirst,dRmaxrot)
          Vcirc=sqrt(Vmaxrot/AEXPN)*6.58e-5 !=10^-5*sqrt(2e33/3.08e24)
c        If(abs(Xm(i)-144.25).lt.0.20.and.abs(Ym(i)-125.33).lt.0.20
c     &          .and.abs(Zm(i)-125.145).lt.0.20)Then
c         write(30,'(2i8,3f9.4,2g11.4)')ifirst,irad
c         write(30,'(16f11.3)')(Overdens(ir),ir=1,16)
c         write(30,'(16i11)')(Nbc(ir,i),ir=1,16)
c         write(30,'(16g11.3)')(Vrmc(ir,i),ir=1,16)
c         write(30,'(16g11.3)')(Vrad(ir),ir=1,16)
c         write(30,'(16g11.3)')(Vrada(ir,i)*Vscale,ir=1,16)
c         write(30,'(16g11.3)')
c     &      (sqrt(Mass(ir)/Rad(ir)/Xscale/AEXPN)*6.58e-5,ir=1,16)
c        EndIf 

         EndIf
         EndDo       ! end iter energy
         EndIf
         Amc(i)     =Mhalo ! store results
         Rmc(i)     = Rhalo 
         RmaxV(i)   = Rmaxrot
         Vrm(i)     =Vcirc
         iRadius(i) =irad

       Enddo   ! end i (Ncentr)
c       write(30,*) ' Clean:',Ncentr,Xscale
c      Do i=1,Ncentr          
c         Xc =Xm(i)  *Xscale*0.7 !          Scale all to physical units: comoving
c         Yc =Ym(i)  *Xscale*0.7
c         Zc =Zm(i)  *Xscale*0.7
c                   write (30,'(2i6,3f8.3,5g11.3)')i,iRadius(i),
c     &     Xm(i)*xscale,Ym(i)*xscale,Zm(i)*xscale,
c     &     Amc(i), Rmc(i), Vrm(i)
c      EndDo 

       Return
       End
C------------------------------------------------------------- 
C                                       Find  halo radius and mass
      SUBROUTINE Decide(Mass,Rad,Overdens,Nrad,irad,Mhalo,Rhalo,
     +                  Vrad,Vmaxrot,Rmaxrot,ifirst,dRmaxrotA )
C--------------------------------------------------------------
      COMMON /OVER/ Ovlim,Ovcnt
      Real   Rad(Nrad),Overdens(Nrad),Mass(Nrad),Mhalo,
     &       Vrad(Nrad)
      DATA  slope/0.999/
         Vmaxrot =0
         Rmaxrot =0. 
         Mhalo=0.
         Rhalo=0.   
      If(ifirst.eq.0)Then
         irad =0         ! too few particles
         RETURN
      EndIf

      If(MAX(Overdens(ifirst),Overdens(2)).LT.Ovcnt)Then
         irad =0         ! max overdensity is less then Ovcnt
         RETURN
      EndIf
      iov  =0
      Do ir=Nrad,ifirst,-1
         If(Overdens(ir).gt.Ovlim .and. iov.eq.0)Then
            iov  =ir
         EndIf
      EndDo
      If(iov.eq.0)Then
         irad =0         ! max overdensity is less then Ovlim
         RETURN
      EndIf
      iovv =0
      Do ir=ifirst,iov
         irleft =max(ir-1,1)
         irr1   =min(ir+1,Nrad)
         irr2   =min(ir+2,Nrad)
         Overright =(2.*Overdens(irr1)+Overdens(irr2))/3.
         Overleft  =(2.*Overdens(ir)+Overdens(irleft))/3.
         If(Overright.gt.0.97*Overleft.or.Vrad(ir).gt.0.5)Then !stop:rising profile
            If(ir.eq.1)Then ! softer conditions for the first bin
               If(Overright.gt.1.2*Overleft)Then
                   iovv =ir
                   goto 20
               EndIf
            else 
                iovv =ir
               goto 20
            EndIf
         EndIf
      EndDo
 20   If(iovv.ne.0)iov =iovv
         irad =iov
            If(irad.lt.ifirst)write (*,*)' err: irad=',irad,ifirst


 30      irleft  =max(irad-2,1)  ! check if overdens is steep
         irright =min(irad+1,Nrad)
         gradlf  =Overdens(irleft) * (Rad(irleft))**slope
         gradrt  =Overdens(irright)* (Rad(irright))**slope
         If(gradlf.lt.gradrt)Then  ! not a steep gradient
            irad =irad -1        ! reduce radius
            iovv =irad
            If(irad.lt.ifirst)Then   ! zero radius -> quit
               irad =0
               RETURN
            Else
               goto 30             ! iterate
            EndIf
         EndIf
c            If(irad.lt.ifirst)write (*,*)' err2: irad=',irad,ifirst
         If(irad.lt.Nrad)Then    ! interpolate to Overlim
            If(Ovlim.gt.Overdens(irad+1) .and.
     .         Overdens(irad).gt.Ovlim)Then       
              Mhalo=(Mass(irad)*(Ovlim-Overdens(irad+1))+
     .             Mass(irad+1)*(Overdens(irad)-Ovlim)) /
     .             (Overdens(irad)-Overdens(irad+1))
              Rhalo=(Rad(irad)*(Ovlim-Overdens(irad+1))+
     .             Rad(irad+1)*(Overdens(irad)-Ovlim)) /
     .             (Overdens(irad)-Overdens(irad+1))
            Else                 ! Overdens is > Overlim 
              Mhalo =Mass(irad)
              Rhalo =Rad(irad)
            EndIf
         Else                 ! Overdens is > Overlim 
            Mhalo =Mass(irad)
            Rhalo =Rad(irad)
         EndIf

         imax  = 0
         Do ir =1,irad  ! find max rotational velocity
            V2    = Mass(ir)/Rad(ir)
            If (V2.gt.Vmaxrot)Then
               Vmaxrot =V2               
               Rmaxrot =Rad(ir)
               imax = ir
            EndIf
         EndDo
                                ! improve Rmaxrot value by finding
                              ! radius containing F(1)/F(2.15)=0.415 fraction
                              ! of mass at Vmax. Then Rmax=2.15*R
         imax = max(imax,1)
         aMrs = 0.415*Mass(imax)
c         write (*,*)  ' == imax=',imax,' mass=',aMrs
c         write (*,'(1P,10G11.3)') Mass
         Do ir=1,irad
            If(Mass(ir).gt.aMrs)Then
               if(ir.eq.1)Then
                  aM0 =0.
                  r0     =0.
               Else
                  aM0 =Mass(ir-1)
                  r0     =Rad(ir-1)
               EndIf
               Rmaxrot2 =2.15*(r0+(Rad(ir)-r0)*(aMrs-aM0)/
     &                          max(Mass(ir)-aM0,1.e-3))
               Rmaxrot =min(Rmaxrot,Rmaxrot2,Rhalo/dRmaxrotA)

               return
            EndIf 
          EndDo 
 50       write (*,'(4i4,1P,10G11.3)') ir,irad,ifirst,iovv,
     &             Mass(ir),Ovlim,   Rmaxrot,
     &                Rmaxrot2,Rad(ir),r0
          write (*,'(1P,10g11.3)')Overdens  

      Return
      End
C------------------------------------------------
C                                  Remove unbound particles for i-th halo
C                                  irad   = radius of the halo
C                                 Vmaxrot = max rotational Velocity
C                                 Rmaxrot = radius at max rot.velocity
      SUBROUTINE RemEng(i,irad,Vmaxrot,Rmaxrot)
C-----------------------------------------------
      Use Structures

      DIMENSION VXxc(Nrad), VYyc(Nrad), VZzc(Nrad),
     &             Vescape2(Nrad)            
      DIMENSION  Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad),
     .      Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad),Varad(Nrad)

         Xc  =Xm(i)  ! store center of mass and velocity      
         Yc  =Ym(i)     
         Zc  =Zm(i)     
            irr = 1  !Max(irad,2)         ! don't take V from the first bin
 10         if(Nbc(irr,i).lt.nmin.and.irr.lt.iCentr)Then
               irr =irr+1
               goto 10
            EndIf 
            VXc =Wxc(irr,i)
            VYc =Wyc(irr,i)
            VZc =Wzc(irr,i)

         Do ir =1,Nrad ! store mean velocity
            VXxc(ir) =Wxc(ir,i) 
            VYyc(ir) =Wyc(ir,i) 
            VZzc(ir) =Wzc(ir,i) 
                                         ! escape velocity for NFW profile
            Vescape2(ir) = (2.15*Vmaxrot)**2 * 
     &           max(log(1.0+2.*Rad(ir)/Rmaxrot)/(Rad(ir)/Rmaxrot),
     &                   0.25) 
               Wxc(ir,i) = 0.
               Wyc(ir,i) = 0.
               Wzc(ir,i) = 0.
               Nbc(ir,i) = 0
               NbcS(ir,i) = 0
               NbcM(ir,i) = 0
               NbcG(ir,i) = 0
               Rrc(ir,i) = 0.
               Vrmc(ir,i)= 0.
               Vrada(ir,i)= 0.
         EndDo
            Call SrejNeib(VXxc,VYyc,VZzc,Xc,Yc,Zc,Vescape2,
     &           Vx1,Vy1,Vz1,Nob1,NobS,NobM,NobG,Rob1,Vxx
     &           ,VXc,VYc,VZc,Varad)
            Do ir=1,Nrad  !  Store Results 
               Wxc(ir,i) =Vx1(ir)
               Wyc(ir,i) =Vy1(ir)
               Wzc(ir,i) =Vz1(ir)
               Nbc(ir,i) =Nob1(ir)
               NbcS(ir,i) =NobS(ir)
               NbcM(ir,i) =NobM(ir)
               NbcG(ir,i) =NobG(ir)
               Rrc(ir,i)  =Rob1(ir)
               Vrmc(ir,i) =Vxx(ir)
               Vrada(ir,i)=Varad(ir)
            EndDo
      Return
      End
C--------------------------------------------------------------
C                                Accumilate statistics only
C                                if relative velocity is not large
C                                Xc,Yc,Zc - center; 
C                               
      SUBROUTINE SrejNeib(VXxc,VYyc,VZzc,Xc,Yc,Zc,Vescape2,
     &          Vx1,Vy1,Vz1,Nob1,NobS,NobM,NobG,Rob1,Vxx
     &           ,VXc,VYc,VZc,Varad)
C--------------------------------------------------------------
      Use Structures

      DIMENSION VXxc(Nrad), VYyc(Nrad), VZzc(Nrad),Vescape2(Nrad)
      DIMENSION  Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad),
     .       Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad),Varad(Nrad)

      iFlag = 0
         Radmax =Rad(Nrad)/FracSearch
         d0      = Rad(1)

 50      dmax = Radmax**2
         Ninside = aNvir*Radmax**3
         iw1 = 1
         iw8 = 8
         iw64 = 64

         Do i=1,Nrad
              Nob1(i) =0
              NobS(i) =0
              NobM(i) =0
              NobG(i) =0
              Rob1(i) =0.
              Vx1(i)  =0.
              Vy1(i)  =0.
              Vz1(i)  =0.
              Vxx(i)  =0.
              Varad(i) =0.
         EndDo
c                                           limits for Label
         i1=INT((Xc-Radmax)/Cell)
         j1=INT((Yc-Radmax)/Cell)
         k1=INT((Zc-Radmax)/Cell)
            i1=MIN(MAX(Nmx,i1),Nbx)
            j1=MIN(MAX(Nmy,j1),Nby)
            k1=MIN(MAX(Nmz,k1),Nbz)
         i2=INT((Xc+Radmax)/Cell)
         j2=INT((Yc+Radmax)/Cell)
         k2=INT((Zc+Radmax)/Cell)
            i2=MIN(MAX(Nmx,i2),Nbx)
            j2=MIN(MAX(Nmy,j2),Nby)
            k2=MIN(MAX(Nmz,k2),Nbz)
C                                        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)Then
                dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
             If(dd.lt.dmax)Then
                  d1 =sqrt(dd)
                    iov =INT(   max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1
                    iov = min(max(1,iov),Nrad)
c                    dv =  (VXxc(iov)-VXp(jp))**2+
c     +                     (VYyc(iov)-VYp(jp))**2+
c     +                     (VZzc(iov)-VZp(jp))**2
                    dv =  (VXc-VXp(jp))**2+
     +                    (VYc-VYp(jp))**2+
     +                    (VZc-VZp(jp))**2
                If(dv .lt. Vescape2(iov))Then
                    ww =iWeight(jp)
                    iww =INT(ww/weightSmall+0.005)
                    Nob1(iov)= Nob1(iov)+ iww
                    If(iww.eq.iw1)NobS(iov)= NobS(iov)+ 1
                    If(iww.eq.iw8)NobM(iov)= NobM(iov)+ 1
                    If(iww.eq.iw64)NobG(iov)= NobG(iov)+ 1
                    Rob1(iov)= Rob1(iov)+ d1*iww
                    Vxx(iov) = Vxx(iov) +
     +                          iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2)
                    Vx1(iov) = Vx1(iov) +iww*VXp(jp)  
                    Vy1(iov) = Vy1(iov) +iww*VYp(jp)
                    Vz1(iov) = Vz1(iov) +iww*VZp(jp)
                    Varad(iov) =Varad(iov) + iww*
     +                   ( (VXp(jp)-VXc)*(Xp(jp)-Xc)+
     +                     (VYp(jp)-VYc)*(Yp(jp)-Yc)+
     +                     (VZp(jp)-VZc)*(Zp(jp)-Zc) )
     &                   /max(d1,1.e-8)
c               write (*,'(i7,i3,2i5,g11.4,2x,3g11.3,2x,6g11.4)') 
c     &            jp,iov,Nob1(iov),iww,Varad(ir),
c     &             Xp(jp),Yp(jp),Zp(jp),VXp(jp),VYp(jp),VZp(jp),
c     &              VXc,VYc,VZc
 
                EndIf
             EndIf
           jp =Lst(jp)
           GoTo10
           EndIf
         EndDo
         EndDo
         EndDo

         Do ir =1,Nrad
     	      Varad(ir) = Varad(ir)/max(Nob1(ir),1)
         EndDo
         Do ir =2,Nrad-1
            If(Nob1(ir).eq.0)Varad(ir)=(Varad(ir-1)+Varad(ir+1))/2.
         EndDo

         Do ir =2,Nrad  ! make intergal quantities
            Nob1(ir)= Nob1(ir)   + Nob1(ir-1)
            NobS(ir)= NobS(ir)   + NobS(ir-1)
            NobM(ir)= NobM(ir)   + NobM(ir-1)
            NobG(ir)= NobG(ir)   + NobG(ir-1)
            Rob1(ir)= Rob1(ir)   + Rob1(ir-1)
            Vxx(ir) = Vxx(ir)    + Vxx(ir-1)
            Vx1(ir) = Vx1(ir)    + Vx1(ir-1)
            Vy1(ir) = Vy1(ir)    + Vy1(ir-1) 
            Vz1(ir) = Vz1(ir)    + Vz1(ir-1) 
         EndDo

        Do ir=1,Nrad
           Vmn =0.
     	      If(Nob1(ir).GT.0)Then
               Nobj =Nob1(ir)
               Vx1(ir) = Vx1(ir)/Nobj
               Vy1(ir) = Vy1(ir)/Nobj
               Vz1(ir) = Vz1(ir)/Nobj
              Vmn     = Vx1(ir)**2+Vy1(ir)**2+Vz1(ir)**2
               Vxx(ir) = SQRT(ABS(Vxx(ir)/Nobj -Vmn))
	         EndIf
         EndDo
         If(Nob1(Nrad).gt.Ninside.and.iFlag.eq.0)Then
             Radmax =Rad(Nrad)
             Iflag = 1
             goto50
          EndIf 

      Return
      End
C--------------------------------------------------
C                             Set Weights of particles for
C                             fast access
      SUBROUTINE SetWeights
C--------------------------------------------------
      Use Structures

      Wpr =NGRID**3/FLOAT(NROW)**3
      If(Nspecies.eq.0)Then  ! old  constant weights
         N_particles =lspecies(1)
         lspecies(1) = NROW**3
         N_particles =lspecies(1)
         Nspecies = 1
         wspecies(1) = 1.
         Ww  = 1.
         Do i=1,N_particles
            iWeight(i) =Ww
         EndDo
      Else
         N_particles =min(lspecies(Nspecies),Nmaxpart)
         jstart =1
         If(iDebug.eq.1)
     +    write (*,*) ' N_Particles=',N_particles
         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)
            Do k=jstart ,jend
               iWeight(k) =wspecies(j)/Wpr
            EndDo
            jstart =jend
         EndDo
      EndIf
      If(iDebug.eq.1)Then
      write (*,*) ' Set Weights for ',N_particles,' particles:',
     &       ' w_first=' ,iWeight(1),' w_last=',  iWeight(N_particles)
      write (*,*) ' Nspecies=',Nspecies,' Weights=',wspecies(1),
     &            wspecies(2), wspecies(3)
      EndIf
      weightSmall = iWeight(1)
      RETURN
      END
C-----------------------------------------------------------
C                        Open file with answers for rerunning
C                        Print help and description
      SUBROUTINE Help
      Use Structures
      

      Open(25,file='InHalos.dat',Status='UNKNOWN')
      Write (*,100)
 100  format(//,55('-'),/
     +1x,'You run BDM - Bound Density Maxima - code'/
     +5x,'Code reads data from two files: PMcrd.DAT',
     +                 ' and PMcrs0.DAT',/
     +5x,'It writes output to two files  :', 
     +   ' Catshort.DAT and Catalog.DAT',/
     +5x,'Code will ask you few questions to set parameters',
     +   ' of the search',/
     +5x,'Your answers will be written into file InHalos.dat')
 




      RETURN
      END


c     ------------------------
      real function seconds ()
c     ------------------------
c
c     purpose: returns elapsed time in seconds
c     uses   : xlf90 utility function timef()
c              (in this case - to be compiled with xlf90 or with -lxlf90)
c              or dtime (Exemplar) or etime (Power Challenge)
c      real*8 timef
      real*8 dummy
      real tarray(2)

      dummy = 0.0

c.... for IBM SP
c      dummy   = timef ()
c      seconds = dummy / 1000

c.... for exemplar
c      dummy = dtime(tarray) 
c      seconds = dtime(tarray)

c.... for hitachi
c      dummy = dtime(tarray)
c       seconds = second()

c.... for power challenge
c      dummy = etime(tarray)
      seconds = etime(tarray)

      return
      end

