!------------------------------------------------- ! ! MODULE setHaloParticles PUBLIC Integer, PARAMETER :: NROW =2048, & NGRID =256, & 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 = 20000000, & ! max # of particles Nlinker = 200, & ! linker-list dimension FracSearch = 5., & ! Nhmax = 150, & ! max # of selected particles in a halo Srms = 0.2, & ! fraction of rms velocity for bound particles nwant = 10000, & ! pool of preselected particles Kdynamic= 100 ! OMP chunk Real :: INPUT,Mhalo,Ovlim,Ovcnt Real, DIMENSION(Nrad) :: Rad(Nrad)=0. Real, DIMENSION(Nmaxpart) :: Xp,Yp,Zp,VXp,VYp,Vzp,iWeight Integer, DIMENSION(Nmaxpart) ::Lagrang Real :: Box=0.,Xscale=0.,Vscale=0.,Dscale=0.,Cell=0., weightSmall Integer :: Np=0,Ncp=0,Ncenter=0, & Nmx,Nbx,Nmy,Nby,Nmz,Nbz Real, ALLOCATABLE, DIMENSION (:) :: Xpp,Ypp,Zpp,VXpp,VYpp,Vzpp,Wpp Integer, ALLOCATABLE, DIMENSION (:) :: LagrangP COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, & Ocurv,extras(100) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER Real :: RECDAT(NRECL),wspecies(10) Integer :: lspecies(10) Integer, ALLOCATABLE, DIMENSION (:,:,:) :: Label Integer, ALLOCATABLE, DIMENSION (:) :: Lst,nBins,nSelect Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc Real, ALLOCATABLE, DIMENSION (:,:) :: aMr,Vvrms,GravPot Integer, ALLOCATABLE, DIMENSION (:,:) :: nPr,nBound EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/) Contains !-------------------------------------------------------------- ! Make linker lists of particles in each cell SUBROUTINE List !-------------------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncp Lst(i)=-1 EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ncp i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) k=INT(Zp(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) k=MIN(MAX(Nmz,k),Nbz) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo write(*,*) ' Made first list' ALLOCATE(Xpp(Ncp),Ypp(Ncp),Zpp(Ncp),Stat=iStat) If(iStat.ne.0)Then write(*,*)'Not enough memory for Xpp in List. No Rearrangment' return EndIf ALLOCATE(VXpp(Ncp),VYpp(Ncp),VZpp(Ncp),Stat=iStat) If(iStat.ne.0)Then write(*,*)'Not enough memory for Vpp in List. No Rearrangment' DEALLOCATE (Xpp,Ypp,Zpp) return EndIf ALLOCATE(LagrangP(Ncp),Wpp(Ncp),Stat=iStat) If(iStat.ne.0)Then write(*,*)'Not enough memory for LagrangP in List. No Rearrangment' DEALLOCATE (Xpp,Ypp,Zpp,VXpp,VYpp,VZpp) return EndIf nn = 0 ! main particles should be first Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx jp =Label(i,j,k) Do while(jp.ne.0) If(jp.le.Np)Then nn =nn +1 LagrangP(nn) =Lagrang(jp) Wpp(nn) =iWeight(jp) Xpp(nn) =Xp(jp) Ypp(nn) =Yp(jp) Zpp(nn) =Zp(jp) VXpp(nn) =VXp(jp) VYpp(nn) =VYp(jp) VZpp(nn) =VZp(jp) EndIf jp =Lst(jp) End DO EndDo EndDo EndDo If(nn.ne.Np)Then write(*,*) ' error in rearranging particles in List' write(*,*) ' New number of particles is=',nn write(*,*) ' It must be equal to old =',Np Stop EndIf ! Restore particles !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Np ! (1,Ncp) changed ! to Np for parallelization Lagrang(i) =LagrangP(i) iWeight(i) =Wpp(i) Xp(i) =Xpp(i) Yp(i) =Ypp(i) Zp(i) =Zpp(i) VXp(i) =VXpp(i) VYp(i) =VYpp(i) VZp(i) =VZpp(i) EndDo DEALLOCATE (Xpp,Ypp,Zpp,VXpp,VYpp,VZpp,LagrangP,Wpp) write(*,*) ' re-arranged particles' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncp Lst(i)=-1 EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ncp i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) k=INT(Zp(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) k=MIN(MAX(Nmz,k),Nbz) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo write(*,*) ' Made second list' Return End SUBROUTINE List !-------------------------------------------------------------- ! Find all neibhours for a center ! Xc,Yc,Zc - center; ! SUBROUTINE SNeib(iHalo) !-------------------------------------------------------------- Real, DIMENSION(nwant) :: Rpart,Ek,Et,Ec Integer, DIMENSION(nwant) :: Lpart,ind Integer :: Labpart(1) x = Xc(iHalo) ; y = Yc(iHalo); z = Zc(iHalo) Vx1=VXc(iHalo) ; Vy1= VYc(iHalo); Vz1= VZc(iHalo) Lpart = 0 Labpart = 0 Ek = 1.e10 Et = 1.e10 nn =0 mm =0 kbin =Nbins(iHalo) RhaloMax =Rad(kbin) Do i=1,kbin nn = nPr(i,iHalo) mm = i If(Rad(i).gt.Rmc(iHalo).or.nn.ge.nwant)exit EndDo If(nn.eq.0)Then Nbound(1,iHalo) =0 ! no particles selected return EndIf Radmax = Rad(mm) dmax = Radmax**2 d0 = Rad(1) ! write(*,'(3i6,3g12.3)')iHalo,nn,mm,Radmax,Cell ! write(*,'(10x,2(3x,3f8.2))') x,y,z,Vx1,Vy1,Vz1 nn = 0 ! counter for pre-selected particles ! limits for Label i1=INT((x-Radmax)/Cell) j1=INT((y-Radmax)/Cell) k1=INT((z-Radmax)/Cell) i1=MIN(MAX(Nmx,i1),Nbx) j1=MIN(MAX(Nmy,j1),Nby) k1=MIN(MAX(Nmz,k1),Nbz) i2=INT((x+Radmax)/Cell) j2=INT((y+Radmax)/Cell) k2=INT((z+Radmax)/Cell) i2=MIN(MAX(Nmx,i2),Nbx) j2=MIN(MAX(Nmy,j2),Nby) k2=MIN(MAX(Nmz,k2),Nbz) ! Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0.and.nn10)STOP Return End SUBROUTINE SNeib !-------------------------------------------------------------- ! Add points around the Box to take into account periodical conditions ! It replicates points periodically and keeps those which are at distance ! less than Radius from the Box. ! N = number of points inside the Box ! Ncp = total number of points SUBROUTINE Points(Radius) !-------------------------------------------------------------- Logical Inside B1 =-Radius B2 =Box+Radius B3 =Box-Radius Ncp =Np Do i=1,Np ! Add dark matter particles x0 =xp(i) y0 =yp(i) z0 =zp(i) Inside = (x0.gt.Radius).and.(x0.lt.B3) Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3)) Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3)) If(.not.Inside)THEN ux =VXp(i) uy =VYp(i) uz =VZp(i) ww=iWeight(i) Do k1 = -1,1 zr =z0+k1*Box If(zr.GT.B1.AND.zr.LT.B2)Then kk =k1**2 Do j1 = -1,1 yr =y0+j1*Box If(yr.GT.B1.AND.yr.LT.B2)Then jj = kk +j1**2 Do i1 = -1,1 xr =x0+i1*Box ii =jj +i1**2 If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)Then Ncp =Ncp +1 If(Ncp.le.Np)Then xp(Ncp) =xr yp(Ncp) =yr zp(Ncp) =zr VXp(Ncp) =ux VYp(Ncp) =uy VZp(Ncp) =uz iWeight(Ncp) =ww Lagrang(Ncp) =i EndIf EndIf EndDo EndIf EndDo EndIf EndDo EndIf EndDo If(Ncp.Gt.Nmaxpart)Then write(*,*)' Too many points:',Ncp,' Increase Nmaxpart =', Nmaxpart STOP EndIf Return End SUBROUTINE Points !-------------------------------- Update Statistics of pairs ! SUBROUTINE SPairs !--------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP SCHEDULE (DYNAMIC,kDynamic) Do i=1,Ncenter Call Sneib(i) EndDo Return End SUBROUTINE SPairs !-------------------------------- Manipulate halo profiles ! SUBROUTINE SmHalo(iHalo) !--------------------------------------------------- Real, Dimension(Nrad) :: tmp kbin =Nbins(iHalo) Vv = Vrms(iHalo) ! smooth rms profile Do i=1,kbin if(Vvrms(i,iHalo)<=1.e-5)Vvrms(i,iHalo)=Vv EndDo Vvrms(1,iHalo) = (Vvrms(1,iHalo)+Vvrms(2,iHalo))/2. Do i=2,kbin-1 tmp(i) =(Vvrms(i-1,iHalo)+Vvrms(i-1,iHalo)+Vvrms(i+1,iHalo))/3. EndDo Vvrms(2:kbin-1,iHalo) = tmp(2:kbin-1) ! Gravitational Potential spherical density ! Units : (km/s)**2 GravPot(kbin,iHalo) = -4.333e-9*aMr(kbin,iHalo)/Rad(kbin) Do k=kbin-1,1,-1 GravPot(k,iHalo) =GravPot(k+1,iHalo) & -4.333e-9*(aMr(k+1,iHalo)-aMr(k,iHalo)) & /(Rad(k)+Rad(k+1))*2. EndDo Return End SUBROUTINE SmHalo !-------------------------------- Manipulate Halos ! SUBROUTINE SmoothHalos !--------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP SCHEDULE (DYNAMIC,kDynamic) Do i=1,Ncenter Call SmHalo(i) EndDo Return End SUBROUTINE SmoothHalos !-------------------------------- Update Statistics of pairs ! SUBROUTINE WriteData !--------------------------------------------------- Do i=1,Ncenter write(2,10)i,Xc(i),Yc(i),Zc(i),VXc(i),VYc(i),VZc(i), & Amc(i),Rmc(i)*1000.,Vrms(i),Vcirc(i), & nSelect(i),(nBound(j,i),j=1,nSelect(i)) EndDo 10 format(i7,3F9.4,3F8.1,g11.3,f8.2,2f7.1,I6,500i8) Return End SUBROUTINE WriteData !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag) !-------------------------------------------------------------- Character :: Line*80, Txt*8='R(kpc/h)', Txt2*8,txt1*4 Do i=1,2 Read(3,'(a)')Line If(iFlag==0)write(*,'(a)') Line EndDo Read(3,'(a4,i4)')txt1,ii If(ii.ne.ISTEP)Then write(*,*)' Error: Snapshot and Catalog are for different time-steps' write(*,'(2(a,i5))') ' Catalog is for ',ii,' Particles are for ',ISTEP stop EndIf Do i=1,11 Read(3,'(a)')Line EndDo Read(3,*,err=10) Rad Rad =Rad /1000. ! convert to Mpch If(iFlag==0)Then write(*,*) '----- Radial bins(Mpch) ----' write(*,'(10f8.3)') Rad EndIf Do read(3,'(a)',err=20)Txt2 If(Txt2==Txt)Exit EndDo Do read (3,*,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,Npart,Rmax,Conc,Nb If(iStat.ne.0)Exit Ncenter = Ncenter +1 If(iFlag==1)Then Xc(Ncenter) = X0 Yc(Ncenter) = Y0 Zc(Ncenter) = Z0 VXc(Ncenter) = VvX VYc(Ncenter) = VvY VZc(Ncenter) = VvZ Amc(Ncenter) = aM Rmc(Ncenter) = rr/1000. Vcirc(Ncenter) = vcrc Vrms(Ncenter) = vrmss Nbins(Ncenter) = Nb EndIf ! write(*,'(a,i7,3g12.3,i4)') ' N=',Ncenter,X0, Y0 , Z0,Nbins Do j=1,Nb read (3,*,err=30,end=30) r2,nn,amass,odens,vv,vc,densd ! write(*,'(10x,i7,3g12.3,i4)')j,r2,amass,vc,nn If(iFlag==1)Then nPr(j,Ncenter) = nn aMr(j,Ncenter) = amass Vvrms(j,Ncenter) = vv EndIf EndDo Enddo If(iFlag==0)write(*,*) ' Read Halos=',Ncenter rewind(3) Return 10 write(*,*) ' Could not read radial bins' stop 20 write(*,*) ' Could not read text R(kpc/h)' stop 30 write(*,*) ' Unexpected End of file' stop End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read particles SUBROUTINE ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! N= number of particles ! x(min/max)s = min/max coordinates of small particles !---------------------------------------------------------- Np =0 xmin =1.e10 ; xmax =-1.e10 ymin =1.e10 ; ymax =-1.e10 zmin =1.e10 ; zmax =-1.e10 xmins =1.e10 ; xmaxs =-1.e10 ymins =1.e10 ; ymaxs =-1.e10 zmins =1.e10 ; zmaxs =-1.e10 Lsmall =lspecies(1) If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IROW,In_page,iL,IN) SCHEDULE(DYNAMIC) DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = NPAGE*(IROW-1) if(mod(IROW,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 Lagrang(i) = i xmin =min(Xp(i),xmin) xmax =max(Xp(i),xmax) ymin =min(Yp(i),ymin) ymax =max(Yp(i),ymax) zmin =min(Zp(i),zmin) zmax =max(Zp(i),zmax) If(i<=Lsmall)Then xmins =min(Xp(i),xmins) xmaxs =max(Xp(i),xmaxs) ymins =min(Yp(i),ymins) ymaxs =max(Yp(i),ymaxs) zmins =min(Zp(i),zmins) zmaxs =max(Zp(i),zmaxs) EndIf ENDDO write(*,*) ' All Particles :',Np write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write(*,*) ' Small Particles:',Lsmall write(*,*) ' x:',xmins,xmaxs write(*,*) ' y:',ymins,ymaxs write(*,*) ' z:',zmins,zmaxs write (*,*) ' Nmaxpart=',Nmaxpart do i=1,1000,200 write (*,'(i8," Test Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i) EndDo Close(21) Return End Subroutine ReadPnt !-------------------------------------------------- ! Set Weights of particles for ! fast access SUBROUTINE SetWeights !-------------------------------------------------- Wpr =NGRID**3/FLOAT(NROW)**3 N_particles =min(lspecies(Nspecies),Nmaxpart) jstart =1 If(N_particles.le.0)Then write (*,*) ' Wrong number of particles:',N_particles write (*,*) ' Nspecies=',Nspecies,' N=',lspecies STop endif Do j=1,Nspecies jend =min(lspecies(j),Nmaxpart) ww = wspecies(j)/Wpr !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (k) Do k=jstart ,jend iWeight(k) =ww EndDo jstart =jend EndDo weightSmall = iWeight(1) RETURN END SUBROUTINE SetWeights ! ------------------------ real function seconds () ! ------------------------ ! ! purpose: returns elapsed time in seconds Integer, SAVE :: first=0,rate=0 Real , SAVE :: t0=0. !real*8 dummy !real tarray(2) If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 t0 =float(i)/float(rate) seconds = 0. Else ! seconds = etime(tarray) CALL SYSTEM_CLOCK(i) seconds = float(i)/float(rate)-t0 EndIf return end function seconds !--------------------------------------------------- ! Read current data from disk/tape, ! Open files ! Nrecl is the number of values in a record ! Npage is the number of particles in a page SUBROUTINE RDTAPE(file1,file2) !--------------------------------------------------- ! Open file on a tape Character (LEN=*) :: file1,file2 OPEN(UNIT=9,FILE=file1, & FORM='UNFORMATTED', STATUS='UNKNOWN') ! Read control information ! and see whether it has proper ! structure READ (9,err=20,end=20) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras WRITE (*,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN,EKIN1,EKIN2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv WRITE (2,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN,EKIN1,EKIN2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv 100 FORMAT(1X,'Header=>',A45,/ & 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ & 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/ & 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I9,/ & 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/& 1x,' Omega_curvature=',F7.3) IF(NROW.NE.NROWC) THEN WRITE (*,*) & ' NROW in PARAMETER and in TAPE-FILE are different' write (*,*) ' Nrow=',NROW,' NROWC=',NROWC ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) & ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF ! Open work file on disk write (*,*) ' start openning PMcrs0' 10 NBYTE = NRECL*4 NACCES= NBYTE / nbyteword OPEN(UNIT=21,FILE=file2,ACCESS='DIRECT',STATUS='UNKNOWN',RECL=NACCES) REWIND 9 write (*,*) ' done openning PMcrs0' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END SUBROUTINE RDTAPE !-------------------------------------------- ! Read current PAGE of particles from disk ! NRECL - length of ROW block in words SUBROUTINE GETROW(IROW,Ifile) !-------------------------------------------- READ (20+Ifile,REC=IROW) RECDAT RETURN END SUBROUTINE GETROW end module SetHaloParticles !------------------------------------------------- ! ! PROGRAM HaloParticles use SetHaloParticles Character*120 :: file1,file2,catname,listname,txt*8 ! ------------------------------------------ Read data write(*,'(a,$)')' Enter snapshot=' read(*,'(a)')txt file1 ='FILES/PMcrd_'//TRIM(txt)//'.DAT' file2 ='FILES/PMcrs0_'//TRIM(txt)//'.DAT' catname='CATALOGS/Catalog.'//TRIM(txt)//'.DAT' listname='Halolist.'//TRIM(txt)//'.DAT' Open( 2,file=listname,Status='UNKNOWN')! short list of halos Open( 3,file=catname,Status='UNKNOWN')! short list of halos tt = seconds() CALL RDTAPE(file1,file2) CALL SetWeights If(extras(100).gt.0.)Then Box =extras(100) Else write(*,'(a,$)')' Enter box size in comoving Mpc/h =' read(*,*) Box EndIf WRITE (*,120) & AEXPN,ISTEP,NROW,NGRID write (*,*) ' Number of mass Species =',Nspecies write (*,*) ' Total number of particles=',lspecies(Nspecies) write (*,*) ' Number of small particles=',lspecies(1) write (*,*) ' Box size (Mpc/h) = ',Box 120 FORMAT(1X, & 5x,'Expansion Parameter= ',F8.3,' Step= ',i5,/ & 5x,'Nrow= ',I4,' Ngrid= ',I4) hsmall= hubble ! Hubble/100. Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= 100.*Xscale/AEXPN ! Scale for velocities Dscale= 2.746e+11*Om0*(Box/NROW)**3 *weightSmall ! mass scale write (*,'(6x,"Mass of smallest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') Dscale write(*,'(6x,"first species=",2g11.3)') wspecies(1),weightSmall If(Nspecies.gt.1)Then write (*,'(6x,"Mass of largest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') & Dscale*wspecies(Nspecies)/(NGRID**3/FLOAT(NROW)**3) EndIf ! ------------------------------------------ Read data CALL ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! particles CALL ReadHalos(0) ! count halos Rmax = Rad(Nrad) CALL Points(Rmax) ! Make more points periodically write (*,*) 'Points made' !write(*,*) '******* time=',seconds() ! set limits for linker list Celloptimal = Rad(Nrad)/FracSearch Cell = max(Box/Nlinker,Celloptimal) Nmx = INT(xmins/Cell)-Nbf ! boundaries for linker list Nbx = INT(xmaxs/Cell)+Nbf Nmy = INT(ymins/Cell)-Nbf Nby = INT(ymaxs/Cell)+Nbf Nmz = INT(zmins/Cell)-Nbf Nbz = INT(zmaxs/Cell)+Nbf write(*,'(2(a,f8.4))') ' Cell for linker list(Mpch)=',Cell,' Optimal=',Celloptimal write(*,'(a,3(2x,2i4))') ' Dimensions for the linker list=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz ww = 4.*((Nbx-Nmx)*(Nby-Nmy)*(Nbz-Nmz)+Ncp)/(1024.**2) write (*,'(a,f9.2)') ' Requested memory for linker-list(Mb)=',ww ALLOCATE (Lst(Ncp),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz),Stat=iStat) If(iStat.ne.0)Then write (*,*) ' Not enough memory to allocate linker list tables' write (*,'(a,8i4)') ' requested limits:', & Nmx,Nbx,Nmy,Nby,Nmz,Nbz ww = -4.*(Nmx-Nbx)*(Nmy-Nby)*(Nmz-Nbz)/1024.**2 write (*,'(a,f9.2)') ' Requested memory(Mb)=',ww Stop EndIf !write(*,*) '******* time=',seconds() Call List ! make linker List, Label !write(*,*) '******* time=',seconds() ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter), Stat =iStat0) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter),Nbins(Ncenter), Stat =iStat1) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter), Stat =iStat2) If(iStat0**2+iStat1**2+iStat2**2.ne.0)Then write(*,*) ' Not enough memory for halo headers.',iStat0,iStat1,iStat2 stop EndIf ALLOCATE (aMr(Nrad,Ncenter),Vvrms(Nrad,Ncenter), Stat =iStat0) ALLOCATE (nPr(Nrad,Ncenter),GravPot(Nrad,Ncenter), Stat =iStat1) If(iStat0**2+iStat1**2.ne.0)Then write(*,*) ' Not enough memory for halo profiles.',iStat0,iStat1 stop EndIf Ncenter = 0 ! initialize halo counter CALL ReadHalos(1) ! read halos Close(3) CALL SmoothHalos write (*,20) Np,Ncenter,Box,Cell write (2,20) Np,Ncenter,Box,Cell 20 Format(5x,'Statistics for ',i9,' points. Maxima=',i8, & 3x,'Box(Mpc)=',f6.1,' Cell for linker list(Mpch)=',f9.5) ALLOCATE (nBound(Nhmax,Ncenter),nSelect(Ncenter),Stat=iStat) If(iStat.ne.0)Then write(*,*) ' Not enough memory for halo profiles.',iStat0,iStat1 stop EndIf !write(*,*) '******* time=',seconds() CALL SPairs ! Accumulate statistics !write(*,*) '******* time=',seconds() CALL WriteData write(*,*) '******* time(sec)=',seconds() End PROGRAM HaloParticles