! ______________________________Find coordinates of particles ! provided by a list of particle IDs ! April 2008 ! A. Klypin (NMSU) ! Module Structures Integer, PARAMETER :: & NROW =1024, & ! maximum number of particles in 1D NGRID =256, & ! zero-level mesh in 1D ! must be equal to ng in ART Nmaxpart = 1024**3, & ! max number of particles for this run ! Nmaxpart must be less than NROW**3 nbyteword = 4, & ! defines length of direct-access:1or4 NPAGE = NROW**2, & ! number of particles in a record NMAX=NGRID/2, & NRECL= NPAGE*6 ! number of words in a record REAL, PARAMETER :: pi=3.14159265 Real :: Box Integer :: List(Nmaxpart) 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) CHARACTER*45 HEADER COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) DIMENSION RECDAT(NRECL),wspecies(10),lspecies(10) EQUIVALENCE (RECDAT(1),XPAR(1)) & ,(wspecies(1),extras(1)), & (lspecies(1),extras(11)) Contains !-------------------------------------------------------------------------- ! Read particles, find those, which are ! in the list SUBROUTINE ListMaxima !--------------------------------------- Character*120 Line Lines =13 ! skip this many lines of the header !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1) DO M1=1,Nmaxpart List(M1) = 0 END DO Do i=1,Lines read(17,'(a)')Line !write(16,'(a)')Line EndDo iList =0 Do read(17,*,iostat=iflag)np,x,y,z if(iflag /= 0)exit List(np) =1 iList =iList +1 EndDo write(*,*) ' Total lines read=',iList End Subroutine ListMaxima !-------------------------------------------------------------------------- ! Read particles, find those, which are ! in the list SUBROUTINE PartCoordinates !--------------------------------------- N =0 xmin =1.e+5 xmax =-1.e+5 ymin =1.e+5 ymax =-1.e+5 zmin =1.e+5 zmax =-1.e+5 Scale =Box/NGRID Vscale= 100.*hubble*Scale/AEXPN ! Scale for velocities ! Loop over particles 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),a,g12.5)') ' Nspecies=',Nspecies, & ' Nparticles=',lspecies(Nspecies), & ' Last weight=',wspecies(Nspecies) Nsp_current =1 WPAR = wspecies(1) ! *8. Lspec_next=lspecies(1)+1 DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last If(IROW/50*50.eq.IROW) & write (*,*) ' Read page=',IROW,' N=',N, & ' Npages=',Npages,In_page iL = NPAGE*(IROW-1) ! Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 If(N.ge.Lspec_next)Then ! next specie write (*,*) ' Another specie. Current=',Nsp_current write (*,*) ' Particle=',N,WPAR, & ' lspecies=',lspecies(Nsp_current), & ' next=',lspecies(Nsp_current+1) Nsp_current=Nsp_current+1 Lspec_next =lspecies(Nsp_current)+1 WPAR = wspecies(Nsp_current) EndIf If(List(N) ==1)Then ! found particle in the List X=(XPAR(IN)-1.)*Scale Y=(YPAR(IN)-1.)*Scale Z=(ZPAR(IN)-1.)*Scale xmin =MIN(xmin,X) ymin =MIN(ymin,Y) zmin =MIN(zmin,Z) xmax =MAX(xmax,X) ymax =MAX(ymax,Y) zmax =MAX(zmax,Z) write(16,'(i10,3f11.5,2x,3g12.4)')N,X,Y,Z, & VX(IN)*Vscale,VY(IN)*Vscale,VZ(IN)*Vscale EndIf ENDDO EndDo write (*,*) ' Number of particles read=', N write (*,*) ' X:',xmin,xmax write (*,*) ' Y:',ymin,ymax write (*,*) ' Z:',zmin,zmax RETURN END Subroutine PartCoordinates !-------------------------------------------------------------------------------------------- ! 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) Close (9) write (*,*) ' done openning PMcrd ----' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END Subroutine RDTAPE !-------------------------------------------------------------------------------------------- Real FUNCTION Input(Line) character(*) :: Line write(*,'(a,$)')Line read(*,*) Input End Function Input End Module Structures !-------------------------------------------------------------------------------------------------- ! Program MaximaPosition use Structures ! ! Read data and open files Open(16,file='ListMaximaNew.dat') Open(17,file='ListMaxima.dat') CALL RDTAPE write (*,*) ' RDTAPE is done' Box = extras(100) If(Box.le.1.e-3)Then Box =INPUT(' Enter box size in Mpc/h =') EndIf ! Call ListMaxima ! read list of particles at dens maxima CALL PartCoordinates ! find their current coordinates and Vs STOP END Program MaximaPosition