
      program read_snapshot

   !------------------------------------------------------------------
   ! This Fortran 90 routine shows how a snapshot of GADGET
   ! which is distributed into several files may be read in.
   ! In this example, coodinates of particles type 1 are read.
   !
   ! Each dump consists of several files, each with a variable 
   ! number of particles. However, these files all have the same format. 
   ! The sequence of information in the files is as follows:
   !
   ! Header
   ! Positions
   ! Velocities
   ! Particle ID's
   !  
   ! Masses (for particle types with vanishing mass entry in the header)
   ! Internal Energy (Gas)
   ! Density (Gas)
   !
   !-------------------------------------------------------------------


      implicit none
      integer snapshot, files,nmax,k

      parameter(nmax = 2*128**3)
c      parameter(nmax = 2*32**3)
c      parameter(nmax = 276498)
      parameter(snapshot=0)
      parameter (files =1)
c      integer, parameter :: SNAPSHOT = 0 ! number of dump
c   integer, parameter :: FILES = 8          ! number of files per snapshot

      character*200 infile 
c     parameter :: path='/afs/mpa/data/volker/lcdm_cr/'

      character*200 filename, snumber, fnumber

      integer*4 npart(0:5), nall(0:5)
      real*8    massarr(0:5)
      real*8    a
      real*8    redshift
      integer*4 fn,i,nstart,flag_sfr,flag_fb,flag_cooling
      integer*4 N,Ntot,NumFiles

c      integer*4 unused(24) 
      character*96 unused 
      real*8 Boxsize
      real*8 OmegaMatter
      real*8 OmegaLambda,HubbleParam


      real*4  pos(3,nmax),vel(3,nmax)
      integer*4 id(nmax)
      real*4 PartPos(3,nmax)
      real*4 uint(nmax/2)
   

      write(*,*) 'Max npart= ',nmax
      write(*,*) 'insert filename'
      read(*,'(a)') filename

     
c      print *,'opening...  '//filename(1:len_trim(filename))
      print *,'opening...  ',filename

   ! now, read in the header

      open (1, file=filename, form='unformatted')

      read(1) npart,massarr,a,redshift,flag_sfr,flag_fb,
     $     nall,flag_cooling,NumFiles,BoxSize,OmegaMatter,OmegaLambda,
     $     HubbleParam,unused
      
      close (1)
   

   ! Now we just read in the coordinates, and only keep those of type=1
   ! The array PartPos(1:3,1:Ntot) will contain all these particles


      Ntot= nall(1)
c      allocate(PartPos(1:3, 1:Ntot))

      nstart=0
      do fn=0, FILES-1 
c         write(snumber,'(I3)') SNAPSHOT
c         write(fnumber,'(I3)') fn
c         do i=1,3 
c            if (snumber(i:i) .eq. ' ') then 
c               snumber(i:i)='0'
c            end if
c         end do
c         filename= path(1:len_trim(path)) // '/snapshot_' // snumber(1:len_trim(snumber)) // '.' // fnumber(verify(fnumber,' '):3)
c         print *,'reading...  '//filename(1:len_trim(filename))


! now, read in the file

         open (1, file=filename, form='unformatted')
c         read (1) npart, massarr, a, redshift, flag_sfr,
c     $        flag_fb, nall, unused

         read(1) npart,massarr,a,redshift,flag_sfr,flag_fb,
     $     nall,flag_cooling,NumFiles,BoxSize,OmegaMatter,OmegaLambda,
     $     HubbleParam,unused


         write(*,*) npart,massarr,a,redshift,flag_sfr,flag_fb,
     $     nall,flag_cooling,NumFiles,BoxSize,OmegaMatter,OmegaLambda,
     $     HubbleParam,unused


         
c         N=sum(npart)
c      allocate(pos(1:3,1:N))

         do i=0,5
            print *,'Type=',i,'    Particles=', nall(i),massarr(i)
         end do

         ntot =0 
         do k=0,5
            ntot = ntot + npart(k)
         enddo
         write(*,*) 'ntot ',ntot
         write(*,*) 'Box Size (Mpc/h) ',BoxSize
         write(*,*) 'Time (aexp) = ',a
         write(*,*) 'Redshift = ',redshift
         write(*,*) 'OmegaM, OmegaL ',OmegaMatter, OmegaLambda
         write(*,*) 'Hubble param ',HubbleParam

         read(1) ((pos(i,k),i=1,3),k=1,ntot)
         read(1) ((vel(i,k),i=1,3),k=1,ntot)
C         read(1) (mass(k),k=1,ntot)
         read(1) (id(k),k=1,ntot)
         if(npart(0).gt.0) then
            read(1) (uint(k),k=1,npart(0))
                         ! Note: If only the positions are needed,
! only those have to be read in !
         endif
         close (1)
      end do
      print *,'Done with reading.'

      write(*,*) 'writing only particles type 0 and 1'

      open(20,file='type0.dat',status='unknown')
      open(21,file='type1.dat',status='unknown')

      do k=1,npart(0)
         write(20,*) pos(1,k),pos(2,k),pos(3,k),uint(k),vel(1,k),id(k)
      enddo

      nstart = 0 
      nstart = npart(0)

      do k=nstart+1,npart(1)+nstart
         write(21,*) pos(1,k),pos(2,k),pos(3,k),id(k)
      enddo
c      nstart = npart(1) + nstart
c      do k=nstart+1,npart(2)+nstart
c         write(22,*) pos(1,k),pos(2,k),pos(3,k)
c      enddo
c      nstart = npart(2) + nstart
c      do k=nstart+1,npart(3)+nstart
c         write(23,*) pos(1,k),pos(2,k),pos(3,k)
c      enddo

         
 

! At this point, the coordinates of all particles of type 1 will
! be stored in the array PartPos(,)

      stop
      end 




