c$Author: fredfox $
c$Date: 2002-04-17 20:56:39 $
c$Revision: 1.4 $
c$Source: /weru/cvs/weps/weps.src/hydro/albedo.for,v $
      real function albedo (bcrlai, snwc, bsfalw,bsfald)

c     + + + purpose + + +
c     this subroutine estimates the net radiation for a given area
c     using known solar radiation, air temperature, and vapor pressure
c     according to wright's (1982) modified version of penman's (1948)
c     relationship.

c     + + + key words + + +
c     radiation, hydrology, weps

c     + + + argument declarations + + +
      real bcrlai
      real snwc
      real bsfalw
      real bsfald

c     + + + argument definitions + + +
c     bsfald   - dry albedo
c     bsfalw   - wet albedo
c     bcrlai   - plant leaf area index
c     snwc   - water content of snow, mm

c     + + + COMMON BLOCK + + +
      include 'p1const.inc'
      include 'p1werm.inc'
      include 'hydro/htheta.inc'

c     + + + local variables + + +
      real albs,albsn,albp
      real sci
      real snci
      real pci

c     + + + local definitions + + +
c     albedo - overall albedo value
c     albs   - soil albedo

c     + + + parameters + + +

c
      parameter   (
     *         albp = 0.23,						!albedo of plants
     *         albsn = 0.60						!albedo of snow
     *         )

c     albp   - plant albedo

c     + + + data initialization + + +

c     + + + end specifications + + +

c     estimate the surface albedo
      if ( snwc .ge. 5.0 ) then
        albedo = albsn                          !snow covers surface & plants
      else
        snci = snwc / 5.0                       !coverage factor for snow
        pci = min(bcrlai/3, 1.0)                !coverage factor for plants based upon leaf area index
        if (pci + snci .gt. 1.0) pci = 1.0 - snci  !make sure factors are sum to 1
        sci = 1.0 - (pci + snci)                !soil albedo factor is what is left over
        if (sci.gt.0.0) then					!need to calc soil albedo
          albs = bsfald + (bsfalw - bsfald) *
     *      (theta(0) - thetaw(1)) / (thetaf(1) - thetaw(1))
          albs = max(albs, bsfalw)   !no less than wet  (wet is less than dry)
          albs = min(albs, bsfald)   !no greater than dry
          albedo = snci*albsn + pci * albp + sci * albs
        else
          albedo = snci*albsn + pci * albp
        endif
      endif

      return
      end