module modis

  use constants
  use hdf_read_module
  
  implicit none

  private

  public :: modis_get_dims, modis_read, modis_close_all,modis_read_bytes
  public :: subset_2d,subset_3d
  public :: MAX_RANK_HDF ! from module HDF
           
  ! specification of a 2d array subset
  type :: subset_2d
     integer(kind=int4), dimension(2) :: start
     integer(kind=int4), dimension(2) :: stride
     integer(kind=int4), dimension(2) :: edge
  end type subset_2d
  type :: subset_3d
     integer(kind=int4), dimension(3) :: start
     integer(kind=int4), dimension(3) :: stride
     integer(kind=int4), dimension(3) :: edge
  end type subset_3d
     

  type :: file_struct
      character(len=256) :: filename
      integer (kind=int4) :: id
  end type file_struct

  ! Keeping the last mod04 file and mod06 file open.
  ! FUTURE: Allow mutiple files of each type to be left open by keeping a vector of 
  ! open files. Could also cache data.
  type(file_struct), target, save :: last_mxd04_file
  type(file_struct), target, save :: last_mxd06_file
  logical, target, save :: mxd04_open = .false.
  logical, target, save :: mxd06_open = .false.

contains

  subroutine modis_get_dims(source, directory, sds_name, time, rank, dims, stat)
    integer(kind=int4), intent(in) :: source
    character(len=*), intent(in) :: directory
    character(len=*), intent(in) :: sds_name 
    type(time_struct), intent(in) :: time
    integer(kind=int4), intent(out) :: rank
    integer(kind=int4), dimension(MAX_RANK_HDF), intent(out) :: dims
    integer(kind=int4), intent(out) :: stat

    type(file_struct), pointer :: file
    integer(kind=int4) :: hdf_status

    stat = STATUS_OK

    call modis_open(source, directory, time, file, stat)
    if(stat /= STATUS_OK) return

    hdf_status = hdf_sds_dimensions_reader(file%id, sds_name, rank, dims)

    if(hdf_status /= sym%SUCCESS) then 
       stat = STATUS_ERR_MISSING_DATA
       return
    end if

  end subroutine modis_get_dims

  ! Data must be pre-allocated to at least the required size. 
  ! Warning: assumes Data has already been initialized to missing values
subroutine modis_read(source, directory, sds_name, subset, time, data, stat)
  integer(kind=int4), intent(in) :: source
  character(len=*), intent(in) :: directory
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  type(time_struct), intent(in) :: time
  real(kind=real8), dimension(:,:), intent(out) :: data
  integer(kind=int4), intent(out) :: stat

  type(file_struct), pointer :: file
  integer(kind=int4) :: sds_type
  integer(kind=int4) :: hdf_stat

  stat = STATUS_OK

  call modis_open(source, directory, time, file, stat)
  if(stat /= STATUS_OK) return

  ! get the SDS type
  hdf_stat = hdf_sds_type_reader(file%id, sds_name, sds_type)
  if(hdf_stat /= sym%SUCCESS) then 
     stat = STATUS_ERR_MISSING_DATA
     return
  end if
  
  ! call the right read routine for the SDS data type
  select case (sds_type)
  case (DFNT_INT8)
     call read_sds_int8(file, sds_name, subset, data, stat)
  case (DFNT_INT16)
     call read_sds_int16(file, sds_name, subset, data, stat)
  case (DFNT_FLOAT32)
     call read_sds_float32(file, sds_name, subset, data, stat)
  case (DFNT_FLOAT64)
     call read_sds_float64(file, sds_name, subset, data, stat)

     ! FUTURE: add cases for other types here

  case default
     stat = STATUS_ERR_DATA_TYPE_INVALID
     return
  end select
  
  if(stat /= STATUS_OK) return

end subroutine modis_read
subroutine modis_read_bytes(source, directory, sds_name, subset2d, time, idata, Fill_Value,rank,stat)
  integer(kind=int4), intent(in) :: source
  character(len=*), intent(in) :: directory
  character(len=*), intent(in) :: sds_name 
  type(subset_3d)  :: subset3d
  type(subset_2d), intent(in) :: subset2d
  type(time_struct), intent(in) :: time
  integer(kind=int1), dimension(:,:,:), intent(out) :: idata
  integer(kind=int1) Fill_Value
  integer(kind=int2), dimension(:,:), allocatable :: idatai2
  integer(kind=int4), intent(out) :: stat

  type(file_struct), pointer :: file
  integer(kind=int4) :: sds_type
  integer(kind=int4) :: hdf_stat
  integer(kind=int4) rank,shapein(3)

  stat = STATUS_OK

  call modis_open(source, directory, time, file, stat)
  if(stat /= STATUS_OK) return

  ! get the SDS type
  hdf_stat = hdf_sds_type_reader(file%id, sds_name, sds_type)
  if(hdf_stat /= sym%SUCCESS) then 
     stat = STATUS_ERR_MISSING_DATA
     return
  end if
  
  ! call the right read routine for the SDS data type
  select case (sds_type)
  case (DFNT_INT8)
     if(rank.eq.2)then
       call read_sds_int8_int8_2d(file, sds_name, subset2d, idata(1,:,:), Fill_Value,stat)
     else
       !subset3d%start(1)=1
   ! maybe should be 0?
       subset3d%start(1)=0
       shapein=shape(idata)
       subset3d%edge(1)=shapein(1) ! edge i number of elements
       subset3d%stride(1)=1
       subset3d%stride(2:3)=subset2d%stride(1:2)
       subset3d%start(2)=subset2d%start(1)
       subset3d%start(3)=subset2d%start(2)
       subset3d%edge(2)=subset2d%edge(1)
       subset3d%edge(3)=subset2d%edge(2)

       call read_sds_int8_int8_3d(file, sds_name, subset3d, idata, Fill_Value,stat)
     endif
  case (DFNT_INT16)
     call read_sds_int16_int16(file, sds_name, subset2d, idatai2, stat)

     ! FUTURE: add cases for other types here

  case default
     stat = STATUS_ERR_DATA_TYPE_INVALID
     return
  end select
  
  if(stat /= STATUS_OK) return

end subroutine modis_read_bytes

subroutine modis_close_all()

  if( mxd04_open ) then
     call modis_close( last_mxd04_file )
     mxd04_open = .false.
  end if

  if( mxd06_open ) then
     call modis_close( last_mxd06_file )
     mxd06_open = .false.
  end if
     
end subroutine modis_close_all

subroutine get_modis_filename(source, directory, time, filename, stat)
  integer(kind=int4), intent(in) :: source
  character(len=*), intent(in) :: directory
  type(time_struct), intent(in) :: time
  character(len=256), intent(out) :: filename
  integer(kind=int4), intent(out) :: stat
  
  character (len=10)  :: prefix
  character (len=4)   :: year_string
  character (len=3)   :: julian_day_string
  character (len=2)   :: hour_string
  character (len=2)   :: minute_string

  stat = STATUS_OK

  if(    source == SOURCE_MODIS_MOD04) then
     prefix = 'MOD04_L2.A'
  elseif(source == SOURCE_MODIS_MYD04) then
     prefix = 'MYD04_L2.A'
  elseif(source == SOURCE_MODIS_MOD06) then
     prefix = 'MOD06_L2.A'
  elseif(source == SOURCE_MODIS_MYD06) then
     prefix = 'MYD06_L2.A'
  else
     stat = STATUS_ERR_SOURCE_INVALID
     return
  end if

  write (year_string, fmt="(I4.4)") time%year
  write (julian_day_string, fmt="(I3.3)") time%jday
  write (hour_string, fmt="(I2.2)") time%hour
  write (minute_string, fmt="(I2.2)") time%minute
  
  filename = trim(directory) // '/' // trim(prefix) // year_string // julian_day_string // '.' // hour_string // minute_string // '.hdf'
  
!  print *, 'filename=', filename

end subroutine get_modis_filename

subroutine modis_close(file)
  type(file_struct), intent(in) :: file
    
  call close_file_hdf_read(file%id, file%filename)
  
end subroutine modis_close

! open file (or if already open, return its file struct)
subroutine modis_open(source, directory, time, file, stat)
  integer(kind=int4), intent(in) :: source
  character(len=*), intent(in) :: directory
  type(time_struct), intent(in) :: time
  type(file_struct), pointer, intent(out) :: file
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4) :: hdf_stat
  logical, pointer :: file_open
  character(len=256) :: filename

  stat = STATUS_OK

  if(source == SOURCE_MODIS_MOD04 .or. source == SOURCE_MODIS_MYD04) then
     file => last_mxd04_file
     file_open => mxd04_open
  elseif(source == SOURCE_MODIS_MOD06 .or. source == SOURCE_MODIS_MYD06) then
     file => last_mxd06_file
     file_open => mxd06_open
  else
     stat = STATUS_ERR_SOURCE_INVALID
     return
  end if

  call get_modis_filename(source, directory, time, filename, stat)
  if(stat /= STATUS_OK) return

  ! if the wrong file is open, close it
  if(file_open) then
     if(file%filename /= filename) then
        call modis_close(file)
        file_open = .false.
     end if
  end if

  ! open file if needed
  if(.not. file_open) then

     file%filename = filename
     hdf_stat = open_file_hdf_read(file%filename, file%id)

     if (hdf_stat == sym%FAILURE) then
        stat = STATUS_ERR_MISSING_DATA 
        return
     end if

     file_open = .true.

  end if

end subroutine modis_open

subroutine read_sds_int8(file, sds_name, subset, data, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  real(kind=real8), dimension(:,:), intent(out) :: data
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  integer(kind=int1), dimension(:,:), allocatable :: buf_int1

  ! SDS attributes
  integer(kind=int1) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
  integer(kind=int1), dimension(2) :: valid_range
  integer(kind=int1), dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start = subset%start
  stride = subset%stride
  edge = subset%edge

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_int1(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Int1)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_int1)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_int1)
     return
  end if

  ! if data is valid, unscale data. We assume Data has already been initialized to missing values  
  where(buf_int1 >= valid_range(1) .and. buf_int1 <= valid_range(2) .and. buf_int1 /= fill_value)
     data(1:edge(1), 1:edge(2)) = (real(buf_int1(:,:), kind=real8) - add_offset) * scale_factor
  end where

  deallocate(buf_int1)

end subroutine read_sds_int8
subroutine read_sds_int8_int8_2d(file, sds_name, subset, idata, fill_value, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  integer(kind=int1), dimension(:,:), intent(out) :: idata
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  integer(kind=int1), dimension(:,:), allocatable :: buf_int1

  ! SDS attributes
  integer(kind=int1) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
!  integer, parameter :: int8 = selected_int_kind(10)
  integer(kind=int1), dimension(2) :: valid_range
!  integer(kind=int4), dimension(:), allocatable :: valid_range_alloc
   integer(kind=int1), dimension(:), allocatable :: valid_range_alloc
!  integer(Kind=int8), dimension(:), allocatable :: valid_range_alloc
!  integer*1, dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat
  integer(kind=int4) :: shapein(2)
  integer(kind=int4) :: i
  shapein=shape(idata)

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start(1:2) = subset%start
  stride(1:2) = subset%stride
  edge(1:2) = subset%edge
!  write(6,*)'start',start,'edge',edge,'stride',stride,'name',trim(sds_name)

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_int1(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Int1)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_int1)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_int1)
     return
  end if
  ! problem is that top of valid range can be all 1 bits which looks like -1
  if ( valid_range(1) .eq. 0 .and. valid_range(2) .eq. -1 )then
    ! since all bits ok just make sure not fill_value
    if(add_offset.eq.0)then
      ! if fill_value is one of the valid_range values just return all bytes read
      if(fill_value == valid_range(1) .or. fill_value == valid_range(2) ) THEN
        ! just return all bytes read
        idata(1:edge(1), 1:edge(2)) = buf_int1(:,:) 
      else
        where( buf_int1 /= fill_value)
         idata(1:edge(1), 1:edge(2)) = buf_int1(:,:) 
        end where
      endif
    else
      where( buf_int1 /= fill_value)
       idata(1:edge(1), 1:edge(2)) = buf_int1(:,:) - int(add_offset)
      end where
    endif
  else
  ! if data is valid, add offset if there. We assume Data has already been initialized to missing values  
  where(buf_int1 >= valid_range(1) .and. buf_int1 <= valid_range(2) .and. buf_int1 /= fill_value)
     idata(1:edge(1), 1:edge(2)) = buf_int1(:,:) - int(add_offset)
  end where
  endif

  deallocate(buf_int1)

end subroutine read_sds_int8_int8_2d
subroutine read_sds_int8_int8_3d(file, sds_name, subset, idata, fill_value,stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_3d), intent(in) :: subset
  integer(kind=int1), dimension(:,:,:), intent(out) :: idata
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(3) :: start
  integer(kind=int4), dimension(3) :: stride
  integer(kind=int4), dimension(3) :: edge
  integer(kind=int1), dimension(:,:,:), allocatable :: buf_int1

  ! SDS attributes
  integer(kind=int1) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
!  integer, parameter :: int8 = selected_int_kind(10)
  integer(kind=int1), dimension(2) :: valid_range
!  integer(kind=int4), dimension(:), allocatable :: valid_range_alloc
   integer(kind=int1), dimension(:), allocatable :: valid_range_alloc
!  integer(Kind=int8), dimension(:), allocatable :: valid_range_alloc
!  integer*1, dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat
  integer(kind=int4) :: shapein(3)
  integer(kind=int4) :: j
  shapein=shape(idata)

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  ! hdf starts at 0
  start(1)=0
  stride(1)=1
  edge(1)=shapein(1)
  start(2:3) = subset%start(2:3)
  stride(2:3) = subset%stride(2:3)
  edge(2:3) = subset%edge(2:3)

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_int1(edge(1), edge(2),edge(3)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Int1)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_int1)
  if (hdf_stat == sym%FAILURE .or. edge(2) /= subset%edge(2) .or. edge(3) /= subset%edge(3)) then
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_int1)
     return
  end if
  ! problem is that top of valid range is all 1 bits which looks like -1
  if ( valid_range(1) .eq. 0 .and. valid_range(2) .eq. -1 )then
    ! since all bits ok just make sure not fill_value
    if(add_offset.eq.0)then
      if(fill_value == valid_range(1) .or. fill_value == valid_range(2) ) THEN
        ! if fill_value is one of the valid_range values just return all bytes read
         idata(1:edge(1), 1:edge(2),1:edge(3)) = buf_int1(:,:,:) 
      else
        where( buf_int1 /= fill_value)
         idata(1:edge(1), 1:edge(2),1:edge(3)) = buf_int1(:,:,:) 
        end where
      endif
    else
      where( buf_int1 /= fill_value)
       idata(1:edge(1), 1:edge(2),1:edge(3)) = buf_int1(:,:,:) - int(add_offset)
      end where
    endif
  else
  ! if data is valid, add offset if there. We assume Data has already been initialized to missing values  
  where(buf_int1 >= valid_range(1) .and. buf_int1 <= valid_range(2) .and. buf_int1 /= fill_value)
     idata(1:edge(1), 1:edge(2),1:edge(3)) = buf_int1(:,:,:) - int(add_offset)
  end where
  endif

  deallocate(buf_int1)

end subroutine read_sds_int8_int8_3d


subroutine read_sds_int16(file, sds_name, subset, data, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  real(kind=real8), dimension(:,:), intent(out) :: data
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  integer(kind=int2), dimension(:,:), allocatable :: buf_int2

  ! SDS attributes
  integer(kind=int2) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
  integer(kind=int2), dimension(2) :: valid_range
  integer(kind=int2), dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start = subset%start
  stride = subset%stride
  edge = subset%edge

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_int2(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Int2)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_int2)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_int2)
     return
  end if

  ! if data is valid, unscale data. We assume Data has already been initialized to missing values  
  where(buf_int2 >= valid_range(1) .and. buf_int2 <= valid_range(2) .and. buf_int2 /= fill_value)
     data(1:edge(1), 1:edge(2)) = (real(buf_int2(:,:), kind=real8) - add_offset) * scale_factor
  end where

  deallocate(buf_int2)

end subroutine read_sds_int16
subroutine read_sds_int16_int16(file, sds_name, subset, idata, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  integer(kind=int2), dimension(:,:), intent(out) :: idata
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  integer(kind=int2), dimension(:,:), allocatable :: buf_int2

  ! SDS attributes
  integer(kind=int2) :: fill_value
  real(kind=real8) :: scale_factor
  integer(kind=int2) :: add_offset
  integer(kind=int2), dimension(2) :: valid_range
  integer(kind=int2), dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start = subset%start
  stride = subset%stride
  edge = subset%edge

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_int2(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Int2)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_int2)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_int2)
     return
  end if

  ! if data is valid, don't unscale data. We assume Data has already been initialized to missing values  
  where(buf_int2 >= valid_range(1) .and. buf_int2 <= valid_range(2) .and. buf_int2 /= fill_value)
     idata(1:edge(1), 1:edge(2)) = buf_int2(:,:) - add_offset
  end where

  deallocate(buf_int2)

end subroutine read_sds_int16_int16

subroutine read_sds_float32(file, sds_name, subset, data, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  real(kind=real8), dimension(:,:), intent(out) :: data
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  real(kind=real4), dimension(:,:), allocatable :: buf_real4

  ! SDS attributes
  real(kind=real4) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
  real(kind=real4), dimension(2) :: valid_range
  real(kind=real4), dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start = subset%start
  stride = subset%stride
  edge = subset%edge

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_real4(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Real4)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_real4)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_real4)
     return
  end if

  ! if data is valid, unscale data. We assume Data has already been initialized to missing values  
  where(buf_real4 >= valid_range(1) .and. buf_real4 <= valid_range(2) .and. buf_real4 /= fill_value)
     data(1:edge(1), 1:edge(2)) = (real(buf_real4(:,:), kind=real8) - add_offset) * scale_factor
  end where

  deallocate(buf_real4)

end subroutine read_sds_float32

subroutine read_sds_float64(file, sds_name, subset, data, stat)
  type(file_struct), intent(in) :: file
  character(len=*), intent(in) :: sds_name 
  type(subset_2d), intent(in) :: subset
  real(kind=real8), dimension(:,:), intent(out) :: data
  integer(kind=int4), intent(out) :: stat

  integer(kind=int4), dimension(2) :: start
  integer(kind=int4), dimension(2) :: stride
  integer(kind=int4), dimension(2) :: edge
  real(kind=real8), dimension(:,:), allocatable :: buf_real8

  ! SDS attributes
  real(kind=real8) :: fill_value
  real(kind=real8) :: scale_factor
  real(kind=real8) :: add_offset
  real(kind=real8), dimension(2) :: valid_range
  real(kind=real8), dimension(:), allocatable :: valid_range_alloc

  integer(kind=int4) :: hdf_stat
  integer(kind=int4) :: alloc_stat

  stat = STATUS_OK

  ! Subsetting parameters can be modified in HDF read routine, so need to copy to local variables
  start = subset%start
  stride = subset%stride
  edge = subset%edge

  ! Read attributes
  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "scale_factor", &
       scale_factor)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "add_offset", &
       add_offset)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "_FillValue", &
       fill_value)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     return
  end if

  ! for some reason the vector attr read methods expect an allocatable array
  allocate(valid_range_alloc(2), stat=alloc_stat)
  if(alloc_stat /= 0) then 
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

  hdf_stat = hdf_sds_attribute_reader(file%id, &
       trim(sds_name), &
       "valid_range", &
       valid_range_alloc)
  if (hdf_stat == sym%FAILURE)  then
     stat = STATUS_ERR_MISSING_ATTRIBUTE
     deallocate(valid_range_alloc)
     return
  end if

  valid_range = valid_range_alloc
  deallocate(valid_range_alloc)

  ! allocate int2 buffer
  allocate( buf_real8(edge(1), edge(2)), stat=alloc_stat)
  if (alloc_stat /= 0) then
     stat = STATUS_ERR_ALLOCATE_FAILED
     return
  end if

!  call set_missing(Buf_Real8)
      
  ! Read data
  hdf_stat = hdf_sds_reader(file%id, &
       trim(sds_name), &
       start, &
       stride, &
       edge, &
       buf_real8)
  if (hdf_stat == sym%FAILURE .or. edge(1) /= subset%edge(1) .or. edge(2) /= subset%edge(2)) then
!     print *, 'Error: Edge=', Edge  
     stat = STATUS_ERR_MISSING_DATA
     deallocate(buf_real8)
     return
  end if

  ! if data is valid, unscale data. We assume Data has already been initialized to missing values  
  where(buf_real8 >= valid_range(1) .and. buf_real8 <= valid_range(2) .and. buf_real8 /= fill_value)
     data(1:edge(1), 1:edge(2)) = (buf_real8(:,:) - add_offset) * scale_factor
  end where

  deallocate(buf_real8)

end subroutine read_sds_float64


end module modis
