!---------------------------------------------------------------------------- ! ! BDM halo finder A.Klypin 2010 ! ! Contains: Module Structures Integer(kind=4), PARAMETER :: & Nrad = 700 ! Number of shells for halo mass profile Real(kind=4) :: & dR, & ! size of radius binning for profiles MaxMemory=800, & ! limit on memory for the run MassMin , & ! minimum halo mass alpha , & ! environment correction factor Rext , & ! scale for resolution correction of vir radius Environ , & ! factor for radius of environment correction OvdensMin ! factor for reducing max density below Overdensity Integer(kind=4) :: & iVirial ! switch: 1 =Virial, 0=200 overdensity crit, 2=200rho_matter Real(kind=4) :: TotalMemory=0.2, t0 ! current memory Real(kind=4) :: Xleft,Xright,Yleft,Yright,Zleft,Zright,dBuffer ! boundaries of domain Character*120 :: CatshortName,CatalogName, & ! names of halo catalogs outputName ! dump file Real(kind=4) :: Om0,Ovdens ! cosmology Real(kind=4) :: MassOne ! current simulation Integer(kind=8) :: Np ! Number of particles in domain real(kind=4), Allocatable :: Vf(:,:,:,:) ! Velocity field Integer(kind=4) :: mesh ! ------------------------- Maxima ---------------------- Integer(kind=8) :: Nmaxima Real(kind=4) :: RadMax Real(kind=4), ALLOCATABLE, DIMENSION(:) :: & ! maxima of density Mvir,Rvir,Mtotal,VmaxM,RmaxM, & xMaxx,yMaxx,zMaxx, & VxMaxx,VyMaxx,VzMaxx, & Xoff, & Dist, MassNext Integer(kind=8), ALLOCATABLE, DIMENSION(:) :: & LstMax, & ! linker-list of maxima MaxIndex, & ! =0 for distinct, =-2 for sub IndexDist, IndexLoc Real(kind=4), ALLOCATABLE, DIMENSION(:,:) :: & MassProf, & ! mass profile DensMax,DistSub ! ------------------------- Halos ---------------------- Integer(kind=4) :: Nhalo Real(kind=4), ALLOCATABLE, DIMENSION(:,:) :: & RadH1,MassH1,VrmsH1,VradH1,VrmsrH1, & RadH2,MassH2,VrmsH2,VradH2,VrmsrH2 Integer(kind=4), ALLOCATABLE, DIMENSION(:,:) :: & NbinH1,NbinH2 !------------------------- Main linker-list -------------- Integer(kind=8), ALLOCATABLE :: Lst(:) ! Linker list Integer(kind=8), ALLOCATABLE :: Label(:,:,:) ! Head of zero-level LL Integer(kind=4) :: Nmx,Nmy,Nmz,Nbx,Nby,Nbz ! limits for linker-list Real(kind=4) :: Cell,Roptimal ! Size in Mpch of a linker-list cell end Module Structures !---------------------------------------------------------------------------- ! Module LinkerList use Structures use Tools Contains !---------------------------------------------------------- ! ! SUBROUTINE BDM(mDENSIT) !---------------------------------------------------------- Use Density Integer(kind=8) :: idummy,ip,NpPM Integer(kind=4) :: jdummy integer(kind=4) :: OMP_GET_MAX_THREADS,OMP_GET_THREAD_NUM character*80 :: Path logical :: op NpPM = Nparticles Np = Nparticles t0 = seconds() Call TimingMain(5,-1) Call ReadParameters(ISTEP) iThreads = OMP_GET_MAX_THREADS() write (*,'(a,i4)') ' Number of threads =',iThreads write (*,'(a,i4)') ' mDENSIT =',mDENSIT write (*,'(a,i4)') ' iVirial =',iVirial Path ='' tstart = seconds() If(mDENSIT==1)Call DENSIT ! density on original Ng mesh tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for Density =',tfinish-tstart,tfinish-t0 Call FindMaxima tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for Maxima =',tfinish-tstart,tfinish-t0 myMemory =Memory(-1_8*NGRID*NGRID*NGRID) DeAllocate (FI) write(*,'(a,i11)') ' Go to RescaleCoords : ',Nparticles Call SetParameters Call RescaleCoords(1) tfinish = seconds() write (*,'(a,i4)') ' iVirial =',iVirial write(*,'(a,i11)') ' Go to AddBuffer : ',Nparticles write(13,'(10x,a,T50,2f10.2)') ' time for Dens+Maxima =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for Dens+Maxima =',tfinish-tstart,tfinish-t0 tstart = seconds() Call AddBuffer tfinish = seconds() write(13,'(10x,a,T50,2f10.2)') ' time for AddBuffer =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for AddBuffer =',tfinish-tstart,tfinish-t0 Call SizeList ALLOCATE (Lst(Np),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) myMemory= Memory(1_8*(Np+(Nbx-Nmx+1_8)*(Nby-Nmy+1_8)*(Nbz-Nmz+1_8))) tstart = seconds() Call List tfinish = seconds() write(13,'(10x,a,T50,2f10.2)') ' time for List =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for List =',tfinish-tstart,tfinish-t0 write(*,*) ' Go to FindDistinctCandidates: ',Nparticles Call FindDistinctCandidates ! get estimates of Mvir and Rvir for all maxima write(*,'(a,i11)') ' Go to ParametersDistinct : ',Nparticles tstart = seconds() Call ParametersDistinct tfinish = seconds() write(13,'(10x,a,T50,2f10.2)') ' time for ParametersDistinct =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for ParametersDistinct =',tfinish-tstart,tfinish-t0 ! Call FindSubs ! Call ParametersSubs write(*,*) ' goto RemoveCloseMaxima' tfinish = seconds() DEALLOCATE (Lst,Label) myMemory= Memory(-1_8*(Np+(Nbx-Nmx+1_8)*(Nby-Nmy+1_8)*(Nbz-Nmz+1_8))) Call SizeListMaxima ALLOCATE (Lst(Nmaxima),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) myMemory= Memory(1_8*(Nmaxima+(Nbx-Nmx+1_8)*(Nby-Nmy+1_8)*(Nbz-Nmz+1_8))) Call ListMaxima tstart = seconds() Call RemoveDuplicates tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for RemoveDuplicates =',tfinish-tstart,tfinish-t0 tstart = seconds() if(abs(alpha)>1.e-3)Then Call VelocityField Call CorrectClose end if tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for CorrectClose =',tfinish-tstart,tfinish-t0 write(*,*) ' Go to Write Catshort: ',Nparticles tstart = seconds() Call WriteFiles tfinish = seconds() write(13,'(10x,a,T50,2f10.2)') ' time for WriteFiles =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for WriteFiles =',tfinish-tstart,tfinish-t0 close(12) !-------- deallocate temporary arrays DEALLOCATE (Lst,Label) myMemory= Memory(-1_8*(Nmaxima+(Nbx-Nmx+1_8)*(Nby-Nmy+1_8)*(Nbz-Nmz+1_8))) DEALLOCATE (Mtotal,Mvir,Rvir,Xoff) myMemory =Memory(-4_8*Nmaxima ) DEALLOCATE (xMaxx,yMaxx,zMaxx) myMemory =Memory(-3_8*Nmaxima ) DEALLOCATE (VxMaxx,VyMaxx,VzMaxx) myMemory =Memory(-3_8*Nmaxima ) DEALLOCATE(LstMax) myMemory =Memory(-1_8*Nmaxima ) DEALLOCATE(VmaxM) myMemory =Memory(-1_8*Nmaxima ) !-------- restore PM structure Call RemoveBuffer(NpPM) myMemory =Memory(1_8*NGRID*NGRID*NGRID) Allocate (FI(NGRID,NGRID,NGRID)) Call TimingMain(5,1) end SUBROUTINE BDM !---------------------------------------------------------- ! Read/create configuration file BDM.config ! SUBROUTINE ReadParameters(jStep) !---------------------------------------------------------- Character :: txt*80,str*10,eqsign*1,Line*120,Catlabel*10 logical :: FileExists ! default set of parameters iVirial = 1 dR = 0.1 dBuffer = 5. alpha = 0.20 Rext = 0.15 MassMin = 1.e12 Environ = 1.5 OvdensMin= 1.5 Inquire(file='BDM.config', Exist = FileExists) If(FileExists)Then open(1,file='BDM.config') rewind(1) Read(1,'(a)')txt write(*,'(a)')txt Read(1,'(a)')txt write(*,'(a)')txt Do Read(1,'(a)',iostat=io) Line If(io /= 0)exit !write(*,'(a)')Line i0 = INDEX(Line,'=') i1 = INDEX(Line,'!') !write(*,*) ' position of = and "!" in current line =',i0,i1 If(i1.gt.i0.and.i0>0)Then backspace(1) Read(1,*)str !write(*,'(a,3x,a)')str backspace(1) If(TRIM(str)=='iVirial'.or.TRIM(str)=='IVIRIAL'.or.TRIM(str)=='ivirial')Then Read(1,*)str,eqsign,iVirial Else If(TRIM(str)=='dR'.or.TRIM(str)=='DR'.or.TRIM(str)=='dr')Then Read(1,*)str,eqsign,dR Else If(TRIM(str)=='Rext'.or.TRIM(str)=='Rextr'.or.TRIM(str)=='RextR')Then Read(1,*)str,eqsign,Rext Else If(TRIM(str)=='alpha'.or.TRIM(str)=='Alpha'.or.TRIM(str)=='ALPHA')Then Read(1,*)str,eqsign,alpha Else If(TRIM(str)=='MassMin'.or.TRIM(str)=='MinMass'.or.TRIM(str)=='massmin')Then Read(1,*)str,eqsign,MassMin Else If(TRIM(str)=='Environ'.or.TRIM(str)=='Envir'.or.TRIM(str)=='Env')Then Read(1,*)str,eqsign,Environ Else If(TRIM(str)=='OvdensMin'.or.TRIM(str)=='OvdMin'.or.TRIM(str)=='Ovd')Then Read(1,*)str,eqsign,OvdensMin Else !write(*,'(3a)') 'Unrecognized parameter: ',str,' Check spelling. I igonore it' Read(1,*)str End If end If EndDo Else open(1,file='BDM.config') ! ------- create default config file write(1,'(a)')'! ------------- BDM configuration file. Spaces between entries are significant' write(1,'(a)')'! These are main parameters:' write(1,10)'iVirial', iVirial, '! switch: 1 =Virial, 0=200 overdensity' write(1,20)'Rext', Rext, '! scale for resolution correction of vir radius' write(1,20)'MassMin', MassMin, '! Minimum halo mass' write(1,20)'dR', dR, '! size of radius binning for profiles' write(1,20)'OvdensMin', OvdensMin, '! factor for reducing max density below Overdensity' write(1,20)'Environ', Environ, '! factor for radius of environment correction' write(1,20)'alpha', alpha, '! environment correction factor' 10 format(10x,a,T20,' = ',i6,T40,a) 20 format(10x,a,T20,' = ',1p,g10.3,T40,a) EndIf close(1) If(MassMin <= 0.)MassMin = 1.e12 If(OvdensMin <= 0.)OvdensMin = 1. If(alpha <= -1.)alpha = 0. write(*,'(/a)') ' ------ current set of parameters:' write(*,'(a)')'! These are main parameters:' write(*,10)'iVirial', iVirial, '! switch: 1 =Virial, 0=200 overdensity' write(*,20)'Rext', Rext, '! scale for resolution correction of vir radius' write(*,20)'MassMin', MassMin, '! Minimum halo mass' write(*,20)'dR', dR, '! size of radius binning for profiles' write(*,20)'OvdensMin', OvdensMin, '! factor for reducing max density below Overdensity' write(*,20)'Environ', Environ, '! factor for radius of environment correction' write(*,20)'alpha', alpha, '! environment correction factor' moment = jstep write(outputName,'(a,i4.4,a)')'CATALOGS/outputB.',jStep,'.dat' !open(13,file=TRIM(outputName),buffered='NO') open(13,file=TRIM(outputName)) !---------------------------- Open files --------------------- SELECT CASE (iVirial) CASE (0) CatLabel = 'W.' ! 200\rho_critical CASE (1) CatLabel = 'V.' ! virial overdensity CASE (2) CatLabel = 'M.' ! 200 matter overdensity CASE (3) CatLabel = 'A.' ! Abacus overdensity: corrected virial end SELECT write(outputName,'(2a,2(i4.4,a))')'CATALOGS/Catshort',TRIM(CatLabel),jStep,'.',Nrealization,'.DAT' open(12,file=TRIM(outputName)) ! write(outputName,'(2a,i4.4,a)')'CATALOGS/Catalog',TRIM(CatLabel),jStep,'.DAT' ! open(20,file=TRIM(outputName),form='unformatted',status ='UNKNOWN') end SUBROUTINE ReadParameters !-------------------------------------------------------------- ! virial overdensity for cosmological model ! at different expansion parameter AEXPN Function OverdenVir() !-------------------------------------------------------------- xx =-(1.-Om0)*AEXPN**3/(Om0+(1.-Om0)*AEXPN**3) OverdenVir =(178.+82.*xx-39.*xx**2)/(1.+xx) write (*,*) ' Overdensity Delta =',OverdenVir,Om0,AEXPN END Function OverdenVir !-------------------------------------------------------------- ! Abacus overdensity for cosmological model ! at different expansion parameter AEXPN Function OverdenAbacus() !-------------------------------------------------------------- xx =-(1.-Om0)*AEXPN**3/(Om0+(1.-Om0)*AEXPN**3) OverdenAbacus =(178.+82.*xx-39.*xx**2)/(1.+xx)*(200./178.) !write (*,*) ' Overdensity Delta =',OverdenVir END Function OverdenAbacus !---------------------------------------------------------- ! ! SUBROUTINE RescaleCoords(iFlag) !---------------------------------------------------------- Integer(kind=8) :: ip Xscale = Box/NGRID ! Scale for comoving coordinates Vscale = 100.*Xscale/AEXPN ! Scale for velocities !Dscale = 2.774e+11*(Box/NROW)**3 ! mass scale !MassOne= Om0*Dscale ! mass of the smallest particle write(*,*) ' Inside rescale coords:', Xscale,Vscale If(iFlag==1)Then ! scale to Mpc and Msun !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) Do ip =1,Np Xpar(ip) = (Xpar(ip)-1.)*Xscale Ypar(ip) = (Ypar(ip)-1.)*Xscale Zpar(ip) = (Zpar(ip)-1.)*Xscale VX(ip) = VX(ip)*Vscale VY(ip) = VY(ip)*Vscale VZ(ip) = VZ(ip)*Vscale end Do else !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) Do ip =1,Np Xpar(ip) = Xpar(ip)/Xscale+1. Ypar(ip) = Ypar(ip)/Xscale+1. Zpar(ip) = Zpar(ip)/Xscale+1. If(Xpar(ip).ge.NGRID)Xpar(ip) =Xpar(ip)-NGRID If(Xpar(ip).lt.1) Xpar(ip) =Xpar(ip)+NGRID If(Ypar(ip).ge.NGRID)Ypar(ip) =Ypar(ip)-NGRID If(Ypar(ip).lt.1) Ypar(ip) =Ypar(ip)+NGRID If(Zpar(ip).ge.NGRID)Zpar(ip) =Zpar(ip)-NGRID If(Zpar(ip).lt.1) Zpar(ip) =Zpar(ip)+NGRID VX(ip) = VX(ip)/Vscale VY(ip) = VY(ip)/Vscale VZ(ip) = VZ(ip)/Vscale end Do end If end SUBROUTINE RescaleCoords !--------------------------------------------------------------------------- ! ! SUBROUTINE WriteFiles integer(kind=8) :: ic,ip,Ncount iHalo = 0 MassMin = max(MassMin,20.*MassOne) Ncount =0 Ncount1 = 0 Do ip=1,Nmaxima x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) If(x>=Xleft.and.x=Yleft.and.y=Zleft.and.zMassmin=' ,Ncount, & ' Nwritten= ',iHalo,' Mass_min=',Massmin end SUBROUTINE WriteFiles !--------------------------------------------------------------------------- ! ! SUBROUTINE WriteProfiles integer(kind=8) :: ic,ip iHalo = 0 Do ip=1,Nmaxima If(Mvir(ip)>10*MassOne)Then iHalo = iHalo +1 x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) If(Mvir(ip)>100.*MassOne)Then !Vrms = sqrt(EkinM(ip)/Mvir(ip)*2.) Cvir = Concentration(Mvir(ip),1.e3*Rvir(ip),VmaxM(ip)) If(Cvir < 0.)Cvir = Rvir(ip)/RmaxM(ip)*2.15 iHalo= ih iStart = 0 Do i=-Nrad+1,0 If(NbinH1(i,ih)>0.and.MassH1(i,ih)>5.*MassOne)Then iStart = i ; exit EndIf EndDo Nlines = 0 ! total lines of profile Do i=-Nrad+1,0 If(NbinH1(i,ih)>0.and.MassH1(i,ih)>5.*MassOne)Then Nlines = Nlines +1 EndIf EndDo write(20) & x,y,z,VxMaxx(ip),VyMaxx(ip),VzMaxx(ip), & Mvir(ip),Mtotal(ip),1.e3*Rvir(ip), VmaxM(ip), & iHalo,Cvir,Mvir(ip)/MassOne,MaxIndex(ip),Xoff(ip),Nlines Radius = 2.*Rvir(ip) Do i=iStart,0 R = Radius*10.**(i*dR) Rin = Radius*10.**((i-1)*dR) Vcirc1 = 6.582e-5*sqrt(MassH1(i,ih)/R)/sqrt(AEXPN) Vcirc2 = 6.582e-5*sqrt(MassH2(i,ih)/R)/sqrt(AEXPN) Volume = 4.1888*(R**3-Rin**3) DensH1 = (MassH1(i,ih)-MassH1(i-1,ih))/Volume*1.e-9 !density Msunh/kpch**3 comoving DensH2 = (MassH2(i,ih)-MassH2(i-1,ih))/Volume*1.e-9 If(NbinH1(i,ih)/= 0.and.MassH1(i,ih)>5.*MassOne)Then write(20) R*1.e3, & NbinH1(i,ih),RadH1(i,ih),MassH1(i,ih),Vcirc1, & DensH1,VrmsH1(i,ih),VradH1(i,ih),VrmsrH1(i,ih), & NbinH2(i,ih),RadH2(i,ih),MassH2(i,ih),Vcirc2, & DensH2,VrmsH2(i,ih),VradH2(i,ih),VrmsrH2(i,ih) End If EndDo End If end If end do close (20) end SUBROUTINE WriteProfiles !--------------------------------------------------------------------------- ! Find profile of each halo and subhalos ! SUBROUTINE GetProfiles integer(kind=8) :: ic,ip,i,Ncount1,Ncount2 Nhalo = Nmaxima write(*,*) ' GetProfiles. Nhalo=',Nhalo write(13,*) ' GetProfiles. Nhalo=',Nhalo iHalo = 0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Nhalo RadH1(:,i) =0. ; MassH1(:,i) =0. ; VrmsH1(:,i) =0. VradH1(:,i) =0. ; VrmsrH1(:,i) =0. ; NbinH1(:,i) =0 RadH2(:,i) =0. ; MassH2(:,i) =0. ; VrmsH2(:,i) =0. VradH2(:,i) =0. ; VrmsrH2(:,i) =0. ; NbinH2(:,i) =0 EndDo Ncount1 = 0. Ncount2 = 0. Do ip=1,Nmaxima If(Mvir(ip)>100.*MassOne)Ncount1 = Ncount1 +1 If(Mtotal(ip)>100.*MassOne)Ncount2 = Ncount2 +1 End Do write(*,*) ' Nmvir >100partilces = ',Ncount1 write(*,*) ' Ntotal>100partilces = ',Ncount2 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ih,ip,ic,x,y,z,aR) Do ip=1,Nmaxima If(Mvir(ip)>100.*MassOne)Then x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) aR = Rvir(ip) if(mod(ih,1000)==0) & write(13,'(i8,5f9.4)') ip,x,y,z,aR Call HaloProfile(x,y,z,aR,ip) EndIf EndDo ! i end SUBROUTINE GetProfiles !--------------------------------------------------------------------------- ! Get profile of a halo SUBROUTINE HaloProfile(x,y,z,aR,ip) !--------------------------------------------------------------------------- Real(kind=4), PARAMETER :: fiScale = 4.333e-9 Real(kind=4) :: Fi(-Nrad:0) Real*8 :: wx,wy,wz integer(kind=8) :: ic,ip,jp Radius = 2.*aR d0 = Radius**2 ! get final statistics of particles factorZ = 100.*sqrt(Om0/AEXPN**3+(1.-Om0)) *AEXPN wx = VxMaxx(ip) ; wy = VyMaxx(ip) ; wz = VzMaxx(ip) Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) ! Get mass profile Do k3 =k1, k2 Do j3 =j1, j2 Do i3 =i1, i2 jp =Label(i3,j3,k3) Do while (jp.ne.0) dd =(x-Xpar(jp))**2+(y-Ypar(jp))**2+(z-Zpar(jp))**2 If(dd< d0) Then r = sqrt(max(dd,1.e-20)) ii = max(min(INT(log10(r/Radius)/dR),0),-Nrad) dx = Xpar(jp) -x dy = Ypar(jp) -y dz = Zpar(jp) -z dvx = VX(jp) - wx +factorZ*dx ! true velocity dvy = VY(jp) - wy +factorZ*dy dvz = VZ(jp) - wz +factorZ*dz vv = dvx**2 + dvy**2 + dvz**2 ! kinetic energy vr = (dvx*dx+dvy*dy+dvz*dz)/r ! radial velocity MassH1(ii,ih) = MassH1(ii,ih) + MassOne RadH1(ii,ih) = RadH1(ii,ih) + r/aR ! radius in virial units VrmsH1(ii,ih) = VrmsH1(ii,ih) + vv VradH1(ii,ih) = VradH1(ii,ih) + vr VrmsrH1(ii,ih) = VrmsrH1(ii,ih) + vr**2 NbinH1(ii,ih) = NbinH1(ii,ih) + 1 EndIf ! dd0.)Ncount = Ncount +1 end Do ! --------------------------- !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (ip,x,y,z,m,jp,dd,k1,k2,k3,j1,j2,j3,i1,i2,i3) Do ip=1,Nmaxima If(Mtotal(ip)>MassMin)Then x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) m = Mtotal(ip) Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) !if(mod(ip,1)==1000)write(*,'(a,12i7)') ' ip=',ip,i1,i2,j1,j2,k1,k2 Do k3 =k1, k2 Do j3 =j1, j2 Do i3 =i1, i2 jp =Label(i3,j3,k3) Do while (jp.ne.0) If(jp/=ip)Then If(Mtotal(jp)>MassMin)Then dd =(x-Xmaxx(jp))**2+(y-Ymaxx(jp))**2+(z-Zmaxx(jp))**2 If((dd.lt.Rvir(jp)**2.or.dd.lt.Rvir(ip)**2))Then If(m.lt.Mtotal(jp))Then Mvir(ip) = 0. Rvir(ip) = 0. Mtotal(ip) = 0. Else Mvir(jp) = 0. Rvir(jp) = 0. Mtotal(jp) = 0. end If End If end If end If jp =Lst(jp) end do end do end do end Do end If End Do ! ip tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for RemoveDuplicates =',tfinish-tstart,tfinish-t0 Do ip=1,Nmaxima If(Mtotal(ip)>MassMin)Ncount1 = Ncount1 +1 end Do write(*,'(2(a,i10))') ' Nhalos with Mvir>MassMin before= ',Ncount,& ' Nhalos after removing dublicates= ',Ncount1 Ncount1 = 0 !! --------------------------- !!!$OMP PARALLEL DO DEFAULT(SHARED) & !!!$OMP PRIVATE (ip,x,y,z,m,jp,dd,k1,k2,k3,j1,j2,j3,i1,i2,i3) ! Do ip=1,Nmaxima ! If(Mtotal(ip)>MassMin)Then ! x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) ! m = Mtotal(ip) ! Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) ! !if(mod(ip,1)==1000)write(*,'(a,12i7)') ' ip=',ip,i1,i2,j1,j2,k1,k2 ! Do k3 =k1, k2 ! Do j3 =j1, j2 ! Do i3 =i1, i2 ! jp =Label(i3,j3,k3) ! Do while (jp.ne.0) ! If(jp/=ip)Then ! dd =(x-Xmaxx(jp))**2+(y-Ymaxx(jp))**2+(z-Zmaxx(jp))**2 ! If((dd.lt.Rvir(jp)**2.or.dd.lt.Rvir(ip)**2).and.Mvir(jp)>MassMin)Then ! !print '(2i9,3f9.2,4es12.4)',ip,jp,sqrt(dd),Rvir(jp),Rvir(ip),Mvir(jp),m ! Mvir(ip) = 0. ! Mtotal(ip) = 0. ! End If ! end If ! jp =Lst(jp) ! end do ! end do ! end do ! end Do ! end If ! End Do ! ip ! tfinish = seconds() ! write(*,'(10x,a,T50,2f10.2)') ' time for RemoveDuplicates =',tfinish-tstart,tfinish-t0 ! Do ip=1,Nmaxima ! If(Mtotal(ip)>MassMin)Ncount1 = Ncount1 +1 ! end Do ! write(*,'(2(a,i10))') ' Nhalos with Mvir>MassMin before= ',Ncount,& ! ' Nhalos after removing dublicates= ',Ncount1 end SUBROUTINE RemoveDuplicates !--------------------------------------------------------------------------- ! ! SUBROUTINE RemoveDuplicatesSimple !--------------------------------------------------------------------------- integer(kind=8) :: ic,ip,i real(kind=4) :: m tstart = seconds() ! --------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ip,x,y,z,m,ic,D2) Do ip=1,Nmaxima If(Mvir(ip)>MassOne)Then x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) m = Mvir(ip) do ic =1,Nmaxima If(ic/=ip)Then D2 = (xMaxx(ic) -x)**2 +(yMaxx(ic) -y)**2 +(zMaxx(ic) -z)**2 If(D2.lt.Rvir(ic)**2.and.m.lt.Mvir(ic))Then Mvir(ip) = 0. End If end If end do end If End Do ! ip tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for RemoveDuplicates =',tfinish-tstart,tfinish-t0 end SUBROUTINE RemoveDuplicatesSimple !--------------------------------------------------------------------------- ! Find parameters of distinct halos ! SUBROUTINE ParametersDistinct !--------------------------------------------------------------------------- integer(kind=8) :: ic,ip,i,Ncount1,Ncount2,Ncount3,Ncount4 tstart = seconds() print *,' ParametersDistinct: Nmaxima= ',Nmaxima ! --------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ip,x,y,z,xv,yv,zv) ! SCHEDULE(DYNAMIC,10000) Do ip=1,Nmaxima x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) xv = VxMaxx(ip); yv = VyMaxx(ip); zv = VzMaxx(ip) Call GetHalo(x,y,z,xv,yv,zv,ip) End Do ! ip tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for ParametersDistinct =',tfinish-tstart,tfinish-t0 Ncount1 = 0 Ncount2 = 0 Ncount3 = 0 Ncount4 = 0 Do ip=1,Nmaxima If(Mvir(ip)>100.*MassOne)Ncount1 = Ncount1 +1 If(Mtotal(ip)>100.*MassOne)Ncount2 = Ncount2 +1 If(Mvir(ip)>MassOne)Ncount3 = Ncount3 +1 If(Mtotal(ip)>MassOne)Ncount4 = Ncount4 +1 End Do write(*,'(2(a,i10))') ' Nmvir >100partilces = ',Ncount1,' Nmvir >1partilce = ',Ncount3 write(*,'(2(a,i10))') ' Ntotal>100partilces = ',Ncount2,' Ntot >1partilce = ',Ncount4 end SUBROUTINE ParametersDistinct !--------------------------------------------------------------------------- ! Get parameters of distinct halos SUBROUTINE GetHalo(x,y,z,xv,yv,zv,ip) !--------------------------------------------------------------------------- integer(kind=8) :: ic,ip,jp,Ncount,Icount Real(kind=4) :: MassP(0:Nrad) Real*8 :: wx,wy,wz,wL,v2, & xcm,ycm,zcm, & ax,ay,az, & x2,y2,z2,xy,xz,yz,r2,Mtot Radius = 2.*Cell Cell2 = Cell**2 factorZ = 100.*sqrt(Om0/AEXPN**3+(1.-Om0)) *AEXPN fact = 1.150e12*Om0 dCell = Box/NGRID dRvmax = 0.1*(Box/NGRID) ! correct Vmax_radius for resolution R0 = Box/NGRID ! one grid size !write(*,'(a,f9.2,a,f9.2,a,f9.2)') ' Overdense = ',Ovdens,' Radius= ',Radius,' Cell= ',Cell !write(*,'(a,7g12.4,i8)') '--- halo=',x,y,z,xv,yv,zv,Mvir(ip),ip 10 MassP = 0. aMa = fact*Ovdens*Radius**3 ! virial mass for outer radial bin d0 = Radius**2 ! get final statistics of bound particles Ncount = 0 Mtot = 0. aM = 0. aR = 0. wx = xv ; wy = yv ; wz = zv Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) ! Get mass profile Ncount =0 Do k3 =k1, k2 Do j3 =j1, j2 Do i3 =i1, i2 jp =Label(i3,j3,k3) Do while (jp.ne.0) dd =(x-Xpar(jp))**2+(y-Ypar(jp))**2+(z-Zpar(jp))**2 If(dd< d0) Then r = sqrt(dd) !ii = max(min(INT(log10(r/Radius)/dLogR),0),-Nrad) ii = min(max(INT(r/dCell/dR),0),Nrad) MassP(ii) = MassP(ii) + MassOne Ncount = Ncount +1 Mtot = Mtot + MassOne EndIf ! dd10.)Stop '---- Too large outer radius in GetHalo' goto 10 End If Do i = Nrad,0,-1 !---- find virial Radius and Mass Rout = dCell*((i+1)*dR) ! outer radius of this bin Rin = dCell*((i)*dR) ! inner radius of this bin dRadius = Rout-Rin If(MassP(i-1).lt.10.*MassOne)exit DeltaIn = MassP(i-1)/(fact*Rin**3) DeltaOut = MassP(i)/(fact*Rout**3) If(DeltaIn > Ovdens.or.i*dR<1.)Then aR = Rout -dRadius*(Ovdens-DeltaOut)/(DeltaIn-DeltaOut) aR = min(max(aR,Rin),Rout) aM = fact*Ovdens*(aR)**3 exit EndIf EndDo Mvir(ip) = aM Rvir(ip) = aR If(aM<10*MassOne)Then Mtotal(ip) =0. Mvir(ip) =0. Rvir(ip) =0. VmaxM(ip) =0. return End If !--- increase radius to correct for resolution !Radius = aR + R0*min(AEXPN*Rext/(aR/R0)**SlopeR/exp(aR/R0/2.5),3.0) Radius = aR + R0*min(sqrt(AEXPN)*Rext/exp(aR/R0/2.5),3.0) !--- re-do binning using virial radius MassP = 0. aMa = fact*Ovdens*Radius**3 ! virial mass for outer radial bin d0 = Radius**2 ! get final statistics of bound particles Ncount = 0 Mtot = 0. wx = 0. ; wy = 0. ; wz = 0. Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) ! Get mass profile and statistics Ncount =0 Do k3 =k1, k2 Do j3 =j1, j2 Do i3 =i1, i2 jp =Label(i3,j3,k3) Do while (jp.ne.0) dd =(x-Xpar(jp))**2+(y-Ypar(jp))**2+(z-Zpar(jp))**2 If(dd< d0) Then r = sqrt(dd) !ii = max(min(INT(log10(r/Radius)/dLogR),0),-Nrad) ii = min(max(INT((r/dCell)/dR),0),Nrad) MassP(ii) = MassP(ii) + MassOne Ncount = Ncount +1 !dx = Xpar(jp) -x !dy = Ypar(jp) -y !dz = Zpar(jp) -z wx = VX(jp) + wx ! average velocity wy = VY(jp) + wy wz = VZ(jp) + wz !ax = VX(jp)**2 + ax ! velocity dispersion !ay = VY(jp)**2 + ay !az = VZ(jp)**2 + az !xcm = xcm + dx ! center of mass !ycm = ycm + dy !zcm = zcm + dz !x2 = x2 + (VX(jp)+factorZ*dx)**2 ! !y2 = y2 + (VY(jp)+factorZ*dy)**2 !z2 = z2 + (VZ(jp)+factorZ*dz)**2 EndIf ! dd0.)Then dd = sqrt(dd) ! distance in kpch !write(*,'(2i10,3x,f9.3,3x,2f8.3,3x,2es12.4,3x,6f9.3)')ip,ic,dd,Rvir(ip),Rvir(ic),Mtotal(ip),Mtotal(ic), & ! x,xMaxx(ic),y,yMaxx(ic),z,zMaxx(ic) if(dd<(Rvir(ip)+Rvir(ic))*Environ)Then Dist(ip) = dd MassNext(ip) = MassNext(ip) + Mtotal(ic) !write(*,'(3x,2i10,f9.4,3x,2f8.3,3x,2es12.4,3x,6f9.3)') ip,ic,dd,Rvir(ip),Rvir(ic),Mtotal(ip),Mtotal(ic), & ! x,xMaxx(ic),y,yMaxx(ic),z,zMaxx(ic) End If end If ic = Lst(ic) end Do end Do end Do end Do End Do ! ip Nclose = 0 dx = Box/mesh Do ip=1, Nmaxima If(Rvir(ip)<0..or.Mtotal(ip)Mtotal(ip))Then Nclose = Nclose +1 If(alpha*Rvir(ip)/Dist(ip)>0.3)write(50,'(i9,2es12.4,3x,3f9.2,3x,3f8.3,3x,3f9.2)')ip,Mtotal(ip),MassNext(ip), & Dist(ip),Rvir(ip),Rvir(ip)/Dist(ip), & xMaxx(ip),yMaxx(ip),zMaxx(ip),VxMaxx(ip),VyMaxx(ip),VzMaxx(ip) Mtotal(ip) = Mtotal(ip)/(1.+alpha*Rvir(ip)/Dist(ip)) Rvir(ip) = Rvir(ip)/(1.+alpha*Rvir(ip)/Dist(ip))**2 x = xMaxx(ip); y = yMaxx(ip); z = zMaxx(ip) i = INT(x/dx+0.5) ! nearest node j = INT(y/dx+0.5) k = INT(z/dx+0.5) sx = VxMaxx(ip) -Vf(1,i,j,k) sy = VyMaxx(ip) -Vf(2,i,j,k) sz = VzMaxx(ip) -Vf(3,i,j,k) VxMaxx(ip) = Vf(1,i,j,k)+ sx*(1.+3.*alpha/(Dist(ip)/Rvir(ip))) VyMaxx(ip) = Vf(2,i,j,k)+ sy*(1.+3.*alpha/(Dist(ip)/Rvir(ip))) VzMaxx(ip) = Vf(3,i,j,k)+ sz*(1.+3.*alpha/(Dist(ip)/Rvir(ip))) !VxMaxx(ip) = VxMaxx(ip)/(1.+alpha*Rvir(ip)/Dist(ip)) !VyMaxx(ip) = VyMaxx(ip)/(1.+alpha*Rvir(ip)/Dist(ip)) !VzMaxx(ip) = VzMaxx(ip)/(1.+alpha*Rvir(ip)/Dist(ip)) end If end Do DEALLOCATE(Dist,MassNext) DEALLOCATE(Vf) write(*,*) ' Number of corrected halos = ',Nclose end SUBROUTINE CorrectClose !--------------------------------------------------------------------------- ! ! SUBROUTINE VelocityField !--------------------------------------------------------------------------- integer(kind=4) :: ic,ip Real(kind=8) :: ww,dd,x,y,z,Rr,R2max,sx,sy,sz mesh =100 dx = Box/mesh write(*,'(a,i4,a,f8.2)') ' velocity mesh =',mesh,' cell size= ',dx ALLOCATE(Vf(4,0:mesh,0:mesh,0:mesh)) Vf(:,:,:,:) = 0. ! !---- Do ip=1, Nmaxima x = xMaxx(ip) ; y = yMaxx(ip) ; z = zMaxx(ip) sx = VxMaxx(ip); sy = VyMaxx(ip); sz = VzMaxx(ip) if(x>Box.or.y>Box.or.z>Box)print *,x,y,z If(Rvir(ip)<0.)Cycle i = min(INT(x/dx),mesh-1) j = min(INT(y/dx),mesh-1) k = min(INT(z/dx),mesh-1) if(k+1>mesh)print *,'Error in Velocity mes: ',x,y,z do i1=i,i+1 do j1=j,j+1 do k1=k,k+1 Vf(4,i1,j1,k1) = Vf(4,i1,j1,k1) +1. Vf(1,i1,j1,k1) = Vf(1,i1,j1,k1) +sx Vf(2,i1,j1,k1) = Vf(2,i1,j1,k1) +sy Vf(3,i1,j1,k1) = Vf(3,i1,j1,k1) +sz end do end do end do end Do Do i=0,mesh Do j=0,mesh Do k=0,mesh if(Vf(4,i,j,k)>0)then ic = Vf(4,i,j,k) Vf(1,i,j,k) = Vf(1,i,j,k)/ic Vf(2,i,j,k) = Vf(2,i,j,k)/ic Vf(3,i,j,k) = Vf(3,i,j,k)/ic end if end Do end Do end Do ! j= 10; k= 20 ! do i=0,25 ! write(50,'(3i3,3x,f8.1,3x,3f9.2)')i,j,k,Vf(4,i,j,k),Vf(1,i,j,k),Vf(2,i,j,k),Vf(3,i,j,k) ! end do ! j= 11; k= 20 ! do i=0,25 ! write(50,'(3i3,3x,f8.1,3x,3f9.2)')i,j,k,Vf(4,i,j,k),Vf(1,i,j,k),Vf(2,i,j,k),Vf(3,i,j,k) ! end do ! j= 12; k= 20 ! do i=0,25 ! write(50,'(3i3,3x,f8.1,3x,3f9.2)')i,j,k,Vf(4,i,j,k),Vf(1,i,j,k),Vf(2,i,j,k),Vf(3,i,j,k) ! end do end SUBROUTINE VelocityField !--------------------------------------------------------------------------- ! get centers and velocities of maxima ! ! SUBROUTINE FindDistinctCandidates !--------------------------------------------------------------------------- integer(kind=8) :: i,im,ic,ip,jp tstart = seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i =1, Nmaxima Mvir(i) = 0. Rvir(i) = 0. VxMaxx(i) = 0. VyMaxx(i) = 0. VzMaxx(i) = 0. EndDo ! --- improve halo positions and halo velocities Do iter =1,1 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (im,x,y,z,xc,yc,zc,xv,yv,zv,nn,i3,j3,k3,i1,i2,j1,j2,k1,k2,jp,dd,Radius,d0) Do im =1, Nmaxima x = xMaxx(im) y = yMaxx(im) z = zMaxx(im) xc = 0. yc = 0. zc = 0. xv = 0. yv = 0. zv = 0. nn = 0 Radius = Cell !!*max(0.5,min(2.,log10(Xoff(im)+10.)/2.)) d0 = Radius**2 Call Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) Do k3 =k1, k2 Do j3 =j1, j2 Do i3 =i1, i2 jp =Label(i3,j3,k3) Do while (jp.ne.0) dd = (x-Xpar(jp))**2+(y-Ypar(jp))**2+(z-Zpar(jp))**2 If(dd< d0) Then xc = xc + Xpar(jp) yc = yc + Ypar(jp) zc = zc + Zpar(jp) xv = xv + VX(jp) yv = yv + VY(jp) zv = zv + VZ(jp) nn = nn + 1 end If jp =Lst(jp) End Do ! jp/= 0 EndDo ! i3 EndDo ! j3 EndDo ! k3 !if(nn.le.1)write(18,'(a,i9,3f9.4,3i9,g12.4)') ' no particles: ',im,x,y,z,nn,k1,k2,Xoff(im) !write(18,'(i9,3f9.4,6i9)') im,x,y,z,nn,k1,k2 If(nn>10)Then xMaxx(im) = xc/nn yMaxx(im) = yc/nn zMaxx(im) = zc/nn If(xMaxx(im).lt.0.) xMaxx(im) = xMaxx(im)+Box If(xMaxx(im).ge.Box)xMaxx(im) = xMaxx(im)-Box If(yMaxx(im).lt.0.) yMaxx(im) = yMaxx(im)+Box If(yMaxx(im).ge.Box)yMaxx(im) = yMaxx(im)-Box If(zMaxx(im).lt.0.) zMaxx(im) = zMaxx(im)+Box If(zMaxx(im).ge.Box)zMaxx(im) = zMaxx(im)-Box VxMaxx(im) = xv/nn VyMaxx(im) = yv/nn VzMaxx(im) = zv/nn Mvir(im) = nn*MassOne Rvir(im) = Radius end If End Do ! write(18,*) ' Iter = ',iter ! Do im=1,Nmaxima ! write(18,'(i9,16f10.4)') im,xMaxx(im),yMaxx(im),zMaxx(im),VxMaxx(im),VyMaxx(im),VzMaxx(im), Xoff(im) ! end Do End Do ! iter 10 tfinish = seconds() write(*,'(10x,a,T50,2f10.2)') ' time for FindDistinctCandidates (secs) =',tfinish-tstart,tfinish-t0 write(13,'(10x,a,T50,2f10.2)') ' time for FindDistinctCandidates (secs) =',tfinish-tstart,tfinish-t0 end SUBROUTINE FindDistinctCandidates !--------------------------------------------------------------------------- ! ! ! ! SUBROUTINE FindMaxima !--------------------------------------------------------------------------- Integer(kind=8) :: i,j,M,Mnow,Nbuff,Nm Integer(kind=4), allocatable, dimension(:) :: Mth Real(kind=4), allocatable, dimension(:,:) :: Xxoff,xxMax,yyMax,zzMax Integer(kind=4) :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM,Nthreads tstart = seconds() Om0 = Om ! Set scales Ovdens = 200./Om0*(Om0+(1.-Om0)*AEXPN**3)/OvdensMin !Ovdens = 200./Om0/1.2 Dmaxim = 1.e7 !$OMP PARALLEL Nthreads =OMP_GET_NUM_THREADS() !$OMP end parallel write(*,*) ' Number of threads =',Nthreads print *,' iVirial=',iVirial write(*,'(a,T40,i11,2f9.4)') ' Inside FindMaxima. Nparticles= ',Nparticles,Ovdens,Om0 !Dmaxim =0. !!$OMP PARALLEL DO DEFAULT(SHARED) & ! ------- find Maximum density !!$OMP PRIVATE (M3,M2,M1) REDUCTION(MAX:Dmaxim) ! DO M3=1,NGRID ! DO M2=1,NGRID ! DO M1=1,NGRID ! Dmaxim =max(FI(M1,M2,M3),Dmaxim) ! end DO ! end DO ! end DO tfinish = seconds() ! write(*,*) ' Maximum density =',Dmaxim ! write(*,*) ' time for Dmax =',tfinish-tstart,tfinish-t0 Nmaxima =0 !$OMP PARALLEL DO DEFAULT(SHARED) & ! ------- count how many maxima we have !$OMP PRIVATE (M3,M2,M1,M3p,M3m,M2p,M2m,M1p,M1m,iMax) REDUCTION(+:Nmaxima) DO M3=1,NGRID M3p =M3+1 ; If(M3p>NGRID)M3p=M3p-NGRID M3m =M3-1 ; If(M3m<1)M3m=M3m+NGRID DO M2=1,NGRID M2p =M2+1 ; If(M2p>NGRID)M2p=M2p-NGRID M2m =M2-1 ; If(M2m<1)M2m=M2m+NGRID DO M1=1,NGRID M1p =M1+1 ; If(M1p>NGRID)M1p=M1p-NGRID M1m =M1-1 ; If(M1m<1)M1m=M1m+NGRID iMax = 1 ! look for all 26 neighbors If(FI(M1,M2,M3).lt.FI(M1m,M2m,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1 ,M2m,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2m,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1, M2,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2p,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1, M2p,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2p,M3m))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2m,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1 ,M2m,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2m,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2p,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1, M2p,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2p,M3))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2m,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1 ,M2m,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2m,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1, M2,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1m,M2p,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1, M2p,M3p))iMax=0 If(FI(M1,M2,M3).lt.FI(M1p,M2p,M3p))iMax=0 If(iMax==1.and.FI(M1,M2,M3)>Ovdens)Then Nmaxima = Nmaxima +1 FI(M1,M2,M3) = FI(M1,M2,M3) +Dmaxim+1. ! assign large number to the maximum end If END DO END DO END DO t1 = seconds() write(*,*) ' time for Find Max =',t1-t0 write(*,*) ' FindMaxima: 1 finished: Nmaxima= ',Nmaxima Nbuff = Nmaxima/Nthreads *7+1000 ! length of the buffer Allocate(Mth(Nbuff),xxMax(Nbuff,Nthreads),yyMax(Nbuff,Nthreads),zzMax(Nbuff,Nthreads)) Allocate(Xxoff(Nbuff,Nthreads)) if(Nmaxima == 0)Stop ' No density maxima found' ALLOCATE (Mtotal(Nmaxima),Mvir(Nmaxima),Rvir(Nmaxima),Xoff(Nmaxima)) myMemory =Memory(4_8*Nmaxima ) ALLOCATE (xMaxx(Nmaxima),yMaxx(Nmaxima),zMaxx(Nmaxima)) myMemory =Memory(3_8*Nmaxima ) ALLOCATE (VxMaxx(Nmaxima),VyMaxx(Nmaxima),VzMaxx(Nmaxima)) myMemory =Memory(3_8*Nmaxima ) ALLOCATE(LstMax(Nmaxima),VmaxM(Nmaxima)) Nmaxx =0 xs = Box/NGRID Mth(:) = 0 Mm = 1 !$OMP PARALLEL DO PRIVATE(M1,M2,M3,Mthread) Firstprivate(Mm) DO M3=1,NGRID !--- make list of maxima DO M2=1,NGRID DO M1=1,NGRID If(FI(M1,M2,M3)>Dmaxim+1.)Then Mthread = OMP_GET_THREAD_NUM()+1 Mth(Mthread) = Mm If(Mm>Nbuff)Then write(*,'(3i6,i11,i4,i11)') M1,M2,M3,Mm,Mthread,Nbuff Stop ' Too many elements for buffer XX' end If xxMax(Mm,Mthread) = (M1-1)*xs yyMax(Mm,Mthread) = (M2-1)*xs zzMax(Mm,Mthread) = (M3-1)*xs Xxoff(Mm,Mthread) = FI(M1,M2,M3) -Dmaxim-1. Mm = Mm +1 end If end DO end DO end DO M = 0 Do i=1,Nthreads M = M + Mth(i) end Do write(*,*) ' Elements on all buffers =',M,Nmaxima i = 0 Do Mthread=1,Nthreads Mnow = Mth(Mthread) Do j = 1,Mnow Xoff(j+i) = Xxoff(j,Mthread) xMaxx(j+i) = xxMax(j,Mthread) yMaxx(j+i) = yyMax(j,Mthread) zMaxx(j+i) = zzMax(j,Mthread) End Do i = i+ Mnow end Do Deallocate(Xxoff,Mth,xxMax,yyMax,zzMax) SELECT CASE (iVirial) CASE (0) Ovdens = 200./Om0*(Om0+(1.-Om0)*AEXPN**3) ! 200\rho_critical CASE (1) Ovdens = OverdenVir() ! virial overdensity CASE (2) Ovdens = 200. ! 200 matter overdens CASE (23) Ovdens = OverdenAbacus() ! Abacus overdens end SELECT t2 = seconds() iOverdens = 0 i300 = 0; i400 =0; i500 = 0; i600 =0; i700 =0 write(*,*) ' time for List Max =',t2-t1 !--- statistics of maxima Do i=1,Nmaxima If(Xoff(i) .ge. Ovdens)iOverdens = iOverdens +1 If(Xoff(i)>300.)Then i300 = i300 +1 if(Xoff(i)>400.)Then i400 = i400 +1 if(Xoff(i)>500.)Then i500 = i500 +1 if(Xoff(i)>600.)Then i600 = i600 +1 if(Xoff(i)>700.)Then i700 = i700 +1 End if end if end if end if end If enddo write(*,'(a,i10,/a,i9,5(a,i8))') ' Number of maxima:',Nmaxima, & ' Above Overdensity =',iOverdens, & ' Above: 300 =',i300,' 400 =',i400,' 500 =',i500,' 600 =',i600,' 700 =',i700 tfinish = seconds() write(13,'(10x,a,T50,2f10.2)') ' time for FindMaxima(secs) =',tfinish-tstart,tfinish-t0 write(*,'(10x,a,T50,2f10.2)') ' time for FindMaxima(secs) =',tfinish-tstart,tfinish-t0 end SUBROUTINE FindMaxima !--------------------------------------------------------------------------- ! Define size and boundaries of linked-list SUBROUTINE SizeListMaxima !--------------------------------------------------------------------------- Real(kind=4) :: MemList,MemMaxList, & fractionMemory = 0.90 ! fraction of memory allocated to Lists Cell = 3.5 !Roptimal k = 0 write(13,'(2(a,g12.4))') ' SizeListMaxima: ',Cell Nmx = (Xleft -dBuffer)/Cell - 1 Nmy = (Yleft -dBuffer)/Cell - 1 Nmz = (Zleft -dBuffer)/Cell - 1 Nbx = (Xright +dBuffer)/Cell + 1 Nby = (Yright +dBuffer)/Cell + 1 Nbz = (Zright +dBuffer)/Cell + 1 write(13,'(a,g12.3,a,6i6)') ' Size List Maxima: Cell=',Cell, & ' Limits=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz end SUBROUTINE SizeListMaxima !--------------------------------------------------------------------------- ! Define size and boundaries of linked-list SUBROUTINE SizeList !--------------------------------------------------------------------------- Real(kind=4) :: MemList,MemMaxList, & fractionMemory = 0.90 ! fraction of memory allocated to Lists !Roptimal = (Nne*MassOne/(Om0*Ovdens))**(1./3.)*9.50e-5 ! Mpch inits !Roptimal = max(Roptimal,Box/NROW) * 1.2 MemMaxList = (MaxMemory - TotalMemory-12.*Np/1024.**3)*fractionMemory If(MemMaxList<0.1)Then STOP ' Stop: Not enough memory !!!' End If Cell = (Box/NGRID*3.) !Roptimal k = 0 write(13,'(2(a,g12.4))') ' SizeList: Roptimal=',Cell,' Allowed Memory=',MemMaxList Do Nmx = (Xleft -dBuffer)/Cell - 1 Nmy = (Yleft -dBuffer)/Cell - 1 Nmz = (Zleft -dBuffer)/Cell - 1 Nbx = (Xright +dBuffer)/Cell + 1 Nby = (Yright +dBuffer)/Cell + 1 Nbz = (Zright +dBuffer)/Cell + 1 MemList = (3.*4.*Real(Nbx-Nmx+1)*Real(Nby-Nmy+1)*Real(Nbz-Nmz+1))/1024.**3 If(MemList < MemMaxList)Exit k = k+1 Cell = Cell *1.1 If(k > 100)STOP ' Stop: too many iterations in SizeList' EndDo write(13,'(a,g12.3,a,6i6)') ' Size List: Cell=',Cell, & ' Limits=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz end SUBROUTINE SizeList !-------------------------------------------------------------- ! Make linker lists of particles in each cell SUBROUTINE List !-------------------------------------------------------------- Integer(kind=8) :: Nm,N0,N1,N2,N3,iMax,jp Integer(kind=4) :: Ndiv(0:1000) Integer(kind=4) :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM,Nthreads !$OMP PARALLEL Nthreads =OMP_GET_NUM_THREADS() !$OMP end parallel Nsplit = Nthreads d = (Nbx-Nmx)/Real(Nsplit) !--- parallelization setup Do i=0,Nsplit Ndiv(i) = Floor(d*i+Nmx+0.5) End Do Ndiv(Nsplit) = Nbx+1 tstart =seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do jp=1,Np Lst(jp)=-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 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (jp,i,j,k,Isplit) Do Isplit =1,Nsplit Do jp=1,Np i=Ceiling(Xpar(jp)/Cell)-1 j=Ceiling(Ypar(jp)/Cell)-1 k=Ceiling(Zpar(jp)/Cell)-1 i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) k=MIN(MAX(Nmz,k),Nbz) If(k.ge.Ndiv(Isplit-1).and.k.lt.Ndiv(Isplit))Then Lst(jp) = Label(i,j,k) Label(i,j,k) = jp End If EndDo end Do t1 = seconds() write(*,*) ' time to make list =',t1-tstart end SUBROUTINE List !-------------------------------------------------------------- ! Make linker lists of maxima in each cell SUBROUTINE ListMaxima !-------------------------------------------------------------- Integer(kind=8) :: Nm,N0,N1,N2,N3,iMax,jp tstart = seconds() !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do jp=1,Nmaxima Lst(jp)=-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,Nmaxima i=Ceiling(Xmaxx(jp)/Cell)-1 j=Ceiling(Ymaxx(jp)/Cell)-1 k=Ceiling(Zmaxx(jp)/Cell)-1 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 t1 = seconds() write(*,*) ' time to make list Max =',t1-tstart end SUBROUTINE ListMaxima !-------------------------------------------------------------- ! eigenvalues and eigenvectors of ! 3x3 symmetric matrix ! Use Reyleigh power iterations to find ! maximum eigenvalue and eigenvector ! then use Trace(A) = sum(eigenvalues) and ! product of eigenvalues = det(A) ! to get other two eigenvalues SUBROUTINE EigenValues(A,x,EigVal) Integer(kind=4), Parameter :: Nmax =200 Real*8, Parameter :: error =2.d-7 Real(kind=4) :: EigVal(3),x(3) Real*8 :: A(3,3),xx(3),y(3),Eig,EigNew,a1,b1,Trace,detA,x1,x2 Integer(kind=4) :: ind(3) xx = 1. ! initial guess Eig = 0. Do kstep=1,Nmax y = MATMUL(A,xx) EigNew = max(y(1),y(2),y(3)) xx = y/EigNew !write(13,'(10x,i4,1p,2g13.5,3x,4g13.5)')kstep,EigNew,Eig,xx If(abs(EigNew-Eig)< error*abs(EigNew))exit Eig = EigNew End Do Eig = EigNew Trace = A(1,1) +A(2,2) +A(3,3) detA = A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) & -A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)) & +A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) ! solve quadratic equation for other eigenvalues a1 = (Trace-Eig)/2. b1 = detA/Eig d = sqrt(max(a1**2-b1,1.d-20)) x1 = a1 + d x2 = a1 - d !EigVal(1) = max(Eig,x1,x2) !EigVal(3) = min(Eig,x1,x2) EigVal = -1.e30 ind = 0 If(Eig.ge.x1.and.Eig.ge.x2)Then ! Eig is the max EigVal(1) = Eig If(x1.ge.x2)Then EigVal(2) = x1 EigVal(3) = x2 Else EigVal(2) = x2 EigVal(3) = x1 EndIf Else If(x1.ge.x2)Then ! x1 is max EigVal(1) = x1 If(Eig.ge.x2)Then EigVal(2) = Eig EigVal(3) = x2 Else EigVal(2) = x2 EigVal(3) = Eig EndIf Else ! x2 is max EigVal(1) = x2 If(Eig.ge.x1)Then EigVal(2) = Eig EigVal(3) = x1 Else EigVal(2) = x1 EigVal(3) = Eig EndIf End If ! If(EigVal(2)>EigVal(1).or.EigVal(3)>EigVal(2))Then ! write(13,'(/5x,1p,10g13.5)')Trace,detA,a1,b1,x1,x2,EigVal ! write(13,'(1p,3g12.4)') A ! write(13,'(a,1p,g13.5)') ' Sum =',EigVal(1)+EigVal(2)+EigVal(3) ! write(13,'(a,1p,g13.5)') ' Prod =',EigVal(1)*EigVal(2)*EigVal(3) ! write(13,'(a,1p,2g13.5)') ' Iter =',EigNew-Eig,Eig ! write(13,'(10x,3i3)') ind ! EndIf x = xx d = sqrt(SUM(x**2)) x = x/d end SUBROUTINE EigenValues !-------------------------------------------------------------- ! ! Find halo concentration using M,R, and Vmax ! M - in Msunh, R - comoving kpch ! Vmax = in km/s real function Concentration(aM,aR,Vmax) ! ------------------------ ! Real*8, parameter :: d0 =1.d0, dCmin =3.d-5, C0= 2.162d0 Real*8 :: M,V,R,A,C,Fc,Vvir,Vratio,dC If(aM<1.e-6.or.aR<1.e-6.or.Vmax <1.e-6)Then Concentration = 0. ; return End If Vvir = 2.076d-3*sqrt(aM/(aR*AEXPN)) Vratio = Vmax/Vvir If(Vratio<1.)Then Concentration = -1. Else dC = 2. C= 2.*C0 A = Vratio**2 !N = 0 Do while (abs(dC)>dCmin) Fc = 0.2162166_8/(log(1.d0+C)/C -1.d0/(1.d0+C)) If(Fc.gt.A.and.C.ge.C0)Then dC = dC/2.d0 C = C - dC Else C = C + dC EndIf !N = N +1 !write(*,'(i6,1p,4G13.6)')N,C,dC,Fc,A End Do Concentration = C End If End function Concentration !--------------------------------------------------------------------------- ! find limits for the linker-list search SUBROUTINE Limits(x,y,z,Radius,i1,i2,j1,j2,k1,k2) !---------------------------------------------------------------------------- i2=Ceiling((x+Radius)/Cell)-1 j2=Ceiling((y+Radius)/Cell)-1 k2=Ceiling((z+Radius)/Cell)-1 i1=Ceiling((x-Radius)/Cell)-1 j1=Ceiling((y-Radius)/Cell)-1 k1=Ceiling((z-Radius)/Cell)-1 i1=MIN(MAX(Nmx,i1),Nbx) j1=MIN(MAX(Nmy,j1),Nby) k1=MIN(MAX(Nmz,k1),Nbz) i2=MIN(MAX(Nmx,i2),Nbx) j2=MIN(MAX(Nmy,j2),Nby) k2=MIN(MAX(Nmz,k2),Nbz) end SUBROUTINE Limits !---------------------------------------------------------------------------- ! Define configuration of files ! Read control info from the first file ! allocate arrays ! Subroutine SetParameters Integer(kind=4) :: jstep Integer :: FileList(3)=(/12,20,30/) logical :: exst character*80 :: Name,CatLabel Character*40 :: txt1,txt2,txt3,txt4,txt5,txt6,txt7,txt8 Character*52 :: txt2b Character*10 :: sxt1,sxt2,sxt3,sxt4 Om0 = Om ! Set scales Xscale = Box/NGRID ! Scale for comoving coordinates Vscale = 100.*Xscale/AEXPN ! Scale for velocities Dscale = 2.774e+11*(Box/NROW)**3 ! mass scale MassOne= Om0*Dscale ! mass of the smallest particle SELECT CASE (iVirial) CASE (0) Ovdens = 200./Om0*(Om0+(1.-Om0)*AEXPN**3) ! 200\rho_critical CASE (1) Ovdens = OverdenVir() ! virial overdensity CASE (2) Ovdens = 200. ! 200 matter overdens CASE (3) Ovdens = OverdenAbacus() ! Abacus overdens end SELECT Xleft = 0. ; Xright = Box Yleft = 0. ; Yright = Box Zleft = 0. ; Zright = Box write(13,*) ' Boundaries(Mpch) =',Xleft,Xright write(13,*) ' Boundaries =',Yleft,Yright write(13,*) ' Boundaries =',Zleft,Zright write(13,*) ' Buffer zone =',dBuffer write(*,*) ' Scales=', Xscale, Vscale Do kf =1,1 kfile = FileList(kf) If(kf==2.or.kf==3)Then WRITE (kfile) HEADER sxt1 =' A =' sxt2 =' Step =' WRITE (kfile) sxt1,AEXPN,sxt2,ASTEP sxt1 =' I =' sxt2 =' Nrow =' sxt3 =' Ngrid=' WRITE (kfile) sxt1,ISTEP,sxt2,NROW,sxt3,NGRID sxt1 =' Omega_0=' sxt2 =' Omega_L=' sxt3 =' hubble =' txt4 =' buffer width (Mpch) =' WRITE (kfile) sxt1,Om0,sxt2,1.-Om0,sxt3,hubble,txt4,dBuffer txt1= ' Number of radial bins =' WRITE (kfile) txt1,Nrad txt1= ' Mass of smallest particle (Msunh) =' WRITE (kfile) txt1,MassOne txt1= ' Overdensity limit =' WRITE (kfile) txt1,Ovdens Txt1 =' XYZ(Mpch) ' Txt2 ='Vxyz(km/s) Mvir/Msunh ' Txt2b='Vxyz(km/s) Mbound Mtot/Msunh ' Txt3 =' Rvir(kpch) Vcirc(km/s)' Txt4 =' Nhalo Nparticles ' Txt5 =' Xoff 2K/Ep-1 Lambda RadRMS/kpch' Txt6 =' b/a c/a MajorAxis: x y z' WRITE (kfile) txt1,txt2,txt3,txt4 Else WRITE (kfile,'(a)') HEADER sxt1 =' A =' sxt2 =' Step =' WRITE (kfile,'(2(a,f8.5))') sxt1,AEXPN,sxt2,ASTEP sxt1 =' I =' sxt2 =' Nrow =' sxt3 =' Ngrid=' WRITE (kfile,'(3(a,i5))') sxt1,ISTEP,sxt2,NROW,sxt3,NGRID sxt1 =' Omega_0=' sxt2 =' Omega_L=' sxt3 =' hubble =' txt4 =' buffer width (Mpch) =' WRITE (kfile,'(6(a,f8.4))') sxt1,Om0,sxt2,1.-Om0,sxt3,hubble,txt4,dBuffer txt1= ' Number of radial bins =' WRITE (kfile,'(a,i4)') txt1,Nrad txt1= ' Mass of smallest particle (Msunh) =' WRITE (kfile,'(a,1p,g12.4)') txt1,MassOne txt1= ' Overdensity limit =' WRITE (kfile,'(a,f8.4)') txt1,Ovdens Txt1 =' XYZ(Mpch) ' Txt2 ='Vxyz(km/s) Mvir/Msunh ' Txt2b='Vxyz(km/s) Mbound Mtot/Msunh ' Txt3 =' Rvir(kpch) Vcirc(km/s)' Txt4 =' Nhalo Nparticles ' Txt5 =' Xoff 2K/Ep-1 Lambda RadRMS/kpch' Txt6 =' b/a c/a MajorAxis: x y z' WRITE (kfile,'(6a)') txt1,txt2,txt3,txt4 end If end Do Txt1= 'R/kpch Npart R/Rvir Mass Vcirc' Txt2= ' Dens/Msun/comvKpch Vrms Vrad Vradrms' Txt3= ' Bound:R/kpch Npart R/Rvir Mass' Txt4= ' Vcirc Dens/Msun/comvKpch Vrms Vrad' Txt5= ' Vradrms' ! WRITE(20,'(3a)') 'R/kpch Npart R/Rvir Mass Vcirc Dens/Msun/comvKpch', & ! ' Vrms Vrad Vradrms Bound:R/kpch Npart R/Rvir Mass', & ! ' Vcirc Dens/Msun/comvKpch Vrms Vrad Vradrms' !100 FORMAT(1X,'Header=>',A45,/ & ! ' A=',F8.3,' Step=',F8.3,/ & ! ' I =',I4,' Nrow=',I4,' Ngrid=',I4,/' Omega_0=',f6.3, & ! ' Omega_L=',f6.3,' h=',f6.3,' buffer width (Mpch) =',f9.4) !110 FORMAT(2(a,T50,f8.3/),4(a,T50,i5/),a,T50,f8.3/,a,T50,g12.4,/4(a,T50,f8.3/)) !120 FORMAT(5x,a,T40,a,T68,a,T81,a,T92,a,T103,a,T121,a,T128,a,T137,a, & ! T152,a,T166,a,T178,a,T187,a,T200,a,T213,a,T222,a,T232,a) write (13,'(6x,"Mass of smallest particle",/9x,"in units M_sun/h is =",3x,g10.3)') MassOne write(*,*) ' Np = ',Np,Nparticles end Subroutine SetParameters !---------------------------------------------------------------------------- ! Add buffer around the computational box ! -- ! -- ! -- move coordinates to new arrays Subroutine AddBuffer Integer(kind=4) :: jstep Integer :: FileList(3)=(/12,20,30/) Real(kind=4), ALLOCATABLE, DIMENSION(:) :: & Xbb, Ybb, Zbb, & ! coords VXbb, VYbb, Vzbb ! velocities Integer(kind=8) :: idummy,ic,ip,iPartMax,i Factor = (1.+4.*dBuffer/Box)**3 iPartMax = Factor*Np ! reserve extra space for total N particles Nparticles = Np ! old number of particles write(*,*) ' Allocate buffers for particles. N=',iPartMax myMemory= Memory(6_8*iPartMax) ALLOCATE(Xbb(iPartMax),Ybb(iPartMax),Zbb(iPartMax)) ALLOCATE(VXbb(iPartMax),VYbb(iPartMax),VZbb(iPartMax)) Xleft = 0. ; Xright = Box Yleft = 0. ; Yright = Box Zleft = 0. ; Zright = Box !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic) Do ic=1,iPartMax VXbb(ic)= 0. VYbb(ic)= 0. VZbb(ic)= 0. Xbb(ic)= 0. Ybb(ic)= 0. Zbb(ic)= 0. End Do ip = 0 Xl = Xleft -dBuffer ; Xr = Xright +dBuffer Yl = Yleft -dBuffer ; Yr = Yright +dBuffer Zl = Zleft -dBuffer ; Zr = Zright +dBuffer write(*,*) ' Replicate coordinates' !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic,xx,yy,zz) Do ic=1,Np !xx = Xp(ic) ; yy =Yp(ic) ; zz =Zp(ic) !Xb(ic) =xx ; Yb(ic) = yy ; Zb(ic) = zz Xbb(ic)= Xpar(ic) Ybb(ic)= Ypar(ic) Zbb(ic)= Zpar(ic) End Do write(*,*) ' Replicate velocities' !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic,xx,yy,zz) Do ic=1,Np VXbb(ic)= VX(ic) VYbb(ic)= VY(ic) VZbb(ic)= VZ(ic) End Do write(*,*) ' add buffer' ip = Np !!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic,xx,yy,zz,x,y,z,k,j,i) Do ic=1,Np xx = Xpar(ic) ; yy =Ypar(ic) ; zz =Zpar(ic) do k = -1,1 z = zz+k*Box if(zZr)cycle do j = -1,1 y = yy+j*Box if(yYr)cycle do i = -1,1 x = xx+i*Box if(xXr)cycle if((xXright) .or. & (yYright) .or. & (zZright))Then ip = ip +1 If(ip>iPartMax)Stop ' Too many buffer particles. Increase Factor in AddBuffer' Xbb(ip) = x Ybb(ip) = y Zbb(ip) = z VXbb(ip)= VX(ic) VYbb(ip)= VY(ic) VZbb(ip)= VZ(ic) end if end do end do end do end Do write(*,*) ' New number of particles = ',ip myMemory=Memory(-6_8*Np) DEALLOCATE(Xpar,Ypar,Zpar) DEALLOCATE(VX,VY,VZ) Np = ip Nparticles = Np myMemory=Memory(6_8*Np) ALLOCATE(Xpar(Np),Ypar(Np),Zpar(Np)) ALLOCATE(VX(Np),VY(Np),VZ(Np)) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) Do ip =1,Np Xpar(ip) = Xbb(ip) ; Ypar(ip) = Ybb(ip) ; Zpar(ip) = Zbb(ip) VX(ip)= VXbb(ip); VY(ip)= VYbb(ip); VZ(ip)= VZbb(ip) EndDo myMemory= Memory(-6_8*iPartMax) DEALLOCATE(Xbb,Ybb,Zbb) DEALLOCATE(VXbb,VYbb,VZbb) xmin = 1.e12; xmax = -1.e12 ymin = 1.e12; ymax = -1.e12 zmin = 1.e12; zmax = -1.e12 !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) & !$OMP REDUCTION(MAX:xmax,ymax,zmax) REDUCTION(MIN:xmin,ymin,zmin) Do ip =1,Np xmin =MIN(xmin,Xpar(ip)) ymin =MIN(ymin,Ypar(ip)) zmin =MIN(zmin,Zpar(ip)) xmax =MAX(xmax,Xpar(ip)) ymax =MAX(ymax,Ypar(ip)) zmax =MAX(zmax,Zpar(ip)) If(Xpar(ip).lt.Xl.or.Ypar(ip).lt.Yl.or.Zpar(ip).lt.Zl)& write(*,'(a,i10,1p,3g14.5)')' Error coord: ',i,Xpar(i),Ypar(i),Zpar(i) EndDo ! write(*,'(a,20es13.5)') ' X= ',(Xpar(i),i=1,10) ! write(*,'(a,20es13.5)') ' Y= ',(Ypar(i),i=1,10) ! write(*,'(a,20es13.5)') ' Z= ',(Zpar(i),i=1,10) ! write(*,'(a,20es13.5)') ' X= ',(Xpar(i),i=Np-9,Np) ! write(*,'(a,20es13.5)') ' Y= ',(Ypar(i),i=Np-9,Np) ! write(*,'(a,20es13.5)') ' Z= ',(Zpar(i),i=Np-9,Np) write(*,'(a,2es13.5)') ' X range= ',xmin,xmax write(*,'(a,2es13.5)') ' Y range= ',ymin,ymax write(*,'(a,2es13.5)') ' Z range= ',zmin,zmax end Subroutine AddBuffer !---------------------------------------------------------------------------- ! Remove buffer around the computational box ! ! Subroutine RemoveBuffer(NpPM) Real(kind=4), ALLOCATABLE, DIMENSION(:) :: & Xbb, Ybb, Zbb, & ! coords VXbb, VYbb, Vzbb ! velocities Integer(kind=8) :: NpPM,ic Xscale = Box/NGRID ! Scale for comoving coordinates Vscale = 100.*Xscale/AEXPN ! Scale for velocities myMemory= Memory(6_8*Nparticles) ALLOCATE(Xbb(Nparticles),Ybb(Nparticles),Zbb(Nparticles)) ALLOCATE(VXbb(Nparticles),VYbb(Nparticles),VZbb(Nparticles)) !--- copy data to buffers write(*,*) ' Copy to new buffer ',NpPM,Xscale !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic) Do ic=1,NpPM Xbb(ic)= Xpar(ic) Ybb(ic)= Ypar(ic) Zbb(ic)= Zpar(ic) VXbb(ic)= VX(ic) VYbb(ic)= VY(ic) VZbb(ic)= VZ(ic) End Do myMemory= Memory(-6_8*Nparticles) DEALLOCATE(Xpar,Ypar,Zpar,VX,VY,VZ) !--- restore data structur, rescale Coords and Veloc myMemory= Memory(6_8*NpPM) ALLOCATE(Xpar(NpPM),Ypar(NpPM),Zpar(NpPM),VX(NpPM),VY(NpPM),VZ(NpPM)) write(*,*) ' Copy to old arrays ',Xscale,Vscale !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic) Do ic=1,NpPM Xpar(ic) =Xbb(ic) Ypar(ic) =Ybb(ic) Zpar(ic) =Zbb(ic) VX(ic) =VXbb(ic) VY(ic) =VYbb(ic) VZ(ic) =VZbb(ic) End Do DEALLOCATE(Xbb,Ybb,Zbb) DEALLOCATE(VXbb,VYbb,VZbb) myMemory= Memory(-6_8*Nparticles) Nparticles = NpPM Np = Nparticles write(*,*) ' Rescale ',Nparticles !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ic) Do ic=1,Nparticles Xpar(ic) = Xpar(ic)/Xscale+1. Ypar(ic) = Ypar(ic)/Xscale+1. Zpar(ic) = Zpar(ic)/Xscale+1. If(Xpar(ic).ge.NGRID)Xpar(ic) =Xpar(ic)-NGRID If(Xpar(ic).lt.1) Xpar(ic) =Xpar(ic)+NGRID If(Ypar(ic).ge.NGRID)Ypar(ic) =Ypar(ic)-NGRID If(Ypar(ic).lt.1) Ypar(ic) =Ypar(ic)+NGRID If(Zpar(ic).ge.NGRID)Zpar(ic) =Zpar(ic)-NGRID If(Zpar(ic).lt.1) Zpar(ic) =Zpar(ic)+NGRID VX(ic) = VX(ic)/Vscale VY(ic) = VY(ic)/Vscale VZ(ic) = VZ(ic)/Vscale End Do end Subroutine RemoveBuffer end Module LinkerList