initInterpMindist.F90 Source File


Source Code

module initInterpMindist

  ! ------------------------------------------------------------------------
  ! This module contains most of the "interesting" logic of initInterp, in terms of
  ! finding the input column (or landunit, patch, etc.) to use as a template for each
  ! output column (etc.).
  !
  ! This is in a separate module to facilitate unit testing, since the full initInterp
  ! involves some awkward dependencies.
  ! ------------------------------------------------------------------------

#include "shr_assert.h"
  use shr_kind_mod   , only: r8 => shr_kind_r8
  use shr_log_mod    , only : errMsg => shr_log_errMsg
  use clm_varctl     , only: iulog
  use abortutils     , only: endrun
  use spmdMod        , only: masterproc
  use clm_varcon     , only: spval, re
  use glcBehaviorMod , only: glc_behavior_type

  implicit none
  private
  save

  ! Public methods

  public :: set_mindist
  public :: set_single_match

  ! Public types

  type, public :: subgrid_special_indices_type
     integer :: ipft_not_vegetated
     integer :: icol_vegetated_or_bare_soil
     integer :: ilun_vegetated_or_bare_soil
     integer :: ilun_crop
     integer :: ilun_landice_multiple_elevation_classes
   contains
     procedure :: is_vegetated_landunit  ! returns true if the given landunit type is natural veg or crop
  end type subgrid_special_indices_type

  type, public :: subgrid_type
     character(len=16) :: name               ! pft, column, landunit, gridcell
     integer , pointer :: ptype(:) => null() ! used for patch type 
     integer , pointer :: ctype(:) => null() ! used for patch or col type
     integer , pointer :: ltype(:) => null() ! used for pft, col or lun type
     real(r8), pointer :: topoglc(:) => null()
     real(r8), pointer :: lat(:)
     real(r8), pointer :: lon(:)
     real(r8), pointer :: coslat(:)
   contains
     procedure :: print_point  ! print info about one point
  end type subgrid_type

  ! Private methods

  private :: set_glcmec_must_be_same_type
  private :: set_icemec_adjustable_type
  private :: do_fill_missing_with_natveg
  private :: is_sametype
  private :: is_baresoil

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

contains

  !-----------------------------------------------------------------------
  subroutine print_point(this, index, unit)
    !
    ! !DESCRIPTION:
    ! Print info about one point in a subgrid_type object
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(subgrid_type), intent(in) :: this
    integer            , intent(in) :: index
    integer            , intent(in) :: unit ! unit to which we should write the info
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'print_point'
    !-----------------------------------------------------------------------

    write(unit,*) 'subgrid level, index = ',&
         this%name, index
    write(unit,*) 'lat, lon = ', this%lat(index), ', ', this%lon(index)
    if (associated(this%ltype)) then
       write(unit,*) 'ltype: ', this%ltype(index)
    end if
    if (associated(this%ctype)) then
       write(unit,*) 'ctype: ', this%ctype(index)
    end if
    if (associated(this%ptype)) then
       write(unit,*) 'ptype: ', this%ptype(index)
    end if
    
  end subroutine print_point


  !=======================================================================

  subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgrido, &
       subgrid_special_indices, glc_behavior, glc_elevclasses_same, &
       fill_missing_with_natveg, mindist_index)

    ! --------------------------------------------------------------------
    ! arguments
    integer            , intent(in)  :: begi, endi 
    integer            , intent(in)  :: bego, endo 
    logical            , intent(in)  :: activei(begi:endi) 
    logical            , intent(in)  :: activeo(bego:endo) 
    type(subgrid_type) , intent(in)  :: subgridi
    type(subgrid_type) , intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    type(glc_behavior_type), intent(in) :: glc_behavior

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    logical            , intent(in)  :: glc_elevclasses_same

    ! If false: if an output type cannot be found in the input, code aborts
    ! If true: if an output type cannot be found in the input, fill with closest natural
    ! veg column (using bare soil for patch-level variables)
    !
    ! NOTE: always treated as true for natural veg and crop landunits/columns/patches in
    ! the output - e.g., if we can't find the right column type to fill crop, we always
    ! use the closest natural veg column, regardless of the value of this flag.
    logical            , intent(in)  :: fill_missing_with_natveg

    integer            , intent(out) :: mindist_index(bego:endo)
    !
    ! local variables
    real(r8) :: dx,dy
    real(r8) :: distmin,dist,hgtdiffmin,hgtdiff    
    integer  :: nsizei, nsizeo
    integer  :: ni,no,nmin,ier,n,noloc
    logical  :: topoglc_present
    logical  :: closest

    ! Whether two glcmec columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glcmec points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical  :: glcmec_must_be_same_type_o(bego:endo)
    ! --------------------------------------------------------------------

    if (associated(subgridi%topoglc) .and. associated(subgrido%topoglc)) then
       topoglc_present = .true.
    else
       topoglc_present = .false.
    end if

    mindist_index(bego:endo) = 0
    distmin = spval

    call set_glcmec_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, &
         glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, &
         glcmec_must_be_same_type_o=glcmec_must_be_same_type_o)

!$OMP PARALLEL DO PRIVATE (ni,no,n,nmin,distmin,dx,dy,dist,closest,hgtdiffmin,hgtdiff)
    do no = bego,endo

       ! Only interpolate onto active points. Otherwise, the mere act of running
       ! init_interp (e.g., of a file onto itself) would lead to changes in a bunch of
       ! inactive points - e.g., going from their cold start initial conditions to some
       ! spunup initial conditions (from the closest active point of that type). This
       ! could potentially lead to different behavior in a transient run, if those points
       ! later became active; that's undesirable.
       if (activeo(no)) then 
          nmin    = 0
          distmin = spval
          hgtdiffmin = spval
          do ni = begi,endi
             if (activei(ni)) then
                if (is_sametype(ni = ni, no = no, &
                     subgridi = subgridi, subgrido = subgrido, &
                     subgrid_special_indices = subgrid_special_indices, &
                     glcmec_must_be_same_type = glcmec_must_be_same_type_o(no), &
                     veg_patch_just_considers_ptype = .true.)) then
                   dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                   dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                        0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                   dist = dx*dx + dy*dy
                   if (topoglc_present) then
                      hgtdiff = abs(subgridi%topoglc(ni) - subgrido%topoglc(no))
                   end if
                   closest = .false.
                   if ( dist < distmin ) then
                      closest = .true.
                      distmin = dist
                      nmin = ni
                      if (topoglc_present) then
                         hgtdiffmin = hgtdiff
                      end if
                   end if
                   if (.not. closest) then
                      ! In *some* cases: For glc_mec points, we first find the closest
                      ! point in lat-lon space, without consideration for column or patch
                      ! type (above). Then, within that closest point, we find the closest
                      ! column in topographic space; this second piece is done here. Note
                      ! that this ordering means that we could choose a column with a very
                      ! different topographic height from the target column, if it is
                      ! closer in lat-lon space than any similar-height columns.
                      if (topoglc_present) then
                         hgtdiff = abs(subgridi%topoglc(ni) - subgrido%topoglc(no))
                         if ((dist == distmin) .and. (hgtdiff < hgtdiffmin)) then
                            closest = .true.
                            hgtdiffmin = hgtdiff
                            distmin = dist
                            nmin = ni
                         end if
                      end if
                   end if
                end if
             end if
          end do
          
          ! If output type is not contained in input dataset, then use closest bare soil,
          ! if this point is one for which we fill missing with natveg.
          if ( distmin == spval .and. &
               do_fill_missing_with_natveg( &
               fill_missing_with_natveg, no, subgrido, subgrid_special_indices)) then
             do ni = begi, endi
                if (activei(ni)) then
                   if ( is_baresoil(ni, subgridi, subgrid_special_indices)) then
                      dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re
                      dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * &
                           0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni))
                      dist = dx*dx + dy*dy
                      if ( dist < distmin )then
                         distmin = dist
                         nmin = ni
                      end if
                   end if
                end if
             end do
          end if

          ! Error conditions
          if ( distmin == spval )then
             write(iulog,*) 'ERROR initInterp set_mindist: &
                  &Cannot find any input points matching output point:'
             call subgrido%print_point(no, iulog)
             write(iulog,*) ' '
             write(iulog,*) 'Consider rerunning with the following in user_nl_clm:'
             write(iulog,*) 'init_interp_fill_missing_with_natveg = .true.'
             write(iulog,*) 'However, note that this will fill all missing types in the output'
             write(iulog,*) 'with the closest natural veg column in the input'
             write(iulog,*) '(using bare soil for patch-level variables).'
             write(iulog,*) 'So, you should consider whether that is what you want.'
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

          mindist_index(no) = nmin

       end if ! end if activeo block
    end do
!$OMP END PARALLEL DO
    
  end subroutine set_mindist

  !-----------------------------------------------------------------------
  subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, &
       subgrid_special_indices, glc_behavior, glc_elevclasses_same, &
       mindist_index)
    !
    ! !DESCRIPTION:
    ! This is an alternative to set_mindist that can work when there is a single matching
    ! input point for each output point (or possibly 0 matching inputs for an inactive
    ! output point). In this method, we only look for points that share the same lat/lon
    ! as the output point.
    !
    ! This method is intended to be used for the case where we want to copy areas from the
    ! input file to the output file (rather than maintaining areas from the surface
    ! dataset); some of the logic here is set up the way it is in order to support this
    ! copying of areas.
    !
    ! !ARGUMENTS:
    integer            , intent(in)  :: begi, endi
    integer            , intent(in)  :: bego, endo
    logical            , intent(in)  :: activeo(bego:endo)
    type(subgrid_type) , intent(in)  :: subgridi
    type(subgrid_type) , intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    type(glc_behavior_type), intent(in) :: glc_behavior

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    !
    ! Must be true for this method! (Otherwise we may not be able to find a single unique
    ! input point for each output point.)
    logical            , intent(in)  :: glc_elevclasses_same

    integer            , intent(out) :: mindist_index(bego:endo)
    !
    ! !LOCAL VARIABLES:
    integer  :: ni, no, ni_match
    logical  :: found
    real(r8) :: dx, dy
    logical  :: ni_sametype

    ! Whether two glcmec columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glcmec points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical  :: glcmec_must_be_same_type_o(bego:endo)

    ! Tolerance in lat/lon for considering a point to be at the same location
    real(r8) :: same_point_tol = 1.e-14_r8

    character(len=*), parameter :: subname = 'set_single_match'
    !-----------------------------------------------------------------------

    if (.not. glc_elevclasses_same) then
       write(iulog,*) subname// &
            ' ERROR: For this init_interp method, the number and bounds of'
       write(iulog,*) 'glacier elevation classes must be the same in input and output'
       call endrun(msg=subname//' ERROR: glc_elevclasses_same must be true for this method')
    end if

    call set_glcmec_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, &
         glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, &
         glcmec_must_be_same_type_o=glcmec_must_be_same_type_o)

!$OMP PARALLEL DO PRIVATE (ni,no,ni_match,found,dx,dy,ni_sametype)
    do no = bego, endo
       ni_match = 0
       found = .false.
       do ni = begi, endi
          dx = abs(subgrido%lon(no)-subgridi%lon(ni))
          dy = abs(subgrido%lat(no)-subgridi%lat(ni))
          if (dx < same_point_tol .and. dy < same_point_tol) then
             ni_sametype = is_sametype(ni = ni, no = no, &
                  subgridi = subgridi, subgrido = subgrido, &
                  subgrid_special_indices = subgrid_special_indices, &
                  glcmec_must_be_same_type = glcmec_must_be_same_type_o(no), &
                  veg_patch_just_considers_ptype = .false.)
             if (ni_sametype) then
                if (found) then
                   write(iulog,*) subname// &
                        ' ERROR: found multiple input points matching output point'
                   call subgrido%print_point(no, iulog)
                   write(iulog,*) ' '
                   write(iulog,*) 'For this init_interp_method: for a given output point,'
                   write(iulog,*) 'we expect to find exactly one input point at the same'
                   write(iulog,*) 'location with the same type.'
                   write(iulog,*) '(This most likely indicates a problem in CTSM, such as'
                   write(iulog,*) 'having multiple columns with the same column type on a'
                   write(iulog,*) 'given gridcell.)'
                   call endrun(msg=subname// &
                        ' ERROR: found multiple input points matching output point')
                else
                   found = .true.
                   ni_match = ni
                end if
             end if
          end if
       end do

       if (found) then
          mindist_index(no) = ni_match
       else
          mindist_index(no) = 0

          if (activeo(no)) then
             write(iulog,*) subname// &
                  ' ERROR: cannot find any input points matching output point'
             call subgrido%print_point(no, iulog)
             write(iulog,*) ' '
             write(iulog,*) 'For this init_interp_method: for a given active output point,'
             write(iulog,*) 'we expect to find exactly one input point at the same'
             write(iulog,*) 'location with the same type.'
             write(iulog,*) 'Note that this requires the input and output grids to be identical.'
             write(iulog,*) '(If you need to interpolate to a different resolution, then'
             write(iulog,*) 'you need to use a different init_interp_method.)'
             call endrun(msg=subname//' ERROR: cannot find any input points matching output point')
          end if
       end if
    end do
!$OMP END PARALLEL DO

  end subroutine set_single_match

  !-----------------------------------------------------------------------
  subroutine set_glcmec_must_be_same_type(bego, endo, dimname, &
       glc_elevclasses_same, glc_behavior, &
       glcmec_must_be_same_type_o)
    !
    ! !DESCRIPTION:
    ! Sets the glcmec_must_be_same_type_o array for each output icemec point
    !
    ! This array will be set to true for icemec output columns/patches for which the
    ! column/patch type must match the input column/patch type to be considered the same
    ! type. This is only valid for icemec points, and is only valid for dimname = 'pft' or
    ! 'column' - for others, the value is undefined.
    !
    ! This assumes that bego and endo match the bounds that are used elsewhere in the
    ! model - i.e., for the decomposition within init_interp to match the model's
    ! decomposition.
    !
    ! !ARGUMENTS:
    integer                 , intent(in)  :: bego    ! beginning index for output points
    integer                 , intent(in)  :: endo    ! ending index for output points
    character(len=*)        , intent(in)  :: dimname ! 'pft', 'column', etc.

    ! True if the number and bounds of glacier elevation classes are the same in the
    ! input and output; false otherwise.
    logical, intent(in) :: glc_elevclasses_same

    type(glc_behavior_type), intent(in)  :: glc_behavior

    logical, intent(out)  :: glcmec_must_be_same_type_o( bego: ) ! see description above

    !
    ! !LOCAL VARIABLES:
    integer :: no

    ! Whether each output icemec point has adjustable column type; only valid for icemec
    ! points, and only valid for subgrid name = pft or column
    logical  :: icemec_adjustable_type_o(bego:endo)

    character(len=*), parameter :: subname = 'set_glcmec_must_be_same_type'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL((ubound(glcmec_must_be_same_type_o) == (/endo/)), errMsg(sourcefile, __LINE__))

    if (.not. glc_elevclasses_same) then
       ! If the number or bounds of the elevation classes differ between input and
       ! output, then we ignore the column and patch types of glcmec points when
       ! looking for the same type - instead using closest topographic height as a
       ! tie-breaker between equidistant columns/patches.
       glcmec_must_be_same_type_o(bego:endo) = .false.
    else
       call set_icemec_adjustable_type(bego=bego, endo=endo, dimname=dimname, &
            glc_behavior=glc_behavior, icemec_adjustable_type_o=icemec_adjustable_type_o)
       do no = bego, endo
          if (icemec_adjustable_type_o(no)) then
             ! If glcmec points in this output cell have adjustable type, then we ignore
             ! the column and patch types of glcmec points when looking for the same
             ! type: we want to find the closest glcmec point in lat-lon space without
             ! regards for column/patch type, because the column/patch types may change
             ! at runtime.
             glcmec_must_be_same_type_o(no) = .false.
          else
             ! Otherwise, we require the column and patch types to be the same between
             ! input and output for this glcmec output point, as is the case for most
             ! other landunits. This is important for a case with interpolation to give
             ! bit-for-bit answers with a case without interpolation (since glcmec
             ! topographic heights can change after initialization, so we can't always
             ! rely on the point with closest topographic height to be the "right" one to
             ! pick as the source for interpolation).
             glcmec_must_be_same_type_o(no) = .true.
          end if
       end do
    end if

  end subroutine set_glcmec_must_be_same_type

  !-----------------------------------------------------------------------
  subroutine set_icemec_adjustable_type(bego, endo, dimname, glc_behavior, &
       icemec_adjustable_type_o)
    !
    ! !DESCRIPTION:
    ! Sets the icemec_adjustable_type_o array for each output icemec point
    !
    ! This array will be set to true for icemec points that have adjustable column type
    ! and false for icemec points that do not. The value will be undefined for non-icemec
    ! points.
    !
    ! This can only be called for the output, not the input!
    !
    ! The output array is only valid for dimname = pft or column; for others, the output
    ! values are undefined.
    !
    ! This assumes that bego and endo match the bounds that are used elsewhere in the
    ! model - i.e., for the decomposition within init_interp to match the model's
    ! decomposition.
    !
    ! !ARGUMENTS:
    integer                 , intent(in)  :: bego    ! beginning index for output points
    integer                 , intent(in)  :: endo    ! ending index for output points
    character(len=*)        , intent(in)  :: dimname ! 'pft', 'column', etc.
    type(glc_behavior_type) , intent(in)  :: glc_behavior
    logical                 , intent(out) :: icemec_adjustable_type_o( bego: ) ! see documentation above
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'set_icemec_adjustable_type'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL((ubound(icemec_adjustable_type_o) == (/endo/)), errMsg(sourcefile, __LINE__))

    select case (dimname)
    case ('pft')
       call glc_behavior%patches_have_dynamic_type_array(bego, endo, &
            icemec_adjustable_type_o(bego:endo))
    case ('column')
       call glc_behavior%cols_have_dynamic_type_array(bego, endo, &
            icemec_adjustable_type_o(bego:endo))
    case ('landunit', 'gridcell')
       ! Do nothing: icemec_adjustable_type_o will be left undefined
    case default
       call endrun(subname//' ERROR: unexpected dimname: '//trim(dimname)//&
            errMsg(sourcefile, __LINE__))
    end select

  end subroutine set_icemec_adjustable_type

  !-----------------------------------------------------------------------
  function do_fill_missing_with_natveg(fill_missing_with_natveg, &
       no, subgrido, subgrid_special_indices)
    !
    ! !DESCRIPTION:
    ! Returns true if the given output point, if missing, should be filled with the
    ! closest natural veg point.
    !
    ! !ARGUMENTS:
    logical :: do_fill_missing_with_natveg  ! function result

    ! whether we should fill ALL missing points with natveg
    logical, intent(in) :: fill_missing_with_natveg

    integer           , intent(in)  :: no
    type(subgrid_type), intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'do_fill_missing_with_natveg'
    !-----------------------------------------------------------------------

    if (subgrido%name == 'gridcell') then
       ! It makes no sense to try to fill missing with natveg for gridcell-level values
       do_fill_missing_with_natveg = .false.
    else if (fill_missing_with_natveg) then
       ! User has asked for all missing points to be filled with natveg
       do_fill_missing_with_natveg = .true.
    else if (subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then
       ! Even if user hasn't asked for it, we fill missing vegetated points (natural veg
       ! and crop) with the closest natveg point. This is mainly to support the common
       ! use case of interpolating non-crop to crop, but also supports adding a new PFT
       ! type.
       do_fill_missing_with_natveg = .true.
    else
       do_fill_missing_with_natveg = .false.
    end if

  end function do_fill_missing_with_natveg


  !=======================================================================

  logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indices, &
       glcmec_must_be_same_type, veg_patch_just_considers_ptype)

    ! --------------------------------------------------------------------
    ! arguments
    integer           , intent(in)  :: ni 
    integer           , intent(in)  :: no 
    type(subgrid_type), intent(in)  :: subgridi
    type(subgrid_type), intent(in)  :: subgrido
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices

    ! Whether two glcmec columns/patches must be the same column/patch type to be
    ! considered the same type. This is only valid for glcmec points, and is only valid
    ! for subgrid name = 'pft' or 'column'.
    logical, intent(in) :: glcmec_must_be_same_type

    ! For vegetated patches (natural veg or crop): if veg_patch_just_considers_ptype is
    ! true, then we consider two vegetated patches to be the same type if they have the
    ! same ptype, without regard to their column and landunit types (as long as both the
    ! input and the output are on either the natural veg or crop landunit). This is
    ! needed to handle the generic crop properly when interpolating from non-crop to
    ! crop, or vice versa.
    !
    ! If false, then they need to have the same column and landunit types, too (as is the
    ! general case).
    logical, intent(in) :: veg_patch_just_considers_ptype
    ! --------------------------------------------------------------------

    is_sametype = .false.

    if (trim(subgridi%name) == 'pft' .and. trim(subgrido%name) == 'pft') then
       if ( .not. glcmec_must_be_same_type .and. &
            subgridi%ltype(ni) == subgrid_special_indices%ilun_landice_multiple_elevation_classes .and. &
            subgrido%ltype(no) == subgrid_special_indices%ilun_landice_multiple_elevation_classes) then
          is_sametype = .true.
       else if (veg_patch_just_considers_ptype .and. &
            subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then
          ! See comment attached to the declaration of veg_patch_just_considers_ptype for
          ! rationale for this logic.
          !
          ! TODO(wjs, 2015-09-15) If we ever allow the same PFT to appear on multiple
          ! columns within a given grid cell, then this logic will need to be made
          ! somewhat more complex: e.g., preferably take something from the same column
          ! type, but if we can't find anything from the same column type, then ignore
          ! column type.

          if (subgrid_special_indices%is_vegetated_landunit(subgridi%ltype(ni)) .and. &
               subgridi%ptype(ni) == subgrido%ptype(no)) then
             is_sametype = .true.
          end if
       else if (subgridi%ptype(ni) == subgrido%ptype(no) .and. &
                subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'column' .and. trim(subgrido%name) == 'column') then
       if ( .not. glcmec_must_be_same_type .and. &
            subgridi%ltype(ni) == subgrid_special_indices%ilun_landice_multiple_elevation_classes  .and. &
            subgrido%ltype(no) == subgrid_special_indices%ilun_landice_multiple_elevation_classes ) then
          is_sametype = .true.
       else if (subgridi%ctype(ni) == subgrido%ctype(no) .and. &
                subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'landunit' .and. trim(subgrido%name) == 'landunit') then
       if (subgridi%ltype(ni) == subgrido%ltype(no)) then
          is_sametype = .true.
       end if
    else if (trim(subgridi%name) == 'gridcell' .and. trim(subgrido%name) == 'gridcell') then
       is_sametype = .true.
    else 
       if (masterproc) then
          write(iulog,*)'ERROR interpinic: is_sametype check on input and output type not supported'
          write(iulog,*)'typei = ',trim(subgridi%name)
          write(iulog,*)'typeo = ',trim(subgrido%name)
       end if
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

  end function is_sametype

  !=======================================================================

  logical function is_baresoil (n, subgrid, subgrid_special_indices)

    ! --------------------------------------------------------------------
    ! arguments
    integer           , intent(in)  :: n 
    type(subgrid_type), intent(in)  :: subgrid
    type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices
    ! --------------------------------------------------------------------

    is_baresoil = .false.

    if (subgrid%name == 'pft') then
       if (subgrid%ptype(n) == subgrid_special_indices%ipft_not_vegetated .and. &
            subgrid%ctype(n) == subgrid_special_indices%icol_vegetated_or_bare_soil .and. &
            subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else if (subgrid%name == 'column') then
       if (subgrid%ctype(n) == subgrid_special_indices%icol_vegetated_or_bare_soil .and. &
            subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else if (subgrid%name == 'landunit') then
       if (subgrid%ltype(n) == subgrid_special_indices%ilun_vegetated_or_bare_soil) then
          is_baresoil = .true.
       end if
    else 
       if (masterproc) then
          write(iulog,*)'ERROR interpinic: is_baresoil subgrid type ',subgrid%name,' not supported'
       end if
       call endrun(msg=errMsg(sourcefile, __LINE__))
    end if

  end function is_baresoil

  !-----------------------------------------------------------------------
  function is_vegetated_landunit(this, ltype)
    !
    ! !DESCRIPTION:
    ! Returns true if the given landunit type is vegetated: either natural veg or crop
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: is_vegetated_landunit  ! function result
    class(subgrid_special_indices_type), intent(in) :: this
    integer, intent(in) :: ltype   ! landunit type of interest
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'is_vegetated_landunit'
    !-----------------------------------------------------------------------

    if (ltype == this%ilun_vegetated_or_bare_soil .or. &
         ltype == this%ilun_crop) then
       is_vegetated_landunit = .true.
    else
       is_vegetated_landunit = .false.
    end if

  end function is_vegetated_landunit


end module initInterpMindist