!------------------------------------------------- ! ! MODULE setHaloParticles PUBLIC Integer, PARAMETER :: NROW =1024, & NGRID =256, & Nmaxpart =512**3, & nbyteword = 4, & Nc =200000, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6, & Nproc = 4, & Nrad =120 Real :: INPUT,Mhalo,Ovlim,Ovcnt Real, DIMENSION(Nrad) :: Overdens,Mass,MnMass, & Radsc Real, DIMENSION(Nmaxpart) :: Xp,Yp,Zp,VXp,VYp,VZp,iWeight Real, DIMENSION(1) :: Xpp,Ypp,Zpp,VXpp,VYpp,VZpp,iWeightp 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 DIMENSION RECDAT(NRECL),wspecies(10),lspecies(10) Real, ALLOCATABLE, DIMENSION(:,:) :: RECDATp EQUIVALENCE (RECDAT(1),XPAR(1)) & ,(wspecies(1),extras(1)), & (lspecies(1),extras(11)) Contains !--------------------------------------------------------- ! Read particles SUBROUTINE ReadPnt(N) !-------------------------------------------------------------- N =0 Iz=1 Iy=1 Ff=1. xshift = 0. yshift = 0. zshift = 0. xmin =1.e10 xmax =-1.e10 ymin =1.e10 ymax =-1.e10 zmin =1.e10 zmax =-1.e10 Nsmall =0 Lsmall =lspecies(1) If(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(*,'(2(a,i9))')' Npage=', NPAGE,' Npages=',Npages write(*,'(2(a,i9))')' Np=',N_particles,' Nspecies=',Nspecies DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = NPAGE*(IROW-1) write (*,*) ' page=',IROW,' Np=',N ! Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 Xpp(N) =XPAR(IN)-1. Ypp(N) =YPAR(IN)-1. Zpp(N) =ZPAR(IN)-1. If(Xpp(N).gt.NGRID)Xpp(N) =Xpp(N)-NGRID If(Ypp(N).gt.NGRID)Ypp(N) =Ypp(N)-NGRID If(Zpp(N).gt.NGRID)Zpp(N) =Zpp(N)-NGRID If(Xpp(N).le.0.)Xpp(N) =Xpp(N)+NGRID If(Ypp(N).le.0.)Ypp(N) =Ypp(N)+NGRID If(Zpp(N).le.0.)Zpp(N) =Zpp(N)+NGRID If(Xpp(N).ge.NGRID)Xpp(N) =Xpp(N)-1.e-4 If(Ypp(N).ge.NGRID)Ypp(N) =Ypp(N)-1.e-4 If(Zpp(N).ge.NGRID)Zpp(N) =Zpp(N)-1.e-4 if(Xpp(N).le.0.)write (*,'(" X:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N if(Ypp(N).le.0.)write (*,'(" Y:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N if(Zpp(N).le.0.)write (*,'(" Z:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N if(Xpp(N).ge.NGRID)write(*,'(" X:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N if(Ypp(N).ge.NGRID)write(*,'(" Y:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N if(Zpp(N).ge.NGRID)write(*,'(" Z:",3g11.4,i10)')Xpp(N),Ypp(N),Zpp(N),N If(N.le.Lsmall)Then Nsmall =Nsmall + 1 xmin =min(Xpp(N),xmin) xmax =max(Xpp(N),xmax) ymin =min(Ypp(N),ymin) ymax =max(Ypp(N),ymax) zmin =min(Zpp(N),zmin) zmax =max(Zpp(N),zmax) EndIf VXpp(N)=VX(IN) VYpp(N)=VY(IN) VZpp(N)=VZ(IN) ENDDO ENDDO write(*,*) ' Small Particles:',Nsmall write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write (*,*) ' Nmaxpart=',Nmaxpart,Np do i=1,10000,200 write (*,'(i8," Velocities=",3g11.4)')i,VXpp(i),VYpp(i),VZpp(i) EndDo Return End SUBROUTINE ReadPnt !--------------------------------------------------------- ! Read particles SUBROUTINE ReadPntP(N) !-------------------------------------------------------------- REAL :: xmin =1.e10,xmax =-1.e10,ymin =1.e10,ymax =-1.e10, & zmin =1.e10,zmax =-1.e10 ALLOCATE (RECDATp(NRECL,Nproc)) Nsmall =0 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 DO IROW=1, Npages, Nproc ! Loop over particle pages !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (jROW,iproc,iL,j,IN) Do iproc =1,Nproc jROW = IROW+iproc-1 In_page =NPAGE If(jROW.eq.Npages)In_page =N_in_last iL = NPAGE*(jROW-1) ! offset for particle number ! write (*,*) ' page=',jROW,' Np=',iL,Npages ! Loop over particles If(jROW.le.Npages)Then READ (21,REC=jROW) (RECDATp(j,iproc),j=1,NRECL) DO IN=1,In_page Xp(iL+IN) =RECDATp(IN,iproc) Yp(iL+IN) =RECDATp(IN+NPAGE,iproc) Zp(iL+IN) =RECDATp(IN+2*NPAGE,iproc) VXp(iL+IN) =RECDATp(IN+3*NPAGE,iproc) VYp(iL+IN) =RECDATp(IN+4*NPAGE,iproc) VZp(iL+IN) =RECDATp(IN+5*NPAGE,iproc) ! write(*,'(5i9,7g12.4)')IROW,iproc,jROW,iL,IN,Xp(iL+IN), & ! Yp(iL+IN), Zp(iL+IN) EndDo EndIf ENDDO EndDo N = lspecies(Nspecies) write(*,*) ' Done' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP REDUCTION(MAX:xmax,ymax,zmax) & !$OMP REDUCTION(MIN:xmin,ymin,zmin) Do i =1,N Xp(i) =Xp(i)-1. Yp(i) =Yp(i)-1. Zp(i) =Zp(i)-1. If(Xp(i).gt.NGRID)Xp(i) =Xp(i)-NGRID If(Yp(i).gt.NGRID)Yp(i) =Yp(i)-NGRID If(Zp(i).gt.NGRID)Zp(i) =Zp(i)-NGRID If(Xp(i).le.0.)Xp(i) =Xp(i)+NGRID If(Yp(i).le.0.)Yp(i) =Yp(i)+NGRID If(Zp(i).le.0.)Zp(i) =Zp(i)+NGRID If(Xp(i).ge.NGRID)Xp(i) =Xp(i)-1.e-4 If(Yp(i).ge.NGRID)Yp(i) =Yp(i)-1.e-4 If(Zp(i).ge.NGRID)Zp(i) =Zp(i)-1.e-4 if(Xp(i).le.0.)write (*,'(" X:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i if(Yp(i).le.0.)write (*,'(" Y:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i if(Zp(i).le.0.)write (*,'(" Z:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i if(Xp(i).ge.NGRID)write(*,'(" X:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i if(Yp(i).ge.NGRID)write(*,'(" Y:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i if(Zp(i).ge.NGRID)write(*,'(" Z:",3g11.4,i10)')Xp(i),Yp(i),Zp(i),i If(i.le.Lsmall)Then xmin =min(Xp(i),xmin) xmax =max(Xp(i),xmax) ymin =min(Yp(i),ymin) ymax =max(Yp(i),ymax) zmin =min(Zp(i),zmin) zmax =max(Zp(i),zmax) EndIf ENDDO write(*,*) ' Small Particles:',lspecies(1) write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write (*,*) ' Nmaxpart=',Nmaxpart,Np ! do i=1,10000,200 ! write (*,'(i8," Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i) ! EndDo DEALLOCATE (RECDATp) Return End SUBROUTINE ReadPntP !--------------------------------------------------- ! 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 !--------------------------------------------------- ! Open file on a tape OPEN(UNIT=9,FILE='PMcrd.DAT', & 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 (16,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='PMcrs0.DAT',ACCESS='DIRECT', & STATUS='UNKNOWN',RECL=NACCES) REWIND 9 write (*,*) ' done openning PMcrs0. Nspecies=',Nspecies write (*,*) ' Nparticles=',lspecies(Nspecies) 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 CALL RDTAPE CALL CPU_TIME(t0) ! CALL ReadPnt(N) CALL CPU_TIME(t1) CALL ReadPntP(N) CALL CPU_TIME(t2) write (*,'(a,3g12.4)') ' cpu=',t1-t0,t2-t0 !Do i=1,N ! if(abs(Xpp(i)-Xp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' x:',i,Xpp(i),Xp(i) ! if(abs(Ypp(i)-Yp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' y:',i,Ypp(i),Yp(i) ! if(abs(Zpp(i)-Zp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' z:',i,Zpp(i),Zp(i) ! if(abs(VXpp(i)-VXp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' vx:',i,VXpp(i),VXp(i) ! if(abs(VYpp(i)-VYp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' vy:',i,VYpp(i),VYp(i) ! if(abs(VZpp(i)-VZp(i)).gt.1.e-5) & ! write(*,'(a,i9,9g12.5)') ' vz:',i,VZpp(i),VZp(i) !EndDo End PROGRAM HaloParticles