!------------------------------------------------- ! Get profiles of halos. ! Reads a halo catalog and particles from a snapshot. ! Accumulates statistics for each halo. No removal of unbound particles ! MODULE setHaloParticles PUBLIC Integer, PARAMETER :: NROW =512, & NGRID =512, & nbyteword = 4, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6, & Nrad = 100, & ! # of radial bins Nbf = 10, & ! width of buffer for linker list Nmaxpart = 512**3+50e6, & ! max # of particles Nlinker = 200, & ! linker-list dimension Kdynamic= 100 ! OMP chunk Real, PARAMETER :: FracSearch = 5. Real :: INPUT,Mhalo,Ovlim,Ovcnt,Vclimit Real, DIMENSION(Nrad) :: Rad=0. Real, DIMENSION(Nmaxpart) :: Xp,Yp,Zp,VXp,VYp,Vzp,iWeight Real :: Box=0.,Xscale=0.,Vscale=0.,Dscale=0.,Cell=0., weightSmall Integer :: Np=0,Ncp=0,Ncenter=0, & Nmx,Nbx,Nmy,Nby,Nmz,Nbz COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, & Ocurv,extras(100) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER Real :: RECDAT(NRECL),wspecies(10) Integer :: lspecies(10) Integer, ALLOCATABLE, DIMENSION (:,:,:) :: Label Integer, ALLOCATABLE, DIMENSION (:) :: Lst,Nparth Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc,Rmaxh,Conc Real, ALLOCATABLE, DIMENSION (:,:) :: aMr,Vvrms,Vrad,RadPr,Vrad2 Integer, ALLOCATABLE, DIMENSION (:,:) :: nPr EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/) Contains !-------------------------------------------------------------- ! Make linker lists of particles in each cell SUBROUTINE List !-------------------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncp Lst(i)=-1 EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ncp i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) k=INT(Zp(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) k=MIN(MAX(Nmz,k),Nbz) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo Return End SUBROUTINE List !-------------------------------------------------------------- ! Find all neibhours for a center ! Xc,Yc,Zc - center; ! SUBROUTINE SNeib(iHalo) !-------------------------------------------------------------- x = Xc(iHalo) ; y = Yc(iHalo); z = Zc(iHalo) Vx1=VXc(iHalo) ; Vy1= VYc(iHalo); Vz1= VZc(iHalo) Radmax = Rad(Nrad) dmax = Radmax**2 d0 = Rad(1) ! write(*,'(3i6,3g12.3)')iHalo,nn,mm,Radmax,Cell ! write(*,'(10x,2(3x,3f8.2))') x,y,z,Vx1,Vy1,Vz1 nn = 0 ! counter for pre-selected particles ! limits for Label i1=INT((x-Radmax)/Cell) j1=INT((y-Radmax)/Cell) k1=INT((z-Radmax)/Cell) i1=MIN(MAX(Nmx,i1),Nbx) j1=MIN(MAX(Nmy,j1),Nby) k1=MIN(MAX(Nmz,k1),Nbz) i2=INT((x+Radmax)/Cell) j2=INT((y+Radmax)/Cell) k2=INT((z+Radmax)/Cell) i2=MIN(MAX(Nmx,i2),Nbx) j2=MIN(MAX(Nmy,j2),Nby) k2=MIN(MAX(Nmz,k2),Nbz) ! Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0)Then dd =(x-Xp(jp))**2+(y-Yp(jp))**2+(z-Zp(jp))**2 ! write(*,'(10x,1p7g12.3)')Xp(jp),Yp(jp),Zp(jp),sqrt(dd),sqrt(dmax) If(dd.lt.dmax)Then d1 =sqrt(dd) ir =INT(max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1 ! bin for this particle ir = min(max(1,ir),Nrad) vv =(Vx1-VXp(jp))**2+(Vy1-VYp(jp))**2+(Vz1-VZp(jp))**2 vrr=((Vx1-VXp(jp))*(x-Xp(jp))+ & (Vy1-VYp(jp))*(y-Yp(jp))+ & (Vz1-VZp(jp))*(z-Zp(jp))) /max(d1,1.e-8) ww =iWeight(jp) RadPr(ir,iHalo) =RadPr(ir,iHalo) +d1*ww nPr(ir,iHalo) =nPr(ir,iHalo) +1 aMr(ir,iHalo) =aMr(ir,iHalo) +ww Vrad(ir,iHalo) =Vrad(ir,iHalo) +vrr*ww Vrad2(ir,iHalo) =Vrad2(ir,iHalo) +vrr**2*ww Vvrms(ir,iHalo) =Vvrms(ir,iHalo) +vv*ww EndIf jp =Lst(jp) GoTo 10 EndIf EndDo EndDo EndDo Return End SUBROUTINE SNeib !-------------------------------------------------------------- ! Add points around the Box to take into account periodical conditions ! It replicates points periodically and keeps those which are at distance ! less than Radius from the Box. ! N = number of points inside the Box ! Ncp = total number of points SUBROUTINE Points(Radius) !-------------------------------------------------------------- Logical Inside B1 =-Radius B2 =Box+Radius B3 =Box-Radius Ncp =Np Do i=1,Np ! Add dark matter particles x0 =xp(i) y0 =yp(i) z0 =zp(i) Inside = (x0.gt.Radius).and.(x0.lt.B3) Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3)) Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3)) If(.not.Inside)THEN ux =VXp(i) uy =VYp(i) uz =VZp(i) ww=iWeight(i) Do k1 = -1,1 zr =z0+k1*Box If(zr.GT.B1.AND.zr.LT.B2)Then kk =k1**2 Do j1 = -1,1 yr =y0+j1*Box If(yr.GT.B1.AND.yr.LT.B2)Then jj = kk +j1**2 Do i1 = -1,1 xr =x0+i1*Box ii =jj +i1**2 If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)Then Ncp =Ncp +1 If(Ncp.le.Np)Then xp(Ncp) =xr yp(Ncp) =yr zp(Ncp) =zr VXp(Ncp) =ux VYp(Ncp) =uy VZp(Ncp) =uz iWeight(Ncp) =ww EndIf EndIf EndDo EndIf EndDo EndIf EndDo EndIf EndDo If(Ncp.Gt.Nmaxpart)Then write(*,*)' Too many points:',Ncp,' Increase Nmaxpart =', Nmaxpart STOP EndIf Return End SUBROUTINE Points !-------------------------------- Update Statistics of pairs ! SUBROUTINE SPairs !--------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP SCHEDULE (DYNAMIC,kDynamic) Do i=1,Ncenter Call Sneib(i) EndDo Return End SUBROUTINE SPairs !-------------------------------- Write results ! SUBROUTINE WriteData !--------------------------------------------------- Ovscale = float(NROW)**3/Box**3 Vcscale = 6.582e-5 ha = AEXPN*100.*sqrt(Om0/AEXPN**3+Oml0) ! ha at aexpn write (2,13) write (2,14) 13 Format(6x,'Coordinates Mpc',9x,'Velocity km/s',7x, & 'Mass Radius Vrms(3D) Vcirc Npart(Vclimit)Then Ncenter = Ncenter +1 If(iFlag==1)Then Xc(Ncenter) = X0 Yc(Ncenter) = Y0 Zc(Ncenter) = Z0 VXc(Ncenter) = VvX VYc(Ncenter) = VvY VZc(Ncenter) = VvZ Amc(Ncenter) = aM Nparth(Ncenter) = Npart Rmaxh(Ncenter) = Rmax Conc(Ncenter) = Cnc Rmc(Ncenter) = rr/1000. Vcirc(Ncenter) = vcrc Vrms(Ncenter) = vrmss EndIf EndIf ! write(*,'(a,i7,3g12.3,i4)') ' N=',Ncenter,X0, Y0 , Z0,Nbins ! Do j=1,Nb ! read (3,*,err=30,end=30) r2,nn,amass,odens,vv,vc,densd ! ! write(*,'(10x,i7,3g12.3,i4)')j,r2,amass,vc,nn ! If(iFlag==1)Then ! nPr(j,Ncenter) = nn ! aMr(j,Ncenter) = amass ! Vvrms(j,Ncenter) = vv ! EndIf ! EndDo Enddo If(iFlag==0)write(*,*) ' Read Halos=',Ncenter rewind(3) Return 10 write(*,*) ' Could not read radial bins' stop 20 write(*,*) ' Could not read text R(kpc/h)' stop 30 write(*,*) ' Unexpected End of file' stop End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read particles SUBROUTINE ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! N= number of particles ! x(min/max)s = min/max coordinates of small particles !---------------------------------------------------------- Np =0 xmin =1.e10 ; xmax =-1.e10 ymin =1.e10 ; ymax =-1.e10 zmin =1.e10 ; zmax =-1.e10 xmins =1.e10 ; xmaxs =-1.e10 ymins =1.e10 ; ymaxs =-1.e10 zmins =1.e10 ; zmaxs =-1.e10 Lsmall =lspecies(1) If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IROW,In_page,iL,IN) SCHEDULE(DYNAMIC) DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = NPAGE*(IROW-1) if(mod(IROW,1).eq.0)write (*,*)' Read page=',IROW,' Np=',iL ! Loop over particles READ (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ If(iL+In_page > Nmaxpart)Stop ' Too many particles. Increase Nmaxpart' DO IN=1,In_page Xp(iL+IN) =(XPAR(IN)-1.)*Xscale Yp(iL+IN) =(YPAR(IN)-1.)*Xscale Zp(iL+IN) =(ZPAR(IN)-1.)*Xscale VXp(iL+IN) =VX(IN) *Vscale VYp(iL+IN) =VY(IN) *Vscale VZp(iL+IN) =VZ(IN) *Vscale EndDo EndDo Np = lspecies(Nspecies) write(*,*) ' Done Reading Points. N=',Np !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP REDUCTION(MAX:xmax,ymax,zmax,xmaxs,ymaxs,zmaxs) & !$OMP REDUCTION(MIN:xmin,ymin,zmin,xmins,ymins,zmins) Do i =1,Np xmin =min(Xp(i),xmin) xmax =max(Xp(i),xmax) ymin =min(Yp(i),ymin) ymax =max(Yp(i),ymax) zmin =min(Zp(i),zmin) zmax =max(Zp(i),zmax) If(i<=Lsmall)Then xmins =min(Xp(i),xmins) xmaxs =max(Xp(i),xmaxs) ymins =min(Yp(i),ymins) ymaxs =max(Yp(i),ymaxs) zmins =min(Zp(i),zmins) zmaxs =max(Zp(i),zmaxs) EndIf ENDDO write(*,*) ' All Particles :',Np write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write(*,*) ' Small Particles:',Lsmall write(*,*) ' x:',xmins,xmaxs write(*,*) ' y:',ymins,ymaxs write(*,*) ' z:',zmins,zmaxs write (*,*) ' Nmaxpart=',Nmaxpart do i=1,1000,200 write (*,'(i8," Test Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i) EndDo Close(21) Return End Subroutine ReadPnt !-------------------------------------------------- ! Set Weights of particles for ! fast access SUBROUTINE SetWeights !-------------------------------------------------- Wpr =NGRID**3/FLOAT(NROW)**3 N_particles =min(lspecies(Nspecies),Nmaxpart) jstart =1 If(N_particles.le.0)Then write (*,*) ' Wrong number of particles:',N_particles write (*,*) ' Nspecies=',Nspecies,' N=',lspecies STop endif Do j=1,Nspecies jend =min(lspecies(j),Nmaxpart) ww = wspecies(j)/Wpr !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (k) Do k=jstart ,jend iWeight(k) =ww EndDo jstart =jend EndDo weightSmall = iWeight(1) RETURN END SUBROUTINE SetWeights ! ------------------------ real function seconds () ! ------------------------ ! ! purpose: returns elapsed time in seconds Integer, SAVE :: first=0,rate=0 Real , SAVE :: t0=0. !real*8 dummy !real tarray(2) If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 t0 =float(i)/float(rate) seconds = 0. Else ! seconds = etime(tarray) CALL SYSTEM_CLOCK(i) seconds = float(i)/float(rate)-t0 EndIf return end function seconds !--------------------------------------------------- ! Read current data from disk/tape, ! Open files ! Nrecl is the number of values in a record ! Npage is the number of particles in a page SUBROUTINE RDTAPE(file1,file2) !--------------------------------------------------- ! Open file on a tape Character (LEN=*) :: file1,file2 OPEN(UNIT=9,FILE=file1, & FORM='UNFORMATTED', STATUS='UNKNOWN') ! Read control information ! and see whether it has proper ! structure READ (9,err=20,end=20) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras WRITE (*,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN,EKIN1,EKIN2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv WRITE (2,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN,EKIN1,EKIN2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv 100 FORMAT(1X,'Header=>',A45,/ & 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ & 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/ & 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ & 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/& 1x,' Omega_curvature=',F7.3) IF(NROW.NE.NROWC) THEN WRITE (*,*) & ' NROW in PARAMETER and in TAPE-FILE are different' write (*,*) ' Nrow=',NROW,' NROWC=',NROWC ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) & ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF ! Open work file on disk write (*,*) ' start openning PMcrs0' 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE=file2,ACCESS='DIRECT',STATUS='UNKNOWN',RECL=NACCES) REWIND 9 write (*,*) ' done openning PMcrs0' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END SUBROUTINE RDTAPE !-------------------------------------------- ! Read current PAGE of particles from disk ! NRECL - length of ROW block in words SUBROUTINE GETROW(IROW,Ifile) !-------------------------------------------- READ (20+Ifile,REC=IROW) RECDAT RETURN END SUBROUTINE GETROW end module SetHaloParticles !------------------------------------------------- ! ! PROGRAM HaloParticles use SetHaloParticles Character*120 :: file1,file2,catname,listname,txt*8 ! ------------------------------------------ Read data write(*,'(a,$)')' Enter xxxxx for snapshot (PMcrs0xxxxx.DAT)=' read(*,'(a)')txt file1 ='PMcrd'//TRIM(txt)//'.DAT' file2 ='PMcrs0'//TRIM(txt)//'.DAT' ! catname='Catshort.999.DAT' ! catname='CatshortB.490.DAT' ! catname='CatshortA.3303.DAT' ! catname='CatshortA486.DAT' listname='Halolist.'//TRIM(txt)//'.DAT' Open( 2,file=listname,Status='UNKNOWN')! output: list of halos write(*,'(a,$)')' Enter xxxxx for catalog (Catshortxxxxx.DAT)=' read(*,'(a)')txt catname ='Catshort'//TRIM(txt)//'.DAT' write(*,'(a,$)')' Enter limit on circular velocity of halos (km/s)=' read(*,*) Vclimit write(*,'(a,$)')' Enter radius of the first bin (kpch) =' read(*,*) R0 Open( 3,file=catname,Status='UNKNOWN')! short list of halos tt = seconds() CALL RDTAPE(file1,file2) CALL SetWeights If(extras(100).gt.0.)Then Box =extras(100) Else write(*,'(a,$)')' Enter box size in comoving Mpc/h =' read(*,*) Box EndIf WRITE (*,120) & AEXPN,ISTEP,NROW,NGRID write (*,*) ' Number of mass Species =',Nspecies write (*,*) ' Total number of particles=',lspecies(Nspecies) write (*,*) ' Number of small particles=',lspecies(1) write (*,*) ' Box size (Mpc/h) = ',Box 120 FORMAT(1X, & 5x,'Expansion Parameter= ',F8.3,' Step= ',i5,/ & 5x,'Nrow= ',I4,' Ngrid= ',I4) hsmall= hubble ! Hubble/100. Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= 100.*Xscale/AEXPN ! Scale for velocities Dscale= 2.746e+11*Om0*(Box/NROW)**3 *weightSmall ! mass scale write (*,'(6x,"Mass of smallest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') Dscale write(*,'(6x,"first species=",2g11.3)') wspecies(1),weightSmall If(Nspecies.gt.1)Then write (*,'(6x,"Mass of largest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') & Dscale*wspecies(Nspecies)/(NGRID**3/FLOAT(NROW)**3) EndIf Dscale= Dscale/weightSmall ! scale factor for particles R0 =R0/1.e3 ! scale to Mpch Do ir=1,Nrad ! all scaled to hubble, not h^(-1) Rad(ir) = R0*((ir-1)/4.+1.)**2! scale radii to comoving Mpc EndDo write (*,*) ' Min and Max Radii =',Rad(1),Rad(Nrad) ! ------------------------------------------ Read data CALL ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! particles CALL ReadHalos(0) ! count halos Rmax = Rad(Nrad) CALL Points(Rmax) ! Make more points periodically write (*,*) 'Points made' !write(*,*) '******* time=',seconds() ! set limits for linker list Celloptimal = Rad(Nrad)/FracSearch Cell = max(Box/Nlinker,Celloptimal) Nmx = INT(xmins/Cell)-Nbf ! boundaries for linker list Nbx = INT(xmaxs/Cell)+Nbf Nmy = INT(ymins/Cell)-Nbf Nby = INT(ymaxs/Cell)+Nbf Nmz = INT(zmins/Cell)-Nbf Nbz = INT(zmaxs/Cell)+Nbf write(*,'(2(a,f8.4))') ' Cell for linker list(Mpch)=',Cell,' Optimal=',Celloptimal write(*,'(a,3(2x,2i4))') ' Dimensions for the linker list=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz ww = 4.*((Nbx-Nmx)*(Nby-Nmy)*(Nbz-Nmz)+Ncp)/(1024.**2) write (*,'(a,f9.2)') ' Requested memory for linker-list(Mb)=',ww ALLOCATE (Lst(Ncp),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz),Stat=iStat) If(iStat.ne.0)Then write (*,*) ' Not enough memory to allocate linker list tables' write (*,'(a,8i4)') ' requested limits:', & Nmx,Nbx,Nmy,Nby,Nmz,Nbz ww = -4.*(Nmx-Nbx)*(Nmy-Nby)*(Nmz-Nbz)/1024.**2 write (*,'(a,f9.2)') ' Requested memory(Mb)=',ww Stop EndIf !write(*,*) '******* time=',seconds() Call List ! make linker List, Label !write(*,*) '******* time=',seconds() ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter), Stat =iStat0) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter), Stat =iStat1) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter), Stat =iStat2) ALLOCATE (Nparth(Ncenter),Rmaxh(Ncenter),Conc(Ncenter), Stat =iStat3) If(iStat0**2+iStat1**2+iStat2**2+iStat3**2.ne.0)Then write(*,*) ' Not enough memory for halo headers.',iStat0,iStat1,iStat2 stop EndIf ALLOCATE (Vvrms(Nrad,Ncenter), Stat =iStat0) ALLOCATE (nPr(Nrad,Ncenter),aMr(Nrad,Ncenter), Stat =iStat1) ALLOCATE (Vrad(Nrad,Ncenter),Vrad2(Nrad,Ncenter), Stat =iStat2) ALLOCATE (RadPr(Nrad,Ncenter), Stat =iStat3) If(iStat0**2+iStat1**2+iStat2**2+iStat3**2.ne.0)Then write(*,*) ' Not enough memory for halo profiles.',iStat0,iStat1 stop EndIf Ncenter = 0 ! initialize halo counter CALL ReadHalos(1) ! read halos Close(3) aMr =0. ! initialize profiles Vvrms =0. Vrad =0. nPr =0. RadPr =0. write (*,20) Np,Ncenter,Box,Vclimit write (2,20) Np,Ncenter,Box,Vclimit 20 Format(5x,'Statistics for ',i9,' points. Maxima=',i8, & 3x,'Box(Mpc)=',f6.1,' Vc_limit_halos(kms)=',f9.5) write(*,*) '--------- time=',seconds(),' Get statistics' CALL SPairs ! Accumulate statistics write(*,*) '--------- time=',seconds(),' Start Write Results' CALL WriteData write(*,*) '-------- time(sec)=',seconds() End PROGRAM HaloParticles