C ______________________________ Selects and Scales PM particles C For setting refinement in initial conditions C December 1998 A. Klypin (aklypin@nmsu.edu) C Input: 1) catalog of halos, which define what regions should be C selected for mass refinement C 2) PMcrd.DAT and PMcrs0.DAT -- C lowest resolution run with one specie c of particles covering all the computational Box C C Output: File with lagrangian coordinates of selected particles. C These coordinates are used by PMstartM to make mass refinement C Format of halo catalog: C First 16 lines -- any text (ignored) C Next line(s): one line per halo (=center): C x,y,z,vx,vy,vz,Mass,Radius,v1,v2 C Only coordinates (x,y,z) and Radius are used. C (x,y,z) = center of sphere of radius Radius C Units: (x,y,z) are in Mpc/h; Radius is in kpc/h C Within the sphere(s) all particles are selected. In the lagrangian C space occupied by the particles mass resolution will be increased. C Spheres may overlap in any possible way INCLUDE 'PMparameters.h' REAL INPUT COMMON /TABLE1/Gtable(3,NROW,NROW,NROW) COMMON /TABLE2/Itable(NROW,NROW,NROW) Character*70 FileASCII,Line,FileCatalog C................................................................... C 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') .... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE ( k,j,i) Do k=1,NROW Do j=1,NROW Do i=1,NROW Gtable(1,i,j,k) =0. Gtable(2,i,j,k) =0. Gtable(3,i,j,k) =0. Itable(i,j,k) =0 EndDo EndDo EndDo CALL RDTAPE CALL SetWeight write (*,*) ' RDTAPE is done' WRITE (17,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + 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=',E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3, + 1x,' Omega_curvature=',F7.3,/ + 1x,' 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 =INPUT(' 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 Do i=1,17 ! skip header lines read(18,'(A)') Line write(*,'(A)') Line EndDo 200 read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr, + Vvx0,Vvy0,Vvz0,Am,Radius,Vrms,Vcirc c 200 read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr, c + Am,Vcir,Radius Radius =Radius/1000. ! Scale radius to Mpc/h rad2 = Radius**2 If(Radius.gt.0.)Then ind = 1 else ind = 0 EndIf xxmax =-1.e+9 xxmin = 1.e+9 vmax =-1.e+9 vmin = 1.e+9 xrand =0. 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 c write (*,*)' Move page=',IROW,' file=',ifile,' N=',In_page iL = NPAGE*(IROW-1) CALL GETROW(IROW,ifile) ! read in a page of particles 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.) ! periodical conditions dX =min(ABS(X-Xcentr),Box-ABS(X-Xcentr)) dY =min(ABS(Y-Ycentr),Box-ABS(Y-Ycentr)) dZ =min(ABS(Z-Zcentr),Box-ABS(Z-Zcentr)) dd =dx**2 +dY**2 +dZ**2 IF(dd.lt.rad2)Then ! get the particle Ipart =IN+iL ! current particle number WPAR =iWeight(Ipart) ! particles weight 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 c 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 Gtable(1,i1,j1,k1) =X Gtable(2,i1,j1,k1) =Y Gtable(3,i1,j1,k1) =Z Itable(i1,j1,k1) = ind c write (17,250) Ipart,X,Y,Z,Vxs,Vys,Vzs,i1,j1,k1 c & ! write to a file EndIf 250 Format(i9,3F9.4,3F9.2,3i5,3i9) ENDDO ENDDO Do k=1,NROW Do j=1,NROW Do i=1,NROW If(Itable(i,j,k).eq.1)write (17,250) & 0,(Gtable(m,i,j,k),m=1,3),0., 0.,0.,i,j,k EndDo EndDo EndDo write (*,*) ' Scaled Coordinates were written to:', & FileASCII write (*,*) ' Number of particles =',Icount write (*,*) ' Center: ', Xcentr,Ycentr,Zcentr,' R=',Radius write (*,*) ' Min/Max of coordinates(Mpc/h)=',xxmin,xxmax write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax goto 200 ! get next center 300 write (*,*) ' <======= the End ========> ' END