modis_albedo.f90

Go to the documentation of this file.

module modis_albedo
 
! HI
!****************************************************************************
  ! !F90
  !
  ! !Description:
  !    This module contains the subroutines used to provide albedo and ecosystem
  !     maps for processing a MOD06 granule.
  !    There is only one callable routine, getAlbedoEco, which returns albedo 
  !     values and IGBP ecosystem classification for each pixel of the 
  !     MOD06 granule,
  !     and for the wavelengths specified.  Also returned are values of 
  !     snow albedos
  !     by ecosystem class for the wavelengths specified.
  !  
  ! !Callable routines:
  !    getAlbedoEco()
  ! 
  ! !Revision History:
  !
  ! Revision 1.0  2003/12/18  12:43:43  EGMoody
  ! Initial revision.
  !
  ! !Team-Unique Header:
  !   Cloud Retrieval Group, NASA Goddard Space Flight Center
  !
  ! !References and Credits:
  !   Written by
  !    Eric Moody
  !    Climate and Radiation Branch, Code 913
  !    NASA/GSFC
  !    Greenbelt MD 20771
  !    Eric.Moody@gsfc.nasa.gov
  !
  ! !Design Notes:
  !
  ! !END
  !*****************************************************************************
 
 
  !Dependencies:
  use generalauxtype
  use nonscience_parameters
  use core_arrays, only : cloudmask
  use mod06albedoecomodule
  use mod06_run_settings, only : set_albedo_bands
  implicit none
 
  private
 
  public :: getalbedoeco, init_snow_stats, lat_start, lat_end
 
  integer, parameter :: MAX_VAR_DIMS = 10
  integer :: lat_start, lat_end
 
  !local variables:
  !counters:
  integer  (kind = integer_fourbyte)  :: i,j,k,l,m,n
 
  !HDF variables:
  integer  (kind = integer_fourbyte)  :: hdfstatus
  integer  (kind = integer_fourbyte  ),  &
             dimension(MAX_VAR_DIMS)     :: hdfstart, hdfstride, hdfedge
  integer  (kind = integer_fourbyte)     :: sds_id, sds_index
 
  integer  (kind = integer_fourbyte), parameter  :: numsnowalbstatwaves = 10
  integer  (kind = integer_fourbyte), parameter  :: numecosystems = 18
  integer  (kind = integer_fourbyte), parameter  ::  numsnowtypes = 2
 
  real :: snowalbedostats (1:numsnowalbstatwaves, 0:numecosystems-1, &
          1:numsnowtypes)
  real :: processsnowalbstats (1:set_albedo_bands,  0:numecosystems-1, &
          1:numsnowtypes)
 
  integer  (kind = integer_fourbyte), parameter  :: numseaicealbwaves = 8
  real     (kind = single), dimension(NumSeaIceAlbWaves) :: meltseaicealbedos, &
            dryseaicealbedos
  real (kind = single), dimension(set_albedo_bands) :: &
    processmeltseaicealbedos, processdryseaicealbedos,  transitionalbedo, & 
    slope, intercept
 
  real   (kind = single), parameter :: percentageofdry = 0.800d0
  real :: wavelength(10)
 
contains
 
 
  subroutine init_snow_stats(SnowAlbedoFN, WavelengthText) 
    
  character (len = *), intent(in ) :: snowalbedofn
  character (len = *),  dimension(:), intent(in ) :: wavelengthtext
        
  real (kind = single), dimension(NumSnowAlbStatWaves) :: possstatwavelengths
  real (kind = single), dimension(NumSeaIceAlbWaves) :: possseaicewavelengths
 
  integer :: errorlevel, i, numalbwavesprocess
  logical :: foundwave
 
  numalbwavesprocess = set_albedo_bands
 
  !*****************************************************************************
  ! Read in the snow albedo statistics from the corresponding file:
  !*****************************************************************************
  call readsnowalbstats  ( snowalbedofn, numsnowtypes, numsnowalbstatwaves,  &
    numecosystems, snowalbedostats, errorlevel )
 
  !*****************************************************************************
  ! For each wavelength to be processed, store the albedo snow statistic in the
  !  cooresponding index in a temp array. Note that 8-10 are broad band, and 
  !  not included.
  !*****************************************************************************
  !Load in the order of the statistical wavelengths from the file:
  possstatwavelengths(1)  = 0.470000
  possstatwavelengths(2)  = 0.555000
  possstatwavelengths(3)  = 0.659000
  possstatwavelengths(4)  = 0.858000
  possstatwavelengths(5)  = 1.240000
  possstatwavelengths(6)  = 1.640000
  possstatwavelengths(7)  = 2.130000
  possstatwavelengths(8)  = 0.0 
  possstatwavelengths(9)  = 0.0
  possstatwavelengths(10) = 0.0
    
  !Loop over the wavelengths to be processed, storing the cooresponding
  ! data from the complete snowAlbedoStats array.  Do not process 3.7,
  ! the last index as there are no cooresponding values, instead compute
  ! from the 2.1µm at the end:
    
  do i = 1, numalbwavesprocess
    read( wavelengthtext(i), * )  wavelength(i) 
  end do
   
  processsnowalbstats= 0.
  do i = 1, numalbwavesprocess-1
    foundwave = .false.
 
    do j = 1, numsnowalbstatwaves-3
      if ( mod(wavelength(i),possstatwavelengths(j)) < 0.05 ) then
        processsnowalbstats(i,:,:) = snowalbedostats(j,:,:)
        foundwave = .true.
        goto 100
      end if
    end do
 100 continue
    if (.not. foundwave) then
      return
    end if
  end do
    
  !compute the 3.7µm (assume 6 = 3.6 and 5 = 2.1µm):
  processsnowalbstats(6,:,:) = processsnowalbstats(5,:,:) * 0.5000
                        
  !*****************************************************************************
  ! For each wavelength to be processed, store the albedo sea ica albedo in the
  !  cooresponding index in a temp array. 
  !*****************************************************************************
  !Load in the order of the statistical wavelengths from the file:
  possseaicewavelengths(1) = 0.470000
  possseaicewavelengths(2) = 0.555000
  possseaicewavelengths(3) = 0.659000
  possseaicewavelengths(4) = 0.858000
  possseaicewavelengths(5) = 1.240000
  possseaicewavelengths(6) = 1.640000
  possseaicewavelengths(7) = 2.130000
  possseaicewavelengths(8) = 3.700000
    
  ! Set the sea ice values to be equal to the dry perm snow/ice 
  ! values or percentage of them:
  dryseaicealbedos(1:7)  = snowalbedostats(1:7,15,1)
  meltseaicealbedos(1:7) = snowalbedostats(1:7,15,1) * percentageofdry
  !3.7 is half of 2.1:
  dryseaicealbedos(8)  =  dryseaicealbedos(7) * 0.5000d0
  meltseaicealbedos(8) = meltseaicealbedos(7) * 0.5000d0
                         
  !Loop over the wavelengths to be processed, storing the cooresponding
  ! data from the complete snowAlbedoStats array.  Process 3.7µm for this
  ! sea ice case
  do i = 1, numalbwavesprocess
    foundwave = .false.
    do j = 1, numseaicealbwaves
       if ( mod(wavelength(i),possseaicewavelengths(j)) < 0.05 ) then
          processdryseaicealbedos(i)  = dryseaicealbedos(j)
          processmeltseaicealbedos(i) = meltseaicealbedos(j)
          foundwave = .true.
          goto 101
       end if
    end do
  101 continue
    if (.not. foundwave) then
!       status = 19
!       return
    end if
  end do
 
  end subroutine init_snow_stats
 
! --------------------------------------------------------------------------
! --------------------------------------------------------------------------
  subroutine getalbedoeco (Lat, Long, JulianDay, EcosystemFN,  &
    emissivity_name, Debug, WavelengthText, ProcessLandAlbedo,  &
    icec, Albedo, emissivity, Ecosystem, Status, sfc_info )
 
  ! WDR access c2_alt_snowice and albedo
  use ch_xfr, only: c2_alt_snowice, c2_sfc_albedo
 
  ! EM -- June 17, 2005
  ! Apparently ncdump does not handle the metadata very well, and there 
  ! was some binary conversion 
  ! issue such that the metadata values are mixed up.  If you use something 
  ! like HDFlook, they are fine.
  ! Here is how ncdump prints them out:
  !  byte IGBP_Land_Cover_Type(NumLatPoints, NumLongPoints) ;
  !    IGBP_Land_Cover_Type:long_name = "IGBP Land Cover Classification Type" ;
  !    IGBP_Land_Cover_Type:units = "class number" ;
  !    IGBP_Land_Cover_Type:valid_range = '\0', '\376' ;
  !    IGBP_Land_Cover_Type:_FillValue = '\377' ;
  !    IGBP_Land_Cover_Type:water = '\0' ;
  !    IGBP_Land_Cover_Type:evergreen needleleaf forest = '\1' ;
  !    IGBP_Land_Cover_Type:evergreen broadleaf forest = '\2' ;
  !    IGBP_Land_Cover_Type:deciduous needleleaf forest = '\3' ;
  !    IGBP_Land_Cover_Type:deciduous broadleaf forest = '\4' ;
  !    IGBP_Land_Cover_Type:mixed forests = '\5' ;
  !    IGBP_Land_Cover_Type:closed shrubland = '\6' ;
  !    IGBP_Land_Cover_Type:open shrublands = '\7' ;
  !    IGBP_Land_Cover_Type:woody savannas = '\10' ;
  !    IGBP_Land_Cover_Type:savannas = '\11' ;
  !    IGBP_Land_Cover_Type:grasslands = '\12' ;
  !    IGBP_Land_Cover_Type:permanent wetlands = '\13' ;
  !    IGBP_Land_Cover_Type:croplands = '\14' ;
  !    IGBP_Land_Cover_Type:urban and built-up = '\15' ;
  !    IGBP_Land_Cover_Type:cropland/natural vegetation mosaic = '\16' ;
  !    IGBP_Land_Cover_Type:snow and ice = '\17' ;
  !    IGBP_Land_Cover_Type:barren or sparsely vegetated = '\20' ;
  !    IGBP_Land_Cover_Type:unclassified = '\376' ;
  ! If you notice it skips 7 and goes to 10 (open shrub to woody savanna) 
  ! and then 17 to 20 (snow/ice to barren).
  !
  ! Here is how my code actually assigns them:
  !
  !*****************************************************************************
  ! Create the Ecosystem_Order array, used in HDF initialization:
  !*****************************************************************************
  !    Ecosystem_Order(0)  = 'Ocean / Water'
  !    Ecosystem_Order(1)  = 'Evergreen Needle Forest'
  !    Ecosystem_Order(2)  = 'Evergreen Broad Forest'
  !    Ecosystem_Order(3)  = 'Deciduous Needle Forest'
  !    Ecosystem_Order(4)  = 'Deciduous Broad Forest'
  !    Ecosystem_Order(5)  = 'Mixed Forest'
  !    Ecosystem_Order(6)  = 'Closed Shrubs'
  !    Ecosystem_Order(7)  = 'Open / Shrubs'
  !    Ecosystem_Order(8)  = 'Woody Savanna'
  !    Ecosystem_Order(9)  = 'Savanna'
  !    Ecosystem_Order(10) = 'Grassland'
  !    Ecosystem_Order(11) = 'Wetland'
  !    Ecosystem_Order(12) = 'Cropland'
  !    Ecosystem_Order(13) = 'Urban'
  !    Ecosystem_Order(14) = 'Crop Mosaic'
  !    Ecosystem_Order(15) = 'Antarctic / Permanant Snow'
  !    Ecosystem_Order(16) = 'Barren / Desert'
  !    Ecosystem_Order(17) = 'Tundra'   
  !
  !Ncdump is quite deficient in this respect and I get asked this all the time.
  !But it is not a coding issue. -- E. Moody
 
  real      (kind = single ), dimension(:,:), intent(in ) :: lat, long
  character (len = *), intent(in ) :: ecosystemfn
  character (len = *), intent(in ) :: emissivity_name
  integer   (kind = integer_fourbyte ),  intent(in ) :: julianday
  logical,  intent(in ) :: debug
  character (len = *), dimension(:),  intent(in ) :: wavelengthtext
  logical,  dimension(:), intent(in ) :: processlandalbedo
  real      (kind = single), dimension(:,:),  intent(in ) :: icec 
  real (kind = single ), dimension(:,:,:),  intent(out) :: emissivity
  integer*2, dimension(:,:,:), intent(out) :: albedo 
  integer*1, dimension(:,:), intent(out) :: ecosystem, sfc_info
  integer   (kind = integer_fourbyte), intent(out) :: status
 
  ! !Description:
  !    This returns albedo values and IGBP ecosystem classification for each 
  !    pixel 
  !    of the MOD06 granule, and for the wavelengths specified.  Also 
  !    returned are 
  !    values of snow albedos by ecosystem class for the wavelengths specified.
  !  
  ! !Input Parameters:
  !    Lat, Long       : Latitude and Longitude arrays from granule.
  !    WavelengthFN    : The albedo filenames for each specified wavelength 
  !                      to process.
  !    EcosystemFN     : The filename for the IGBP Ecosystem Classification map.
  !    Julian Day      : Julian Day
  !    Debug           : Logical Flag for printing during debugging.
  !    WavelengthText  : Text array of wavelengths to be processed.
  !
  ! !Output Parameters:
  !    Wavelength      : Real values of the Text array.
  !    Albedo          : 3D array containing the albedo values for each 
  !                     lat/lon point
  !                      and for each wavelength specified.
  !    SnowAlbedo      : Snow albedos by ecosystem class for each 
  !                       wavelength specified.
  !    Ecosystem       : IGBP ecosystem classification for each lat/lon point.
  !    Status          : Flag specifying if this routine worked correctly.
  !
  ! !Revision History:
  !   See Module revision history at the beginning of the file.
  !
  ! !Team-Unique Header:
  !   Cloud Retrieval Group, NASA Goddard Space Flight Center
  !
  ! !References and Credits:
  !   Written by
  !    Eric Moody
  !    Climate and Radiation Branch, Code 913
  !    NASA/GSFC
  !    Greenbelt MD 20771
  !    Eric.Moody@gsfc.nasa.gov
  !
  ! !Design Notes:
  !  Albedos values for each lat/lon point are provided by finding 
  !  the nearest neighbor
  !    from the ancillary albedo maps.
  !  Wavelengths to be processed are specified through the WavelengthText 
  !  array and
  !    corresponding WavelengthFN array.
  !  There is no ancillary files for 3.7 micron, so this is assumed to be 
  !  half the 2.1 value.
  !  The 3.7 micron always has to be last in the WavelengthText array and 
  !   WavelengthFN array,
  !    however the WavelengthFN for 3.7 is not used, so it can be fill.
  !  As the albedo ancillary data is for land only, water values are 
  !  assumed to be 5% for
  !    all wavelengths.
  !  Sea ice is approximated with perm snow/ice values for a winter time, 
  !  and a percentage 
  !    during a melt season.  A transition period is approximated between 
  !  winter and melt.
  !  
  !  The following status error values are used:
  !  1 = Array Sizes not valid
  !  2 = Invalid min/max long/lat
  !  3 = Invalid global map start/stop x/y
  !  4 = Could not gain access to Albedo HDF file
  !  5 = Could not gain access to Albedo SDS
  !  6 = Could not read Albedo SDS
  !  7 = Could not end access to Albedo SDS
  !  8 = Could not end access to Albedo HDF file
  !  9 = Could not gain access to Ecosystem HDF file
  !  10 = Could not gain access to Ecosystem SDS
  !  11 = Could not read Ecosystem SDS
  !  12 = Could not end access to Ecosystem SDS
  !  13 = Could not end access to Ecosystem HDF file
  !  14 = Trouble with Albedo Snow Statistics
  !  15 = Trouble with Reading NISE data
  !  16 = Pixel's computed global map index invalid
  !  17 = Pixel's geolocation invalid
  !  18 = Wavelengths do not match for snow statistics
  !  19 = Wavelengths do not match for sea ice values
  !
  !
  ! !END
 
  !Local Variables
 
  !
  !Size variables:
  !
  !Number of Albedo Wavelengths to be processed.
  integer  (kind = integer_fourbyte) :: numalbwavesprocess
  !Number of Columns and Rows of the granule to be processed
  integer  (kind = integer_fourbyte) :: numberofcols, numberofrows
 
  !
  !Min/Max Lat/Lon:
  !
  real     (kind = single) :: minlat, maxlat, minlong, maxlong
 
  ! 
  ! Albedo/Ecosystem map variables:
  !
  ! Global map dims and maximum number of albedo wavelengths:
  integer  (kind = integer_fourbyte)  :: nummapcols, nummaprows
  integer  (kind = integer_fourbyte)  :: nummapcols_alb, nummaprows_alb
                                                                 
  ! Stores the resolution (in degrees) of the global map
  real     (kind = single)    :: map_resolution
  real     (kind = single)  :: map_resolution_alb
  ! Stores the upper right corner global map coords, used to compute rest 
  ! of map coords.
  real  (kind = single)   :: westernmostlongitude, northernmostlatitude
  real  (kind = single)   :: westernmostlongitude_alb, northernmostlatitude_alb
  ! Indices for section of global map that is needed for supplying values 
  ! for granule:
  integer  (kind = integer_fourbyte)  :: startx, stopx, starty, stopy, &
    numberofxpoints,  numberofypoints  
  integer  (kind = integer_fourbyte) :: startx_alb, stopx_alb,  &
    starty_alb, stopy_alb
  real (kind = single) :: startxval, stopxval, startyval, stopyval 
  ! Stores the real and integer values of the section of the global albedo map.
  ! Used to scale the albedo values read in as integers from the HDF files.
  real  (kind = single ), parameter :: scalefactor = 0.00100
  ! Store the NISE snow type:
  INTEGER(integer_fourbyte), allocatable, DIMENSION(:,:)    :: snowicetype
  ! Store the section of the global ecosystem map needed for supplying 
  ! `values for granule:
  integer (kind = integer_onebyte), allocatable, dimension(:,:) :: ecosystemmap 
  ! Dummy index variables:
  integer  (kind = integer_fourbyte) :: xindex, yindex 
  integer  (kind = integer_fourbyte) :: xindex_alb, yindex_alb 
  ! HDF variables:
  ! File IDs
  integer  (kind = integer_fourbyte) :: ecomapfid, albmapfid
  ! SDS IDs
  integer  (kind = integer_fourbyte) :: ecomapsdsid, albmapsdsid    
  ! Stores the SDS name
  character(len =  MAX_SDS_NAME_LEN) :: sdsname
 
  ! Error handling:
  INTEGER(integer_fourbyte)    :: errorlevel
 
  ! Wavelength descriptions used in snow statistics:
  !Variables used for sea ice albedos:
  LOGICAL :: arcticmeltingseason, wintertomelttransition, melttowintertransition
  integer  (kind = integer_fourbyte), parameter :: transitiondays = 10, &
    nhmeltbeg = 152, nhmeltend = 244, shmeltbeg = 334, shmeltend = 61
    
  !Water albedo values:
  real  (kind = single), parameter  :: wateralbedo = 0.050000
  real, dimension(:), allocatable :: albedo_real4
  real, parameter :: emissivity_fill = 0.985 ! according to Wisconsin
    
  ! emissivity map variables
  real, parameter :: emissivity_res = 0.05
  integer, parameter :: num_emiss_cols = 7200
  integer, parameter :: num_emiss_rows = 3600
  real, dimension(:,:,:), allocatable :: emissivity_map
  integer*2, dimension(:,:), allocatable :: emissivity_map_int
  integer :: emiss_x, emiss_y
  real :: emiss_west, emiss_north
    
  integer :: file_id, var_id, err_code, start(2), stride(2), edge(2)
    
  ! WDR definitions needed in replacing ftn hdf i/o package
  integer c_sfstart
  external c_sfstart
  integer c_sfselect
  external c_sfselect
  integer c_sfn2index
  external c_sfn2index
  integer c_sfginfo
  external c_sfginfo
  integer c_sfrdata
  external c_sfrdata
  integer c_sfendacc
  external c_sfendacc
  integer c_sfend
  external c_sfend
  integer :: dfacc_read = 1 
  integer :: fail = -1
 
  ! WDR this will support the re-read only when lat, long range falls outside 
  ! the previous range of min, max (lat, lon)
  real, dimension(2) :: min_geo_sav = (/ -9000., -9000. /), &
                         max_geo_sav = (/ -9000., -9000. /)
  real :: geo_margin = 5.  ! 5 degree extra margin to allocate for
 
  !*****************************************************************************
  ! Determine the size of the arrays:  
  !*****************************************************************************
  !Determine dimensions:
  numberofcols      = size( albedo(:,:,:),   dim = 1 )
  numberofrows      = size( albedo(:,:,:),   dim = 2 )
  numalbwavesprocess = size( albedo(:,:,:),   dim = 3 )  
    
  if( .not. allocated( albedo_real4) ) &
    allocate(albedo_real4(numalbwavesprocess))
    
    !print*,'c,R ',NumberOfCols,NumberOfRows
    !Error checking:
    if ( (numberofcols < 1)      .or. &
      (numberofrows < 1)      .or. &
      (numalbwavesprocess < 2) .or. &
      (numecosystems < 1)          ) then
      status = error
      call local_message_handler( &
        'Problem with albedo array size, contact SDST',status,'getAlbedoEco')
      return     
    end if
 
    !************************************************************************
    ! Determine the min/max of lat/lon, which will be used to determine what 
    ! portion
    !  of the albedo/eco maps need to be read in. 
    !************************************************************************
    ! Determine min and max:
    if (maxval(lat) == -999. .or. maxval(long) == -999.) then
      ! there is no gelocation data in the scan.
      albedo(:,:,:) = -999.  
       emissivity(:,:,:) = -999.
      return
    else
      minlat  = minval( lat,  &
        mask = ( (lat  >= -90.0000 .and. lat  <= 90.00000) ) )
      maxlat  = maxval( lat,  &
        mask = ( (lat  >= -90.0000 .and. lat  <= 90.00000) ) )
      minlong = minval( long, &
        mask = ( (long >= -180.000 .and. long <= 180.0000) ) )
      maxlong = maxval( long, &
        mask = ( (long >= -180.000 .and. long <= 180.0000) ) )
    endif
 
    ! WDR see if we need to re-read the grids for geo bounds outside current
    !PRINT*, 'Decide whether to read the albedo/emiss tables'
    if( ( minlat < min_geo_sav(1) ) .or. ( maxlat > max_geo_sav(1) ) .or. &
        ( minlong < min_geo_sav(2) ) .or. ( maxlong > max_geo_sav(2) ) ) then
      !PRINT*, 'Reading albedo/emiss tables'
      !PRINT*, 'minLat, maxLat, minLong, maxLong ', minLat, maxLat, minLong, maxLong
      !PRINT*, 'min_geo_sav, max_geo_sav ', min_geo_sav, max_geo_sav
 
      ! set the new margin-enhanced range and re-read from files
      minlat = minlat - geo_margin
      maxlat = maxlat + geo_margin
      minlong = minlong - geo_margin
      maxlong = maxlong + geo_margin
 
      if( minlat < -90. ) minlat = -90.
      if( maxlat > 90. ) maxlat = 90.
      if( minlong < -180. ) minlong = -180.
      if( maxlong > 180. ) maxlong = 180.
 
      min_geo_sav(1) = minlat
      min_geo_sav(2) = minlong
      max_geo_sav(1) = maxlat
      max_geo_sav(2) = maxlong
      !PRINT*, 'NEW min_geo_sav, max_geo_sav ', min_geo_sav, max_geo_sav
 
  
      !**********************************************************************
      ! Set up the global albedo/ecosystem map resolution and determine the 
      ! start/stop
      !  points for the portion of the global map to be read in.  This 
      ! is done by
      !  computing the map indices from the min/max lat/lon values.
      !**********************************************************************
      ! 9-30-11 For Collection 4 (Moody) surface albedo and ecosystem maps 
      !  were identical resolution,
      !  the resolution of the C5 SA dataset however is double C4, thus 
      ! appropriate modifications  
      !   have been made below to read the higher resolution Albedo map 
      ! data. - GTA
      
      nummapcols = 21600
      nummaprows = 10800
        
      map_resolution = 360.000d0 / float(nummapcols)
      
      ! Compute the WesternMostLongitude and NothernMostLatitude of the maps:
      westernmostlongitude = map_resolution / 2.0 - 180.0
      northernmostlatitude = 90.0 - map_resolution / 2.0
      !print*,'mapres=',Map_Resolution,westernMostLongitude,northernMostLatitude
      !print*,'minLong=',minLong,maxLong,minLat,maxLat
     
      ! Because the albedo map res. is different than the ecosystem map, 
      ! save the albedo map
      !  parameters for later determination surface albedo for each 
      ! granule pixel  - GTA
        
      ! Compute the start/stop x/y indices based on min/max lat/lon values:
      ! use anint so that it rounds to nearest index.
      startx = anint( (minlong-westernmostlongitude) / (map_resolution) ) + 1
      stopx  = anint( (maxlong-westernmostlongitude) / (map_resolution) ) + 1
      starty = anint( (northernmostlatitude-maxlat ) / (map_resolution) ) + 1
      stopy  = anint( (northernmostlatitude-minlat ) / (map_resolution) ) + 1
      
      ! Due to roundoff, geocoordinates of -90, 90, -180, 180 can result in 
      ! indices of one beyond the map bounds, so correct that here:
      if (startx == 0) startx = 1
      if (stopx  == nummapcols+1 ) stopx = nummapcols
      if (starty == 0) starty = 1
      if (stopy  == nummaprows+1 ) stopy = nummaprows
      
      ! Compute the number of points:
      numberofxpoints = stopx - startx + 1
      numberofypoints = stopy - starty + 1
  
      ! Error Checking:
      if ( (startx < 1) .or. (stopx > nummapcols) .or. &
        (starty < 1) .or. (stopy > nummaprows) .or. &
        (startx > stopx) .or. (starty > stopy)      ) then
        ! status = error
        ! call local_message_handler('Problem with loop index, contact SDST', &
        !  status,'getAlbedoEco')
        ! return  
         if (stopx > nummapcols) stopx = nummapcols
        if (stopy > nummaprows) stopy = nummaprows
        if (startx > stopx) startx = stopx-1
        if (starty > stopy) starty = stopy-1
        if (startx < 1) startx = 1
        if (starty < 1) starty = 1
        
      end if
  
      !***********************************************************************
      ! Allocate the global albedo and ecosytem map storage arrays based 
      ! upon these
      !  start and stop values.  Note that the arrays will be indexed using 
      ! start/stop
      !  variables, instead of 1-NumberOfZPoints, for ease of selecting 
      ! the nearest
      !  neighbor to the inputted lat/lon grid.
      !************************************************************************
      ! WDR re-allocate here
      if( allocated(ecosystemmap) ) deallocate( ecosystemmap )
      allocate( ecosystemmap(startx:stopx, starty:stopy ) )
                    
      !*************************************************************************
      ! Set up the HDF read variables, making note that HDF arrays start at 0
      !  instead of F90's 1:
      !************************************************************************
      ! Define the starting point
      hdfstart(:) = 0
      hdfstart(1) = startx - 1
      hdfstart(2) = starty - 1
        
      ! Define the stride = 1 so that it doesn't skip any data.
      hdfstride(:) = 1
        
      ! Define the Number of Points to read:
      hdfedge(:) = 1
      hdfedge(1) = numberofxpoints
      hdfedge(2) = numberofypoints
          
      !************************************************************************
      ! Read in the portion of the ecosystem map from the corresponding file:
      !***********************************************************************
      ! Define the SDS Name:
      sdsname = "IGBP_Land_Cover_Type"
      
      ! Open the HDF file for reading:
      ecomapfid = c_sfstart(trim(ecosystemfn), dfacc_read)
      if (ecomapfid == fail) then
        status = failure
        call local_message_handler('Problem detected ecosysyem file sfstart', &
          status,'getAlbedoEco')
        return
      end if
        
      ! Obtain the data set ID's:
      ecomapsdsid   = c_sfselect(ecomapfid, c_sfn2index(ecomapfid, trim(sdsname)))
      if (ecomapsdsid == fail) then
        status = failure
        call local_message_handler('Problem detected ecosystem file sfselect', &
          status,'getAlbedoEco')
        return
      end if
      
      ! Read in the data:
      hdfstatus = c_sfrdata(ecomapsdsid, hdfstart, hdfstride, hdfedge, &
        ecosystemmap)
      !Error checking:
      if (hdfstatus == fail) then
        status = failure
        call local_message_handler('Problem detected ecosystem file sfrdata', &
           status,'getAlbedoEco')
        return
      end if
          
      ! End Access to the datasets and files:
      hdfstatus = c_sfendacc( ecomapsdsid )
      if (hdfstatus == fail) then
        status = failure
        call local_message_handler('Problem detected ecosystem file sfendacc', &
          status,'getAlbedoEco')
        return
      end if
      hdfstatus = c_sfend( ecomapfid   )   
      if (hdfstatus == fail) then
        status = failure
        call local_message_handler('Problem detected ecossystem file sfend', &
          status,'getAlbedoEco')
        return
      end if
      !************************************************************************
      ! Read in the portion of the emissivity map from the corresponding file:
      !***********************************************************************
            
      ! Compute the WesternMostLongitude and NothernMostLatitude of the maps:
      emiss_west = emissivity_res / 2.0 - 180.0
      emiss_north = 90.0 - emissivity_res / 2.0
        
      ! Compute the start/stop x/y indices based on min/max lat/lon values:
      !use anint so that it rounds to nearest index.
      startx = anint( (minlong-emiss_west) / (emissivity_res) ) + 1
      stopx  = anint( (maxlong-emiss_west) / (emissivity_res) ) + 1
      starty = anint( (emiss_north-maxlat ) / (emissivity_res) ) + 1
      stopy  = anint( (emiss_north-minlat ) / (emissivity_res) ) + 1
        
      ! Due to roundoff, geocoordinates of -90, 90, -180, 180 can result in 
      ! indices of one beyond the map bounds, so correct that here:
      if (startx == 0) startx = 1
      if (stopx  == num_emiss_cols+1 ) stopx = num_emiss_cols
      if (starty == 0) starty = 1
      if (stopy  == num_emiss_rows+1 ) stopy = num_emiss_rows
        
      ! Compute the number of points:
      numberofxpoints = stopx - startx + 1
      numberofypoints = stopy - starty + 1
  
      ! Error Checking:
      if ( (startx < 1) .or. (stopx > num_emiss_cols) .or. &
        (starty < 1) .or. (stopy > num_emiss_rows) .or. &
        (startx > stopx) .or. (starty > stopy)      ) then
    
            if (stopx > nummapcols) stopx = nummapcols
            if (stopy > nummaprows) stopy = nummaprows
            if (startx > stopx) startx = stopx-1
            if (starty > stopy) starty = stopy-1
            if (startx < 1) startx = 1
            if (starty < 1) starty = 1
    
        ! status = error
        ! call local_message_handler('Problem with loop index, contact SDST', &
        ! status,'getAlbedoEco')
        ! return  
      end if
        
      ! Define the starting point
      hdfstart(:) = 0
      hdfstart(1) = startx - 1
      hdfstart(2) = starty - 1
        
      ! Define the stride = 1 so that it doesn't skip any data.
      hdfstride(:) = 1
        
      ! Define the Number of Points to read:
      hdfedge(:) = 1
      hdfedge(1) = numberofxpoints
      hdfedge(2) = numberofypoints  
        
      ! WDR re-allocate here instead
      if( allocated( emissivity_map ) ) then
        deallocate( emissivity_map )
      end if
      allocate( emissivity_map(startx:stopx, starty:stopy, 2) )
      allocate( emissivity_map_int(startx:stopx, starty:stopy ) )
    
      ecomapfid = c_sfstart(trim(emissivity_name), dfacc_read)
    
      do i=1, 2
    
        if (i==1) sdsname = "emiss7"
        if (i==2) sdsname = "emiss14"
        
        ecomapsdsid = c_sfselect(ecomapfid, c_sfn2index(ecomapfid, &
          trim(sdsname)))
            
        hdfstatus = c_sfrdata(ecomapsdsid, hdfstart, hdfstride, hdfedge, &
          emissivity_map_int)
        emissivity_map(:,:,i) = emissivity_map_int * scalefactor
            
        hdfstatus = c_sfendacc(ecomapsdsid)
        
      end do
    
      hdfstatus = c_sfend(ecomapfid)
    
      deallocate(emissivity_map_int)
      ! print*,'to snowalb'
  
    end if
  
    !*****************************************************************************
    ! Read in the NISE snow type from the corresponding file:
    !*****************************************************************************
    allocate( snowicetype(1:numberofcols, 1:numberofrows   ) ) 
    snowicetype = 0  ! WDR-UIV
  
    if (errorlevel == fail) then
      status = failure
      call local_message_handler('Problem reported from GetNISEType', &
        status,'getAlbedoEco')
      return
    end if
                                       
    ! print*,'tp latlon loop',NumberOfCols,NumberOfRows
  
    ! WDR we insert the l2 snow-ice array here
    snowicetype = c2_alt_snowice
    !*****************************************************************************
    ! Loop over each pixel in the lat/lon grid, compute the global coordinates and
    !  set the pixel's albedo and ecosystem values to the corresponding global
    !  albedo and ecosystem map values.  This is essentially a nearest neighbor 
    !  calculation.
    !*****************************************************************************
    do i = 1, numberofcols
      do j = 1, numberofrows
        
        sfc_info(i,j) = 0
        ! Check to see if the geolocation is valid, also check we have a cloud. 
        if ( (lat(i,j) >= -90.0000) .and. &
          (lat(i,j) <=  90.0000) .and. &
          (long(i,j) >= -180.000) .and. &
          (long(i,j) <=  180.000) ) then 
  
        ! put the emissivity values where they belong 
        emiss_x = anint( (long(i,j)-emiss_west) / (emissivity_res) ) + 1
       emiss_y = anint( (emiss_north-lat(i,j)) / (emissivity_res) ) + 1
  
        if (emiss_x == 0) emiss_x = 1
        if (emiss_x  == num_emiss_cols+1 ) emiss_x = num_emiss_cols
        if (emiss_y == 0) emiss_y = 1
        if (emiss_y  == num_emiss_rows+1 ) emiss_y = num_emiss_rows
  
  
       emissivity(i,j, 1) = emissivity_map(emiss_x, emiss_y, 1)
        emissivity(i,j, 2) = emissivity_map(emiss_x, emiss_y, 2)
            
        if (emissivity(i,j,1) < 0.) emissivity(i,j,1) = emissivity_fill
        if (emissivity(i,j,2) < 0.) emissivity(i,j,2) = emissivity_fill
            
        ! end emissivity additions
  
        ! print*,'rows=',j,i,Map_Resolution,Long(i,j),westernMostLongitude, &
        !  northernMostLatitude,Lat (i,j)
           
        ! Valid geolocation, so proceed.
            
        ! Compute the ecosystem map x,y indices from this pixel's lat/long values.
        xindex = anint( (long(i,j)-westernmostlongitude) / (map_resolution) ) + 1
        yindex = anint( (northernmostlatitude-lat(i,j)) / (map_resolution) ) + 1
 
        ! Due to roundoff, geocoordinates of -90, 90, -180, 180 can result in 
        ! indices of one beyond the map bounds, so correct that here:
        if (xindex == 0) xindex = 1
          if (xindex  == nummapcols+1 ) xindex = nummapcols
          if (yindex == 0) yindex = 1
          if (yindex  == nummaprows+1 ) yindex = nummaprows
          ! print*,'xIndex, NumMapCols,yIndex, NumMapRows=',xIndex, &
          ! NumMapCols,yIndex, NumMapRows
          ! Checks to make sure indices are within bounds of the map:
          if (      (xindex < 1) .or. (xindex > nummapcols) &
            .or. (yindex < 1) .or. (yindex > nummaprows) ) then
            status = error
            call local_message_handler( &
       'Problem detected, bad array indexes (main loop) ',status,'getAlbedoEco')
              return
          end if
  
          ! Implement Rich Frey (UWisc) interpolation scheme in which he 
          ! sets the map
          ! to 200 if it should be set to snow/ice.  For sea ice, he also 
          ! uses the icec data:
#ifndef GEOS5_SFC
          if (lat(i,j) >= 40.0 .or. lat(i,j) <= -40.0) then
            if (ecosystemmap(xindex,yindex) == 0) then
              ! This is water, so check if the interpolation flagged as 200, 
              ! or the icec data
              ! says this is ice, or the NISE has misidentified it as Perm Snow:
              if (     (icec(i,j) > 0.5 .and. icec(i,j) <= 1.0) &
                .or. (snowicetype(i,j) == 101) &
                .or. (snowicetype(i,j) == 200) ) then
                ! This is water, so this will be sea ice.  Set snowIceType to 95%
                ! so that the sea ice flag is triggered
                snowicetype(i,j) = 95
              end if
            else
              if (     (snowicetype(i,j) >= 1 .and. snowicetype(i,j) <= 100) &
                .or. (snowicetype(i,j) == 200) ) then
                ! This is snow covered land, so set to dry snow:
                snowicetype(i,j) = 103
              end if
            end if
          end if
#else
          if (ecosystemmap(xindex,yindex) == 0) then 
            if ( snowicetype(i,j) >= 30 .and. snowicetype(i,j) <= 100) then 
              snowicetype(i,j) = snowicetype(i,j) 
             else
                snowicetype(i,j) = 0
             endif
#ifdef SSFR_INST
            if ( icec(i,j) > 0. .and. icec(i,j) <= 1.0 ) &
                snowicetype(i,j) = nint(icec(i,j)*100)
#else
            if ( icec(i,j) > 0.5 .and. icec(i,j) <= 1.0 ) &
                snowicetype(i,j) = nint(icec(i,j)*100)
#endif
          endif
  
#endif
  
          ! Store the ecosystem value:
          ecosystem(i,j) = ecosystemmap(xindex,yindex)
          ! print*,xIndex,yIndex,EcosystemMap(xIndex,yIndex)
          ! WDR we use the sfc albedo from l2gen instead
          albedo_real4(:) = c2_sfc_albedo( i, j, :)
                  
          ! Set the albedo to either land values, water, or snow/ice, depending
          ! on the scenario:
          if (      (snowicetype(i,j) >= 1)   &
            .and. (snowicetype(i,j) <= 100) &
            .and. (ecosystem(i,j) == 0)     ) then
  
            ecosystem(i,j) = 15
  
            !This is sea ice, so set the either wet or dry scenario:
                
             sfc_info(i,j) = 1      
                
            ! Determine arctic melting season and set albedo values:
            if (lat(i,j) >= 0.0) then
            ! Arctic melting season: June 1st (day 152) to September 1st (day 244):
            SELECT CASE (julianday)
              CASE ( (nhmeltbeg-transitiondays):nhmeltbeg-1)
                ! transition period between winter and melt, so linear 
                ! fit between them
                ! to get albedo value:
                slope(:) = (processmeltseaicealbedos(:) - &
                  processdryseaicealbedos(:)) / real(transitiondays) 
                intercept(:) = processdryseaicealbedos(:) - slope(:)  &
                  * real(nhmeltbeg-transitiondays) 
                transitionalbedo(:) = slope(:) * real(julianday) + intercept(:)
                albedo_real4(:) = (transitionalbedo(:) *  &
                  real(snowicetype(i,j))/100.000)         &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
              CASE ( nhmeltbeg:nhmeltend )
                ! Melt Sea Ice Values:
                albedo_real4(:) = (processmeltseaicealbedos(:) * &
                  real(snowicetype(i,j))/100.000) &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
              CASE ( nhmeltend+1:(nhmeltend+transitiondays) )
                ! transition period between melt and winter, so linear 
                ! fit between them
                ! to get albedo value:
                slope(:) = (processdryseaicealbedos(:)- &
                  processmeltseaicealbedos(:)) / real(transitiondays) 
                intercept(:) = processmeltseaicealbedos(:) - slope(:)  &
                  * real(nhmeltend) 
                transitionalbedo(:) = slope(:) * real(julianday) + intercept(:)    
                albedo_real4(:) = (transitionalbedo(:) * &
                  real(snowicetype(i,j))/100.000)         &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
              CASE DEFAULT
                ! Dry Sea Ice values::
                albedo_real4(:) = (processdryseaicealbedos(:) *  &
                  real(snowicetype(i,j))/100.000) &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
            END SELECT
    
          else           
            ! Determine southern hemisphere melting season:
            ! Antarctic melting season: December 1st (day 334) to 
            ! March 1st (day 61):
            SELECT CASE (julianday)
              CASE ( (shmeltbeg-transitiondays):shmeltbeg-1)
                ! transition period between winter and melt, so linear 
                ! fit between them
                ! to get albedo value:
                slope(:) = (processmeltseaicealbedos(:)- &
                  processdryseaicealbedos(:))   &
                  / real(transitiondays) 
                intercept(:) = processdryseaicealbedos(:) - slope(:)  &
                  * real(shmeltbeg-transitiondays) 
                transitionalbedo(:) = slope(:) * real(julianday) + intercept(:)    
                albedo_real4(:) = (transitionalbedo(:) * &
                  real(snowicetype(i,j))/100.000)         &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
                      
              CASE ( shmeltbeg:366 )
                ! Melt Sea Ice Values:
                albedo_real4(:) = (processmeltseaicealbedos(:) * &
                  real(snowicetype(i,j))/100.000) &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
              CASE ( 1:shmeltend )
                ! Melt Sea Ice Values:
                albedo_real4(:) = (processmeltseaicealbedos(:) * &
                  real(snowicetype(i,j))/100.000) &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
              CASE (shmeltend+1:(shmeltend+transitiondays) )
                ! transition period between melt and winter, so linear fit 
                !  between them
                ! to get albedo value:
                slope(:) = (processdryseaicealbedos(:)- &
                  processmeltseaicealbedos(:)) / real(transitiondays) 
                intercept(:) = processmeltseaicealbedos(:) - slope(:) * &
                  real(shmeltend) 
              
                transitionalbedo(:) = slope(:) * real(julianday) + intercept(:) 
                albedo_real4(:) = (transitionalbedo(:) * &
                  real(snowicetype(i,j))/100.000)         &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
    
              CASE DEFAULT
                ! Dry Sea Ice values::
                albedo_real4(:) = (processdryseaicealbedos(:) * &
                  real(snowicetype(i,j))/100.000) &
                  + (wateralbedo * real(100-snowicetype(i,j))/100.000)
            END SELECT
                
          end if
 
          albedo(i,j,:) = nint(albedo_real4(:)*1000.)
 
      else if ( ecosystem(i,j) == 15 .or. snowicetype(i,j) == 101 ) then
        ! This is permanent snow/ice, so set to either map or stat value:
        ! WDR new logic with c2_sfc_albedo
        if ((ecosystem(i,j) == 15) .and. &
          (      c2_sfc_albedo(i,j,1) > 0 &
          .and. c2_sfc_albedo(i,j,1) < 1.) ) then
 
          ! This pixel is perm snow class and has a valid map value,
          ! so set to the map's value:
          albedo(i,j,:) = nint( (1.0 - emissivity(i,j,1))*1000.)
 
          sfc_info(i,j) = 3
 
        else
          ! This pixel isn't flag as perm ice eco, so set it to the
          ! statistical value for perm ice:
          ecosystem(i,j) = 15
          albedo(i,j,:) = nint(processsnowalbstats(:, ecosystem(i,j), 1)*1000.)
          albedo(i,j,6) = nint( (1.0 - emissivity(i,j,1))*1000.)
        
          sfc_info(i,j) = 3
        end if
#ifndef GEOS5_SFC
        else if ( snowicetype(i,j) == 103 ) then
          !This is dry snow, so set the pixel to dry snow stats:
          ! if (Ecosystem(i,j) == 0) Ecosystem(i,j) = 15
          albedo(i,j,:) = nint( processsnowalbstats(:, ecosystem(i,j), 1) * 1000.)
          albedo(i,j,6) = nint( (1.0 - emissivity(i,j,1))*1000.)
          ecosystem(i,j)= 15
  
          sfc_info(i,j) = 3
 
#else
        ! use the snow fraction as true fraction 
        else if (snowicetype(i,j) >= 30 .and. snowicetype(i,j) <= 100 .and. &
          ecosystem(i,j) /= 0) then 
          albedo(i,j,:) = nint( (c2_sfc_albedo(i,j,:)* &
            (1.-snowicetype(i,j)*0.01) + &
            snowicetype(i,j)*0.01* &
            processsnowalbstats(:, ecosystem(i,j), 1)) * 1000.)
          albedo(i,j,6) = nint( (1.0 - emissivity(i,j,1))*1000.)
          ecosystem(i,j)= 15
  
          sfc_info(i,j) = 3
 
#endif
 
        else if ( snowicetype(i,j) == 104 ) then
          ! This is wet snow, so set the pixel to wet snow stats:
          ! if (Ecosystem(i,j) == 0) Ecosystem(i,j) = 15
          albedo(i,j,:) = nint(processsnowalbstats(:, ecosystem(i,j), 2) * 1000.)
          albedo(i,j,6) = nint((1.0 - emissivity(i,j,1))*1000.)
          ecosystem(i,j)= 15
  
          sfc_info(i,j) = 3
  
  
        else if (ecosystem(i,j) == 0) then
          ! Water, so set all wavelengths to the water albedo value, 5%:
          albedo(i,j,:) = nint(wateralbedo*1000.)
  
          sfc_info(i,j) = 0
        else
  
          ! This is a snow-free, ocean free, and not perm snow/ice pixel,
          ! so just store the value from the albedo map:
          albedo(i,j,:) = c2_sfc_albedo(i,j,:) * 1000.
          ! use the 1-emissivity for the 3.7um albedo map
          albedo(i,j,6) = nint((1.0 - emissivity(i,j,1))*1000.)
 
          sfc_info(i,j) = 2
        end if
 
      else
    
        ! Bad geolocation, so set albedos and emissivity to fill value
        albedo(i,j,:) = -999
        emissivity(i,j,:) = -999.
   
      end if
    end do
  end do
  !*****************************************************************************
  ! Deallocate arrays:
  !*****************************************************************************
  ! WDR we'll be keeping some of these arrays
  !WDRdeallocate( AlbedoMap        )
  !WDRdeallocate( EcosystemMap     ) 
  !WDRdeallocate( emissivity_map)
  deallocate( snowicetype      )
  !WDRdeallocate(albedo_real4)
 
  if (status > 0 )then
    status = error
    call local_message_handler('Undetermined problem detected in getAlbedoEco',&
      status,'getAlbedoEco')
  endif 
    
  end subroutine getalbedoeco
 
end module modis_albedo