!---------------------------------------------------------------------------------------- ! Find properties of large distinct halos ! Do not remove unbound particles ! Input: Halo Catalog (short list) ! Dark matter particles ! Output: Halo catalog ! Parameters: ! outer radius in units of virial radius ! number of particles for potential energy ! MODULE setHaloParticles PUBLIC Integer, PARAMETER :: & NROW = 4096, & NGRID = 256, & nbyteword = 4, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6, & Nrad = 100, & ! max # of radial bins Nbf = 10, & ! width of buffer for linker list Nlinker = 300, & ! linker-list size Nhmax = 20000, & ! max # of selected particles for Epotential Nmax_s = 50, & ! max # of selected particles for Bound nwant = 20000, & ! pool of preselected particles for Bound kDynamic = 1000 ! OMP chunk for parallelization Integer*8, PARAMETER :: & Nmaxpart = 2048_8**3*1.3 ! max # of particles Real, PARAMETER :: & FracSearch = 5., & ! Rout = 1., & ! outer radius in units halo virial radius Risolate = 1., & ! isolation: no halos inside Risolate*Rvir Srms = 0.2, & ! fraction of rms velocity for bound particles Unbound = 10., & ! remove particles V>Unbound*Vcirc aMmin_dist = 5.e10, & ! min mass of distinct halos dLog = 0.0667, & ! bin size in log(r) epsilon = 2., & ! force softening for Epot in kpc units eps2 = (epsilon/1000.)**2 ! Real :: INPUT,Mhalo,Ovlim,Ovcnt,Rdmax Real, DIMENSION(Nhmax) :: dXe,dYe,dZe,we Real, DIMENSION(Nmaxpart) :: Xp,Yp,Zp,VXp,VYp,Vzp,iWeight Integer*8, DIMENSION(Nmaxpart) ::Lagrang Real :: Box=0.,Xscale=0.,Vscale=0.,Dscale=0.,Cell=0., weightSmall Integer*8 :: Np=0,Ncp=0 Integer :: Ncenter=0, & Nmx,Nbx,Nmy,Nby,Nmz,Nbz,Nbin_max,Nbin_min COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EK,EK1,EK2,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*8 :: lspecies(10) Integer*8, ALLOCATABLE, DIMENSION (:,:,:) :: Label Integer*8, ALLOCATABLE, DIMENSION (:,:) :: nBound Integer*8, ALLOCATABLE, DIMENSION (:) :: Lst Integer*4, ALLOCATABLE, DIMENSION (:) :: iSeed,LDistinct,nSelect Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc, & Amh,dRh,Ekin,Epot, & Anhx,Anhy,Anhz, & Vhrad Real, ALLOCATABLE, DIMENSION (:,:) :: aMr,Vvrms,Rbin,aJx,aJy,aJz,Vrad,Vrad2 Integer, ALLOCATABLE, DIMENSION (:,:) :: nPr EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/,dXe,dYe,dZe,we) Contains !------------------------------------------------------------------------- ! Make linker lists of particles in each cell ! Start with sorting particles by cells: ! particles in the same cell are in contiguous memory ! Need additional memory for particles: Dummy, iDummy arrays SUBROUTINE List !------------------------------------------------------------------------ !Real, ALLOCATABLE, DIMENSION (:) :: Dummy !Integer*8, ALLOCATABLE, DIMENSION (:) :: iDummy Integer*8 :: i,jp,nn !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncp Lst(i)= -1 EndDo write(*,*) ' Lst initialized' !$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 write(*,*) ' Label initialized' !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (jp,i,j,k,k1) ! Do k1 =Nmz,Nbz ! write(*,*)' k1=',k1 Do jp=1,Ncp k=INT(Zp(jp)/Cell) k=MIN(MAX(Nmz,k),Nbz) ! If(k.eq.k1)Then i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp ! EndIf EndDo ! EndDo write(*,*) ' Made first list' ! return ! ALLOCATE(iDummy(Ncp)) ! pointers: old position -> new position in arrays ! write(*,*) ' --- working on rearranging particles so that particles are sorted by cells' ! nn = 0 ! Do k=Nmz,Nbz ! !write (*,*) ' re-arrange k=',k ! Do j=Nmy,Nby ! Do i=Nmx,Nbx ! jp =Label(i,j,k) ! Do while(jp.ne.0) ! If(jp.le.Ncp)Then ! nn =nn +1 ! iDummy(jp) = nn ! assign the pointer ! EndIf ! jp =Lst(jp) ! End DO ! EndDo ! EndDo ! EndDo ! If(nn.ne.Ncp)Then ! write(*,*) ' error in rearranging particles in List' ! write(*,*) ' New number of particles is=',nn ! write(*,*) ' It must be equal to old =',Ncp ! Stop ! EndIf ! write(*,*) ' --- sort particles: move to a buffer, then place back ' ! ALLOCATE(Dummy(Ncp)) ! buffer for other arrays ! X ------------------------------------- !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = Xp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Xp(i) = Dummy(i) ! EndDo ! Y !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = Yp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Yp(i) = Dummy(i) ! EndDo ! Z !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = Zp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Zp(i) = Dummy(i) ! EndDo ! VX------------------------------------- !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = VXp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! VXp(i) = Dummy(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) =VYp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! VYp(i) = Dummy(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = VZp(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! VZp(i) = Dummy(i) ! EndDo ! ! Weight ------------------------------------- !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! Dummy(iDummy(i)) = iWeight(i) ! EndDo !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (i) ! Do i=1,Ncp ! iWeight(i) = Dummy(i) ! EndDo ! DEALLOCATE (Dummy, iDummy) ! write(*,*) ' Finished re-arranging particles for data locality' ! write(*,*) ' --- Now start making final linker 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 second list: particles are sorted and linked' ! Return End SUBROUTINE List !------------------------------------------------------------------------- ! Make linker list of Halos SUBROUTINE ListHalos !------------------------------------------------------------------------ Integer*8 :: i,jp !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncenter 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,Ncenter i=INT(Xc(jp)/Cell) j=INT(Yc(jp)/Cell) k=INT(Zc(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 ListHalos !!------------------------------------------------------------------------- !! Estimate potential energy for selected particles !! Input: nn - number of particles in the list !! dXe,dYe,dZe, we = coordinates and weights !! epsilon = force softening !! Output: sum of pairs Sum(i=1,N-1; j=i+1,N){mi*mj/} ! Function GetEpot(nn) !!------------------------------------------------------------------------! !Real*8 :: ee ! ee =0. ! Do i=1,nn-1 ! w = we(i) ! x =dXe(i) ; y =dYe(i) ; z =dZe(i) ! Do j=i+1,nn ! ee = ee + w*we(j)/( & ! (dXe(j)-x)**2+(dYe(j)-y)**2+(dZe(j)-z)**2+eps2) ! EndDo ! j ! EndDo ! i ! GetEpot = ee ! Return ! End Function GetEpot !-------------------------------------------------------------- ! Make list of distinct halos ! Input: Ncenter - numer of halos, {Xc ...} = halo data ! Cell = size of linker-list cell ! Radmax = maximum halo radius ! Output: LDistinct(i) = 1 for distinct halo ! SUBROUTINE HaloDistinct !-------------------------------------------------------------- Radmax =Rdmax /Rout *2. ! maximum virial radius Lx = Nbx ! length of periodical box Ly = Nby Lz = Nbz !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncenter LDistinct(i) =1 EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,i,j,k,x,y,z,R0,i1,j1,k1,i2,j2,k2,jp,dd,i3,j3,k3,aM0) Do ic=1,Ncenter x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 =Rmc(ic) aM0 =Amc(ic) 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 k3=k1,k2 k = k3 ! periodical boundaries If(k >Lz)k = k-Lz If(k <0) k = k+Lz If(k>Nbz.or.kLz)j = j-Lz If(j <0) j = j+Lz If(j>Nby.or.jLz)i = i-Lz If(i <0) i = i+Lz If(i>Nbx.or.iic)LDistinct(ic) =0 EndIf ! dd R0 jp =Lst(jp) End Do ! jp /=0 EndDo ! i EndDo ! j EndDo ! k EndDo ! ic nn=0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) REDUCTION(+:nn) Do i=1,Ncenter If(LDistinct(i).eq.1)nn=nn+1 EndDo write(*,*) ' Number of distinct halos =',nn Return End SUBROUTINE HaloDistinct !------------------------------------------------------------------------ ! Get information on a distinct halo iHalo ! Input: Xc,Yc,Zc - center; ! iHalo = number of halo ! Output: updated halo statistics SUBROUTINE SNeib(iHalo) !----------------------------------------------------------------------- Real*8 :: am,xm,ym,zm,Ek,Et,angx,angy,angz,vrr,ee Integer :: omp_get_thread_num Integer*8 :: jp iS = 964771 x = Xc(iHalo) ; y = Yc(iHalo); z = Zc(iHalo) Vx1=VXc(iHalo) ; Vy1= VYc(iHalo); Vz1= VZc(iHalo) vrms2 = (Unbound*Vcirc(iHalo))**2 ! do not take too fast particles R0 = Rmc(iHalo)*Rout ! outer radius Rvir = Rmc(iHalo) ! virial radius dmax = R0**2 iProc = omp_get_thread_num() +1 ! current thread in range 1-Nproc jSeed= iSeed(iProc) ! random seed nPartExpect = Amc(iHalo)/Dscale ! expected number of small ! particles inside virial radius If(nPartExpect 0.999.or.Randd(jseed)10)STOP Return End SUBROUTINE Sbound !----------------------------------------------------------------- ! Inverse cumulative chi-2 distribution with ! 3 degrees of freedom: chi_3(x)=frac. Need it for finding ! velocity, which has given fraction of particles. ! Finds x, ! which has given fraction of particles frac. ! fSigma = v/rms_v = [0-infty] ! frac = [0-1] ! error of approximation is less than ! 7.5 percent for v<4.5sigma ! 15 percent at 6 sigma, ! 30 percent at 10sigma FUNCTION fSigma(frac) real, parameter :: beta =1.,thresh=0.9 real, parameter :: C =log(1.-thresh) If(frac',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) 110 format(1x,' Outer radius in virial radius= ',f6.2,' dLog =',f6.3, & ' Isolation radius (Rvir units)=',f6.2, & /1x,' Unbound factor (Vcirc units)= ',f6.2, & /1x,' Nparticles for direct energy summation=',i7, & ' force smoothing for energy (kpch)=',f6.2) scaleM = Dscale/weightSmall 50 format(a,T19,a, T45,a,T61,a,T71,a,T80,a, & T87,a,T92,a,T103,a,T113,a,T120,a,T128,2a) write(2,50) & ' Npart','Coords(Mpch)','Veloc(km/s)','Mvir', & 'Rvir','Vrms','Vcirc','Mvir_unbnd','Vrms','VirRat', & 'dR/Rvir','AngMom','Lambda' write(40,50) & ' Npart','Coords(Mpch)','Veloc(km/s)','Mvir', & 'Rvir','Vrms','Vcirc','Mvir_unbnd','Vrms','VirRat', & 'dR/Rvir','AngMom','Lambda' write(40,'(T3,a,T12,a,T20,a,T28,a,T40,a,T55,a, & T63,a,T71,a,T90,a,T92,a,T99,a)') & 'Npart','Rout/Rvir','Rbin','M(0)Then rr = Rbin(j,i) Else rr = 10.**((j+0.5)*dlog) EndIf If(rr > 0.01.and.rr<=Rout)Then zbin = 10.**((j+1)*dlog) Vch = sqrt(aMsum/zbin*scaleM/scaleR)*2.08e-3 Vol = 4.188*(zbin**3 -10.**(3.*j*dlog)) write(40,12)nPr(j,i),zbin,Rbin(j,i),aMsum* Dscale/weightSmall, & aMr(j,i)* Dscale/weightSmall/Vol/scaleR**3, & Vch, Vvrms(j,i),Vrad(j,i), & Vrad2(j,i),aJx(j,i),aJy(j,i),aJz(j,i), & sqrt(aJx(j,i)**2 +aJy(j,i)**2 +aJz(j,i)**2) EndIf EndDo EndIf EndDo 10 format(i7,3F9.4,3F8.1,g11.3,f8.2,2f7.1,g11.4,f8.2,f8.3,f8.4,f9.4,3f8.5) 12 format(i8,2F9.4,2g13.5,2x,4f9.2,2x,4f9.2) Return End SUBROUTINE WriteData !----------------------------------------------------------- ! dump bound particles data SUBROUTINE WriteDataBound !----------------------------------------------------------- WRITE (2) & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv WRITE (40,100) & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv write(2) Nmax_s,Ncenter,Rout,dLog,Risolate,Unbound write(40,110) Nmax_s,Ncenter,Rout,dLog,Risolate,Unbound 100 FORMAT(1X, & 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) 110 format(1x,'Max number of selected particles=',i5, & /1x,'Number of halos=',i12,& /1x,' Outer radius in virial radius= ',f6.2,' dLog =',f6.3, & ' Isolation radius (Rvir units)=',f6.2, & /1x,'Unbound factor (Vcirc units)= ',f6.2) Do i=1,Ncenter If(i<1001) & write(40,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)) write(2)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.3,3F8.1,g11.3,f8.2,2f7.1,I6,500i11) Return End SUBROUTINE WriteDataBound !--------------------------------------------------------- ! Read halos catalog ! aMmin_dist -- minimum halo mass ! iFlag = 0 - read and count halos, but do not store them ! = 1 - read and store halos SUBROUTINE ReadHalos(iFlag,aMmin_dist) !-------------------------------------------------------------- Character :: Line*80, Txt*6='Number', Txt2*8,txt1*4 Rdmax =0. 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 read(3,*,err=20)Txt2 read(3,*,err=20)Txt2 Do read(3,*,err=20)Txt2 !write(*,*)TRIM(Txt2) If(TRIM(Txt2)==Txt)Exit EndDo Do read (3,45,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,Npart,Rmax,Conc !,Nb ! write (*,'(i9,3f10.4,f8.2,i10,2f9.2,i6)') Ncenter, X0, Y0 , Z0,vcrc,Npart,Rmax,Conc,Nb 45 Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) If(iStat.ne.0)Exit If(aM>aMmin_dist)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 Rmc(Ncenter) = rr/1000. Vcirc(Ncenter) = vcrc Vrms(Ncenter) = vrmss Rdmax = MAX(rr/1000.,Rdmax) EndIf EndIf ! test mass Enddo 40 Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) If(iFlag==0)write(*,'(a,i8)') ' Read Halos=',Ncenter If(iFlag==1)write(*,'(a,i8,a,f8.3,a)') ' Read Halos=',Ncenter, & ' Max Radius=',Rdmax,'Mpch' 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 !---------------------------------------------------------- Integer*8 :: i,Lsmall,JPAGE,N_particles,iL,IN,N_in_last Np =0 JPAGE = NPAGE 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)/JPAGE +1 N_in_last=N_particles -JPAGE*(Npages-1) EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,IROW,In_page,iL,IN) DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = JPAGE*(IROW-1) if(mod(IROW,1).eq.0)write (*,*)' Read page=',IROW,' Np=',iL ! Loop over particles If(lspecies(Nspecies)>1024**3)Then j = (IROW-1)/64 +1 if(j>8.or.j<1)write(*,*) ' Error in file number: ',j read(20+j,REC=mod(IROW-1,64)+1)xpar,ypar,zpar,vx,vy,vz Else READ (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ EndIf If(iL+In_page > Nmaxpart)Stop ' Too many particles. Increase Nmaxpart' DO IN=1,In_page i = iL +IN Xp(i) =(XPAR(IN)-1.)*Xscale Yp(i) =(YPAR(IN)-1.)*Xscale Zp(i) =(ZPAR(IN)-1.)*Xscale VXp(i) =VX(IN) *Vscale VYp(i) =VY(IN) *Vscale VZp(i) =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) Lagrang(i) = i 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 !-------------------------------------------------- Integer*8 :: jend,N_particles,jstart,k 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(txt) !--------------------------------------------------- ! Open file on a tape Character (LEN=*) :: txt Character*120 :: file1 file1 ='PMFILES/PMcrd.0'//TRIM(txt)//'.DAT' 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,EK,EK1,EK2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras WRITE (*,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & 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 Jfiles =1 If(lspecies(Nspecies)> 1024**3)Jfiles=8 Do i =1,Jfiles write(file1,'(a,i1,a)')'PMFILES/PMcrs',i-1,'.0'//TRIM(txt)//'.DAT' write(*,'(2a,i3)') ' Open file=',file1,i OPEN(UNIT=20+i,FILE=TRIM(file1),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo 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 !------------------------------------------- random number generator FUNCTION RANDd(M) !------------------------------------------------ DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4 & /453815927,2147483648.,2147483647,536870912,131072,256, & 16777216,4,16384,8388608,128/ ML=M/K1*K1 M1=(M-ML)*L1 ML=M/K2*K2 M2=(M-ML)*L2 ML=M/K3*K3 M3=(M-ML)*L3 ML=M/K4*K4 M4=(M-ML)*L4 M5=KI-M IF(M1.GE.M5)M1=M1-KI-1 ML=M+M1 M5=KI-ML IF(M2.GE.M5)M2=M2-KI-1 ML=ML+M2 M5=KI-ML IF(M3.GE.M5)M3=M3-KI-1 ML=ML+M3 M5=KI-ML IF(M4.GE.M5)M4=M4-KI-1 ML=ML+M4 M5=KI-ML IF(LC.GE.M5)ML=ML-KI-1 M=ML+LC RANDd=M/AM RETURN END FUNCTION RANDd end module SetHaloParticles !-------------------------------------------------------------------------------------- ! ! PROGRAM HaloParticles use SetHaloParticles Character*120 :: catname,listSh,listCat,txt*8,haloList,halo_ascii integer :: omp_get_max_threads Logical :: FindDistinct,FindBound FindDistinct = .true. FindBound = .true. ! ------------------------------------------ Read data write(*,'(a,$)')' Enter snapshot=' read(*,'(a)')txt catname='CATALOGS/CatshortA.'//TRIM(txt)//'.DAT' listSh ='CatshortDistinctA.'//TRIM(txt)//'.DAT' listCat ='CatalogDistinctA.'//TRIM(txt)//'.DAT' haloList ='HaloListA.'//TRIM(txt)//'.DAT' halo_ascii='HaloListA_ascii.'//TRIM(txt)//'.DAT' tt = seconds() CALL RDTAPE(txt) ! read control info, open particles file CALL SetWeights ! set iWeight array and weightSmall 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 xmins =0.; xmaxs =250. ymins =0.; ymaxs =250. zmins =0.; zmaxs =250. CALL ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! read particles: Np If(.not.FindDistinct)goto 200 Open( 2,file=listSh,Status='UNKNOWN')! short list of halos Open( 40,file=listCat,Status='UNKNOWN')! short list of halos Open( 3,file=catname,Status='UNKNOWN')! short list of halos ! Find properties of Distinct/Isolated halos CALL ReadHalos(0,aMmin_dist) ! count halos: Ncenter Cell = Box/Nlinker Nmx = INT(xmins/Cell) ! boundaries for linker list Nbx = INT(xmaxs/Cell) Nmy = INT(ymins/Cell) ; Nby = INT(ymaxs/Cell) Nmz = INT(zmins/Cell) ; Nbz = INT(zmaxs/Cell) write (*,'(a,3(3x,2i5))') ' Limits: ',Nmx,Nbx,Nmy,Nby,Nmz,Nbz ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter)) write (*,*) ' Allocated Centers' Ncenter = 0 ! initialize halo counter CALL ReadHalos(1,aMmin_dist) ! read halos Close(3) Rdmax = Rdmax*Rout ! set maximum radius write (*,*) ' Read Halos: ',Ncenter ALLOCATE (Lst(Ncenter),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) write (*,*) ' Allocated List' Call ListHalos ! make linker list for Halos write (*,*) ' ListHalos made: ',Ncenter ALLOCATE (LDistinct(Ncenter)) Call HaloDistinct ! make list of distinct halos write (*,*) ' Found Distinct halos' DEALLOCATE(Lst,Label) Ncp = Np CALL Points(Rdmax) ! Make more points periodically write (*,*) 'Points made' write(*,*) '******* time=',seconds() ! set limits for linker list Celloptimal = Rdmax/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)) write(*,*) '******* time=',seconds() Call List ! make linker List, Label write(*,*) '******* time=',seconds() Nbin_max = INT(log10(Rout)/dlog)+1 Nbin_min = Nbin_max - Nrad-1 ALLOCATE (aMr(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Vvrms(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (nPr(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Rbin(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Vrad(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Vrad2(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (aJx(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (aJy(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (aJz(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Amh(Ncenter),dRh(Ncenter),Ekin(Ncenter),Epot(Ncenter)) ALLOCATE (Anhx(Ncenter),Anhy(Ncenter),Anhz(Ncenter),Vhrad(Ncenter)) !$OMP PARALLEL DO DEFAULT(SHARED) & ! initiate in parallel !$OMP PRIVATE (ic,j) Do ic =1,Ncenter Do j=Nbin_min, Nbin_max aMr(j,ic) = 0. ; Vvrms(j,ic) = 0. ; nPr(j,ic) = 0. Rbin(j,ic) = 0. ; Vrad(j,ic) = 0. ; Vrad2(j,ic) = 0. aJx(j,ic) = 0. ; aJy(j,ic) = 0. ; aJz(j,ic) = 0. EndDo Amh(ic) = 0. ; dRh(ic) = 0. ; Ekin(ic) = 0. Epot(ic) = 0. ; Vhrad(ic) = 0. Anhx(ic) = 0. ; Anhy(ic) = 0. ; Anhz(ic) = 0. EndDo Nproc = omp_get_max_threads() ! Number of cores Nskip = 1e6 ! skip this number of seeds write(*,*) ' Number of cores =',Nproc ALLOCATE(iSeed(Nproc)) ! Set seeds for selection of particles jseed = 9121071 Do i=1,Nproc iSeed(i) = jseed Do j=1,Nskip xx = Randd(jseed) EndDo EndDo write (*,20) Np,Ncenter,Box,Cell write (2,20) Np,Ncenter,Box,Cell 20 Format(5x,'Statistics for ',i10,' points. Maxima=',i9, & 3x,'Box(Mpc)=',f6.1,' Cell for linker list(Mpch)=',f9.5) write(*,*) '******* time=',seconds() CALL SPairs ! Accumulate statistics write(*,*) '******* time=',seconds() CALL WriteData write(*,*) '******* time(sec)=',seconds() DEALLOCATE (aMr,Vvrms,nPr,Rbin,Vrad,Vrad2,aJx,aJy,aJz,Amh,dRh) DEALLOCATE (Ekin,Epot,Anhx,Anhy,Anhz,Vhrad) DEALLOCATE (Lst,Label) DEALLOCATE (Xc,Yc,Zc,VXc,VYc,VZc,Amc,Rmc,Vrms,Vcirc) close (2) close (40) 200 If(.not.FindBound)Stop ! Find Bound particles of all halos ------------------------------ write (*,*) ' Find most bound particles ' Open( 3,file=catname,Status='UNKNOWN')! short list of halos Open( 2,file=haloList,Status='UNKNOWN',form='unformatted')! short list of halos Open( 40,file=halo_ascii,Status='UNKNOWN')! short list of halos Ncenter =0 CALL ReadHalos(0,0.) ! count halos: Ncenter Cell = Box/Nlinker Nmx = INT(xmins/Cell) ! boundaries for linker list Nbx = INT(xmaxs/Cell) Nmy = INT(ymins/Cell) ; Nby = INT(ymaxs/Cell) Nmz = INT(zmins/Cell) ; Nbz = INT(zmaxs/Cell) write (*,'(a,3(3x,2i5))') ' Limits: ',Nmx,Nbx,Nmy,Nby,Nmz,Nbz ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter)) write (*,*) ' Allocated Centers=',Ncenter Ncenter = 0 ! initialize halo counter CALL ReadHalos(1,0.) ! read halos Rdmax = Cell !min(Rdmax/5.,Cell) ! Take Cell size or 1/5 of maximum virial radius Ncp = Np CALL Points(Rdmax) ! Make more points periodically write(*,*) ' *** time after Points =',seconds() ALLOCATE (Lst(Ncp),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) ALLOCATE (nBound(Nmax_s,Ncenter),nSelect(Ncenter)) write(*,*) ' *** time=',seconds() Call List ! make linker List, Label write(*,*) ' *** time after List=',seconds() CALL SHalos ! Accumulate statistics write(*,*) ' *** time after SHalos=',seconds() CALL WriteDataBound write(*,*) ' *** time final (sec)=',seconds() End PROGRAM HaloParticles