C ______________________________ Convert PM format to ASCII C November 1997 A. Klypin (aklypin@nmsu.edu) C C Scales internal PM coordinates and velocities C to coordinates in Mpc/h and km/s C In PM the coordinates are in the range 1 - (NGRID+1) C velocities are P = a_expansion*V_pec/(x_0H_0) C where x_0 = comoving cell_size=Box/Ngrid C H_0 = Hubble at z=0 C C NROW = number of particles in 1D C NGRID= number of cells in 1D MODULE STRUCTURES parameter ( ngrid = 128 ) ! # of 0 lv cells in 1d parameter ( nrow = 8192 ) ! # of particles in 1d parameter ( nspecies = 7 ) ! # of particle species parameter ( nb = 100e6 ) ! # of particles in a buffer parameter ( npage = nrow**2 ) ! # of particles in a page parameter ( nrecl = npage * 6 ) ! length of particle row in words parameter ( nbyteword = 4 ) ! defines length of direct-access record logical :: inside Character :: FileASCII*50,fname*120 Integer :: NACCES ! integer*8 variables INTEGER*8 :: JPAGE,N_particles Integer*4 lspecies(nspecies) ! read buffers COMMON / RUN / AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROWC,NGRIDC,Nspecs,Nseed,Om0,Oml0,hubble,Wp5, & Ocurv,extras(100) DIMENSION wspecies(nspecies) EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) character*45 HEADER common / HEADDR/ HEADER COMMON /ROW/ xpar(NPAGE),ypar(NPAGE),zpar(NPAGE), & vxx(NPAGE),vyy(NPAGE),vzz(NPAGE) REAL, ALLOCATABLE :: Xc(:,:),Yc(:,:),Zc(:,:), & VXc(:,:),VYc(:,:),VZc(:,:) Integer, ALLOCATABLE :: Np_th(:) !$OMP THREADPRIVATE (/ROW/) ! make copies for every thread End Module STRUCTURES !----------------------------------------------------------------- PROGRAM PM_to_ASCIIm USE STRUCTURES REAL :: INPUT character :: moment*120 Integer :: OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS Integer*8 :: Ipart,iL JPAGE = NPAGE NACCES= NRECL*4 / nbyteword C C Read data and open files WRITE (*,'(A,$)') 'Enter Name for output ASCII file = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)')'Enter snapshot (0-for standard files)= ' READ (*,*) moment If(TRIM(moment)=='0'.or.TRIM(moment)=='')Then Open(4,file ='PMcrd.DAT',form ='UNFORMATTED',status ='UNKNOWN') Else write(fname,'(a,a,a)')'PMcrd.',TRIM(moment),'.DAT' Open (4,file =fname,form ='UNFORMATTED',status ='UNKNOWN') EndIf CALL RDdata Jfiles =1 If(lspecies(Nspecies)> 1024**3)Jfiles=8 If(TRIM(moment)=='0'.or.TRIM(moment)=='')Then Do i =1,Jfiles write(fname,'(a,i1,a)')'PMcrs',i-1,'.DAT' OPEN(UNIT=20+i,FILE=TRIM(fname),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo Else Open (4,file =fname,form ='UNFORMATTED',status ='UNKNOWN') Do i =1,Jfiles write(fname,'(a,i1,a,a,a)')'PMcrs',i-1,'.',TRIM(moment),'.DAT' OPEN(UNIT=20+i,FILE=TRIM(fname),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo EndIf write (*,*) ' RDTAPE is done' Fraction =INPUT(' Enter fraction of particles =') If(Fraction.le.0.)Stop Ifraction = INT(1./Fraction+0.01) Radius =INPUT(' Enter radius of sphere or zero for all in box =') Rad2 = 0. If(Radius.gt.0.)Then write (*,'(A,$)') ' Enter Center: x,y,z =' read (*,*) xcen,ycen,zcen xl =xcen -Radius xr =xcen +Radius yl =ycen -Radius yr =ycen +Radius zl =zcen -Radius zr =zcen +Radius Rad2 = Radius**2 EndIf xrand = 0. ncount = 0 write (*,*) ' Every ',Ifraction,' particle will be selected' WRITE (17,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + xcen,ycen,zcen,Fraction,Radius 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=',I9,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Center=',3F9.3,' Fraction =',F7.3,' Radius=',F7.3) If(extras(100).ne.0.)Then Box = extras(100) else Box =INPUT(' Enter box size in Mpc/h =') endif write (*,*) ' Box size is = ',Box,'Mpc/h' BoxV =Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities ScaleC = Box/NGRID ! scale factor for Coordinates N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -JPAGE*(Npages-1) write (*,*) ' Pages =',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last ! Allocate buffers for reading particles nthr = OMP_GET_MAX_THREADS() write(*,*) ' Number of threads=',nthr ALLOCATE (Xc(nb,nthr),Yc(nb,nthr),Zc(nb,nthr)) ALLOCATE (VXc(nb,nthr),VYc(nb,nthr),VZc(nb,nthr)) ALLOCATE (Np_th(nthr)) Np_th = 0 ! initialize counters for each thread xmax =-1.e10; ymax =-1.e10; zmax =-1.e10 vxmax =-1.e10; vymax=-1.e10; vzmax=-1.e10 xmin =1.e10 ; ymin =1.e10 ; zmin =1.e10 vxmin =1.e10 ; vymin=1.e10 ; vzmin=1.e10 C Loop over particles !$OMP PARALLEL DO DEFAULT(SHARED), !$OMP+ PRIVATE (IROW,In_page,iL,in,j,i_n, !$OMP+ X,Y,Z,Vxs,Vys,Vzs,Nseed,inside,xrand,rr), !$OMP+ REDUCTION(MAX:xmax,ymax,zmax), !$OMP+ REDUCTION(MIN:xmin,ymin,zmin), !$OMP+ REDUCTION(MAX:vxmax,vymax,vzmax), !$OMP+ REDUCTION(MIN:vxmin,vymin,vzmin), !$OMP+ SCHEDULE(DYNAMIC,4) DO IROW=1,Npages i_n = OMP_GET_THREAD_NUM()+1 In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last write (*,'(" Read page=",2i5,i10)')IROW,i_n,Np_th(i_n) iL = JPAGE*(IROW-1) If(lspecies(Nspecies)>1024**3)Then j = (IROW-1)/64 +1 if(j>8.or.j<1)write(*,*) ' Error in file number: ',j read(20+j,REC=mod(IROW-1,64)+1)xpar,ypar,zpar,vxx,vyy,vzz Else read(21, rec = IROW ) xpar,ypar,zpar,vxx,vyy,vzz EndIf DO IN=1,In_page X =ScaleC* (XPAR(IN)-1.) Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) Vxs=ScaleV* VXx(IN) Vys=ScaleV* VYy(IN) Vzs=ScaleV* VZz(IN) xmax =MAX(xmax,X) ymax =MAX(ymax,Y) zmax =MAX(zmax,Z) vxmax =MAX(vxmax,Vxs) vymax =MAX(vymax,Vys) vzmax =MAX(vzmax,Vzs) xmin =MIN(xmin,X) ymin =MIN(ymin,Y) zmin =MIN(zmin,Z) vxmin =MIN(vxmin,Vxs) vymin =MIN(vymin,Vys) vzmin =MIN(vzmin,Vzs) Ipart =IN+iL ! current particle number ! Do i =1,Nspecies ! if(Ipart.le.lspecies(i))Then ! W =wspecies(i) ! goto 50 ! endif ! EndDo inside =.true. If(Radius.gt.0.)Then if(X.lt.xl.or.X.gt.xr)inside =.false. if(Y.lt.yl.or.Y.gt.yr)inside =.false. if(Z.lt.zl.or.Z.gt.zr)inside =.false. EndIf If(inside)Then If(Fraction.lt.0.999)xrand =RANDd(Nseed) ! random fraction If(xrand.le.Fraction)Then !rr = (X-xc)**2+(Y-yc)**2+(Z-zc)**2 !If(rr.le.Rad2)Then !write(*,*)xrand,Fraction,inside !write(*,'(25x,4g12.6,i12,i4)')Z,zl,zr,Radius,Ipart,i_n Np_th(i_n) = Np_th(i_n) +1 If(Np_th(i_n).gt.nb)Stop ' Not enough buffer length nb' Xc(Np_th(i_n),i_n) = X Yc(Np_th(i_n),i_n) = Y Zc(Np_th(i_n),i_n) = Z VXc(Np_th(i_n),i_n) = VXs VYc(Np_th(i_n),i_n) = VYs VZc(Np_th(i_n),i_n) = VZs EndIf EndIf ENDDO ENDDO Icount = SUM(Np_th) Do i=1,nthr ! loop over all threads If(Np_th(i).gt.0)Then ! if not empty, dump particles Do j=1, Np_th(i) write (17,'(3f10.5,3f10.2)') Xc(j,i),Yc(j,i),Zc(j,i), & VXc(j,i),VYc(j,i),VZc(j,i) EndDo EndIf EndDo write (*,*) ' Scaled Coordinates were written to:', & FileASCII write (*,*) ' Number of particles total selected =',Icount write (*,*) ' on different threads=',Np_th write (*,*) ' Min/Max of x coordinates(Mpc/h)=',xmin,xmax write (*,*) ' y coordinates =',ymin,ymax write (*,*) ' z coordinates =',zmin,zmax write (*,*) ' Min/Max of x velocities (km/s) =',vxmin,vxmax write (*,*) ' y velocities =',vymin,vymax write (*,*) ' z velocities =',vzmin,vzmax END PROGRAM PM_to_ASCIIm C--------------------------------------------------- C Read current data from disk/tape, C Open files C Nrecl is the number of values in a record C Npage is the number of particles in a page SUBROUTINE RDdata C--------------------------------------------------- USE STRUCTURES READ (4,err=20,end=20) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecs,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,NGRID (PMparamters.h)=',NROW,NGRID write (*,*) ' NROW,NGRID (PMcrd.DAT) =',NROWC,NGRIDC ENDIF IF(NGRID.NE.NGRIDC) THEN WRITE (*,*) + ' NGRID in PARAMETER and in TAPE-FILE are different:' write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC ENDIF write(*,*) ' Number of particles =',lspecies(Nspecies) CLOSE (4) RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END SUBROUTINE RDdata C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End FUNCTION INPUT C------------------------------------------------ C random number generator FUNCTION RANDd(M) C------------------------------------------------ DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4 + /453815927,2147483648.,2147483647,536870912,131072,256, + 16777216,4,16384,8388608,128/ ML=M/K1*K1 M1=(M-ML)*L1 ML=M/K2*K2 M2=(M-ML)*L2 ML=M/K3*K3 M3=(M-ML)*L3 ML=M/K4*K4 M4=(M-ML)*L4 M5=KI-M IF(M1.GE.M5)M1=M1-KI-1 ML=M+M1 M5=KI-ML IF(M2.GE.M5)M2=M2-KI-1 ML=ML+M2 M5=KI-ML IF(M3.GE.M5)M3=M3-KI-1 ML=ML+M3 M5=KI-ML IF(M4.GE.M5)M4=M4-KI-1 ML=ML+M4 M5=KI-ML IF(LC.GE.M5)ML=ML-KI-1 M=ML+LC RANDd=M/AM RETURN END FUNCTION Randd