module subgridAveMod

!-----------------------------------------------------------------------
!BOP
!
! !MODULE: subgridAveMod
!
! !DESCRIPTION:
! Utilities to perfrom subgrid averaging
!
! !USES:
  use shr_kind_mod, only: r8 => shr_kind_r8
  use clmtype
  use clm_varcon, only : spval, isturb,  icol_roof, icol_sunwall, icol_shadewall, &
                         icol_road_perv, icol_road_imperv
  use clm_varctl, only : iulog
  use abortutils, only : endrun

! !PUBLIC TYPES:
  implicit none
  save
!
! !PUBLIC MEMBER FUNCTIONS:
  public :: p2c   ! Perfrom an average from pfts to columns
  public :: p2l   ! Perfrom an average from pfts to landunits
  public :: p2g   ! Perfrom an average from pfts to gridcells
  public :: c2l   ! Perfrom an average from columns to landunits
  public :: c2g   ! Perfrom an average from columns to gridcells
  public :: l2g   ! Perfrom an average from landunits to gridcells

  interface p2c
     module procedure p2c_1d
     module procedure p2c_2d
     module procedure p2c_1d_filter
     module procedure p2c_2d_filter
  end interface
  interface p2l
     module procedure p2l_1d
     module procedure p2l_2d
  end interface
  interface p2g
     module procedure p2g_1d
     module procedure p2g_2d
  end interface
  interface c2l
     module procedure c2l_1d
     module procedure c2l_2d
  end interface
  interface c2g
     module procedure c2g_1d
     module procedure c2g_2d
  end interface
  interface l2g
     module procedure l2g_1d
     module procedure l2g_2d
  end interface
!
! !PRIVATE MEMBER FUNCTIONS:
  private :: build_scale_l2g
  private :: create_scale_l2g_lookup
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!EOP

! WJS (10-14-11): TODO:
! 
! - I believe that scale_p2c, scale_c2l and scale_l2g should be included in the sumwt
! accumulations (e.g., sumwt = sumwt + wtgcell * scale_p2c * scale_c2l * scale_l2g), but
! that requires some more thought to (1) make sure that is correct, and (2) make sure it
! doesn't break the urban scaling. (See also my notes in create_scale_l2g_lookup.)
!   - Once that is done, you could use a scale of 0, avoiding the need for the use of
!   spval and the special checks that requires.
!
! - Currently, there is a lot of repeated code to calculate scale_c2l. This should be
! cleaned up.
!   - At a minimum, should collect the repeated code into a subroutine to eliminate this
!   repitition
!   - The best thing might be to use a lookup array, as is done for scale_l2g
! -----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2c_1d
!
! !INTERFACE:
  subroutine p2c_1d (lbp, ubp, lbc, ubc, parr, carr, p2c_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from pfts to columns.
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_col
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp              ! beginning and ending pft
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column
    real(r8), intent(in)  :: parr(lbp:ubp)         ! pft array
    real(r8), intent(out) :: carr(lbc:ubc)         ! column array
    character(len=*), intent(in) :: p2c_scale_type ! scale type
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: pi,p,c,index           ! indices
    real(r8) :: scale_p2c(lbp:ubp)     ! scale factor for column->landunit mapping
    logical  :: found                  ! temporary for error check
    real(r8) :: sumwt(lbc:ubc)         ! sum of weights
    real(r8), pointer :: wtcol(:)      ! weight of pft relative to column
    integer , pointer :: pcolumn(:)    ! column index of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in column
    integer , pointer :: pfti(:)       ! initial pft index in column
!------------------------------------------------------------------------

    wtcol    => clm3%g%l%c%p%wtcol
    pcolumn  => clm3%g%l%c%p%column
    npfts    => clm3%g%l%c%npfts
    pfti     => clm3%g%l%c%pfti

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2c_1d error: scale type ',p2c_scale_type,' not supported'
       call endrun()
    end if

    carr(lbc:ubc) = spval
    sumwt(lbc:ubc) = 0._r8
    do p = lbp,ubp
       if (wtcol(p) /= 0._r8) then
          if (parr(p) /= spval) then
             c = pcolumn(p)
             if (sumwt(c) == 0._r8) carr(c) = 0._r8
             carr(c) = carr(c) + parr(p) * scale_p2c(p) * wtcol(p)
             sumwt(c) = sumwt(c) + wtcol(p)
          end if
       end if
    end do
    found = .false.
    do c = lbc,ubc
       if (sumwt(c) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = c
       else if (sumwt(c) /= 0._r8) then
          carr(c) = carr(c)/sumwt(c)
       end if
    end do
    if (found) then
       write(iulog,*)'p2c error: sumwt is greater than 1.0 at c= ',index
       call endrun()
    end if

  end subroutine p2c_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2c_2d
!
! !INTERFACE:
  subroutine p2c_2d (lbp, ubp, lbc, ubc, num2d, parr, carr, p2c_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from landunits to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_col
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp              ! beginning and ending pft
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(r8), intent(in)  :: parr(lbp:ubp,num2d)   ! pft array
    real(r8), intent(out) :: carr(lbc:ubc,num2d)   ! column array
    character(len=*), intent(in) :: p2c_scale_type ! scale type
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,pi,p,c,index         ! indices
    real(r8) :: scale_p2c(lbp:ubp)     ! scale factor for column->landunit mapping
    logical  :: found                  ! temporary for error check
    real(r8) :: sumwt(lbc:ubc)         ! sum of weights
    real(r8), pointer :: wtcol(:)      ! weight of pft relative to column
    integer , pointer :: pcolumn(:)    ! column index of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in column
    integer , pointer :: pfti(:)       ! initial pft index in column
!------------------------------------------------------------------------

    wtcol    => clm3%g%l%c%p%wtcol
    pcolumn  => clm3%g%l%c%p%column
    npfts    => clm3%g%l%c%npfts
    pfti     => clm3%g%l%c%pfti

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2c_2d error: scale type ',p2c_scale_type,' not supported'
       call endrun()
    end if

    carr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do p = lbp,ubp
          if (wtcol(p) /= 0._r8) then
             if (parr(p,j) /= spval) then
                c = pcolumn(p)
                if (sumwt(c) == 0._r8) carr(c,j) = 0._r8
                carr(c,j) = carr(c,j) + parr(p,j) * scale_p2c(p) * wtcol(p)
                sumwt(c) = sumwt(c) + wtcol(p)
             end if
          end if
       end do
       found = .false.
       do c = lbc,ubc
          if (sumwt(c) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index = c
          else if (sumwt(c) /= 0._r8) then
             carr(c,j) = carr(c,j)/sumwt(c)
          end if
       end do
       if (found) then
          write(iulog,*)'p2c_2d error: sumwt is greater than 1.0 at c= ',index,' lev= ',j
          call endrun()
       end if
    end do 
  end subroutine p2c_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2c_1d_filter
!
! !INTERFACE:
  subroutine p2c_1d_filter (numfc, filterc, pftarr, colarr)
!
! !DESCRIPTION:
! perform pft to column averaging for single level pft arrays
!
! !USES:
    use clm_varpar, only : max_pft_per_col
    use clm_varcon, only : istice_mec
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: numfc
    integer , intent(in)  :: filterc(numfc)
    real(r8), pointer     :: pftarr(:)
    real(r8), pointer     :: colarr(:)
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer :: fc,c,pi,p,l         ! indices
    integer , pointer :: npfts(:)
    integer , pointer :: pfti(:)
    integer , pointer :: pftf(:)
    integer , pointer :: clandunit(:)
    integer , pointer :: ltype(:)
    real(r8), pointer :: wtcol(:)
    real(r8), pointer :: wtgcell(:)
!-----------------------------------------------------------------------

    npfts     => clm3%g%l%c%npfts
    pfti      => clm3%g%l%c%pfti
    pftf      => clm3%g%l%c%pftf
    clandunit => clm3%g%l%c%landunit
    ltype     => clm3%g%l%itype
    wtcol     => clm3%g%l%c%p%wtcol
    wtgcell   => clm3%g%l%c%p%wtgcell

    do fc = 1,numfc
       c = filterc(fc)
       l = clandunit(c)
       colarr(c) = 0._r8
       do p = pfti(c), pftf(c)
          ! Note: some glacier_mec pfts may have zero weight
          if (wtgcell(p) > 0._r8 .or. ltype(l)==istice_mec) colarr(c) = colarr(c) + pftarr(p) * wtcol(p)
       end do
    end do

  end subroutine p2c_1d_filter

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2c_2d_filter
!
! !INTERFACE:
  subroutine p2c_2d_filter (lev, numfc, filterc, pftarr, colarr)
!
! !DESCRIPTION:
! perform pft to column averaging for multi level pft arrays
!
! !USES:
    use clm_varpar, only : max_pft_per_col

! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lev
    integer , intent(in)  :: numfc
    integer , intent(in)  :: filterc(numfc)
    real(r8), pointer     :: pftarr(:,:)
    real(r8), pointer     :: colarr(:,:)
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer :: fc,c,pi,p,j    ! indices
    integer , pointer :: npfts(:)
    integer , pointer :: pfti(:)
    integer , pointer :: pftf(:)
    real(r8), pointer :: wtcol(:)
!-----------------------------------------------------------------------

    npfts => clm3%g%l%c%npfts
    pfti  => clm3%g%l%c%pfti
    pftf  => clm3%g%l%c%pftf
    wtcol => clm3%g%l%c%p%wtcol

    do j = 1,lev
       do fc = 1,numfc
          c = filterc(fc)
          colarr(c,j) = 0._r8
          do p = pfti(c), pftf(c)
             colarr(c,j) = colarr(c,j) + pftarr(p,j) * wtcol(p)
          end do
       end do
    end do

  end subroutine p2c_2d_filter

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2l_1d
!
! !INTERFACE:
  subroutine p2l_1d (lbp, ubp, lbc, ubc, lbl, ubl, parr, larr, &
       p2c_scale_type, c2l_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from pfts to landunits
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_lu
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp              ! beginning and ending pft indices
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    real(r8), intent(in)  :: parr(lbp:ubp)         ! input column array
    real(r8), intent(out) :: larr(lbl:ubl)         ! output landunit array
    character(len=*), intent(in) :: p2c_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: pi,p,c,l,index         ! indices
    logical  :: found                  ! temporary for error check
    real(r8) :: sumwt(lbl:ubl)         ! sum of weights
    real(r8) :: scale_p2c(lbc:ubc)     ! scale factor for pft->column mapping
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor for column->landunit mapping
    real(r8), pointer :: wtlunit(:)    ! weight of pft relative to landunit
    integer , pointer :: pcolumn(:)    ! column of corresponding pft
    integer , pointer :: plandunit(:)  ! landunit of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in landunit
    integer , pointer :: pfti(:)       ! initial pft index in landunit
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    canyon_hwr => clm3%g%l%canyon_hwr
    ltype      => clm3%g%l%itype
    ctype      => clm3%g%l%c%itype
    clandunit  => clm3%g%l%c%landunit
    wtlunit    => clm3%g%l%c%p%wtlunit
    pcolumn    => clm3%g%l%c%p%column
    plandunit  => clm3%g%l%c%p%landunit
    npfts      => clm3%g%l%npfts
    pfti       => clm3%g%l%pfti

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'p2l_1d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2l_1d error: scale type ',p2c_scale_type,' not supported'
       call endrun()
    end if

    larr(:) = spval
    sumwt(:) = 0._r8
    do p = lbp,ubp
       if (wtlunit(p) /= 0._r8) then
          c = pcolumn(p)
          if (parr(p) /= spval .and. scale_c2l(c) /= spval) then
             l = plandunit(p)
             if (sumwt(l) == 0._r8) larr(l) = 0._r8
             larr(l) = larr(l) + parr(p) * scale_p2c(p) * scale_c2l(c) * wtlunit(p)
             sumwt(l) = sumwt(l) + wtlunit(p)
          end if
       end if
    end do
    found = .false.
    do l = lbl,ubl
       if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = l
       else if (sumwt(l) /= 0._r8) then
          larr(l) = larr(l)/sumwt(l)
       end if
    end do
    if (found) then
       write(iulog,*)'p2l_1d error: sumwt is greater than 1.0 at l= ',index
       call endrun()
    end if

  end subroutine p2l_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2l_2d
!
! !INTERFACE:
  subroutine p2l_2d(lbp, ubp, lbc, ubc, lbl, ubl, num2d, parr, larr, &
       p2c_scale_type, c2l_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from pfts to landunits
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_lu
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp              ! beginning and ending pft indices
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(r8), intent(in)  :: parr(lbp:ubp,num2d)   ! input pft array
    real(r8), intent(out) :: larr(lbl:ubl,num2d)   ! output gridcell array
    character(len=*), intent(in) :: p2c_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,pi,p,c,l,index       ! indices
    logical  :: found                  ! temporary for error check
    real(r8) :: sumwt(lbl:ubl)         ! sum of weights
    real(r8) :: scale_p2c(lbc:ubc)     ! scale factor for pft->column mapping
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor for column->landunit mapping
    real(r8), pointer :: wtlunit(:)    ! weight of pft relative to landunit
    integer , pointer :: pcolumn(:)    ! column of corresponding pft
    integer , pointer :: plandunit(:)  ! landunit of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in landunit
    integer , pointer :: pfti(:)       ! initial pft index in landunit
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    canyon_hwr => clm3%g%l%canyon_hwr
    ltype      => clm3%g%l%itype
    clandunit  => clm3%g%l%c%landunit
    ctype      => clm3%g%l%c%itype
    wtlunit   => clm3%g%l%c%p%wtlunit
    pcolumn   => clm3%g%l%c%p%column
    plandunit => clm3%g%l%c%p%landunit
    npfts     => clm3%g%l%npfts
    pfti      => clm3%g%l%pfti

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'p2l_2d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2l_2d error: scale type ',p2c_scale_type,' not supported'
       call endrun()
    end if

    larr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do p = lbp,ubp
          if (wtlunit(p) /= 0._r8) then
             c = pcolumn(p)
             if (parr(p,j) /= spval .and. scale_c2l(c) /= spval) then
                l = plandunit(p)
                if (sumwt(l) == 0._r8) larr(l,j) = 0._r8
                larr(l,j) = larr(l,j) + parr(p,j) * scale_p2c(p) * scale_c2l(c) * wtlunit(p)
                sumwt(l) = sumwt(l) + wtlunit(p)
             end if
          end if
       end do
       found = .false.
       do l = lbl,ubl
          if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index = l
          else if (sumwt(l) /= 0._r8) then
             larr(l,j) = larr(l,j)/sumwt(l)
          end if
       end do
       if (found) then
          write(iulog,*)'p2l_2d error: sumwt is greater than 1.0 at l= ',index,' j= ',j
          call endrun()
       end if
    end do

  end subroutine p2l_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2g_1d
!
! !INTERFACE:
  subroutine p2g_1d(lbp, ubp, lbc, ubc, lbl, ubl, lbg, ubg, parr, garr, &
       p2c_scale_type, c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from pfts to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_gcell
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp            ! beginning and ending pft indices
    integer , intent(in)  :: lbc, ubc            ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl            ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg            ! beginning and ending gridcell indices
    real(r8), intent(in)  :: parr(lbp:ubp)       ! input pft array
    real(r8), intent(out) :: garr(lbg:ubg)       ! output gridcell array
    character(len=*), intent(in) :: p2c_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!  !LOCAL VARIABLES:
!EOP
    integer  :: pi,p,c,l,g,index       ! indices
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_p2c(lbp:ubp)     ! scale factor
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of pfts relative to gridcells
    integer , pointer :: pcolumn(:)    ! column of corresponding pft
    integer , pointer :: plandunit(:)  ! landunit of corresponding pft
    integer , pointer :: pgridcell(:)  ! gridcell of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in gridcell
    integer , pointer :: pfti(:)       ! initial pft index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    canyon_hwr => clm3%g%l%canyon_hwr
    ltype      => clm3%g%l%itype
    clandunit  => clm3%g%l%c%landunit
    ctype      => clm3%g%l%c%itype
    wtgcell   => clm3%g%l%c%p%wtgcell
    pcolumn   => clm3%g%l%c%p%column
    pgridcell => clm3%g%l%c%p%gridcell
    plandunit => clm3%g%l%c%p%landunit
    npfts     => clm3%g%npfts
    pfti      => clm3%g%pfti

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'p2g_1d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2g_1d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    garr(:) = spval
    sumwt(:) = 0._r8
    do p = lbp,ubp
       if (wtgcell(p) /= 0._r8) then
          c = pcolumn(p)
          l = plandunit(p)
          if (parr(p) /= spval .and. scale_c2l(c) /= spval .and. scale_l2g(l) /= spval) then
             g = pgridcell(p)
             if (sumwt(g) == 0._r8) garr(g) = 0._r8
             garr(g) = garr(g) + parr(p) * scale_p2c(p) * scale_c2l(c) * scale_l2g(l) * wtgcell(p)
             sumwt(g) = sumwt(g) + wtgcell(p)
          end if
       end if
    end do
    found = .false.
    do g = lbg, ubg
       if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = g
       else if (sumwt(g) /= 0._r8) then
          garr(g) = garr(g)/sumwt(g)
       end if
    end do
    if (found) then
       write(iulog,*)'p2g_1d error: sumwt is greater than 1.0 at g= ',index
       call endrun()
    end if

  end subroutine p2g_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: p2g_2d
!
! !INTERFACE:
  subroutine p2g_2d(lbp, ubp, lbc, ubc, lbl, ubl, lbg, ubg, num2d, &
       parr, garr, p2c_scale_type, c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from pfts to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !USES:
    use clm_varpar, only : max_pft_per_gcell
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbp, ubp              ! beginning and ending pft indices
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending gridcell indices
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(r8), intent(in)  :: parr(lbp:ubp,num2d)   ! input pft array
    real(r8), intent(out) :: garr(lbg:ubg,num2d)   ! output gridcell array
    character(len=*), intent(in) :: p2c_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,pi,p,c,l,g,index     ! indices
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_p2c(lbp:ubp)     ! scale factor
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of pfts relative to gridcells
    integer , pointer :: pcolumn(:)    ! column of corresponding pft
    integer , pointer :: plandunit(:)  ! landunit of corresponding pft
    integer , pointer :: pgridcell(:)  ! gridcell of corresponding pft
    integer , pointer :: npfts(:)      ! number of pfts in gridcell
    integer , pointer :: pfti(:)       ! initial pft index in gridcell
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    canyon_hwr   => clm3%g%l%canyon_hwr
    ltype        => clm3%g%l%itype
    clandunit    => clm3%g%l%c%landunit
    ctype        => clm3%g%l%c%itype
    wtgcell      => clm3%g%l%c%p%wtgcell
    pcolumn      => clm3%g%l%c%p%column
    pgridcell    => clm3%g%l%c%p%gridcell
    plandunit    => clm3%g%l%c%p%landunit
    npfts        => clm3%g%npfts
    pfti         => clm3%g%pfti

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'p2g_2d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    if (p2c_scale_type == 'unity') then
       do p = lbp,ubp
          scale_p2c(p) = 1.0_r8
       end do
    else
       write(iulog,*)'p2g_2d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    garr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do p = lbp,ubp
          if (wtgcell(p) /= 0._r8) then
             c = pcolumn(p)
             l = plandunit(p)
             if (parr(p,j) /= spval .and. scale_c2l(c) /= spval .and. scale_l2g(l) /= spval) then
                g = pgridcell(p)
                if (sumwt(g) == 0._r8) garr(g,j) = 0._r8
                garr(g,j) = garr(g,j) + parr(p,j) * scale_p2c(p) * scale_c2l(c) * scale_l2g(l) * wtgcell(p)
                sumwt(g) = sumwt(g) + wtgcell(p)
             end if
          end if
       end do
       found = .false.
       do g = lbg, ubg
          if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index = g
          else if (sumwt(g) /= 0._r8) then
             garr(g,j) = garr(g,j)/sumwt(g)
          end if
       end do
       if (found) then
          write(iulog,*)'p2g_2d error: sumwt gt 1.0 at g/sumwt = ',index,sumwt(index)
          call endrun()
       end if
    end do

  end subroutine p2g_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: c2l_1d
!
! !INTERFACE:
  subroutine c2l_1d (lbc, ubc, lbl, ubl, carr, larr, c2l_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to landunits
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc      ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl      ! beginning and ending landunit indices
    real(r8), intent(in)  :: carr(lbc:ubc) ! input column array
    real(r8), intent(out) :: larr(lbl:ubl) ! output landunit array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: ci,c,l,index           ! indices
    integer  :: max_col_per_lu         ! max columns per landunit; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor for column->landunit mapping
    real(r8) :: sumwt(lbl:ubl)         ! sum of weights
    real(r8), pointer :: wtlunit(:)    ! weight of landunits relative to gridcells
    integer , pointer :: clandunit(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in landunit
    integer , pointer :: coli(:)       ! initial column index in landunit
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtlunit    => clm3%g%l%c%wtlunit
    clandunit  => clm3%g%l%c%landunit
    ncolumns   => clm3%g%l%ncolumns
    coli       => clm3%g%l%coli

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'c2l_1d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    larr(:) = spval
    sumwt(:) = 0._r8
    do c = lbc,ubc
       if (wtlunit(c) /= 0._r8) then
          if (carr(c) /= spval .and. scale_c2l(c) /= spval) then
             l = clandunit(c)
             if (sumwt(l) == 0._r8) larr(l) = 0._r8
             larr(l) = larr(l) + carr(c) * scale_c2l(c) * wtlunit(c)
             sumwt(l) = sumwt(l) + wtlunit(c)
          end if
       end if
    end do
    found = .false.
    do l = lbl,ubl
       if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = l
       else if (sumwt(l) /= 0._r8) then
          larr(l) = larr(l)/sumwt(l)
       end if
    end do
    if (found) then
       write(iulog,*)'c2l_1d error: sumwt is greater than 1.0 at l= ',index
       call endrun()
    end if

  end subroutine c2l_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: c2l_2d
!
! !INTERFACE:
  subroutine c2l_2d (lbc, ubc, lbl, ubl, num2d, carr, larr, c2l_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to landunits
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc            ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl            ! beginning and ending landunit indices
    integer , intent(in)  :: num2d               ! size of second dimension
    real(r8), intent(in)  :: carr(lbc:ubc,num2d) ! input column array
    real(r8), intent(out) :: larr(lbl:ubl,num2d) ! output landunit array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,l,ci,c,index         ! indices
    integer  :: max_col_per_lu         ! max columns per landunit; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_c2l(lbc:ubc)        ! scale factor for column->landunit mapping
    real(r8) :: sumwt(lbl:ubl)         ! sum of weights
    real(r8), pointer :: wtlunit(:)    ! weight of column relative to landunit
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in landunit
    integer , pointer :: coli(:)       ! initial column index in landunit
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtlunit    => clm3%g%l%c%wtlunit
    clandunit  => clm3%g%l%c%landunit
    ncolumns   => clm3%g%l%ncolumns
    coli       => clm3%g%l%coli

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'c2l_2d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    larr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do c = lbc,ubc
          if (wtlunit(c) /= 0._r8) then
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval) then
                l = clandunit(c)
                if (sumwt(l) == 0._r8) larr(l,j) = 0._r8
                larr(l,j) = larr(l,j) + carr(c,j) * scale_c2l(c) * wtlunit(c)
                sumwt(l) = sumwt(l) + wtlunit(c)
             end if
          end if
       end do
       found = .false.
       do l = lbl,ubl
          if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index = l
          else if (sumwt(l) /= 0._r8) then
             larr(l,j) = larr(l,j)/sumwt(l)
          end if
       end do
       if (found) then
          write(iulog,*)'c2l_2d error: sumwt is greater than 1.0 at l= ',index,' lev= ',j
          call endrun()
       end if
    end do

  end subroutine c2l_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: c2g_1d
!
! !INTERFACE:
  subroutine c2g_1d(lbc, ubc, lbl, ubl, lbg, ubg, carr, garr, &
       c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending landunit indices
    real(r8), intent(in)  :: carr(lbc:ubc)         ! input column array
    real(r8), intent(out) :: garr(lbg:ubg)         ! output gridcell array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: ci,c,l,g,index         ! indices
    integer  :: max_col_per_gcell      ! max columns per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of columns relative to gridcells
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: cgridcell(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in gridcell
    integer , pointer :: coli(:)       ! initial column index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtgcell    => clm3%g%l%c%wtgcell
    clandunit  => clm3%g%l%c%landunit
    cgridcell  => clm3%g%l%c%gridcell
    ncolumns   => clm3%g%ncolumns
    coli       => clm3%g%coli

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'c2l_1d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    garr(:) = spval
    sumwt(:) = 0._r8
    do c = lbc,ubc
       if ( wtgcell(c) /= 0._r8) then
          l = clandunit(c)
          if (carr(c) /= spval .and. scale_c2l(c) /= spval .and. scale_l2g(l) /= spval) then
             g = cgridcell(c)
             if (sumwt(g) == 0._r8) garr(g) = 0._r8
             garr(g) = garr(g) + carr(c) * scale_c2l(c) * scale_l2g(l) * wtgcell(c)
             sumwt(g) = sumwt(g) + wtgcell(c)
          end if
       end if
    end do
    found = .false.
    do g = lbg, ubg
       if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = g
       else if (sumwt(g) /= 0._r8) then
          garr(g) = garr(g)/sumwt(g)
       end if
    end do
    if (found) then
       write(iulog,*)'c2g_1d error: sumwt is greater than 1.0 at g= ',index
       call endrun()
    end if

  end subroutine c2g_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: c2g_2d
!
! !INTERFACE:
  subroutine c2g_2d(lbc, ubc, lbl, ubl, lbg, ubg, num2d, carr, garr, &
       c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending gridcell indices
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(r8), intent(in)  :: carr(lbc:ubc,num2d)   ! input column array
    real(r8), intent(out) :: garr(lbg:ubg,num2d)   ! output gridcell array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,ci,c,g,l,index       ! indices
    integer  :: max_col_per_gcell      ! max columns per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_c2l(lbc:ubc)     ! scale factor
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of columns relative to gridcells
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: cgridcell(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in gridcell
    integer , pointer :: coli(:)       ! initial column index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(r8), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtgcell    => clm3%g%l%c%wtgcell
    clandunit  => clm3%g%l%c%landunit
    cgridcell  => clm3%g%l%c%gridcell
    ncolumns   => clm3%g%ncolumns
    coli       => clm3%g%coli

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0_r8
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l) 
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0_r8
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c) 
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0_r8
             end if
          else
             scale_c2l(c) = 1.0_r8
          end if
       end do
    else
       write(iulog,*)'c2g_2d error: scale type ',c2l_scale_type,' not supported'
       call endrun()
    end if

    garr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do c = lbc,ubc
          if (wtgcell(c) /= 0._r8) then
             l = clandunit(c)
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval .and. scale_l2g(l) /= spval) then
                g = cgridcell(c)
                if (sumwt(g) == 0._r8) garr(g,j) = 0._r8
                garr(g,j) = garr(g,j) + carr(c,j) * scale_c2l(c) * scale_l2g(l) * wtgcell(c)
                sumwt(g) = sumwt(g) + wtgcell(c)
             end if
          end if
       end do
       found = .false.
       do g = lbg, ubg
          if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index = g
          else if (sumwt(g) /= 0._r8) then
             garr(g,j) = garr(g,j)/sumwt(g)
          end if
       end do
       if (found) then
          write(iulog,*)'c2g_2d error: sumwt is greater than 1.0 at g= ',index
          call endrun()
       end if
    end do

  end subroutine c2g_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: l2g_1d
!
! !INTERFACE:
  subroutine l2g_1d(lbl, ubl, lbg, ubg, larr, garr, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from landunits to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbl, ubl       ! beginning and ending sub landunit indices
    integer , intent(in)  :: lbg, ubg       ! beginning and ending gridcell indices
    real(r8), intent(in)  :: larr(lbl:ubl)  ! input landunit array
    real(r8), intent(out) :: garr(lbg:ubg)  ! output gridcell array
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: li,l,g,index           ! indices
    integer  :: max_lu_per_gcell       ! max landunits per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of landunits relative to gridcells
    integer , pointer :: lgridcell(:)  ! gridcell of corresponding landunit
    integer , pointer :: nlandunits(:) ! number of landunits in gridcell
    integer , pointer :: luni(:)       ! initial landunit index in gridcell
!------------------------------------------------------------------------

    wtgcell    => clm3%g%l%wtgcell
    lgridcell  => clm3%g%l%gridcell
    nlandunits => clm3%g%nlandunits
    luni       => clm3%g%luni

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    garr(:) = spval
    sumwt(:) = 0._r8
    do l = lbl,ubl
       if (wtgcell(l) /= 0._r8) then
          if (larr(l) /= spval .and. scale_l2g(l) /= spval) then
             g = lgridcell(l)
             if (sumwt(g) == 0._r8) garr(g) = 0._r8
             garr(g) = garr(g) + larr(l) * scale_l2g(l) * wtgcell(l)
             sumwt(g) = sumwt(g) + wtgcell(l)
          end if
       end if
    end do
    found = .false.
    do g = lbg, ubg
       if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
          found = .true.
          index = g
       else if (sumwt(g) /= 0._r8) then
          garr(g) = garr(g)/sumwt(g)
       end if
    end do
    if (found) then
       write(iulog,*)'l2g_1d error: sumwt is greater than 1.0 at g= ',index
       call endrun()
    end if

  end subroutine l2g_1d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: l2g_2d
!
! !INTERFACE:
  subroutine l2g_2d(lbl, ubl, lbg, ubg, num2d, larr, garr, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from landunits to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbl, ubl             ! beginning and ending column indices
    integer , intent(in)  :: lbg, ubg             ! beginning and ending gridcell indices
    integer , intent(in)  :: num2d                ! size of second dimension
    real(r8), intent(in)  :: larr(lbl:ubl,num2d)  ! input landunit array
    real(r8), intent(out) :: garr(lbg:ubg,num2d)  ! output gridcell array
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,g,li,l,index         ! indices
    integer  :: max_lu_per_gcell       ! max landunits per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(r8) :: scale_l2g(lbl:ubl)     ! scale factor
    real(r8) :: sumwt(lbg:ubg)         ! sum of weights
    real(r8), pointer :: wtgcell(:)    ! weight of landunits relative to gridcells
    integer , pointer :: lgridcell(:)  ! gridcell of corresponding landunit
    integer , pointer :: nlandunits(:) ! number of landunits in gridcell
    integer , pointer :: luni(:)       ! initial landunit index in gridcell
!------------------------------------------------------------------------

    wtgcell   => clm3%g%l%wtgcell
    lgridcell => clm3%g%l%gridcell
    nlandunits => clm3%g%nlandunits
    luni       => clm3%g%luni

    call build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)

    garr(:,:) = spval
    do j = 1,num2d
       sumwt(:) = 0._r8
       do l = lbl,ubl
          if (wtgcell(l) /= 0._r8) then
             if (larr(l,j) /= spval .and. scale_l2g(l) /= spval) then
                g = lgridcell(l)
                if (sumwt(g) == 0._r8) garr(g,j) = 0._r8
                garr(g,j) = garr(g,j) + larr(l,j) * scale_l2g(l) * wtgcell(l)
                sumwt(g) = sumwt(g) + wtgcell(l)
             end if
          end if
       end do
       found = .false.
       do g = lbg,ubg
          if (sumwt(g) > 1.0_r8 + 1.e-6_r8) then
             found = .true.
             index= g
          else if (sumwt(g) /= 0._r8) then
             garr(g,j) = garr(g,j)/sumwt(g)
          end if
       end do
       if (found) then
          write(iulog,*)'l2g_2d error: sumwt is greater than 1.0 at g= ',index,' lev= ',j
          call endrun()
       end if
    end do

  end subroutine l2g_2d

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: build_scale_l2g
!
! !INTERFACE:
  subroutine build_scale_l2g(l2g_scale_type, lbl, ubl, scale_l2g)
!
! !DESCRIPTION:
! Fill the scale_l2g(lbl:ubl) array with appropriate values for the given l2g_scale_type.
! This array can later be used to scale each landunit in forming grid cell averages.
!
! !USES:
     use clm_varcon, only : max_lunit
!
! !ARGUMENTS:
     implicit none
     character(len=*), intent(in)  :: l2g_scale_type     ! scale factor type for averaging
     integer         , intent(in)  :: lbl, ubl           ! beginning and ending column indices
     real(r8)        , intent(out) :: scale_l2g(lbl:ubl) ! scale factor 
!
! !REVISION HISTORY:
! Created by Bill Sacks 10/11
!
!
! !LOCAL VARIABLES:
!EOP
     real(r8) :: scale_lookup(max_lunit) ! scale factor for each landunit type
     integer  :: l                       ! index
     integer , pointer :: ltype(:)       ! landunit type
!-----------------------------------------------------------------------

     ltype      => clm3%g%l%itype
     
     call create_scale_l2g_lookup(l2g_scale_type, scale_lookup)

     do l = lbl,ubl
        scale_l2g(l) = scale_lookup(ltype(l))
     end do

  end subroutine build_scale_l2g

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: create_scale_l2g_lookup
!
! !INTERFACE:
  subroutine create_scale_l2g_lookup(l2g_scale_type, scale_lookup)
! 
! DESCRIPTION:
! Create a lookup array, scale_lookup(1..max_lunit), which gives the scale factor for
! each landunit type depending on l2g_scale_type
!
! !USES:
     use clm_varcon, only : istsoil, istice, istdlak, istslak, istwet, isturb, istice_mec,&
                            istcrop, max_lunit, spval
!
! !ARGUMENTS:
     implicit none
     character(len=*), intent(in)  :: l2g_scale_type           ! scale factor type for averaging
     real(r8)        , intent(out) :: scale_lookup(max_lunit)  ! scale factor for each landunit type
!
! !REVISION HISTORY:
! Created by Bill Sacks 10/11
!
!EOP
!-----------------------------------------------------------------------

     ! ------------ WJS (10-14-11): IMPORTANT GENERAL NOTES ------------
     !
     ! Since scale_l2g is not currently included in the sumwt accumulations, you need to
     ! be careful about the scale values you use. Values of 1 and spval are safe
     ! (including having multiple landunits with value 1), but only use other values if
     ! you know what you are doing! For example, using a value of 0 is NOT the correct way
     ! to exclude a landunit from the average, because the normalization will be done
     ! incorrectly in this case: instead, use spval to exclude a landunit from the
     ! average. Similarly, using a value of 2 is NOT the correct way to give a landunit
     ! double relative weight in general, because the normalization won't be done
     ! correctly in this case, either.
     !
     ! In the longer-term, I believe that the correct solution to this problem is to
     ! include scale_l2g (and the other scale factors) in the sumwt accumulations
     ! (e.g., sumwt = sumwt + wtgcell * scale_p2c * scale_c2l * scale_l2g), but that
     ! requires some more thought to (1) make sure that is correct, and (2) make sure it
     ! doesn't break the urban scaling.
     !
     ! -----------------------------------------------------------------


     ! Initialize scale_lookup to spval for all landunits. Thus, any landunit that keeps
     ! the default value will be excluded from grid cell averages.
     scale_lookup(:) = spval

     if (l2g_scale_type == 'unity') then
        scale_lookup(:) = 1.0_r8
     else if (l2g_scale_type == 'veg') then
        scale_lookup(istsoil) = 1.0_r8
        scale_lookup(istcrop) = 1.0_r8
     else if (l2g_scale_type == 'ice') then
        scale_lookup(istice) = 1.0_r8
        scale_lookup(istice_mec) = 1.0_r8
     else if (l2g_scale_type == 'nonurb') then
        scale_lookup(:) = 1.0_r8
        scale_lookup(isturb) = spval
     else
        write(iulog,*)'scale_l2g_lookup_array error: scale type ',l2g_scale_type,' not supported'
        call endrun()
     end if

  end subroutine create_scale_l2g_lookup

end module subgridAveMod