MODULE Atmosphere

  USE Constants
  USE Pixel_Common
  USE Area_Read, ONLY: Areastr
  USE Modis
  
  IMPLICIT NONE

  PRIVATE

  PUBLIC :: Get_Atmosphere_Data
  
  ! FIXME: temporary. need to make this configurable somehow
  CHARACTER(Len=256), PARAMETER :: MODIS_L2_DIRECTORY = './modis_data/'

  CHARACTER(Len=9), PARAMETER :: MODIS_SDS_NAME_LON = 'Longitude'
  CHARACTER(Len=8), PARAMETER :: MODIS_SDS_NAME_LAT = 'Latitude'

CONTAINS
   
  !===============================================================================
  !
  ! Get_Atmosphere_Data
  ! 
  ! Get atmospheric data corresponding to the current segment being processed. 
  ! Currently functionality is limited to retrieving 2-dimensional MODIS L2 data 
  ! (MOD04, MYD04, MOD06 and MYD06), and the current instrument must be MODIS.
  !
  ! Data that is at a resolution lower than the L1 resolution will be returned
  ! at native (low-res) resolution. A limitation of this experimental version 
  ! is that the mapping of low-res data to the current
  ! segment is approximate; in this version it is the caller's 
  ! responsibility to compare the low-res and L1 longitudes and latitudes. In a
  ! future version the low res data will be better matched spatially to the 
  ! L1 data, possibly by applying a resampling algorithm.
  !
  ! The native resolution longitude and latitude can be obtained via the optional
  ! arguments Native_Lon and Native_Lat. The longitude and latitude come from the
  ! L2 file if they are available at the same resolution as the SDS being read, 
  ! otherwise from the L1 file. 
  ! 
  ! The Atmosphere_Data array must have the POINTER attribute. If it has not been
  ! allocated, it will be allocated in Get_Atmosphere_Data and the caller will be 
  ! responsible for deallocation when it is no longer needed. To minimize memory
  ! allocation overhead, it is suggested that the Atmosphere_Data array be 
  ! declared with the SAVE attribute and re-used across successive segments, and 
  ! finally be deallocated at the end of processing the last segment.
  !
  ! It is recommended that the optional Stat argument is supplied and checked to  
  ! detect error conditions. Data should be assumed to be invalid if an error has
  ! occurred.
  !
  ! Inputs: 
  !
  !  Source                 data source, currently only MODIS MXD04 and 06: 
  !                         SOURCE_MODIS_MOD04, SOURCE_MODIS_MYD04, 
  !                         SOURCE_MODIS_MOD06, SOURCE_MODIS_MYD06
  !  Field Name             name of the desired field, i.e. the SDS name in the 
  !                         file
  !
  ! Outputs:
  !
  !  Atmosphere_Data        array of data corresponding to the current segment,
  !                         at L1 or native resolution
  !  Native_Lon (optional)  array of longitude at the native resolution of the SDS
  !  Native_Lat (optional)  array of latitude at the native resolution of the SDS
  !  Stat (optional)        error status: STATUS_OK means no errors or warnings
  !
  !===============================================================================
  ! FUTURE: add a scalar version and create an interface
  SUBROUTINE Get_Atmosphere_Data( &
       Source, &
       Field_Name, &
       Atmosphere_Data, &
       Native_Lon, &
       Native_Lat, &
!       Interp_Type, &
       Lines_Add, &
       Stat )

    ! inputs and outputs
    INTEGER(Kind=Int4), INTENT(IN) :: Source
    CHARACTER(Len=*), INTENT(IN) :: Field_Name
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, INTENT(INOUT) :: Atmosphere_Data 
!    integer(kind=int4), intent(in), optional :: Interp_Type
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, OPTIONAL, INTENT(INOUT) :: Native_Lon
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, OPTIONAL, INTENT(INOUT) :: Native_Lat
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat
    INTEGER(Kind=Int4), OPTIONAL :: Lines_Add

    ! local clones and pointers to optionals
!    integer(kind=int4) :: opt_Interp_Type
    INTEGER(Kind=Int4), POINTER :: Opt_Stat
    INTEGER(Kind=Int4), TARGET :: Fs_Stat

    ! lists of valid options
    INTEGER(Kind=Int4), DIMENSION(4), PARAMETER :: VALID_SOURCES = (/SOURCE_MODIS_MOD04, SOURCE_MODIS_MYD04, SOURCE_MODIS_MOD06, SOURCE_MODIS_MYD06/)
!    integer(kind=int4), dimension(2), parameter :: VALID_INTERP_TYPES = (/INTERP_TYPE_NONE, INTERP_TYPE_SCALE/)
    
    INTEGER(Kind=Int4) :: Sds_Rank
    INTEGER(Kind=Int4), DIMENSION(MAX_RANK_HDF) :: Sds_Dims
    INTEGER(Kind=Int4) :: L2_Lon_Rank
    INTEGER(Kind=Int4), DIMENSION(MAX_RANK_HDF) :: L2_Lon_Dims
    INTEGER(Kind=Int4) :: L2_Lat_Rank
    INTEGER(Kind=Int4), DIMENSION(MAX_RANK_HDF) :: L2_Lat_Dims
    INTEGER(Kind=Int4), DIMENSION(2) :: Atmosphere_Data_Res
    TYPE(Subset_2D) :: Subset
    INTEGER(Kind=Int4), DIMENSION(2) :: Scale_Factor
    logical mod04
    integer(kind=Int4) maxend,maxedge,maxadd

    ! handle optional inputs
!!$    if(present(interp_type)) then
!!$       opt_interp_type = interp_type
!!$    else
!!$       opt_interp_type = INTERP_TYPE_SCALE
!!$    end if

    IF(PRESENT(Stat)) THEN
       Opt_Stat => Stat
    ELSE
       Opt_Stat => Fs_Stat
    END IF
          
    Opt_Stat = STATUS_OK
    
    ! validate inputs
    IF (.NOT. Valid_Argument(Source, VALID_SOURCES)) THEN
       Opt_Stat = STATUS_ERR_SOURCE_INVALID
       RETURN
    END IF
    mod04=.false.
    if(Source.eq.VALID_SOURCES(1).or.Source.eq.VALID_SOURCES(2))then
      mod04=.true.
    endif
    
!!$    if (.not. valid_argument(opt_interp_type, VALID_INTERP_TYPES)) then
!!$       opt_stat = STATUS_ERR_INTERP_TYPE_INVALID
!!$       return
!!$    end if

    ! Determine the SDS resolution 
    CALL Modis_Get_Dims(Source, &
         TRIM(MODIS_L2_DIRECTORY), &
         TRIM(Field_Name), &
         Timestr, &
         Sds_Rank, &
         Sds_Dims, &
         Opt_Stat)
    IF(Opt_Stat /= STATUS_OK) RETURN

    ! SDS must be 2d
    IF(Sds_Rank /= 2) THEN
       Opt_Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
       RETURN
    END IF

    ! create the HDF subset struct based on current segment and SDS resolution
    CALL Get_Subset((/Sds_Dims(1), Sds_Dims(2)/), &
         Subset, &
         Scale_Factor, &
         Opt_Stat)
    IF(Opt_Stat /= STATUS_OK) RETURN

!    if( opt_interp_type == INTERP_TYPE_NONE) then
    if(mod04.and.present(Lines_Add))then
!     increment sizes of direction 2
!     add to end first
!     determine maximum end in j (line) direction that is within original grid
      maxend=min(sds_dims(2)-1,Subset%start(2)+Subset%edge(2)-1+Lines_Add)
!     determine max edge from maxend
      maxedge=maxend-Subset%start(2)+1
      Subset%edge(2)=maxedge
!     now add to beginning if possible
!     determine maximum can add to start
      maxadd=min(Lines_Add,Subset%start(2))
!     decrease start and increment edge
      
      Subset%start(2)=Subset%start(2)-maxadd
      Subset%edge(2)=Subset%edge(2)+maxadd
    endif
    Atmosphere_Data_Res = Subset%Edge

    ! allocate output array if needed, and clear it
    CALL Allocate_If_Needed(Atmosphere_Data_Res, Atmosphere_Data, Opt_Stat)
    IF( Opt_Stat /= STATUS_OK) RETURN

    ! read data from SDS
    CALL Modis_Read(Source, &
         TRIM(MODIS_L2_DIRECTORY), &
         TRIM(Field_Name), &
         Subset, &
         Timestr, &
         Atmosphere_Data, &
         Opt_Stat)
    IF(Opt_Stat /= STATUS_OK) RETURN

    ! upscale data if needed
!    if(opt_interp_type == INTERP_TYPE_SCALE) call upscale_data(subset%edge, scale_factor, Atmosphere_Data, Stat)

    ! get native lon if present
    IF( PRESENT(Native_Lon)) THEN

       ! allocate lon array if needed, same res as data
       CALL Allocate_If_Needed(Atmosphere_Data_Res, Native_Lon, Opt_Stat)
       IF( Opt_Stat /= STATUS_OK) RETURN

       ! first choice is the longitude variable in the L2 file. Get the dims.
       CALL Modis_Get_Dims(Source, &
            TRIM(MODIS_L2_DIRECTORY), &
            TRIM(MODIS_SDS_NAME_LON), &
            Timestr, &
            L2_Lon_Rank, &
            L2_Lon_Dims, &
            Opt_Stat)
       IF(Opt_Stat /= STATUS_OK) RETURN

       ! SDS must be 2d
       IF(L2_Lon_Rank /= 2) THEN
          Opt_Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
          RETURN
       END IF

       ! if the L2 lon is the same shape as the SDS, get it from the file
       IF(Sds_Dims(1) == L2_Lon_Dims(1) .AND. Sds_Dims(2) == L2_Lon_Dims(2)) THEN
          
          CALL Modis_Read(Source, &
               TRIM(MODIS_L2_DIRECTORY), &
               TRIM(MODIS_SDS_NAME_LON), &
               Subset, &
               Timestr, &
               Native_Lon, &
               Opt_Stat)
          IF(Opt_Stat /= STATUS_OK) RETURN

       ! otherwise, if the L1 lon is the same shape, use it
       ELSEIF (Sds_Dims(1) == Areastr%Num_Elem .AND. Sds_Dims(2) == Areastr%Num_Line) THEN
          
          Native_Lon = Sat%Lon
          
          ! else no matching lon is available: error
       ELSE
          
          Opt_Stat = STATUS_ERR_MISSING_DATA
          RETURN
          
       ENDIF

    END IF

    ! get native lat if present
    IF( PRESENT(Native_Lat)) THEN

       ! allocate lat array if needed, same res as data
       CALL Allocate_If_Needed(Atmosphere_Data_Res, Native_Lat, Opt_Stat)
       IF( Opt_Stat /= STATUS_OK) RETURN

       ! first choice is the latitude variable in the L2 file. Get the dims.
       CALL Modis_Get_Dims(Source, &
            TRIM(MODIS_L2_DIRECTORY), &
            TRIM(MODIS_SDS_NAME_LAT), &
            Timestr, &
            L2_Lat_Rank, &
            L2_Lat_Dims, &
            Opt_Stat)
       IF(Opt_Stat /= STATUS_OK) RETURN

       ! SDS must be 2d
       IF(L2_Lat_Rank /= 2) THEN
          Opt_Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
          RETURN
       END IF

       ! if the L2 lat is the same shape as the SDS, get it from the file
       IF(Sds_Dims(1) == L2_Lat_Dims(1) .AND. Sds_Dims(2) == L2_Lat_Dims(2)) THEN
          
          CALL Modis_Read(Source, &
               TRIM(MODIS_L2_DIRECTORY), &
               TRIM(MODIS_SDS_NAME_LAT), &
               Subset, &
               Timestr, &
               Native_Lat, &
               Opt_Stat)
          IF(Opt_Stat /= STATUS_OK) RETURN

       ! otherwise, if the L1 lat is the same shape, use it
       ELSEIF (Sds_Dims(1) == Areastr%Num_Elem .AND. Sds_Dims(2) == Areastr%Num_Line) THEN          
          Native_Lat = Sat%Lat
          
          ! else no matching lat is available: error
       ELSE
          Opt_Stat = STATUS_ERR_MISSING_DATA
          RETURN
       ENDIF

    END IF

    
  END SUBROUTINE Get_Atmosphere_Data

  ! Get the subset in the SDS that corresponds to the current segment, based on
  ! the resolution of the SDS. 
  SUBROUTINE Get_Subset(Sds_Resolution, Subset, Scale_Factor, Stat)
    INTEGER(Kind=Int4), DIMENSION(2), INTENT(IN) :: Sds_Resolution
    TYPE(Subset_2D), INTENT(OUT) :: Subset
    INTEGER(Kind=Int4), DIMENSION(2), INTENT(OUT) :: Scale_Factor
    INTEGER(Kind=Int4), INTENT(OUT) :: Stat

    INTEGER(Kind=Int4), DIMENSION(2) :: L1_Resolution

    Stat = STATUS_OK

    L1_Resolution = (/Areastr%Num_Elem, Areastr%Num_Line/)
    
    ! Currently SDS res higher than L1 res is not supported. If needed, should be
    ! easy to add
    IF(ANY(Sds_Resolution > L1_Resolution)) THEN
       Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
       RETURN
    END IF

    Scale_Factor = L1_Resolution / Sds_Resolution  ! integer division
    
    ! aspect ratios must be the same
    IF(Scale_Factor(1) /= Scale_Factor(2)) THEN
       Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
       RETURN
    END IF
    
    IF(MODULO(Sat%Xstart - 1, Scale_Factor(1)) /= 0 .OR. &
         MODULO(Sat%Ystart - 1, Scale_Factor(2)) /= 0) THEN 
       Stat = STATUS_ERR_SEG_START_INVALID
       RETURN
    END IF
    
    IF(MODULO(Sat%Nscans_Per_Segment, Scale_Factor(2)) /= 0)THEN
       Stat = STATUS_ERR_SEG_SIZE_INVALID
       RETURN
    END IF

    ! Set HDF subsetting parameters
!!$    subset%Start(1) = sat%Xstart - 1
!!$    subset%Start(2) = sat%Ystart - 1
!!$    subset%Stride(1) = sat%Xstride
!!$    subset%Stride(2) = 1
!!$    subset%Edge(1) = sat%Nx
!!$    subset%Edge(2) = sat%Ny

    Subset%Start(1) = (Sat%Xstart - 1) / Scale_Factor(1)
    Subset%Start(2) = (Sat%Ystart - 1) / Scale_Factor(2)
    Subset%Stride(1) = Sat%Xstride
    Subset%Stride(2) = 1
    Subset%Edge(1) = Sat%Nx / Scale_Factor(1)
    Subset%Edge(2) = Sat%Ny / Scale_Factor(2)

  END SUBROUTINE Get_Subset

  ! assumes scale_factor >= 1
  SUBROUTINE Upscale_Data(Unscaled_Res, Scale_Factor, DATA, Stat)
        INTEGER(Kind=Int4), DIMENSION(2), INTENT(IN) :: Unscaled_Res
        INTEGER(Kind=Int4), DIMENSION(2), INTENT(IN) :: Scale_Factor
        REAL(Kind=Real8), DIMENSION(:,:), INTENT(INOUT) :: DATA 
        INTEGER(Kind=Int4), INTENT(OUT) :: Stat
        
        INTEGER(Kind=Int4) :: i, j

        Stat = STATUS_OK

        ! check data array is large enough
        IF( ANY((Unscaled_Res * Scale_Factor) > SHAPE(DATA))) THEN
           Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
           RETURN
        END IF

        DO j = Unscaled_Res(2), 1, -1
           DO i = Unscaled_Res(1), 1, -1
              DATA((i-1) * Scale_Factor(1) + 1 : i * Scale_Factor(1), (j-1) * Scale_Factor(2) + 1 : j * Scale_Factor(2)) = DATA(i,j)
           END DO
        END DO

  END SUBROUTINE Upscale_Data


  ! FUTURE: this routine is also in surface.f90. Consider promoting to a common location (interface_utils.f90?)
  ! return true if the argument is in the list of valid arguments
  FUNCTION Valid_Argument(Arg, Valid_Args) RESULT(Is_Valid)
    INTEGER(Kind=Int4), INTENT(IN) :: Arg
    INTEGER(Kind=Int4), DIMENSION(:), INTENT(IN) :: Valid_Args
    LOGICAL :: Is_Valid
    Is_Valid = ANY(Arg == Valid_Args)
  END FUNCTION Valid_Argument
 
  ! FUTURE: this routine could also be promoted to a common location, and surface.f90 could be modified to call it.
  SUBROUTINE Allocate_If_Needed(Dims, DATA, Stat)
    INTEGER(Kind=Int4), DIMENSION(2) :: Dims
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, INTENT(INOUT) :: DATA 
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat

    INTEGER(Kind=Int4) :: Alloc_Status

    Stat = STATUS_OK

    ! allocate output array if needed, and clear it
!   check size and deallocate if not correct size
    IF(ASSOCIATED(DATA)) THEN   
     IF(ANY(SHAPE(DATA) .ne. Dims)) THEN
       deallocate(data)
     endif
    END IF 
    IF(ASSOCIATED(DATA)) THEN   
       CALL Set_Missing(DATA)
       IF(ANY(SHAPE(DATA) < Dims)) THEN
          Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
          RETURN
       END IF
    ELSE

!       print *, 'allocating data array...'

       ALLOCATE(DATA(Dims(1), Dims(2)), Stat=Alloc_Status)
       IF(Alloc_Status /= 0) THEN 
          Stat = STATUS_ERR_ALLOCATE_FAILED
          RETURN
       END IF
       CALL Set_Missing(DATA)
    END IF

  END SUBROUTINE Allocate_If_Needed

END MODULE Atmosphere