!-------------------------------------------------------------------------------- ! Read Files ! ! (nx,ny,nz) - number of boxes in each direction ! Rhalo = width of buffer in Mpch around each region Module SetArrs Integer, PARAMETER :: & NGRID = 256, & nbyteword = 4, & NPAGE = 1024**2 ! # particles in a record Integer*8,PARAMETER ::Np=2048_8**3 ! max number of particles Real :: AEXPN,ASTEP,PARTW,Om0,Oml0,hubble,NGRIDC Integer :: jStep,Nspecies,Nseed,ISTEP,Nx,Ny,Nz,iFormat CHARACTER*45 HEADER Real*4, DIMENSION(NPAGE) :: Xp,Yp,Zp,VXp,VYp,Vzp,Wpar Integer*8, DIMENSION(NPAGE) :: id_part Real :: Box,Rad,xc,yc,zc,Xscale, & Vscale,dR Integer*8 :: iCount Contains !-------------------------------------------------------------------------------- ! Define configuration of files ! Read control info from the first file ! allocate arrays ! Subroutine ReadConfig character*80 :: Name logical :: ext,found write(*,'(a,$)')' Enter snapshot number =' read(*,*)jStep ! read header information from the first existing found = .false. Do jFile = 1,1024 !write(Name,'(a,i4.4,a,i4.4,a)') 'PMss.',jFile,'.',jStep,'.DAT' write(Name,'(a,i4.4,a,i4.4,a)') 'PMss.',jStep,'.',jFile,'.DAT' Inquire(file=Name,exist =ext) If(ext)Then open(1,file=Name,form='unformatted') Read(1)HEADER Read(1)AEXPN,ASTEP,ISTEP,NROWC,NGRIDC,Nspecies, & Nseed,Om0,Oml0,hubble,Box Read(1) k,Nx,Ny,Nz,dR,nBuffer Read(1) xL,xR,yL,yR,zL,zR write(*,'(a)')HEADER write(*,'(a,f9.4,5x,a,2i5)')' aexpn=',AEXPN,' Step=',ISTEP,jStep if(ISTEP .ne. jStep)Stop ' -- wrong time step --' write(*,'(a,3i5)')' Configuration of files: in each direction=',Nx,Ny,Nz Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= Box*100./Ngrid/AEXPN ! Scale for comoving coordinates found = .true. iFormat = 0 read(1,iostat=ierr)iNfile read(1,iostat=ierr)Nrecord write(*,*) ' In this file N=',iNfile write(*,*) ' Nrecord=',Nrecord If(ierr/=0)Stop 'Error reading particle file in SetArrs ' Read(1,iostat=ierr) (Xp(m),m=1,Nrecord), & (Yp(m),m=1,Nrecord), & (Zp(m),m=1,Nrecord) Xscale= 1. ! Scale for comoving coordinates Vscale= 1. ! Scale for comoving coordinates write(*,*) ' io err =',ierr,' Must be zero for no errors' If(ierr/=0)Then backspace(1) read(1,iostat=ierr)(Xp(j), j=1,Nrecord) If(ierr/=0)Stop 'Error Xp reading' read(1,iostat=ierr)(Yp(j), j=1,Nrecord) If(ierr/=0)Stop 'Error Yp reading' read(1,iostat=ierr)(Zp(j), j=1,Nrecord) If(ierr/=0)Stop 'Error Zp reading' read(1,iostat=ierr)(VXp(j),j=1,Nrecord) If(ierr/=0)Stop 'Error VXp reading' read(1,iostat=ierr)(VYp(j),j=1,Nrecord) If(ierr/=0)Stop 'Error VYp reading' read(1,iostat=ierr)(VZp(j),j=1,Nrecord) If(ierr/=0)Stop 'Error VZp reading' read(1,iostat=ierr)(id_part(j),j=1,Nrecord) If(ierr/=0)Stop 'Error ID reading' read(1,iostat=ierr)(Wpar(j),j=1,Nrecord) If(ierr/=0)Stop 'Error Wpar reading' iFormat =1 END If exit End If End Do close(1) If(.not.found)Stop ' No files for this snapshot are found' write(*,*) ' You are reading iFormat =',iFormat write(Name,'(a,i4.4,a)') 'Part.',jStep,'.DAT' open(2,file=Name,form='formatted') write(2,'(a)')HEADER write(2,'(a,f9.4,5x,a,i5,a,f8.2)')' aexpn=',AEXPN,' Step=',ISTEP,' Box=',Box end Subroutine ReadConfig end Module SetArrs !-------------------------------------------------------------------------------- ! Program ReadFiles use SetArrs Logical :: inside character*120 :: Name Call ReadConfig ! set parameters ! write(*,'(/a,$)')' Enter Radius and center (radius <0 for all) =' read(*,*)Rad,xc,yc,zc write(Name,'(a,i4.4,4(a,i4.4),a)') 'Part.',jStep,'.', & INT(Rad),'.',INT(xc),'.',INT(yc),'.',INT(zc), & '.DAT' open(2,file=TRIM(Name),form='formatted') write(2,'(a)')HEADER write(2,'(a,f9.4,5x,a,i5,a,f8.2,a,4f8.2)')' aexpn=',AEXPN,' Step=',ISTEP,' Box=',Box, & ' Radius and center= ',Rad,xc,yc,zc If(Rad.gt.0.)Then xxL = (xc-Rad)/Xscale ; xxR = (xc+Rad)/Xscale ! limits in internal units yyL = (yc-Rad)/Xscale ; yyR = (yc+Rad)/Xscale zzL = (zc-Rad)/Xscale ; zzR = (zc+Rad)/Xscale xxL = (xc-Rad)/Xscale ; xxR = (xc+Rad)/Xscale ! limits in internal units yyL = (yc-Rad)/Xscale ; yyR = (yc+Rad)/Xscale zzL = (zc-Rad)/Xscale ; zzR = (zc+Rad)/Xscale Else xxL = 0. ; xxR = Box !Ngrid +1. ! limits in internal units yyL = 0. ; yyR = Box !Ngrid +1. zzL = 0. ; zzR = Box !Ngrid +1. End If write(*,'(a,1p,6g12.4)') ' Limits =',xxL,xxR,yyL,yyR,zzL,zzR write(*,'(/a,$)')' Enter fraction of particles to write =' read(*,*)fraction iSeed = 102191 ! seed for random numbers iCount = 0 ! total selected particles inside = .false. Do i0=1,Nx ! test whether the region is inside one domain Do j0=1,Ny Do k0=1,Nz node = i0+(j0-1)*Nx +(k0-1)*Nx*Ny Call ReadHeader(i0,j0,k0,xL,xR,yL,yR,zl,zR,iFlag) If(iFlag == 1)Then ! file exists inside = ((xR.ge.xxR).and.(xL.le.xxL)) inside = ((yR.ge.yyR).and.(yL.le.yyL)).and.inside inside = ((zR.ge.zzR).and.(zL.le.zzL)).and.inside If(inside)Then Call ReadParticles(xxL,xxR,yyL,yyR,zzl,zzR,iSeed,fraction) write(*,*) ' Total Number of particles =',iCount Stop ' All particles were in one file ' EndIf close(1) EndIf EndDo ! k EndDo ! j EndDo ! i Do k0=1,Nx ! read all files and test individually Do j0=1,Ny Do i0=1,Nz Call ReadHeader(i0,j0,k0,xL,xR,yL,yR,zl,zR,iFlag) If(iFlag == 1)Then ! file exists xL = xL + dR ; xR =xR -dR ! do not take particles in yL = yL + dR ; yR =yR -dR ! buffer around the region zL = zL + dR ; zR =zR -dR ! to avoid dublication inside = (min(xR,xxR)>max(xL,xxL)) inside = (min(yR,yyR)>max(yL,yyL)).and.inside inside = (min(zR,zzR)>max(zL,zzL)).and.inside If(inside)Then xL = max(xL,xxL) xR = min(xR,xxR) yL = max(yL,yyL) yR = min(yR,yyR) zL = max(zL,zzL) zR = min(zR,zzR) write(*,'(a,6f9.3,a,f9.2)') ' Limits: ',xL,xR,yL,yR,zl,zR,' Buffer width=',dR Call ReadParticles(xL,xR,yL,yR,zl,zR,iSeed,fraction) EndIf close(1) End If EndDo ! k EndDo ! j EndDo ! i write(*,*) ' Total Number of particles =',iCount Stop end Program ReadFiles !--------------------------------------------------------- ! Read particles SUBROUTINE ReadHeader(i0,j0,k0,xL,xR,yL,yR,zl,zR,iFlag) !-------------------------------------------------------------- use SetArrs character*80 :: Name logical :: ext,found node = i0+(j0-1)*Nx +(k0-1)*Nx*Ny iFlag = 0 !write(*,'(a,10i7)') ' node=',node,i0,j0,k0,Nx,Ny,Nz !write(Name,'(a,i4.4,a,i4.4,a)') 'PMss.',node,'.',jStep,'.DAT' write(Name,'(a,i4.4,a,i4.4,a)') 'PMss.',jStep,'.',node,'.DAT' !write(*,'(2a)') ' file=',Name Inquire(file=TRIM(Name),exist =ext) If(ext)Then open(1,file=TRIM(Name),form='unformatted') Read(1)HEADER Read(1)AEXPN,ASTEP,ISTEP,NROWC,NGRIDC,Nspecies, & Nseed,Om0,Oml0,hubble,Box Read(1) k,Nx,Ny,Nz,dR0 Read(1) xL,xR,yL,yR,zl,zR Read(1) iNfile write(*,'(5x,a,i5,3x,a,6f8.2,3x,a,i10)') 'node= ',node, & 'limits= ',xL,xR, & yL,yR, zL,zR, & 'Current selected particles =',iCount iFlag=1 EndIf end SUBROUTINE ReadHeader !--------------------------------------------------------- ! Read particles SUBROUTINE ReadParticles(xL,xR,yL,yR,zl,zR,iSeed,fraction) !-------------------------------------------------------------- use SetArrs Logical :: inside,takeit takeit = .true. !read(1)nn ! total number of particles in this file ! write(*,*) ' Number of particles in this node =',nn Do read(1,iostat=ierr)Nrecord ! number of particles in this record If(ierr /= 0) exit write(*,'(2(3x,a,i10))') 'Number of particles in this record =',Nrecord, & 'Count selected=',iCount If(iFormat == 0)Then Read(1,iostat=ierr) (Xp(m),m=1,Nrecord), & (Yp(m),m=1,Nrecord), & (Zp(m),m=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat =0 ' read(1,iostat=ierr) (VXp(m),m=1,Nrecord), & (VYp(m),m=1,Nrecord), & (VZp(m),m=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat =0 ' read(1,iostat=ierr) (Wpar(m),m=1,Nrecord), & (id_part(m),m=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat =0 ' Else read(1,iostat=ierr)(Xp(j), j=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat=1' read(1)(Yp(j), j=1,Nrecord) read(1)(Zp(j), j=1,Nrecord) read(1)(VXp(j),j=1,Nrecord) read(1)(VYp(j),j=1,Nrecord) read(1)(VZp(j),j=1,Nrecord) read(1,iostat=ierr)(id_part(j),j=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat=1' read(1,iostat=ierr)(Wpar(j),j=1,Nrecord) If(ierr/=0)stop ' Error reading Particles iFormat=1' End If Do i=1,Nrecord inside = ( xL <= Xp(i)).and. (Xp(i)