program readtree

! Example of how to read in a merger tree in HDF5 format

! The idea is to make this calling routine use Fortran-95 only and hide all the
! underlying HDF libraries and conversions to/from C from the user.

  ! We need to load up the HDF5 type definitions
  use hdf5

  ! Load function interfaces
  use mt

  ! This defines the merger tree structures as per the standard
  ! and defines some global vairables.
  use mt_data

  implicit none

  character*(CSIZE) :: file
  character*(CSIZE) :: description
  integer :: i

  ! Read in data file name
  if (command_argument_count().eq.1) then
     call get_command_argument(1,file)
  else
     stop 'Usage: readtree <treefile>'
  end if

  ! Open data file
  ! call mt_open(file,ioflags,description)
  ! Inputs: 
  !    file - file name
  !    ioflags - read/write string;
  !           'r' file must exist and will be opened readonly
  !           'rw' file must exist and will be opened for writing
  !           'wo' file will be created, over-writing any existing file
  !           'w' file will be created, will fail if already exists
  !    description - comment to be written/appended to file
  ! NOTE: description cannot be replaced with a string constant
  call mt_open(file,'r',description)
  print*,'description=',description(1:len_trim(description))

  ! Read snaptable from HDF file
  call mt_read_snaptable()
  print*,'Snap='
  do i=1,int(NSnap)
     print*,Snap(i)
  end do
  print*,'SnapProp='
  do i=1,int(NSnapProp)
     print*,SnapProp(i)
  end do

  ! Read merger tree from HDF file
  call mt_read_mergertree()
  print*,'MergerTree='
  do i=1,int(NHalo)
     print*,mtHaloID(i),mtSnapshot(i),mtDensity(i),mtMass(i),mtFirstSubhaloIndex(i),&
          &mtNeighbourIndex(i),mtHostHaloIndex(i),mtFirstProgenitorIndex(i),&
          &mtNextSiblingIndex(i),mtDescendantIndex(i)
  end do

  ! Close file
  call mt_close()

  ! Tidy up.  All this does is deallocate some arrays.
  call mt_tidy()

end program readtree
