module MixReader
  public openMixFile, closeMixFile, readVariable
  private
  integer fileId 

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine openMixFile(filename)
    integer DFACC_READ,DFACC_RDONLY
    parameter(DFACC_READ   = 1)
    parameter(DFACC_RDONLY = 1)
    integer DF_MAXFNLEN
    parameter(DF_MAXFNLEN = 256)

    character(len=*) :: filename

    integer sfstart

    write(*,*) "Opening file ", filename
    fileId = sfstart(filename, DFACC_RDONLY)
    call errorCheck(fileId, 'sfstart')

    ! Read some basic attributes.
    
  end subroutine openMixFile

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine closeMixFile()
    integer status
    integer sfend
    status = sfend(fileId)
    call errorCheck(status,'sfend')
  end subroutine closeMixFile

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine readVariable(varname,data)
    character(len=*) :: varname
    real, allocatable :: data(:,:)

    integer :: errorId
    integer :: varId
    integer :: selectId
    integer :: datatype
    integer :: nAttributes
    integer :: rank
    integer :: dimensions(32)
    integer :: indexStart(32)
    integer :: indexEnd(32)
    character(len=1024) :: dataName
    integer :: i

    integer sfn2index,sfselect,sfginfo

    integer nPoints

    do i=1,32
       dimensions(i) = 0
       indexStart(i) = 0
       indexEnd(i) = 0
    enddo

    ! SDnametoindex
    varId = sfn2index(fileId, varName)
    call errorCheck(varId,'sfn2index')
    ! sdselect
    selectId = sfselect(fileId, varId)
    call errorCheck(selectId,'sfselect')
    ! sdgetinfo
    errorId = sfginfo(selectId,dataName,rank, dimensions, datatype,nAttributes)
    call errorCheck(errorId,'sfginfo')

    if (rank .ne. 2) then
       STOP "Only works with 2d variables!"
    endif
    do i=0,rank
       indexStart(i+1) = 0
       indexEnd(i+1) = dimensions(i+1)
    enddo

    ! Allocate memory
    if (.not. allocated(data)) then
       allocate(data (dimensions(1), dimensions(2) ))
    endif
    
    errorId = sfrdata(selectId,indexStart,(/1,1/),indexEnd,data)
    call errorCheck(errorId, "sfrdata")

    errorId = sfendacc(selectId)
    call errorCheck(errorId, 'sfendacc')
    
  end subroutine readVariable

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine errorCheck(status,msg)
    integer status
    character(len=*) :: msg
    if (status .eq. -1) then
       write(*,*) msg
       stop "HDF Error"
    endif
  end subroutine errorCheck
end module MixReader