MODULE Atmosphere

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

  PRIVATE

  PUBLIC :: Get_Atmosphere_Data,Get_Atmosphere_Bytes
  
  ! 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. 
  !
  ! Extra lines of data can be requested with the option Pad_Lines. The size of
  ! returned data array and optional longitude and latitude arrays is increased
  ! by the number of pad lines that are available.
  ! 
  ! 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
  !  Pad_Lines              number of extra lines of data requested above and
  !                         below the array. 
  !
  ! 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, &
       Pad_Lines, &
       Native_Lon, &
       Native_Lat, &
!       Interp_Type, &
       Stat )

    ! inputs and outputs
    INTEGER(Kind=Int4), INTENT(IN) :: Source
    CHARACTER(Len=*), INTENT(IN) :: Field_Name
    REAL(Kind=Real8), DIMENSION(:,:), POINTER :: Atmosphere_Data  ! F95 does not allow intent attr for pointer dummy args 
!    integer(kind=int4), intent(in), optional :: Interp_Type
    INTEGER(Kind=Int4), INTENT(IN), OPTIONAL :: Pad_Lines   
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, OPTIONAL :: Native_Lon ! ""
    REAL(Kind=Real8), DIMENSION(:,:), POINTER, OPTIONAL :: Native_Lat ! ""
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat

    ! 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

    ! 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
    
!!$    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, &
         Pad_Lines=Pad_Lines)
    IF(Opt_Stat /= STATUS_OK) RETURN

!    if( opt_interp_type == INTERP_TYPE_NONE) then
    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
  ! adaptation of Get_Atmospere_DATA to read Quality Flag Bytes
  SUBROUTINE Get_Atmosphere_Bytes( &
       Source, &
       Field_Name, &
       Atmosphere_Bytes, &
       Fill_Value, &
       Pad_Lines, &
!       Interp_Type, &
       Stat )

    ! inputs and outputs
    INTEGER(Kind=Int4), INTENT(IN) :: Source
    CHARACTER(Len=*), INTENT(IN) :: Field_Name
    Integer(Kind=int1), DIMENSION(:,:,:), POINTER :: Atmosphere_Bytes  ! F95 does not allow intent attr for pointer dummy args 
!    integer(kind=int4), intent(in), optional :: Interp_Type
    INTEGER(kind=int1) :: Fill_Value
    INTEGER(Kind=Int4), INTENT(IN), OPTIONAL :: Pad_Lines   
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat

    ! 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), DIMENSION(3) :: Atmosphere_Data_Res
    TYPE(Subset_2D) :: Subset
    INTEGER(Kind=Int4), DIMENSION(2) :: Scale_Factor
!    write(6,*)'get_atmosphere_bytes ',trim(Field_name)

    ! 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
    
!!$    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 or 3d
    IF(Sds_Rank /= 3.and.Sds_Rank /=2) THEN
       Opt_Stat = STATUS_ERR_ARRAY_SHAPE_INVALID
       RETURN
    END IF
    if(sds_Rank.eq.2)then
      ! 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, &
         Pad_Lines=Pad_Lines)
         IF(Opt_Stat /= STATUS_OK) RETURN
         Atmosphere_Data_Res(1)=1
         
    else

      ! create the HDF subset struct based on current segment and SDS resolution
      ! Dimesion (1) is number of bytes in Quality Flag. Thus nc and nr are in dims(2) and dims(3)
      CALL Get_Subset((/Sds_Dims(2), Sds_Dims(3)/), &
         Subset, &
         Scale_Factor, &
         Opt_Stat, &
         Pad_Lines=Pad_Lines)
         IF(Opt_Stat /= STATUS_OK) RETURN
         Atmosphere_Data_Res(1)=Sds_Dims(1)
    endif

!    if( opt_interp_type == INTERP_TYPE_NONE) then
    Atmosphere_Data_Res(2:3) = Subset%Edge

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

    ! read data from SDS
    CALL Modis_Read_Bytes(Source, &
         TRIM(MODIS_L2_DIRECTORY), &
         TRIM(Field_Name), &
         Subset, &
         Timestr, &
         Atmosphere_Bytes, Fill_value, Sds_Rank,&
         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)


    
  END SUBROUTINE Get_Atmosphere_Bytes

  ! 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, Pad_Lines)
    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), INTENT(IN), OPTIONAL :: Pad_Lines

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

    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
    
    ! x and y scale factors should be the same
    ! Set x and y to the larger scale factor. Should allow any extra lines or rows in 
    ! the SDS to be ignored.
    Scale_Factor = MAX(Scale_Factor(1), Scale_Factor(2))  
    
    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) / 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)

    ! Adjust for pad lines
    IF(PRESENT(Pad_Lines)) THEN
       ! FIXME: check for pad_lines < 0, return STATUS_ERR_OFFSET_INVALID (must be added 
       ! to constant.f90)
       Yend = MIN(Subset%Start(2) + Subset%Edge(2) + Pad_Lines, Sds_Resolution(2))
       Subset%Start(2) = MAX(Subset%Start(2) - Pad_Lines, 0)
       Subset%Edge(2) = Yend - Subset%Start(2)
    END IF
    
  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, Arr, Stat)
    INTEGER(Kind=Int4), DIMENSION(2) :: Dims
    REAL(Kind=Real8), DIMENSION(:,:), POINTER :: Arr    ! F95 does not allow intent attr for pointer dummy args            
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat

    INTEGER(Kind=Int4) :: Alloc_Status

    Stat = STATUS_OK

    ! allocate output array if needed, and clear it
    IF(ASSOCIATED(Arr)) THEN   
       IF(ANY(SHAPE(Arr) /= Dims)) THEN
!          print *, 'deallocating wrong size array...'           
!          write(6,*)' deallocating wrong size array...',shape(arr),dims
          call flush(6)
          DEALLOCATE(Arr)
       END IF
    END IF

    IF(.NOT. ASSOCIATED(Arr)) THEN   
!       print *, 'allocating array...' 
!          write(6,*)' allocating array...',dims
!          call flush(6)
       ALLOCATE(Arr(Dims(1), Dims(2)), Stat=Alloc_Status)
       IF(Alloc_Status /= 0) THEN 
          Stat = STATUS_ERR_ALLOCATE_FAILED
          RETURN
       END IF
    END IF
    
    CALL Set_Missing(Arr)

  END SUBROUTINE Allocate_If_Needed
  SUBROUTINE Allocate_If_Needed_Bytes(Dims, Arr, Stat)
    USE Constants
    INTEGER(Kind=Int4), DIMENSION(3) :: Dims
    Integer(Kind=int1), DIMENSION(:,:,:), POINTER :: Arr    ! F95 does not allow intent attr for pointer dummy args            
    INTEGER(Kind=Int4), INTENT(OUT), OPTIONAL, TARGET :: Stat

    INTEGER(Kind=Int4) :: Alloc_Status

    Stat = STATUS_OK

    ! allocate output array if needed, and clear it
    IF(ASSOCIATED(Arr)) THEN   
       IF(ANY(SHAPE(Arr) /= Dims)) THEN
!          print *, 'deallocating wrong size array...'           
          DEALLOCATE(Arr)
       END IF
    END IF

    IF(.NOT. ASSOCIATED(Arr)) THEN   
!       print *, 'allocating array...' 
       ALLOCATE(Arr(Dims(1), Dims(2),Dims(3)), Stat=Alloc_Status)
       IF(Alloc_Status /= 0) THEN 
          Stat = STATUS_ERR_ALLOCATE_FAILED
          RETURN
       END IF
    END IF
    
    CALL Set_Missing_array_int1(Arr)

  END SUBROUTINE Allocate_If_Needed_Bytes
subroutine set_missing_array_int1(array)
  integer(kind=int1), dimension(:,:,:), intent(inout) :: array
  array = MISSING_VALUE_INT1
end subroutine set_missing_array_int1

END MODULE Atmosphere
