subroutine mt_read_mergertree()

! Reads a merger tree from an HDF file in the official mt format
! You may edit this routine to read in only those fields/arrays that you are interested in.

! TODO: Currently this routine only reads in data from arrays.
!       It can be adapted to read in records by copying the format of mt_read_snapshot.f90

  use hdf5
  use iso_c_binding
  use mt
  use mt_data

  implicit none

  ! Abbreviations for c_strings needed in function calls
  character(csize) :: dot, int64_t, int32_t, float_t, double_t
  integer(hid_t) :: gid
  integer :: NHaloint         ! To avoid compilation type mismatch

  ! Make the abbreviations
  dot=mt_str_f2c('.')
  int64_t=mt_str_f2c('int64')
  int32_t=mt_str_f2c('int32')
  float_t=mt_str_f2c('float')
  double_t=mt_str_f2c('double')

  ! First check to see whether the MergerTree group exists and either open it or make it
  gid=mt_gopen(mt_fid,mt_str_f2c('MergerTree'))

  ! Read in the number of halos in the merger tree
  call mt_aget_int(gid,dot,mt_str_f2c('NHalo'),Nhaloint);
  NHalo=NHaloint
  print*,'NHalo=',NHalo

  ! You may edit this routine to read in only those arrays that you are interested in.
  ! HaloID
  allocate(mtHaloID(NHalo))
  call mt_dread(gid, mt_str_f2c('HaloID'), int64_t, mtHaloID)
  ! Snapshot
  allocate(mtSnapshot(NHalo))
  call mt_dread(gid, mt_str_f2c('Snapshot'), int32_t, mtSnapshot)
  ! Density
  allocate(mtDensity(NHalo))
  call mt_dread(gid, mt_str_f2c('Density'), float_t, mtDensity)
  ! Mass
  allocate(mtMass(NHalo))
  call mt_dread(gid, mt_str_f2c('Mass'), float_t, mtMass)
  ! FirstSubhaloIndex
  allocate(mtFirstSubhaloIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('FirstSubhaloIndex'), int32_t, mtFirstSubhaloIndex)
  ! NeighbourIndex
  allocate(mtNeighbourIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('NeighbourIndex'), int32_t, mtNeighbourIndex)
  ! HostHaloIndex
  allocate(mtHostHaloIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('HostHaloIndex'), int32_t, mtHostHaloIndex)
  ! FirstProgenitorIndex
  allocate(mtFirstProgenitorIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('FirstProgenitorIndex'), int32_t, mtFirstProgenitorIndex)
  ! NextSiblingIndex
  allocate(mtNextSiblingIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('NextSiblingIndex'), int32_t, mtNextSiblingIndex)
  ! DescendantIndex
  allocate(mtDescendantIndex(NHalo))
  call mt_dread(gid, mt_str_f2c('DescendantIndex'), int32_t, mtDescendantIndex)

  ! Now close the open group
  call mt_gclose(gid)

end subroutine mt_read_mergertree
