program readsnap
      implicit none

! data from snapshot file
      character(len=4) :: label
      integer, dimension(6) :: nparttab,nparttottab,nparttothwtab
      integer :: blocksize,npart,num_files,lpt_scalingfactor
      integer :: flag_sfr,flag_feedback,flag_cooling,flag_stellarage, &
                 flag_metals,flag_entropy_instead_u, &
                 flag_doubleprecision,flag_ic_info
      double precision :: time,redshift,boxsize, &
                          omega0,omegalambda,hubbleparam
      double precision, dimension(6) :: masstab
      integer, dimension(:),   allocatable :: id
      real,    dimension(:),   allocatable :: mass,rho,ein,temp,tmp
      real,    dimension(:,:), allocatable :: pos,vel

! other variables
      character(len=1000) :: filename
      integer :: i,j,k
      logical :: havemass

! read data from snapshot file
      call get_command_argument(1,filename)
      write(0,'("Reading snapshot file ",a,"...")') trim(filename)
      open(1,file=filename,status='old',form='unformatted')
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      read(1) nparttab,masstab,time,redshift,flag_sfr,flag_feedback, &
              nparttottab,flag_cooling,num_files,boxsize, &
              omega0,omegalambda,hubbleparam,flag_stellarage, &
              flag_metals,nparttothwtab,flag_entropy_instead_u, &
              flag_doubleprecision,flag_ic_info,lpt_scalingfactor
      write(0,'("  box size: ",g10.3)') boxsize
      write(0,'("  total:  ",i8," particles")') sum(nparttab)
      write(0,'("  type ",i1,": ",i8," particles, mass",es12.4)') &
                           (i,nparttab(i),masstab(i),i=1,6)
      npart = sum(nparttab)
      write(0,'("Using ",i0," particles.")') npart
      allocate(id(npart),mass(npart),rho(npart),ein(npart),temp(npart), &
               pos(3,npart),vel(3,npart),tmp(npart))
      id   = 0.
      rho  = 0.
      ein  = 0.
      temp = 0.
      pos  = 0.
      vel  = 0.
      j=1
      do i=1,6
        k = j+nparttab(i)
        mass(j:k-1) = masstab(i)
        j = k
      enddo

    1 continue
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      if(label.eq.'ID  ') then
        read(1) id
      elseif(label.eq.'POS ') then
        read(1) pos
      elseif(label.eq.'VEL ') then
        read(1) vel
      elseif(label.eq.'MASS') then
        k = (blocksize-8)/4
        read(1) tmp(1:k)
        j = 1
        do i=1,npart
          if(mass(i).eq.0.) then
            if(j.gt.k) stop 'Not enough masses. Stop.'
            mass(i) = tmp(j)
            j = j+1
          endif
        enddo
      elseif(label.eq.'RHO ') then
        read(1) rho(1:nparttab(1))
      elseif(label.eq.'U   ') then
        read(1) ein(1:nparttab(1))
      elseif(label.eq.'TEMP') then
        read(1) temp(1:nparttab(1))
      else
        read(1) ! just skip the blocks we're not interested in
      endif
      goto 1
    2 continue
      close(1)
    7 format(a4,": ",i8," bytes")

    do i=1,npart
      if(id(i).lt.0) id(i)=id(i)+2147483648
    enddo

! output particle data as text
      write(0,'("Writing particle data...")')
      write(*,'("# n1=",i0,"; n2=",i0,"; n3=",i0,"; n4=",i0,"; n5=",i0,"; n6=",i0)') &
        nparttab
      write(*,'(i12,8es12.4)') &
        (id(i),mass(i),ein(i),pos(:,i),vel(:,i),i=1,npart)

      write(0,'("Done.")')
      end