read_points Subroutine

   subroutine read_points(filename, m)
      character(len=*), intent(in) :: filename
      type(mesh_t), intent(inout) :: m

      integer :: unit_id, ios, i, id
      real(rk) :: x, y, z

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('mesh_io', 'could not open '//trim(filename))

      read(unit_id, *, iostat=ios) m%npoints
      if (ios /= 0 .or. m%npoints <= 0) call fatal_error('mesh_io', 'invalid points header')

      allocate(m%points(3, m%npoints))
      do i = 1, m%npoints
         read(unit_id, *, iostat=ios) id, x, y, z
         if (ios /= 0) call fatal_error('mesh_io', 'failed reading point')
         if (id < 1 .or. id > m%npoints) call fatal_error('mesh_io', 'point id out of range')
         m%points(:, id) = [x, y, z]
      end do

      close(unit_id)
   end subroutine read_points