! ______________________________ Selects and Scales PM particles ! For setting refinement in initial conditions ! December 1998 A. Klypin (aklypin@nmsu.edu) ! Input: 1) catalog of halos, which define what regions should be ! selected for mass refinement ! 2) PMcrd.DAT and PMcrs0.DAT -- ! lowest resolution run with one specie ! of particles covering all the computational Box ! ! Output: File with lagrangian coordinates of selected particles. ! These coordinates are used by PMstartM to make mass refinement ! Format of halo catalog: ! First 16 lines -- any text (ignored) ! Next line(s): one line per halo (=center): ! x,y,z,vx,vy,vz,Mass,Radius,v1,v2 ! Only coordinates (x,y,z) and Radius are used. ! (x,y,z) = center of sphere of radius Radius ! Units: (x,y,z) are in Mpc/h; Radius is in kpc/h ! Within the sphere(s) all particles are selected. In the lagrangian ! space occupied by the particles mass resolution will be increased. ! Spheres may overlap in any possible way Module SetData PARAMETER (NROW =8192) ! maximum number of particles in 1D PARAMETER (NGRID =256) ! zero-level mesh in 1D ! must be equal to ng in ART ! Integer 4 <-> 8 PARAMETER (NhalosMax = 10000) ! maximum number of halos Integer*8,PARAMETER ::Nmaxpart=512**3 ! max number of particles for this run ! Nmaxpart must be less than NROW**3 Integer*8 lspecies(10) Integer :: Nhalos PARAMETER (nbyteword = 4) ! defines length of direct-access:1or4 PARAMETER (NPAGE = min(1024**2,NROW**2) ) ! number of particles in a record PARAMETER (NMAX=NGRID/2) PARAMETER (NRECL= NPAGE*6) ! number of words in a record 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 COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) DIMENSION RECDAT(NRECL),wspecies(10) EQUIVALENCE (RECDAT(1),XPAR(1)) & ,(wspecies(1),extras(1)), & (lspecies(1),extras(11)) REAL INPUT Character*70 FileASCII,Line,FileCatalog real*4, DIMENSION(NhalosMax) :: Xcentr,Ycentr,Zcentr,Rcentr CONTAINS !--------------------------------------------------- ! 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' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END SUBROUTINE RDTAPE !---------------------------------- Read in variables REAL FUNCTION GET(text) !------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x GET =x Return End FUNCTION GET end Module SetData !---------------------------------------------------------------------------------------------------- Program LagrangianCoordinates use SetData ! Read data and open files WRITE (*,'(A,$)') ' Enter Output File Name = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') ' Enter File Name for Input Halo Catalog = ' READ (*,'(A)') FileCatalog OPEN(18,FILE=FileCatalog,STATUS='UNKNOWN') CALL RDTAPE write (*,*) ' RDTAPE is done' WRITE (17,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EKIN, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv 100 FORMAT(' Header=>',A45,/, & ' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/, & ' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/, & ' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I11,/ , & ' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3, & ' Omega_curvature=',F7.3,/, & ' Npart Euler Coordinates (Mpc/h)', & ' Velocities (km/s)', & ' Lagrangian Coords (limits: 1-N_1d)') If(extras(100).gt.0.)Then Box =extras(100) write(*,*) ' Box size =',Box,' Mpc/h' Else Box =GET(' Enter box size in comoving Mpc/h =') EndIf IF(Nspecies.ne.1)Then write (*,*) '--- You can use this program only for runs with', & ' one specie of particles at the lowest mass resolution.' write (*,*)' --- The number of species is ',Nspecies,' STOP' Stop Endif BoxV = Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities ScaleC = Box/NGRID ! scale factor for Coordinates write(*,*) ' Box= ',Box write(*,*) ' Ngrid = ',NGRID,NGRIDC write(*,*) ' ScaleC= ',ScaleC Do i=1,6 ! skip header lines read(18,'(A)') Line write(*,'(A)') Line EndDo Nhalos = 0 factorVir = 3. ! increase region of refinement by this factor Do read(18,*,iostat=ierr ) Xcentr1,Ycentr1,Zcentr1, & Vvx0,Vvy0,Vvz0,dd,Am,Rcentr1,Vrms,Vcirc if(ierr/=0)exit Nhalos = Nhalos +1 If(Nhalos.gt.NhalosMax)stop '----- too many halos. Change NhalosMax ' ! 200 read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr, ! Am,Vcir,Rcentr Rcentr(Nhalos) =factorVir *Rcentr1/1000. ! Scale radius to Mpc/h Xcentr(Nhalos) = Xcentr1 Ycentr(Nhalos) = Ycentr1 Zcentr(Nhalos) = Zcentr1 write (*,*) ' Center: ', Xcentr1,Ycentr1,Zcentr1,' Rad=',Rcentr1 EndDo write(*,*) ' Number of halos read =',Nhalos if(Nhalos==0)stop '===> No halos is read' xxmax =-1.e+9 xxmin = 1.e+9 vmax =-1.e+9 vmin = 1.e+9 Icount =0 ifile =1 N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) N_line=INT((lspecies(1))**0.33333+0.5) N2 =N_line**2 write (*,*) ' Pages =',Npages,' Species =',Nspecies write (*,*) ' Nparticles =',N_particles,' N_in_last=',N_in_last write (*,*) ' Particles_1d =',N_line DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last write (*,*)' Move page=',IROW,' file=',ifile,' N=',In_page iL = NPAGE*(IROW-1) !CALL GETROW(IROW,ifile) ! read in a page of particles READ (20+Ifile,REC=IROW) RECDAT DO IN=1, In_page ! Loop over particles X =ScaleC* (XPAR(IN)-1.) ! scale coordinates Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) Ipart =IN+iL ! current particle number Vxs=ScaleV* VX(IN) ! scale velocities Vys=ScaleV* VY(IN) Vzs=ScaleV* VZ(IN) Icount =Icount +1 ! count particles xxmax =MAX(xxmax,X,Y,Z) xxmin =MIN(xxmin,X,Y,Z) vmax =MAX(vmax,Vxs,Vys,Vzs) vmin =MIN(vmin,Vxs,Vys,Vzs) ! find lagrangian coords ! Jpart =N_particles -Ipart+1 ! reverse the order of particles Jpart =Ipart ! direct order of particles k1 =INT((Jpart-1)/N2)+1 jj =Jpart-(k1-1)*N2 j1 =INT((jj-1)/N_line)+1 i1 =jj-(j1-1)*N_line if(Z<18..and.Z>10.)write(55,'(3f9.3,3x,i10,3i5)')X,Y,Z,Ipart,i1,j1,k1 DO j =1,Nhalos rad2 = Rcentr(j)**2 ! periodical conditions dX =min(ABS(X-Xcentr(j)),Box-ABS(X-Xcentr(j))) dY =min(ABS(Y-Ycentr(j)),Box-ABS(Y-Ycentr(j))) dZ =min(ABS(Z-Zcentr(j)),Box-ABS(Z-Zcentr(j))) dd =dx**2 +dY**2 +dZ**2 IF(dd.lt.rad2)Then ! get the particle write (17,'(i9,3F9.4,3F11.2,3i5,3f10.4)') Ipart,X,Y,Z,Vxs,Vys,Vzs,i1,j1,k1,Rcentr(j)!,sqrt(dd),sqrt(rad2) ! exit ! write to a file EndIf EndDO ENDDO ENDDO write (*,*) ' Scaled Coordinates were written to:', FileASCII write (*,*) ' Number of particles =',Icount write (*,*) ' Min/Max of coordinates(Mpc/h)=',xxmin,xxmax write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax 300 write (*,*) ' <======= the End ========> ' END Program LagrangianCoordinates