!---------------------------------------------------------------------------------------- ! Re-estimate Vcirc of halos ! ! Input: Halo Catalog (short list) ! Dark matter particles ! Output: Halo catalog ! MODULE setHaloParticles PUBLIC Integer, PARAMETER :: & NROW = 4096, & NGRID = 256, & nbyteword = 4, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6, & Nbin = 10, & ! # radial bins Nbf = 10, & ! width of buffer for linker list Nlinker = 300, & ! linker-list size kDynamic = 1000 ! OMP chunk for parallelization Integer*8, PARAMETER :: & Nmaxpart = 2048_8**3*1.3 ! max # of particles Real, PARAMETER :: & Unbound = 2., & ! Limit for removing unbound particles ! V = Vunbound*Vcirc dlog =1./float(Nbin), & ! log bin size FracSearch =5. ! defines cell-size for particles linker list: ! cell = Rdmax/FracSearch Real :: INPUT,Mhalo,Ovlim,Ovcnt,Rdmax Real, DIMENSION(Nmaxpart) :: Xp,Yp,Zp,VXp,VYp,Vzp,iWeight Real :: Box=0.,Xscale=0.,Vscale=0.,Dscale=0.,Cell=0., weightSmall Integer*8 :: Np=0,Ncp=0 Integer :: Ncenter=0, & Nmx,Nbx,Nmy,Nby,Nmz,Nbz,Nbin_max,Nbin_min COMMON / ROW / XPAR(NPAGE),YPAR(NPAGE),ZPAR(NPAGE), & VX(NPAGE),VY(NPAGE),VZ(NPAGE) COMMON / CONTROL/ AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EK,EK1,EK2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5, & Ocurv,extras(100) COMMON / HEADDR/ HEADER CHARACTER*45 HEADER Real :: RECDAT(NRECL),wspecies(10) Integer*8 :: lspecies(10) Integer*8, ALLOCATABLE, DIMENSION (:,:,:) :: Label Integer*4, ALLOCATABLE, DIMENSION (:,:) :: nPr(:,:) Integer*8, ALLOCATABLE, DIMENSION (:) :: Lst Real*4, ALLOCATABLE, DIMENSION (:,:) :: aMr,Vvrms,Rbin Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc, & Vcirc2,Amh EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/) Contains !------------------------------------------------------------------------- ! Make linker lists of particles in each cell ! Start with sorting particles by cells: ! particles in the same cell are in contiguous memory ! Need additional memory for particles: Dummy, iDummy arrays SUBROUTINE List !------------------------------------------------------------------------ !Real, ALLOCATABLE, DIMENSION (:) :: Dummy !Integer*8, ALLOCATABLE, DIMENSION (:) :: iDummy Integer*8 :: i,jp,nn !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncp Lst(i)= -1 EndDo write(*,*) ' Lst initialized' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx Label(i,j,k)=0 EndDo EndDo EndDo write(*,*) ' Label initialized' !!$OMP PARALLEL DO DEFAULT(SHARED) & !!$OMP PRIVATE (jp,i,j,k,k1) ! Do k1 =Nmz,Nbz ! write(*,*)' k1=',k1 Do jp=1,Ncp k=INT(Zp(jp)/Cell) k=MIN(MAX(Nmz,k),Nbz) ! If(k.eq.k1)Then i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp ! EndIf EndDo ! EndDo End SUBROUTINE List !------------------------------------------------------------------------ ! Get information on a halo iHalo ! Input: Xc,Yc,Zc - center; ! iHalo = number of halo ! Output: updated halo statistics SUBROUTINE SNeib(iHalo) !----------------------------------------------------------------------- Real*8 :: am,xm,ym,zm Integer*8 :: jp x = Xc(iHalo) ; y = Yc(iHalo); z = Zc(iHalo) Vx1=VXc(iHalo) ; Vy1= VYc(iHalo); Vz1= VZc(iHalo) vrms2 = (Unbound*Vcirc(iHalo))**2 ! do not take too fast particles Rvir = Rmc(iHalo) ! virial radius R0 = Rvir dmax = Rvir**2 xm =0. ; ym =0.; zm =0. am =0. mm =0 nn = 0 ! counter for pre-selected particles ! limits for Label i1=INT((x-R0)/Cell) j1=INT((y-R0)/Cell) k1=INT((z-R0)/Cell) i1=MIN(MAX(Nmx,i1),Nbx) j1=MIN(MAX(Nmy,j1),Nby) k1=MIN(MAX(Nmz,k1),Nbz) i2=INT((x+R0)/Cell) j2=INT((y+R0)/Cell) k2=INT((z+R0)/Cell) i2=MIN(MAX(Nmx,i2),Nbx) j2=MIN(MAX(Nmy,j2),Nby) k2=MIN(MAX(Nmz,k2),Nbz) ! Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) Do while (jp.ne.0) dd =(x-Xp(jp))**2+(y-Yp(jp))**2+(z-Zp(jp))**2 If(dd.lt.dmax)Then d1 =sqrt(dd) w = iWeight(jp) dx = Xp(jp) - x dy = Yp(jp) - y dz = Zp(jp) - z dvx = VXp(jp) - Vx1 +100.*dx ! correct for hubble flow dvy = VYp(jp) - Vy1 +100.*dy ! get proper velocity dvz = VZp(jp) - Vz1 +100.*dz ee = dvx**2 +dvy**2 +dvz**2 If(ee < vrms2)Then !vr = (dvx*dx+dvy*dy+dvz*dz)/d1 !ax = dy*dvz - dz*dvy !ay = dz*dvx - dx*dvz !az = dx*dvy - dy*dvx dd = d1/Rvir ! radius in units of virial radius ind = Max(Min(INT(log10(dd)/dlog+100)-100,Nbin_max),Nbin_min) aMr(ind,iHalo) = aMr(ind,iHalo) + w Vvrms(ind,iHalo) = Vvrms(ind,iHalo) + ee*w nPr(ind,iHalo) = nPr(ind,iHalo) + 1 Rbin(ind,iHalo) = Rbin(ind,iHalo) + dd*w EndIf ! ee < vrms2 EndIf ! dd < dmax jp =Lst(jp) End Do ! jp /= 0 EndDo EndDo EndDo Amh(iHalo) = am Do ind=Nbin_min,Nbin_max ! Normalize profiles If(nPr(ind,iHalo)==0)Then w = 1. Else w = aMr(ind,iHalo) EndIf Vvrms(ind,iHalo) = sqrt(Vvrms(ind,iHalo)/w) Rbin(ind,iHalo) = Rbin(ind,iHalo)/w EndDo ! ind Do ind=Nbin_min+1,Nbin_max ! Mass profile aMr(ind,iHalo) = aMr(ind-1,iHalo) +aMr(ind,iHalo) EndDo Return End SUBROUTINE SNeib !-------------------------------------------------------------- ! Add points around the Box to take into account periodical conditions ! It replicates points periodically and keeps those which are at distance ! less than Radius from the Box. ! Np = number of points inside the Box ! Ncp = total number of points including periodiacally replicated SUBROUTINE Points(Radius) !-------------------------------------------------------------- Logical Inside Integer*8 :: i,jp B1 =-Radius B2 =Box+Radius B3 =Box-Radius Ncp =Np ! remember primary points !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,k1,j1,i1,x0,y0,z0,Inside,ux,uy,uz,ww,zr,kk,yr,jj,xr,ii) Do i=1,Np ! Add dark matter particles x0 =xp(i) y0 =yp(i) z0 =zp(i) !if(mod(i,1000000)==0)write(*,*)' Points replication: ',i/1.e6 Inside = (x0.gt.Radius).and.(x0.lt.B3) Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3)) Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3)) If(.not.Inside)THEN ux =VXp(i) uy =VYp(i) uz =VZp(i) ww=iWeight(i) Do k1 = -1,1 zr =z0+k1*Box If(zr.GT.B1.AND.zr.LT.B2)Then kk =k1**2 Do j1 = -1,1 yr =y0+j1*Box If(yr.GT.B1.AND.yr.LT.B2)Then jj = kk +j1**2 Do i1 = -1,1 xr =x0+i1*Box ii =jj +i1**2 If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)Then !$OMP critical Ncp =Ncp +1 If(Ncp.le.Nmaxpart )Then xp(Ncp) =xr yp(Ncp) =yr zp(Ncp) =zr VXp(Ncp) =ux VYp(Ncp) =uy VZp(Ncp) =uz iWeight(Ncp) =ww EndIf !$OMP end critical EndIf EndDo EndIf EndDo EndIf EndDo EndIf EndDo If(Ncp.Gt.Nmaxpart)Then write(*,*)' Too many points:',Ncp,' Increase Nmaxpart =', Nmaxpart STOP EndIf Return End SUBROUTINE Points !-------------------------------- Update Statistics for DIstinct halos ! SUBROUTINE SPairs !--------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP SCHEDULE (DYNAMIC,kDynamic) Do i=1,Ncenter Call Sneib(i) EndDo Return End SUBROUTINE SPairs !--------------------------------------------------------- ! dump results ! SUBROUTINE WriteData !--------------------------------------------------------- WRITE (2,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv WRITE (40,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & NROWC,NGRID,NRECL,Om0,Oml0,hubble, & Ocurv write(2,110) dLog,Unbound,Nhmax,epsilon write(40,110) dLog,Unbound,Nhmax,epsilon 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) 110 format(1x,' dLog =',f6.3, & /1x,' Unbound factor (Vcirc units)= ',f6.2) scaleM = Dscale/weightSmall 50 format(a,T19,a, T45,a,T61,a,T71,a,T80,a, & T87,a,T92,a,T103,a,T113,a,T120,a,T128,2a) write(2,50) & ' Npart','Coords(Mpch)','Veloc(km/s)','Mvir', & 'Rvir','Vrms','Vcirc','Mvir_unbnd','Vrms','VirRat', & 'dR/Rvir','AngMom','Lambda' write(40,50) & ' Npart','Coords(Mpch)','Veloc(km/s)','Mvir', & 'Rvir','Vrms','Vcirc','Mvir_unbnd','Vrms','VirRat', & 'dR/Rvir','AngMom','Lambda' write(40,'(T3,a,T12,a,T20,a,T28,a,T40,a,T55,a, & T63,a,T71,a,T90,a,T92,a,T99,a)') & 'Npart','Rout/Rvir','Rbin','M( Nmaxpart)Stop ' Too many particles. Increase Nmaxpart' DO IN=1,In_page i = iL +IN ! Xp(i) =(XPAR(IN)-1.)*Xscale ! Yp(i) =(YPAR(IN)-1.)*Xscale ! Zp(i) =(ZPAR(IN)-1.)*Xscale ! VXp(i) =VX(IN) *Vscale ! VYp(i) =VY(IN) *Vscale ! VZp(i) =VZ(IN) *Vscale EndDo EndDo Np = lspecies(Nspecies) write(*,*) ' Done Reading Points. N=',Np !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) & !$OMP REDUCTION(MAX:xmax,ymax,zmax,xmaxs,ymaxs,zmaxs) & !$OMP REDUCTION(MIN:xmin,ymin,zmin,xmins,ymins,zmins) Do i =1,Np xmin =min(Xp(i),xmin) xmax =max(Xp(i),xmax) ymin =min(Yp(i),ymin) ymax =max(Yp(i),ymax) zmin =min(Zp(i),zmin) zmax =max(Zp(i),zmax) If(i<=Lsmall)Then xmins =min(Xp(i),xmins) xmaxs =max(Xp(i),xmaxs) ymins =min(Yp(i),ymins) ymaxs =max(Yp(i),ymaxs) zmins =min(Zp(i),zmins) zmaxs =max(Zp(i),zmaxs) EndIf ENDDO write(*,*) ' All Particles :',Np write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax write(*,*) ' Small Particles:',Lsmall write(*,*) ' x:',xmins,xmaxs write(*,*) ' y:',ymins,ymaxs write(*,*) ' z:',zmins,zmaxs write (*,*) ' Nmaxpart=',Nmaxpart do i=1,1000,200 write (*,'(i8," Test Velocities=",3g11.4)')i,VXp(i),VYp(i),VZp(i) EndDo Close(21) Return End Subroutine ReadPnt !-------------------------------------------------- ! Set Weights of particles for ! fast access SUBROUTINE SetWeights !-------------------------------------------------- Integer*8 :: jend,N_particles,jstart,k Wpr =NGRID**3/FLOAT(NROW)**3 N_particles =min(lspecies(Nspecies),Nmaxpart) jstart =1 If(N_particles.le.0)Then write (*,*) ' Wrong number of particles:',N_particles write (*,*) ' Nspecies=',Nspecies,' N=',lspecies STop endif Do j=1,Nspecies jend =min(lspecies(j),Nmaxpart) ww = wspecies(j)/Wpr !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (k) Do k=jstart ,jend iWeight(k) =ww EndDo jstart =jend EndDo weightSmall = iWeight(1) RETURN END SUBROUTINE SetWeights ! ------------------------ real function seconds () ! ------------------------ ! ! purpose: returns elapsed time in seconds Integer, SAVE :: first=0,rate=0 Real , SAVE :: t0=0. !real*8 dummy !real tarray(2) If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 t0 =float(i)/float(rate) seconds = 0. Else ! seconds = etime(tarray) CALL SYSTEM_CLOCK(i) seconds = float(i)/float(rate)-t0 EndIf return end function seconds !--------------------------------------------------- ! 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(txt) !--------------------------------------------------- ! Open file on a tape Character (LEN=*) :: txt Character*120 :: file1 file1 ='PMFILES/PMcrd.0'//TRIM(txt)//'.DAT' OPEN(UNIT=9,FILE=file1, & 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,EK,EK1,EK2,AU0,AEU0, & NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 & ,Ocurv,extras WRITE (*,100) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & EK,EK1,EK2, & 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 Jfiles =1 If(lspecies(Nspecies)> 1024**3)Jfiles=8 Do i =1,Jfiles write(file1,'(a,i1,a)')'PMFILES/PMcrs',i-1,'.0'//TRIM(txt)//'.DAT' write(*,'(2a,i3)') ' Open file=',file1,i OPEN(UNIT=20+i,FILE=TRIM(file1),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCES) EndDo REWIND 9 write (*,*) ' done openning PMcrs0' RETURN 20 write (*,*) ' Error reading the header file: Abort' stop END SUBROUTINE RDTAPE !-------------------------------------------- ! Read current PAGE of particles from disk ! NRECL - length of ROW block in words SUBROUTINE GETROW(IROW,Ifile) !-------------------------------------------- READ (20+Ifile,REC=IROW) RECDAT RETURN END SUBROUTINE GETROW end module SetHaloParticles !-------------------------------------------------------------------------------------- ! ! PROGRAM HaloParticles use SetHaloParticles Character*120 :: catname,listSh,listCat,txt*8,haloList,halo_ascii integer :: omp_get_max_threads ! ------------------------------------------ Read data write(*,'(a,$)')' Enter snapshot=' read(*,'(a)')txt catname='CATALOGS/CatshortA.'//TRIM(txt)//'.DAT' listSh ='CatshortC.'//TRIM(txt)//'.DAT' listCat ='CatalogC.'//TRIM(txt)//'.DAT' tt = seconds() CALL RDTAPE(txt) ! read control info, open particles file CALL SetWeights ! set iWeight array and weightSmall If(extras(100).gt.0.)Then Box =extras(100) Else write(*,'(a,$)')' Enter box size in comoving Mpc/h =' read(*,*) Box EndIf WRITE (*,120) & AEXPN,ISTEP,NROW,NGRID write (*,*) ' Number of mass Species =',Nspecies write (*,*) ' Total number of particles=',lspecies(Nspecies) write (*,*) ' Number of small particles=',lspecies(1) write (*,*) ' Box size (Mpc/h) = ',Box 120 FORMAT(1X, & 5x,'Expansion Parameter= ',F8.3,' Step= ',i5,/ & 5x,'Nrow= ',I4,' Ngrid= ',I4) hsmall= hubble ! Hubble/100. Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= 100.*Xscale/AEXPN ! Scale for velocities Dscale= 2.746e+11*Om0*(Box/NROW)**3 *weightSmall ! mass scale write (*,'(6x,"Mass of smallest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') Dscale write(*,'(6x,"first species=",2g11.3)') wspecies(1),weightSmall If(Nspecies.gt.1)Then write (*,'(6x,"Mass of largest particle",/9x, & "in units M_sun/h is =",3x,g10.3)') & Dscale*wspecies(Nspecies)/(NGRID**3/FLOAT(NROW)**3) EndIf ! ------------------------------------------ Read data xmins =0.; xmaxs =250. ymins =0.; ymaxs =250. zmins =0.; zmaxs =250. CALL ReadPnt(xmins,ymins,zmins,xmaxs,ymaxs,zmaxs) ! read particles: Np Open( 2,file=listSh,Status='UNKNOWN')! short list of halos Open( 40,file=listCat,Status='UNKNOWN')! short list of halos Open( 3,file=catname,Status='UNKNOWN')! short list of halos ! Find properties of Distinct/Isolated halos CALL ReadHalos(0) ! count halos: Ncenter ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter)) ALLOCATE (Amh(Ncenter),Vcirc2(Ncenter)) write (*,*) ' Allocated Centers' CALL ReadHalos(1) ! read halos Close(3) write (*,*) ' Read Halos: ',Ncenter Ncp = Np CALL Points(Rdmax) ! Make more points periodically write (*,*) 'Points made' write(*,*) '******* time=',seconds() ! set limits for linker list Celloptimal = Rdmax/FracSearch Cell = max(Box/Nlinker,Celloptimal) Nmx = INT(xmins/Cell)-Nbf ! boundaries for linker list Nbx = INT(xmaxs/Cell)+Nbf Nmy = INT(ymins/Cell)-Nbf Nby = INT(ymaxs/Cell)+Nbf Nmz = INT(zmins/Cell)-Nbf Nbz = INT(zmaxs/Cell)+Nbf write(*,'(2(a,f8.4))') ' Cell for linker list(Mpch)=',Cell,' Optimal=',Celloptimal write(*,'(a,3(2x,2i4))') ' Dimensions for the linker list=',Nmx,Nbx,Nmy,Nby,Nmz,Nbz ww = 4.*((Nbx-Nmx)*(Nby-Nmy)*(Nbz-Nmz)+Ncp)/(1024.**2) write (*,'(a,f9.2)') ' Requested memory for linker-list(Mb)=',ww ALLOCATE (Lst(Ncp),Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) write(*,*) '******* time=',seconds() Call List ! make linker List, Label write(*,*) '******* time=',seconds() Nbin_max = INT(1./dlog+0.5) Nbin_min = 0 ALLOCATE (aMr(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Vvrms(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (nPr(Nbin_min:Nbin_max,Ncenter)) ALLOCATE (Rbin(Nbin_min:Nbin_max,Ncenter)) !$OMP PARALLEL DO DEFAULT(SHARED) & ! initiate in parallel !$OMP PRIVATE (ic,j) Do ic =1,Ncenter Do j=Nbin_min, Nbin_max aMr(j,ic) = 0. ; Vvrms(j,ic) = 0. Rbin(j,ic) = 0. ; nPr(j,ic) = 0. EndDo EndDo write (*,20) Np,Ncenter,Box,Cell write (2,20) Np,Ncenter,Box,Cell 20 Format(5x,'Statistics for ',i10,' points. Maxima=',i9, & 3x,'Box(Mpc)=',f6.1,' Cell for linker list(Mpch)=',f9.5) write(*,*) '******* time=',seconds() CALL SPairs ! Accumulate statistics write(*,*) '******* time=',seconds() CALL WriteData write(*,*) '******* time(sec)=',seconds() End PROGRAM HaloParticles