c     file name: et.for

      subroutine   et (rn, bwudav, elev, bwtdmx, bwtdmn)

c     + + + PURPOSE + + +
c     This subroutine calculates daily potential evapotranspiration
c     using Van Bavel's (1966) revised combination method.
c     DATE:  09/17/93
c     MODIFIED:  10/06/93

c     + + + KEY WORDS + + +
c     evapotranspiration

c     + + + ARGUMENT DECLARATIONS + + +
      real rn
      real bwudav
      real uza
      real elev
      real bwtdmn
      real bwtdmx

c     + + + ARGUMENT DEFINITIONS + + +
c     rn     - Net radiation (mj/m^2/day)
c     bwudav      - Wind speed (m/s) at  10 m
c     uza    - Wind speed (m/s) at  height za
c     elev   - Elevation of the site (m)

c     + + + PARAMETERS + + +

      real a1, a2, a3
      parameter   (a1 = 824.4996, a2 = 35702.8022, a3 = -607945000.)

c     a1,a2,a3 - constants used to compute barametric pressure

      real b1, b2, b3
      parameter   (b1 = 67.5242, b2 = 149.531, b3 = -4859.0665)

c     b1,b2,b3 - constants used in eqn for (svpg0)

      real c1, c2, c3
      parameter   (c1 = 129.07487, c2 = 149.195, c3 = -4160.7968)

c     c1,c2,c3 - constants used to compute saturated vapor pressure

      real d1, d2
      parameter   (d1 = 2.5002773719, d2 = 0.0023644939)

c     d1,d2 - constants used to compute latent heat of vaporization

      real za
      parameter   (za = 2.0)

c     za - height of measurement of meteorological sensors (m)

      real e
      parameter   (e = 0.622)

c     e - water to air molecular weights ratio

      real vk
      parameter   (vk = 0.41)

c     vk - Von Karman's constant

      real z0
      parameter   (z0 = 0.0025)

c     z0 - roughness thickness (m)

c     + + + COMMON BLOCKS + + +
*$noereference
      include 'w1clig.inc'
      include 'h1et.inc'
*$reference

c     + + + LOCAL VARIABLES + + +
      real term1, term2, term3
      real bp
      real svpg
      real svpg0
      real vlh
      real vpa
      real vpd
      real vps
      real vpsmn
      real vpsmx
      real arho
      real ttc
c     real etpr
c     real etpw

c     + + + LOCAL DEFINITIONS + + +
c     term1  - Temporay local variables
c     bp     - Barometric pressure (kpa)
c     svpg   - Ratio of saturation vapor pressure curve slope
c              to the pychrometric constant, adjusted to
c              ambient barometric pressure (unitless)
c     svpg0  - Unadjusted ratio of the saturation vapor pressure
c              curve slope to the pychrometric constant (unitless)
c     vlh    - Latent heat of vaporization (mj/kg)
c     vpa    - Actual vapor pressure (kpa)
c     vpd    - Saturation vapor pressure deficit (kpa)
c     vps    - Saturated vapor pressure (kpa)
c     vpsmn  - Saturated vapor pressure at min air temp (kpa)
c     vpsmx  - Saturated vapor pressure at max air temp (kpa)
c     arho   - Air density (kg/m^3)
c     ttc    - Turbulent transfer coefficient (kg/m^2/kpa/day)
c     etpr   - Potential evapotransp due to radiation (mm/day)
c     etpw   - Potential evapotranspiration due to wind (mm/day)

c     + + + END SPECIFICATIONS + + +

C ***      write(*,*) 'et: bwtdmn, bwtdmx, awzdpt ',
C ***     *  bwtdmn, bwtdmx, awzdpt
      bp = a1*exp(((elev+a2)**2.)/a3)                    !h-15
      svpg0 = b1*exp(((awtdav-b2)**2.)/b3)               !h-14
      svpg = svpg0*(101.325/bp)                          !h-13
      vlh=d1 - (d2*awtdav)                               !h-24
      if (awzdpt .eq. 0.0) then
         ah0drat = 10.0
       else
         ah0drat = rn / (vlh * awzdpt)
      end if
      term1 = svpg * (rn/vlh)                            !h-12(a)
      vpa = c1*exp(((awtdpt-c2)**2.)/c3)                 !h-29
      vpsmn = c1*exp(((bwtdmn-c2)**2.)/c3)               !h-27
      vpsmx = c1*exp(((bwtdmx-c2)**2.)/c3)               !h-28
      vps = (vpsmn + vpsmx)/2.0                          !h-26
      vpd = vps - vpa                                    !h-25
      arho = 1000.*((bp/101.325)*((0.001293)/(1+(0.00367*awtdav))))  !h-32
      uza = bwudav * (za / 10.0) ** 0.142857
      ttc =(arho*e*(vk**2.)*uza*86400.)/(bp*(log(za/z0))**2.)        !h-30
      term2 = vpd*ttc                                    !h-12(b)
      term3 = svpg + 1.0                                 !h-12(c)
      ahzetp = (term1+term2)/term3                       !h-12
      if ( ahzetp  .le. 0.0 )  ahzetp = 0.0
      awrrh = vpa / vps
c
      return
      end