!
!-BEGIN_PROLOG-----------------------------------------------------------------
!
! NAME: mabel_h5_reader.f90
!
! PURPOSE: Read Mabel L1A, L1B or L2A HDF5 file.  This reads
!    some selected parameters from various HDF5 groups.  The parameter
!    list is not all-inclusive, but meant to show examples of how to read
!    different types of data.
!
! FILES ACCESSED: Input HDF5 file
!
! REFERENCES: See Mabel Detailed Design Documents
!
! COMMENTS: None.
!
! TO COMPILE: h5fc mabel_h5_reader.f -o mabel_h5_reader
!
! HISTORY:
!
!   ---------- --------- ---- -------------------------------------------------
!   2011-07-26 jelee1    0007 Create example HDF5 reader
!   2011-07-27 jelee1    0044 gps_sec_offset is bad value in L1B
!   2012-05-24 jelee1    0132 Added L2A support
!   2012-10-01 jelee1    0160 Update sample Fortran reader for Release 008
!   2014-01-29 jelee1    0222 Clean up unused vars (MABEL)
!   2014-06-05 jelee1    0286 Change gps_sec_offset to granule_gps_epoch
!   2014-11-07 jelee1    0286 Adjusted paths for MABEL 2.0.1.
!
!-END_PROLOG-------------------------------------------------------------------
!
program mabel_h5_reader

  use HDF5
  use H5LT

  character(len=255) :: &
    infile                    ! Input filename
!
! HDF5 IDs
!
  integer(HID_T) :: &
    h5_fileid                 ! HDF5 file pointer
!
! Metadata parameters
!
  character(len=255) :: &
    granule_type              ! Type of input granule
!
! NOTE: HDF5 likes to read datasets as arrays, not single elements.
!       So to read a single element, you have to declare it as an array
!       of dimension(1).
!
  double precision, dimension(1) :: &
    granule_gps_epoch         ! GPS seconds epoch for this granule
!
! Flight_param parameters
!
  character(len=255) :: &
    location                  ! Descriptive flight location
  integer, allocatable :: &
    channel_1064(:)           ! List of channel numbers for this granule
!
! Range/Elev parameters
!
  real, allocatable :: &
    delta_time(:), &          ! Delta time (from granule_gps_epoch)
    range(:)                  ! Range or Elevation (depending on granule type)
  double precision, allocatable :: &
    gps_sec(:)                ! GPS seconds (epoch+delta)
!
! Misc variables
!
  integer(HSIZE_T),dimension(1) :: &
    dims                      ! HDF5 dimensions
  integer(SIZE_T) :: i, &
    i_size                    ! HDF5 size
  integer :: &
    rank, &                   ! HDF5 rank
    i_class, &                ! HDF5 class
    i_res                     ! Result code
  character(len=255) :: &
    m_label, m_path           ! Strings used for generating labels
  logical :: &
    l_found                   ! Flag indicating if data were found
!
! Get input file from command line
!
  if (command_argument_count() < 1) then
    print *, "Usage: mabel_h5_reader INPUTFILE"
    call exit(3)
  endif
  call get_command_argument(1,infile)
!
! Open the HDF5 library
!
  call H5open_f(i_res)
!
! Open HDF5 File
!
  call H5Fopen_f(infile, H5F_ACC_RDONLY_F, h5_fileid, i_res)
  if (i_res/=0) then
    print *, "IO error opening: " // trim(infile)
    call exit(3)
  endif
!
! Read a simple dataset from the ancillary group.  granule_gps_epoch is
! important since it needs to be added to the delta_time in order
! to get a full GPS seconds-of-week value.
!
  dims=1
  rank=1
  granule_gps_epoch=0.0d0
  call h5ltread_dataset_f(h5_fileid, "/ancillary_data/granule_gps_epoch", &
    H5T_NATIVE_DOUBLE, granule_gps_epoch, dims, i_res)
  if (i_res == 0) then
    print *, "/ancillary_data/granule_gps_epoch=", granule_gps_epoch(1)
  else
    print *, "Error reading /ancillary_data/granule_gps_epoch"
    call exit(3)
  endif
!
! Lets also read the attribute "identifier_product_type". We will need to
! know the product type when we generate paths for the photon information.
!
  granule_type=""
  call H5LTget_attribute_string_f(h5_fileid, "/", &
    "identifier_product_type", granule_type, i_res)
  print *, "identifier_product_type=",trim(granule_type)
!
! Read a simple string dataset from flight_parameters.
!
  location=""
  call H5LTread_dataset_string_f(h5_fileid, &
    "/flight_parameters/flight_location", location, i_res)
  if (i_res == 0) then
    print *, "/flight_parameters/location=",trim(location)
  else
    print *, "Error reading /flight_parameters/flight_location"
    call exit(3)
  endif
!
! Read an array from the flightparam group. Lets read the values for the
! 1064 channels. To do this, we need to find the size of the array, allocate
! the local variables, and then read the data. Channel numbered "0" are
! not valid channels
!
  call h5ltget_dataset_info_f(h5_fileid, &
    "/flight_parameters/channel_1064", dims, i_class, i_size, i_res)
  if (dims(1) > 0 .and. i_res == 0) then
    if (allocated(channel_1064)) deallocate(channel_1064)
    allocate(channel_1064(dims(1)))
    call h5ltread_dataset_f(h5_fileid, "/flight_parameters/channel_1064", &
      H5T_NATIVE_INTEGER, channel_1064, dims, i_res)
    do i = 1, dims(1)
      if (channel_1064(i)==0) exit
      print *, "/flight_parameters/channel_1064=", channel_1064(i)
    enddo
  else
    print *, "Error reading channel_1064 from flight_parameters group"
    call exit(3)
  endif
!
! Open the first range channel group. This is a little tricky since during
! processing, bad channels may have been deleted.
!
  l_found=.TRUE.
  find_channel : do i = 1, size(channel_1064)

!
! Set the measurement path & label based on the granule_type. If you always
! just wanted to read L2A data, you wouldn't go through this -- you would just
! hardcode the path and parameter name and be done with it.
!
    if (granule_type(1:9) .eq. 'MABEL_L1A') then
      write(m_path,fmt='(a,i3.3)') "/range/channel",channel_1064(i)
      m_label="/range_uncorr"
    elseif (granule_type(1:9) .eq. 'MABEL_L1B') then
      write(m_path,fmt='(a,i3.3)') "/range/channel",channel_1064(i)
      m_label="/calibrated_range_m"
    elseif (granule_type(1:9) .eq. 'MABEL_L2A') then
      write(m_path,fmt='(a,i3.3,a)') "/channel",channel_1064(i),"/photon"
      m_label="/ph_h"
    else
      print *, "Error: Unknown granule type: "//trim(granule_type)
      call exit(3)
    endif
!
! See if the path exists
!
    call h5ltpath_valid_f(h5_fileid,trim(m_path),.FALSE.,l_found, i_res)
    if (l_found) then
      exit find_channel
    else
      print *, trim(m_path) // " must have been skipped during processing."
    endif
  enddo find_channel
!
! Make sure we did find a channel
!
  if (.not. l_found) then
    print *, "Error finding a valid channel"
    call exit(3)
  endif
  print *, "Processing channel=", channel_1064(i)
!
! Read the delta_time from the first channel group. This reads
! the entire granule into memory. A more sophisticated example may read the
! the data as hyperslab chunks in order to conserve memory space.
!
  call h5ltget_dataset_info_f(h5_fileid, trim(m_path)//"/delta_time", &
    dims, i_class, i_size, i_res)
  if (dims(1) > 0 .and. i_res==0) then
    if (allocated(delta_time)) deallocate(delta_time)
    allocate(delta_time(dims(1)))
    call h5ltread_dataset_f(h5_fileid, trim(m_path)//"/delta_time", &
      H5T_NATIVE_REAL, delta_time, dims, i_res)
    if (i_res/=0) then
      print *, "Error reading "//trim(m_path)//"/delta_time."
      call exit(3)
    endif
  else
    print *, "Error reading "//trim(m_path)//"/delta_time."
    call exit(3)
  endif
!
! Now, add the granule_gps_epoch to create a full gps_seconds
!
  if (allocated(gps_sec)) deallocate(gps_sec)
  allocate(gps_sec(dims(1)))
  gps_sec=delta_time+granule_gps_epoch(1)
!
! Read the range or elevation. We used the granule_type earlier to figure
! out what type of file we were reading.
!
  call h5ltget_dataset_info_f(h5_fileid, trim(m_path)//trim(m_label), &
    dims, i_class, i_size, i_res)
  if (dims(1) > 0 .and. i_res==0) then
    if (allocated(range)) deallocate(range)
    allocate(range(dims(1)))
    call h5ltread_dataset_f(h5_fileid, trim(m_path)//trim(m_label), &
      H5T_NATIVE_REAL, range, dims, i_res)
    if (i_res/=0) then
      print *, "Error reading "//trim(m_path)//trim(m_label)
      call exit(3)
    endif
  else
    print *, "Error reading "//trim(m_path)//trim(m_label)
    call exit(3)
  endif
!
! Print the first 10 time/range pairs
!
  do i = 1, min(10,dims(1))
    print *, trim(m_label)//"=", gps_sec(i),range(i)
  enddo
!
! De-allocate allocated data
!
  if (allocated(channel_1064)) deallocate(channel_1064)
  if (allocated(delta_time)) deallocate(delta_time)
  if (allocated(gps_sec)) deallocate(gps_sec)
  if (allocated(range)) deallocate(range)
!
! Close the HDF5 File
!
  call H5Fclose_f(h5_fileid, i_res)
!
! Close the HDF5 library
!
  call H5close_f(i_res)

  end program mabel_h5_reader
