/WPS/geogrid/src/source_data_module.F

Large files files are truncated, but you can click here to view the full file

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Module: source_data_module
!
! Description: 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module source_data_module

   use hash_module
   use list_module
   use module_debug
   use misc_definitions_module
 
   ! Parameters
   integer, parameter :: RETURN_LANDMASK = 0, &
                         RETURN_DOMCAT_LM = 1, &
                         RETURN_DFDX_LM = 2, &
                         RETURN_DFDY_LM = 3, &
                         RETURN_FIELDNAME = 4, &
                         RETURN_DOMCAT = 5, &
                         RETURN_DFDX = 6, &
                         RETURN_DFDY = 7
   integer, parameter :: MAX_LANDMASK_CATEGORIES = 100
 
   ! Module variables
   integer :: num_entries
   integer :: next_field = 1
   integer :: output_field_state = RETURN_LANDMASK 
   integer, pointer, dimension(:) :: source_proj, source_wordsize, source_endian, source_fieldtype, &
                  source_dest_fieldtype, source_priority, source_tile_x, source_tile_y, &
                  source_tile_z, source_tile_z_start, source_tile_z_end, source_tile_bdr, &
                  source_category_min, source_category_max, source_smooth_option, &
                  source_smooth_passes, source_output_stagger, source_row_order
   integer :: source_iswater, source_islake, source_isice, source_isurban, source_isoilwater
   real, pointer, dimension(:) :: source_dx, source_dy, source_known_x, source_known_y, &
                  source_known_lat, source_known_lon, source_masked, source_truelat1, source_truelat2, &
                  source_stdlon, source_scale_factor, source_missing_value, source_fill_missing
   character (len=128), pointer, dimension(:) :: source_fieldname, source_path, source_interp_string, &
                  source_dominant_category, source_dominant_only, source_dfdx, source_dfdy, &
                  source_z_dim_name, source_units, source_descr, source_geotiff_file
   character (len=128) :: source_mminlu
   logical, pointer, dimension(:) :: is_proj, is_wordsize, is_endian, is_fieldtype, &
                  is_dest_fieldtype, is_priority, is_tile_x, is_tile_y, is_tile_z, &
                  is_tile_z_start, is_tile_z_end, is_tile_bdr, is_category_min, &
                  is_category_max, is_masked, &
                  is_dx, is_dy, is_known_x, is_known_y, &
                  is_known_lat, is_known_lon, is_truelat1, is_truelat2, is_stdlon, &
                  is_scale_factor, is_fieldname, is_path, is_dominant_category, &
                  is_dominant_only, is_dfdx, is_dfdy, is_z_dim_name, &
                  is_smooth_option, is_smooth_passes, is_signed, is_missing_value, &
                  is_fill_missing, is_halt_missing, is_output_stagger, is_row_order, &
                  is_units, is_descr, is_subgrid, is_geotiff

   type (list), pointer, dimension(:) :: source_res_path, source_interp_option, source_landmask_land, &
                                         source_landmask_water
   type (hashtable) :: bad_files   ! Track which files produce errors when we try to open them
   type (hashtable) :: duplicate_fnames  ! Remember which output fields we have returned

 
   contains
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_datalist
   !
   ! Purpose:
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine get_datalist()
 
      use gridinfo_module
      use stringutil
  
      implicit none
  
      ! Parameters
      integer, parameter :: BUFSIZE = 256
  
      ! Local variables
      integer :: nparams, idx, eos, ispace, comma, i, j, funit
      logical :: have_specification, is_used
      character (len=128) :: res_string, path_string, interp_string, landmask_string
      character (len=BUFSIZE) :: buffer
  
      nparams = 0
      num_entries = 0
  
      do funit=10,100
         inquire(unit=funit, opened=is_used)
         if (.not. is_used) exit
      end do
      open(funit,file=trim(opt_geogrid_tbl_path)//'GEOGRID.TBL',form='formatted',status='old',err=1000)
  
      !
      ! We will first go through the file to determine how many field
      !   specifications there are
      !
    10 read(funit,'(a)',end=20,err=1001) buffer
      call despace(buffer)
  
      ! Is this line a comment?
      if (buffer(1:1) == '#') then
  
      ! Are we beginning a new field specification?
      else if (index(buffer,'=====') /= 0) then
         if (nparams > 0) num_entries = num_entries + 1  
         nparams = 0
  
      else
         eos = index(buffer,'#')
         if (eos /= 0) buffer(eos:BUFSIZE) = ' ' 
   
         ! Does this line contain at least one parameter specification?
         if (index(buffer,'=') /= 0) then
            nparams = nparams + 1
         end if
      end if
      go to 10
  
    20 rewind(funit)
  
      ! In case the last entry ended without a ======== line
      if (nparams > 0) num_entries = num_entries + 1  
  
      call mprintf(.true.,STDOUT,'Parsed %i entries in GEOGRID.TBL', i1=num_entries)
  
      !
      ! Now that we know how many fields the user has specified, allocate
      !   the properly sized arrays
      !
      allocate(source_wordsize(num_entries))
      allocate(source_endian(num_entries))
      allocate(source_fieldtype(num_entries))
      allocate(source_dest_fieldtype(num_entries))
      allocate(source_proj(num_entries))
      allocate(source_priority(num_entries))
      allocate(source_dx(num_entries))
      allocate(source_dy(num_entries))
      allocate(source_known_x(num_entries))
      allocate(source_known_y(num_entries))
      allocate(source_known_lat(num_entries))
      allocate(source_known_lon(num_entries))
      allocate(source_truelat1(num_entries))
      allocate(source_truelat2(num_entries))
      allocate(source_stdlon(num_entries))
      allocate(source_fieldname(num_entries))
      allocate(source_path(num_entries))
      allocate(source_interp_string(num_entries))
      allocate(source_tile_x(num_entries))
      allocate(source_tile_y(num_entries))
      allocate(source_tile_z(num_entries))
      allocate(source_tile_z_start(num_entries))
      allocate(source_tile_z_end(num_entries))
      allocate(source_category_min(num_entries))
      allocate(source_category_max(num_entries))
      allocate(source_tile_bdr(num_entries))
      allocate(source_masked(num_entries))
      allocate(source_output_stagger(num_entries))
      allocate(source_row_order(num_entries))
      allocate(source_dominant_category(num_entries))
      allocate(source_dominant_only(num_entries))
      allocate(source_dfdx(num_entries))
      allocate(source_dfdy(num_entries))
      allocate(source_scale_factor(num_entries))
      allocate(source_z_dim_name(num_entries))
      allocate(source_units(num_entries))
      allocate(source_descr(num_entries))
      allocate(source_smooth_option(num_entries))
      allocate(source_smooth_passes(num_entries))
      allocate(source_missing_value(num_entries))
      allocate(source_fill_missing(num_entries))
      allocate(source_res_path(num_entries))
      allocate(source_interp_option(num_entries))
      allocate(source_landmask_land(num_entries))
      allocate(source_landmask_water(num_entries))
      allocate(source_geotiff_file(num_entries))
      do i=1,num_entries
         call list_init(source_res_path(i))
         call list_init(source_interp_option(i))
         call list_init(source_landmask_land(i))
         call list_init(source_landmask_water(i))
      end do
  
      allocate(is_wordsize(num_entries))
      allocate(is_endian(num_entries))
      allocate(is_fieldtype(num_entries))
      allocate(is_dest_fieldtype(num_entries))
      allocate(is_proj(num_entries))
      allocate(is_priority(num_entries))
      allocate(is_dx(num_entries))
      allocate(is_dy(num_entries))
      allocate(is_known_x(num_entries))
      allocate(is_known_y(num_entries))
      allocate(is_known_lat(num_entries))
      allocate(is_known_lon(num_entries))
      allocate(is_truelat1(num_entries))
      allocate(is_truelat2(num_entries))
      allocate(is_stdlon(num_entries))
      allocate(is_fieldname(num_entries))
      allocate(is_path(num_entries))
      allocate(is_tile_x(num_entries))
      allocate(is_tile_y(num_entries))
      allocate(is_tile_z(num_entries))
      allocate(is_tile_z_start(num_entries))
      allocate(is_tile_z_end(num_entries))
      allocate(is_category_min(num_entries))
      allocate(is_category_max(num_entries))
      allocate(is_tile_bdr(num_entries))
      allocate(is_masked(num_entries))
      allocate(is_halt_missing(num_entries))
      allocate(is_output_stagger(num_entries))
      allocate(is_row_order(num_entries))
      allocate(is_dominant_category(num_entries))
      allocate(is_dominant_only(num_entries))
      allocate(is_dfdx(num_entries))
      allocate(is_dfdy(num_entries))
      allocate(is_scale_factor(num_entries))
      allocate(is_z_dim_name(num_entries))
      allocate(is_units(num_entries))
      allocate(is_descr(num_entries))
      allocate(is_smooth_option(num_entries))
      allocate(is_smooth_passes(num_entries))
      allocate(is_signed(num_entries))
      allocate(is_missing_value(num_entries))
      allocate(is_fill_missing(num_entries))
      allocate(is_subgrid(num_entries))
      allocate(is_geotiff(num_entries))

      source_wordsize=0
      source_endian=0
      source_fieldtype=0
      source_dest_fieldtype=0
      source_proj=0
      source_priority=0
      source_dx=0
      source_dy=0
      source_known_x=0
      source_known_y=0
      source_known_lat=0
      source_known_lon=0
      source_truelat1=0
      source_truelat2=0
      source_stdlon=0
      source_fieldname=' '
      source_path=' '
      source_interp_string=' '
      source_tile_x=0
      source_tile_y=0
      source_tile_z=0
      source_tile_z_start=0
      source_tile_z_end=0
      source_category_min=0
      source_category_max=0
      source_tile_bdr=0
      source_masked=0
      source_output_stagger=0
      source_row_order=0
      source_dominant_category=' '
      source_dominant_only=' '
      source_dfdx=' '
      source_dfdy=' '
      source_scale_factor=0
      source_z_dim_name=' '
      source_units=' '
      source_descr=' '
      source_smooth_option=0
      source_smooth_passes=0
      source_missing_value=0
      source_fill_missing=0

      is_wordsize=.false.
      is_endian=.false.
      is_fieldtype=.false.
      is_dest_fieldtype=.false.
      is_proj=.false.
      is_priority=.false.
      is_dx=.false.
      is_dy=.false.
      is_known_x=.false.
      is_known_y=.false.
      is_known_lat=.false.
      is_known_lon=.false.
      is_truelat1=.false.
      is_truelat2=.false.
      is_stdlon=.false.
      is_fieldname=.false.
      is_path=.false.
      is_tile_x=.false.
      is_tile_y=.false.
      is_tile_z=.false.
      is_tile_z_start=.false.
      is_tile_z_end=.false.
      is_category_min=.false.
      is_category_max=.false.
      is_tile_bdr=.false.
      is_masked=.false.
      is_halt_missing=.false.
      is_output_stagger=.false.
      is_row_order=.false.
      is_dominant_category=.false.
      is_dominant_only=.false.
      is_dfdx=.false.
      is_dfdy=.false.
      is_scale_factor=.false.
      is_z_dim_name=.false.
      is_units=.false.
      is_descr=.false.
      is_smooth_option=.false.
      is_smooth_passes=.false.
      is_signed=.false.
      is_missing_value=.false.
      is_fill_missing=.false.
      is_subgrid=.false.

      write(source_mminlu,'(a4)') 'USGS'
      source_iswater = 16
      source_islake = -1
      source_isice = 24
      source_isurban = 1
      source_isoilwater = 14
  
      ! 
      ! Actually read and save the specifications
      !
      nparams = 0
      i = 1
    30 buffer = ' '
      read(funit,'(a)',end=40,err=1001) buffer
      call despace(buffer)
  
      ! Is this line a comment?
      if (buffer(1:1) == '#') then
         ! Do nothing.
  
      ! Are we beginning a new field specification?
      else if (index(buffer,'=====') /= 0) then   !{
         if (nparams > 0) i = i + 1  
         nparams = 0
         if (i <= num_entries) then
!BUG: Are these initializations needed now that we've added initializations for
!     all options after their initial allocation above?
            is_path(i) = .false.
            is_masked(i) = .false.
            is_halt_missing(i) = .false.
            is_output_stagger(i) = .false.
            is_dominant_category(i) = .false.
            is_dominant_only(i) = .false.
            is_dfdx(i) = .false.
            is_dfdy(i) = .false.
            is_dest_fieldtype(i) = .false.
            is_priority(i) = .false.
            is_z_dim_name(i) = .false.
            is_smooth_option(i) = .false.
            is_smooth_passes(i) = .false.
            is_fill_missing(i) = .false.
            is_subgrid(i) = .false.
         end if
  
      else
         ! Check whether the current line is a comment
         if (buffer(1:1) /= '#') then
            have_specification = .true. 
         else
            have_specification = .false.
         end if

         ! If only part of the line is a comment, just turn the comment into spaces 
         eos = index(buffer,'#')
         if (eos /= 0) buffer(eos:BUFSIZE) = ' ' 
   
         do while (have_specification)   !{
   
            ! If this line has no semicolon, it may contain a single specification,
            !   so we set have_specification = .false. to prevent the line from being 
            !   processed again and "pretend" that the last character was a semicolon
            eos = index(buffer,';')
            if (eos == 0) then
               have_specification = .false.
               eos = BUFSIZE
            end if
    
            idx = index(buffer(1:eos-1),'=')
    
            if (idx /= 0) then   !{
               nparams = nparams + 1
     
               if (index('name',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_fieldname(i) = .true.
                  source_fieldname(i) = ' '
                  source_fieldname(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
     
               else if (index('priority',trim(buffer(1:idx-1))) /= 0) then
                  is_priority(i) = .true.
                  read(buffer(idx+1:eos-1),'(i10)') source_priority(i)
     
               else if (index('dest_type',trim(buffer(1:idx-1))) /= 0) then
                  if (index('continuous',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_dest_fieldtype(i) = .true.
                     source_dest_fieldtype(i) = CONTINUOUS
                  else if (index('categorical',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_dest_fieldtype(i) = .true.
                     source_dest_fieldtype(i) = CATEGORICAL
                  end if
     
               else if (index('interp_option',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  interp_string = ' '
                  interp_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
                  ispace = index(interp_string,':')
                  if (ispace /= 0) then
                     write(res_string,'(a)') interp_string(1:ispace-1)
                  else
                     res_string = 'default'
                  end if
                  write(interp_string,'(a)') trim(interp_string(ispace+1:128))
                  if (list_search(source_interp_option(i), ckey=res_string, cvalue=interp_string)) then
                     call mprintf(.true., WARN, &
                                  'In GEOGRID.TBL entry %i, multiple interpolation methods are '// &
                                  'given for resolution %s. %s will be used.', &
                                  i1=i, s1=trim(res_string), s2=trim(interp_string))
                  else
                     call list_insert(source_interp_option(i), ckey=res_string, cvalue=interp_string)
                  end if
     
               else if (index('smooth_option',trim(buffer(1:idx-1))) /= 0) then
                  if ((index('1-2-1',trim(buffer(idx+1:eos-1))) /= 0) .and. &
                      (len_trim(buffer(idx+1:eos-1)) == len('1-2-1'))) then
                     is_smooth_option(i) = .true.
                     source_smooth_option(i) = ONETWOONE
                  else if ((index('smth-desmth',trim(buffer(idx+1:eos-1))) /= 0) .and. & 
                      (len_trim(buffer(idx+1:eos-1)) == len('smth-desmth'))) then
                     is_smooth_option(i) = .true.
                     source_smooth_option(i) = SMTHDESMTH
                  else if ((index('smth-desmth_special',trim(buffer(idx+1:eos-1))) /= 0) .and. &
                      (len_trim(buffer(idx+1:eos-1)) == len('smth-desmth_special'))) then
                     is_smooth_option(i) = .true.
                     source_smooth_option(i) = SMTHDESMTH_SPECIAL
                  end if
     
               else if (index('smooth_passes',trim(buffer(1:idx-1))) /= 0) then
                  is_smooth_passes(i) = .true.
                  read(buffer(idx+1:eos-1),'(i10)') source_smooth_passes(i)
       
               else if (index('rel_path',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  path_string = ' '
                  path_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
                  if (path_string(ispace-idx-1:ispace-idx-1) /= '/') &
                     path_string(ispace-idx:ispace-idx) = '/'
                  if (path_string(1:1) == '/') then
                     path_string(1:127) = path_string(2:128)
                     path_string(128:128) = ' '
                  end if
                  ispace = index(path_string,':')
                  if (ispace /= 0) then
                     write(res_string,'(a)') path_string(1:ispace-1)
                  else
                     res_string = 'default'
                  end if
                  write(path_string,'(a)') trim(geog_data_path)//trim(path_string(ispace+1:128))
                  if (list_search(source_res_path(i), ckey=res_string, cvalue=path_string)) then
                     call mprintf(.true., WARN, &
                                  'In GEOGRID.TBL entry %i, multiple paths are given for '// &
                                  'resolution %s. %s will be used.', &
                                  i1=i, s1=trim(res_string), s2=trim(path_string))
                  else
                     call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
                  end if
     
               else if (index('abs_path',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  path_string = ' '
                  path_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
                  if (path_string(ispace-idx-1:ispace-idx-1) /= '/') &
                     path_string(ispace-idx:ispace-idx) = '/'
                  ispace = index(path_string,':')
                  if (ispace /= 0) then
                     write(res_string,'(a)') path_string(1:ispace-1)
                  else
                     res_string = 'default'
                  end if
                  write(path_string,'(a)') trim(path_string(ispace+1:128))
                  if (list_search(source_res_path(i), ckey=res_string, cvalue=path_string)) then
                     call mprintf(.true., WARN, &
                                  'In GEOGRID.TBL entry %i, multiple paths are given for '// &
                                  'resolution %s. %s will be used.', &
                                  i1=i, s1=trim(res_string), s2=trim(path_string))
                  else
                     call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
                  end if
    
               else if (index('output_stagger',trim(buffer(1:idx-1))) /= 0) then
                  if (index('M',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_output_stagger(i) = .true.
                     source_output_stagger(i) = M
                  else if (index('U',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_output_stagger(i) = .true.
                     source_output_stagger(i) = U
                  else if (index('V',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_output_stagger(i) = .true.
                     source_output_stagger(i) = V
                  else if (index('HH',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_output_stagger(i) = .true.
                     source_output_stagger(i) = HH
                  else if (index('VV',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_output_stagger(i) = .true.
                     source_output_stagger(i) = VV
                  end if
     
               else if ((index('landmask_water',trim(buffer(1:idx-1))) /= 0) .and. &
                        (len_trim(buffer(1:idx-1)) == 14)) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  landmask_string = ' '
                  landmask_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
                  ispace = index(landmask_string,':')
                  if (ispace /= 0) then
                     write(res_string,'(a)') landmask_string(1:ispace-1)
                  else
                     res_string = 'default'
                  end if
                  write(landmask_string,'(a)') trim(landmask_string(ispace+1:128))
                  if (list_search(source_landmask_water(i), ckey=res_string, cvalue=landmask_string)) then
                     call mprintf(.true., WARN, &
                                  'In GEOGRID.TBL entry %i, multiple landmask category specifications are given for '// &
                                  'resolution %s. %s will be used.', &
                                  i1=i, s1=trim(res_string), s2=trim(landmask_string))
                  else
                     call list_insert(source_landmask_water(i), ckey=res_string, cvalue=landmask_string)
                  end if

               else if ((index('landmask_land',trim(buffer(1:idx-1))) /= 0) .and. &
                        (len_trim(buffer(1:idx-1)) == 13)) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  landmask_string = ' '
                  landmask_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
                  ispace = index(landmask_string,':')
                  if (ispace /= 0) then
                     write(res_string,'(a)') landmask_string(1:ispace-1)
                  else
                     res_string = 'default'
                  end if
                  write(landmask_string,'(a)') trim(landmask_string(ispace+1:128))
                  if (list_search(source_landmask_land(i), ckey=res_string, cvalue=landmask_string)) then
                     call mprintf(.true., WARN, &
                                  'In GEOGRID.TBL entry %i, multiple landmask category specifications are given for '// &
                                  'resolution %s. %s will be used.', &
                                  i1=i, s1=trim(res_string), s2=trim(landmask_string))
                  else
                     call list_insert(source_landmask_land(i), ckey=res_string, cvalue=landmask_string)
                  end if
     
               else if ((index('masked',trim(buffer(1:idx-1))) /= 0) .and. &
                        (len_trim(buffer(1:idx-1)) == 6)) then
                  if (index('water',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_masked(i) = .true.
                     source_masked(i) = 0.
                  else if (index('land',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_masked(i) = .true.
                     source_masked(i) = 1.
                  end if
     
               else if (index('fill_missing',trim(buffer(1:idx-1))) /= 0) then
                  is_fill_missing(i) = .true.
                  read(buffer(idx+1:eos-1),*) source_fill_missing(i)
       
               else if (index('halt_on_missing',trim(buffer(1:idx-1))) /= 0) then
                  if (index('yes',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_halt_missing(i) = .true.
                  else if (index('no',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_halt_missing(i) = .false.
                  end if
     
               else if (index('dominant_category',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_dominant_category(i) = .true.
                  source_dominant_category(i) = ' '
                  source_dominant_category(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
      
               else if (index('dominant_only',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_dominant_only(i) = .true.
                  source_dominant_only(i) = ' '
                  source_dominant_only(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
     
               else if (index('df_dx',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_dfdx(i) = .true.
                  source_dfdx(i) = ' '
                  source_dfdx(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
     
               else if (index('df_dy',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_dfdy(i) = .true.
                  source_dfdy(i) = ' '
                  source_dfdy(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
     
               else if (index('z_dim_name',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  is_z_dim_name(i) = .true.
                  source_z_dim_name(i) = ' '
                  source_z_dim_name(i)(1:ispace-idx) = buffer(idx+1:ispace-1)

               else if (index('subgrid',trim(buffer(1:idx-1))) /= 0) then
                  ispace = idx+1
                  do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
                     ispace = ispace + 1
                  end do 
                  if (index('yes',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_subgrid(i) = .true.
                  else if (index('no',trim(buffer(idx+1:eos-1))) /= 0) then
                     is_subgrid(i) = .false.
                  end if

               else
                   call mprintf(.true., WARN, 'In GEOGRID.TBL, unrecognized option %s in '// &
                                'entry %i.',i1=idx, s1=buffer(i:idx-1))
       
               end if
    
            end if   !} index(buffer(1:eos-1),'=') /= 0
    
            buffer = buffer(eos+1:BUFSIZE)
         end do   ! while eos /= 0 }
  
      end if   !} index(buffer, '=====') /= 0
      go to 30
  
    40 close(funit)
  
      ! Check the user specifications for gross errors
      if ( .not. check_data_specification() ) then
         call datalist_destroy()
         call mprintf(.true.,ERROR,'Errors were found in either index files or GEOGRID.TBL.')
      end if
  
      call hash_init(bad_files)
  
      return
  
    1000 call mprintf(.true.,ERROR,'Could not open GEOGRID.TBL')
  
    1001 call mprintf(.true.,ERROR,'Encountered error while reading GEOGRID.TBL')
 
   end subroutine get_datalist
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   ! Name: get_source_params
   !
   ! Purpose: For each field, this routine reads in the metadata in the index file 
   !   for the first available resolution of data specified by res_string. Also 
   !   based on res_string, this routine sets the interpolation sequence for the
   !   field. This routine should be called prior to processing a field for each
   !   domain.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   subroutine get_source_params(res_string)
 
      use stringutil
#ifdef _HAS_GEOTIFF
      use geotiff_module, only : open_geotiff,merge_geotiff_header
#endif
 
      implicit none
  
      ! Parameters
      integer, parameter :: BUFSIZE = 256
  
      ! Arguments
      character (len=128), intent(in) :: res_string
  
      ! Local variables
      integer :: idx, i, is, ie, ispace, eos, iquoted, funit
      character (len=128) :: temp_data, temp_interp
      character (len=256) :: test_fname
      character (len=BUFSIZE) :: buffer
      logical :: have_specification, is_used

      ! For each entry in the GEOGRID.TBL file
      do idx=1,num_entries

         ! Initialize metadata
         is_wordsize(idx)  = .false.
         is_endian(idx)    = .false.
         is_row_order(idx) = .false.
         is_fieldtype(idx) = .false.
         is_proj(idx)      = .false.
         is_dx(idx)        = .false.
         is_dy(idx)        = .false.
         is_known_x(idx)   = .false.
         is_known_y(idx)   = .false.
         is_known_lat(idx) = .false.
         is_known_lon(idx) = .false.
         is_truelat1(idx)  = .false.
         is_truelat2(idx)  = .false.
         is_stdlon(idx)    = .false.
         is_tile_x(idx)    = .false.
         is_tile_y(idx)    = .false.
         is_tile_z(idx)    = .false.
         is_tile_z_start(idx) = .false.
         is_tile_z_end(idx)   = .false.
         is_category_min(idx) = .false.
         is_category_max(idx) = .false.
         is_tile_bdr(idx)     = .false.
         is_fieldname(idx)    = .false.
         is_scale_factor(idx) = .false.
         is_units(idx)        = .false.
         is_descr(idx)        = .false.
         is_signed(idx)       = .false.
         is_missing_value(idx) = .false.
         is_geotiff(idx)       = .false.
   
         
         ! Set the interpolator sequence for the field to be the first value in res_string that matches
         !   the resolution keyword for an interp_sequence specification
         is = 1
         ie = index(res_string(is:128),'+') - 1
         if (ie <= 0) ie = 128
         temp_interp = res_string(is:ie)
         do while (.not. list_search(source_interp_option(idx), ckey=temp_interp, cvalue=source_interp_string(idx)) &
                   .and. is <= 128)
            call mprintf(.true., INFORM, 'For %s, couldn''t find interpolator sequence for '// &
                         'resolution %s.', &
                         s1=trim(source_fieldname(idx)), s2=trim(temp_interp))
            is = ie+2
            ie = is + index(res_string(is:128),'+') - 2
            if (ie - is <= 0) ie = 128
            temp_interp = res_string(is:ie)
         end do

         if (is > 128) then
            temp_interp = 'default'
            if (list_search(source_interp_option(idx), ckey=temp_interp, cvalue=source_interp_string(idx))) then
               call mprintf(.true., INFORM, 'Using default interpolator sequence for %s.', &
                            s1=trim(source_fieldname(idx)))
            else
               call mprintf(.true., ERROR, 'Could not find interpolator sequence for requested resolution for %s.'// &
                            ' The sources specified in namelist.wps is not listed in GEOGRID.TBL.', &
                            s1=trim(source_fieldname(idx)))
            end if
         else
            call mprintf(.true., INFORM, 'Using %s interpolator sequence for %s.', &
                         s1=temp_interp, s2=trim(source_fieldname(idx)))
         end if
   
         ! Set the data source for the field to be the first value in res_string that matches
         !   the resolution keyword for an abs_path or rel_path specification
         is = 1
         ie = index(res_string(is:128),'+') - 1
         if (ie <= 0) ie = 128
         temp_data = res_string(is:ie)
         do while (.not. list_search(source_res_path(idx), ckey=temp_data, cvalue=source_path(idx)) &
                   .and. is <= 128)
            call mprintf(.true., INFORM, 'For %s, couldn''t find %s data source.', &
                         s1=trim(source_fieldname(idx)), s2=trim(temp_data))
            is = ie+2
            ie = is + index(res_string(is:128),'+') - 2
            if (ie - is <= 0) ie = 128
            temp_data = res_string(is:ie)
         end do

         if (is > 128) then
            temp_data = 'default'
            if (list_search(source_res_path(idx), ckey=temp_data, cvalue=source_path(idx))) then
               call mprintf(.true., INFORM, 'Using default data source for %s.', &
                            s1=trim(source_fieldname(idx)))
            else
               call mprintf(.true., ERROR, 'Could not find data resolution for requested resolution for %s. '// &
                            'The source specified in namelist.wps is not listed in GEOGRID.TBL.', &
                            s1=trim(source_fieldname(idx)))
            end if
         else
            call mprintf(.true., INFORM, 'Using %s data source for %s.', &
                         s1=temp_data, s2=trim(source_fieldname(idx)))
         end if

         call mprintf(trim(temp_data) /= trim(temp_interp),WARN,'For %s, using %s data source with %s interpolation sequence.', &
                      s1=source_fieldname(idx), s2=temp_data, s3=temp_interp)

         write(test_fname, '(a)') trim(source_path(idx))//'index'
     
         !
         ! Open the index file for the data source for this field, and read in metadata specs
         !
         do funit=10,100
            inquire(unit=funit, opened=is_used)
            if (.not. is_used) exit
         end do
         open(funit,file=trim(test_fname),form='formatted',status='old',err=1000)
   
      30 buffer = ' '
         read(funit,'(a)',end=40, err=1001) buffer
         call despace(buffer)
         
         ! Is this line a comment?
         if (buffer(1:1) == '#') then
            ! Do nothing.
     
         else
            have_specification = .true. 
      
            ! If only part of the line is a comment, just turn the comment into spaces 
            eos = index(buffer,'#')
            if (eos /= 0) buffer(eos:BUFSIZE) = ' ' 
      
            do while (have_specification)   !{
      
               ! If this line has no semicolon, it may contain a single specification,
               !   so we set have_specification = .false. to prevent the line from being 
               !   processed again and pretend that the last character was a semicolon
               eos = index(buffer,';')
               if (eos == 0) then
                  have_specification = .false.
                  eos = BUFSIZE
               end if
       
               i = index(buffer(1:eos-1),'=')
       
               if (i /= 0) then   !{
       
                  if (index('projection',trim(buffer(1:i-1))) /= 0) then
                     if (index('lambert',trim(buffer(i+1:eos-1))) /= 0) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_LC
                     else if (index('polar_wgs84',trim(buffer(i+1:eos-1))) /= 0 .and. &
                              len_trim('polar_wgs84') == len_trim(buffer(i+1:eos-1))) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_PS_WGS84
                     else if (index('albers_nad83',trim(buffer(i+1:eos-1))) /= 0 .and. &
                              len_trim('albers_nad83') == len_trim(buffer(i+1:eos-1))) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_ALBERS_NAD83
                     else if (index('polar',trim(buffer(i+1:eos-1))) /= 0 .and. &
                              len_trim('polar') == len_trim(buffer(i+1:eos-1))) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_PS
                     else if (index('mercator',trim(buffer(i+1:eos-1))) /= 0) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_MERC
                     else if (index('regular_ll',trim(buffer(i+1:eos-1))) /= 0) then
                        is_proj(idx) = .true.
                        source_proj(idx) = PROJ_LATLON
                     end if
        
                  else if (index('type',trim(buffer(1:i-1))) /= 0) then
                     if (index('continuous',trim(buffer(i+1:eos-1))) /= 0) then
                        is_fieldtype(idx) = .true.
                        source_fieldtype(idx) = CONTINUOUS
                     else if (index('categorical',trim(buffer(i+1:eos-1))) /= 0) then
                        is_fieldtype(idx) = .true.
                        source_fieldtype(idx) = CATEGORICAL
                     end if
        
                  else if (index('signed',trim(buffer(1:i-1))) /= 0) then
                     if (index('yes',trim(buffer(i+1:eos-1))) /= 0) then
                        is_signed(idx) = .true.
                     else if (index('no',trim(buffer(i+1:eos-1))) /= 0) then
                        is_signed(idx) = .false.
                     end if
        
                  else if (index('units',trim(buffer(1:i-1))) /= 0) then
                     ispace = i+1
                     iquoted = 0
                     do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
                        if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
                        ispace = ispace + 1
                     end do 
                     is_units(idx) = .true.
                     source_units(idx) = ' '
                     if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
                     if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
                     source_units(idx)(1:ispace-i) = buffer(i+1:ispace-1)
        
                  else if (index('description',trim(buffer(1:i-1))) /= 0) then
                     ispace = i+1
                     iquoted = 0
                     do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
                        if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
                        ispace = ispace + 1
                     end do 
                     is_descr(idx) = .true.
                     source_descr(idx) = ' '
                     if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
                     if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
                     source_descr(idx)(1:ispace-i) = buffer(i+1:ispace-1)
        
                  else if (index('mminlu',trim(buffer(1:i-1))) /= 0) then
                     ispace = i+1
                     iquoted = 0
                     do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
                        if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
                        ispace = ispace + 1
                     end do 
                     if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
                     if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
                     source_mminlu(1:ispace-i) = buffer(i+1:ispace-1)
        
                  else if (index('iswater',trim(buffer(1:i-1))) /= 0) then
                     read(buffer(i+1:eos-1),*) source_iswater
          
                  else if (index('islake',trim(buffer(1:i-1))) /= 0) then
                     read(buffer(i+1:eos-1),*) source_islake
          
                  else if (index('isice',trim(buffer(1:i-1))) /= 0) then
                     read(buffer(i+1:eos-1),*) source_isice
          
                  else if (index('isurban',trim(buffer(1:i-1))) /= 0) then
                     read(buffer(i+1:eos-1),*) source_isurban
          
                  else if (index('isoilwater',trim(buffer(1:i-1))) /= 0) then
                     read(buffer(i+1:eos-1),*) source_isoilwater
          
                  else if (index('dx',trim(buffer(1:i-1))) /= 0) then
                     is_dx(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_dx(idx)
          
                  else if (index('dy',trim(buffer(1:i-1))) /= 0) then
                     is_dy(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_dy(idx)
          
                  else if (index('known_x',trim(buffer(1:i-1))) /= 0) then
                     is_known_x(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_known_x(idx)
          
                  else if (index('known_y',trim(buffer(1:i-1))) /= 0) then
                     is_known_y(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_known_y(idx)
          
                  else if (index('known_lat',trim(buffer(1:i-1))) /= 0) then
                     is_known_lat(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_known_lat(idx)
          
                  else if (index('known_lon',trim(buffer(1:i-1))) /= 0) then
                     is_known_lon(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_known_lon(idx)
          
                  else if (index('stdlon',trim(buffer(1:i-1))) /= 0) then
                     is_stdlon(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_stdlon(idx)
          
                  else if (index('truelat1',trim(buffer(1:i-1))) /= 0) then
                     is_truelat1(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_truelat1(idx)
          
                  else if (index('truelat2',trim(buffer(1:i-1))) /= 0) then
                     is_truelat2(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_truelat2(idx)
          
                  else if (index('wordsize',trim(buffer(1:i-1))) /= 0) then
                     is_wordsize(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_wordsize(idx)
        
                  else if (index('endian',trim(buffer(1:i-1))) /= 0) then
                     if (index('big',trim(buffer(i+1:eos-1))) /= 0) then
                        is_endian(idx) = .true.
                        source_endian(idx) = BIG_ENDIAN
                     else if (index('little',trim(buffer(i+1:eos-1))) /= 0) then
                        is_endian(idx) = .true.
                        source_endian(idx) = LITTLE_ENDIAN
                     else
                        call mprintf(.true.,WARN,'Invalid value for keyword ''endian'' '// &
                                     'specified in index file. BIG_ENDIAN will be used.')
                     end if
          
                  else if (index('row_order',trim(buffer(1:i-1))) /= 0) then
                     if (index('bottom_top',trim(buffer(i+1:eos-1))) /= 0) then
                        is_row_order(idx) = .true.
                        source_row_order(idx) = BOTTOM_TOP
                     else if (index('top_bottom',trim(buffer(i+1:eos-1))) /= 0) then
                        is_row_order(idx) = .true.
                        source_row_order(idx) = TOP_BOTTOM
                     end if
          
                  else if (index('tile_x',trim(buffer(1:i-1))) /= 0) then
                     is_tile_x(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_x(idx)
        
                  else if (index('tile_y',trim(buffer(1:i-1))) /= 0) then
                     is_tile_y(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_y(idx)
        
                  else if (index('tile_z',trim(buffer(1:i-1))) /= 0) then
                     is_tile_z(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_z(idx)
        
                  else if (index('tile_z_start',trim(buffer(1:i-1))) /= 0) then
                     is_tile_z_start(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_z_start(idx)
        
                  else if (index('tile_z_end',trim(buffer(1:i-1))) /= 0) then
                     is_tile_z_end(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_z_end(idx)
        
                  else if (index('category_min',trim(buffer(1:i-1))) /= 0) then
                     is_category_min(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_category_min(idx)
        
                  else if (index('category_max',trim(buffer(1:i-1))) /= 0) then
                     is_category_max(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_category_max(idx)
        
                  else if (index('tile_bdr',trim(buffer(1:i-1))) /= 0) then
                     is_tile_bdr(idx) = .true.
                     read(buffer(i+1:eos-1),'(i10)') source_tile_bdr(idx)
        
                  else if (index('missing_value',trim(buffer(1:i-1))) /= 0) then
                     is_missing_value(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_missing_value(idx)
        
                  else if (index('scale_factor',trim(buffer(1:i-1))) /= 0) then
                     is_scale_factor(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_scale_factor(idx)

                  else if (index('geotiff',trim(buffer(1:i-1))) /= 0) then
                     is_geotiff(idx) = .true.
                     read(buffer(i+1:eos-1),*) source_geotiff_file(idx)
          
                  else
                     call mprintf(.true., WARN, 'In %s, unrecognized option %s in entry %i.', &
                                  s1=trim(test_fname), s2=buffer(1:i-1), i1=i)
          
                  end if
       
               end if   !} index(buffer(1:eos-1),'=') /= 0
       
               buffer = buffer(eos+1:BUFSIZE)
            end do   ! while eos /= 0 }
     
         end if   !} index(buffer, '=====') /= 0

         go to 30

     40 continue
    
#ifdef _HAS_GEOTIFF
         if(is_geotiff(idx)) then
           if(source_geotiff_file(idx)(1:1) .ne. '/')then
             source_geotiff_file(idx)=TRIM(source_path(idx))//TRIM(source_geotiff_file(idx))
           endif
           call open_geotiff(source_geotiff_file(idx))
           call merge_geotiff_header(source_geotiff_file(idx),is_proj(idx),source_proj(idx),&
                        is_fieldtype(idx),source_fieldtype(idx),           &
                        is_units(idx),source_units(idx),                   &
                        is_descr(idx),source_descr(idx),is_dx(idx),        &
                        source_dx(idx),is_dy(idx),source_dy(idx),          &
                        is_known_x(idx),source_known_x(idx),               &
                        is_known_y(idx),source_known_y(idx),               &
                        is_known_lat(idx),source_known_lat(idx),           &
                        is_known_lon(idx),source_known_lon(idx),           &
                        is_stdlon(idx),source_stdlon(idx),                 &
                        is_truelat1(idx),source_truelat1(idx),             &
                        is_truelat2(idx),source_truelat2(idx),             &
                        is_row_order(idx),source_row_order(idx),           &
                        is_tile_x(idx),source_tile_x(idx),                 &
                        is_tile_y(idx),source_tile_y(idx),                 &
                        is_tile_z(idx),source_tile_z(idx),                 &
                        is_tile_z_start(idx),source_tile_z_start(idx),     &
                        is_tile_z_end(idx),  source_tile_z_end(idx),       &
                        is_category_min(idx),source_category_min(idx),     &
                        is_category_max(idx),source_category_max(idx),     &
                        is_tile_bdr(idx),source_tile_bdr(idx),             &
                        is_missing_value(idx),source_missing_value(idx))
           is_wordsize(idx)=.true.
           source_wordsize(idx)=4
           is_scale_factor(idx)=.true.
           source_scale_factor(idx)=1.
           call mprintf(.true.,INFORM,'For geotiff file, %s, using the following parameters:', &
                        s1=TRIM(source_geotiff_file(idx)))
           if(is_descr(idx)) &
             call mprintf(.true.,INFORM,'description=%s',s1=TRIM(sou…