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
      real,    dimension(:,:), allocatable :: pos,vel

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

! read data from snapshot file
      write(0,'("Reading snapshot file...")')
      call get_command_argument(1,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 = nparttab(1) ! we use particle type 1 here (gas particles)
      allocate(id(npart),mass(npart),rho(npart),ein(npart), &
               pos(3,npart),vel(3,npart))
      id   = 0.
      mass = 0.
      rho  = 0.
      ein  = 0.
      pos  = 0.
      vel  = 0.
      havemass = .false.
    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
        read(1) mass
        havemass = .true.
      elseif(label.eq.'RHO ') then
        read(1) rho
      elseif(label.eq.'U   ') then
        read(1) ein
      else
        read(1) ! just skip the blocks we're not interested in
      endif
      goto 1
    2 continue
      close(1)
    7 format(a4,": ",i8," bytes")

      if(.not.havemass) then
        write(0,'("no individual masses present")')
        mass = masstab(1)
      endif

! output particle data as text
      write(0,'("Writing particle data...")')
      write(*,'(2i8,9es12.4)') &
        (i,id(i),mass(i),rho(i),ein(i),pos(:,i),vel(:,i),i=1,npart)

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