!-------------------------------------------------------------------------------- ! read particles ! make 3D density map ! make 2d projection Module SetArrs Integer, parameter :: N=1024, Nz=2048 Real*4, DIMENSION(0:N,0:N,0:Nz) :: dens Real*4, DIMENSION(0:N,0:N) :: Sdens Integer, PARAMETER :: & NROW = 4096, & NGRID = 256, & nbyteword = 4, & NPAGE = NROW**2, & ! # particles in a record NRECL = NPAGE*6 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) EQUIVALENCE (wspecies(1),extras(1)), & (lspecies(1),extras(11)) !$OMP THREADPRIVATE(/ROW/) Contains !-------------------------------------------------------------------------------------------- ! ! Open Files, read control information ! Subroutine OpenPM Character :: fname*120 Open(4,file ='PMcrd.DAT',form ='UNFORMATTED',status ='UNKNOWN') READ (4,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 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) NACCES= NRECL*4 / nbyteword Jfiles =1 If(lspecies(Nspecies)> 1024**3)Jfiles=8 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 RETURN 20 write (*,*) ' Error reading the header file: Abort' stop end Subroutine OpenPM end Module SetArrs !-------------------------------------------------------------------------------- ! Program GetDensity use SetArrs Character :: Line*80 real*8 :: dx,x,y,z,d1,d2,d3,av,sig,fmin,fmax,ww integer*8 :: Lsmall,N_particles,JPAGE,iL,IN,Np real*4 :: distr(-30:20),dd(-30:20) Box =250. ! computational box pBox =31./16. ! size of the plotted box dx = pBox/(N-1) ! step for coordinates ! define border of the box ! xc = 139. ; yc =68.; zc =15. xc = 49. ; yc =166.9; zc =92.2 xmin =xc-pBox/2. ; xmax =xc +pBox/2. ymin =yc-pBox/2. ; ymax =yc +pBox/2. zmin =zc-dx*Nz/2. ; zmax =zc +dx*Nz/2. Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= 100.*Xscale/AEXPN ! Scale for velocities write(*,*) ' Xlimits=',xmin,xmax write(*,*) ' Ylimits=',ymin,ymax write(*,*) ' Zlimits=',zmin,zmax x0 = xmin y0 = ymin z0 = zmin !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(i,j,k) Do k=0,Nz Do j=0,N Do i=0,N dens(i,j,k) =0. Enddo EndDo Enddo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE(i,j) Do j=0,N Do i=0,N sdens(i,j) =0. Enddo EndDo Call OpenPM JPAGE = NPAGE ! make it int*8 Lsmall = lspecies(1) N_particles = lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/JPAGE +1 N_in_last = N_particles -JPAGE*(Npages-1) write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' N_in_last=',N_in_last Np =0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IROW,In_page,iL,IN,x,y,z,vvx,vvy,vvz,i,j,k,i1,j1,k1, & !$OMP d1,d2,d3,w1,w2,w3,w4,w5,w6,w7,w8) & !$OMP REDUCTION(+:Np) DO IROW=16, 300 !Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = JPAGE*(IROW-1) if(mod(IROW,1).eq.0)write (*,*)' Read page=',IROW,' Np=',iL ! Loop over particles 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,vx,vy,vz Else READ (21,REC=IROW) XPAR,YPAR,ZPAR,VX,VY,VZ EndIf DO IN=1,In_page x =(XPAR(IN)-1.)*Xscale ! this is the place where we have y =(YPAR(IN)-1.)*Xscale ! all coordinates available z =(ZPAR(IN)-1.)*Xscale !vvx =VX(IN) *Vscale !vvy =VY(IN) *Vscale !vvz =VZ(IN) *Vscale !write(*,'(3x,3i10,3x,3f9.3)')IROW,IN,iL+IN,x,y,z If(x>xmin .and.xymin .and.yzmin .and.z=dx)& ! write(*,'(3(f8.3,i5,2g11.4),i8,f8.4,1P,16g12.4)')x-x0,i,d1,dx-d1,y-y0,j,d2,dx-d2,z-z0,k,d3,dx-d3,Np,dx,w1,w2,w3,w4 !$OMP atomic dens(i,j,k) = dens(i,j,k) +w1 !$OMP atomic dens(i1,j,k) = dens(i1,j,k) +w2 !$OMP atomic dens(i,j1,k) = dens(i,j1,k) +w3 !$OMP atomic dens(i1,j1,k) = dens(i1,j1,k) +w4 !$OMP atomic dens(i,j,k1) = dens(i,j,k1) +w5 !$OMP atomic dens(i1,j,k1) = dens(i1,j,k1) +w6 !$OMP atomic dens(i,j1,k1) = dens(i,j1,k1) +w7 !$OMP atomic dens(i1,j1,k1) = dens(i1,j1,k1) +w8 EndIf EndDo EndDo write(*,*) ' Done Reading Points. N=',Np, lspecies(Nspecies) close (1) open(2,file='dens2.dat') CALL Statistics print *,' Now write to file' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ww) do j=1,N do i=1,N ww =0. do k=1,Nz ww = max(ww,dens(i,j,k)) ! +dens(i,j,k)**2 Enddo sdens(i,j) =ww/dx**3 !*1.e9 enddo enddo distr = 0. ; dd = 0. dlog = 0.25 w2 =1.e5 w1 = -1.e5 ww =0 do j=1,N do i=1,N w1 =max(w1,sdens(i,j)) w2 =min(w2,sdens(i,j)) ww = ww + sdens(i,j) ind =INT(log10(sdens(i,j))/dlog+1000.)-1000 ind =min(max(-30,ind),20) distr(ind) = distr(ind) +1 dd(ind) = dd(ind) +sdens(i,j) enddo enddo !distr = distr write(*,*) ' Max Sdens=',w1 write(*,*) ' Min Sdens=',w2 write(*,*) ' Av Sdens=',ww/N**2 do i=-30,20 write(*,'(3x,i3,1P,3g12.4)')i,10.**(i*dlog),distr(i),dd(i)/max(1.,distr(i)) enddo do j=1,N write (2,'(1P,16g12.4)')(sdens(i,j),i=1,N) enddo Stop end Program GetDensity !------------------------------------------------- Subroutine Statistics use SetArrs Real*4 :: av, sig, fmin, fmax Integer*8 :: nem nem =0 av =0. sig=0. fmax =-1.e10 fmin =1.e10 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j) REDUCTION(+:nem,av,sig) & !$OMP REDUCTION(MAX:fmax) REDUCTION(MIN:fmin) Do k=1,Nz Do j=1,N Do i=1,N if(dens(i,j,k)<1.e-10)nem=nem+1 av = av +dens(i,j,k) sig = sig +dens(i,j,k)**2 fmax = max(fmax,dens(i,j,k)) fmin = min(fmin,dens(i,j,k)) if(dens(i,j,k)<0.)write(*,'(3x,g12.4,3i6)') dens(i,j,k),i,j,k EndDo EndDo EndDo acells =float(N)**2*Nz write(*,*) ' Number of Cells=',acells, ' Empty=',nem write(*,*) ' Average =',av/acells,' rms =',sqrt(max(sig,1.d-10)/acells) write(*,*) ' Min =',fmin, ' max =',fmax return End Subroutine Statistics