c$Author: fredfox $
c$Date: 2001-07-03 19:05:16 $
c$Revision: 1.16.2.1 $
c$Source: /weru/cvs/weps/weps.src/hydro/transp.for,v $
      subroutine   transp (
     i                     layrsn,
     i                     bszlyd,
     i                     bszlyt,
     i                     prtd,
     o                     wsf)

c     + + + PURPOSE + + +
c     This subroutine determines the actual plant transpiration first
c     by distributing the potential rate of plant transpiration
c     throughout the root zone.  the actual plant transpiration is
c     then obtained by adjusting the potential plant traspiration
c     on the basis of soil water availability.
c     DATE:  09/22/93
c     MODIFIED:  10/06/93
c     MODIFIED:  08/03/95
c     MODIFIED:  05/01/99 - LEW

c     + + + KEYWORDS + + +
c     transpiration

c     + + + ARGUMENT DECLARATIONS + + +
      integer layrsn
      real bszlyd(*)
      real bszlyt(*)
      real prtd
      real wsf

c     + + + ARGUMENT DEFINITIONS + + +
c     layrsn - Number of soil layers used in simulation
c     bszlyd  - Depth to bottom of soil layer from surface (mm)
c     prtd   - Plant root depth (m)
c     wsf    - Plant growth water stress factor (unitless)
c     bszlyt - Layer thickness (mm)

c     + + + PARAMETERS + + +
      real   wud
      parameter   (wud = 3.0650)
      real   wuc
      parameter   (wuc = 1.0)

c     wud - Water use distribution, a depth parameter of 3.065 is
c           used in weps assuming about 30 percent of the total
c           water use comes from the top 10 percent of the root
c           zone
c     wuc - Water use compensation factor, determines how much a plant
c           can draw water from lower soil layers when higher soil layers
c           are dry. wuc = 0 means water will be withdrawn according to
c           the wud distribution without compensation if upper layers are 
c           dry. wuc = 1 means that if water is available in lower layer, 
c           more water than indicated by wud will be withdrawn.

c     + + + FUNCTIONS CALLED + + +

      real availwc
      real acplwu

c     availwc - Available water content (mm/mm)
c     acplwu - Actual water use rate from soil layer (mm/day)

c     + + + COMMON BLOCKS + + +

      include 'p1werm.inc'
      include 'p1unconv.inc'
      include 'h1et.inc'

c     + + + LOCAL COMMON BLOCKS + + +

      include 'hydro/htheta.inc'
C ***      include 'hydro/hlayrs.inc'

c     + + + LOCAL VARIABLES + + +

      integer k
      real depth
      real awcr(mnsz)
      real wua(mnsz)
      real wup(mnsz)
      real twu
      real wup_fac(0:mnsz)
      real temp

c     + + + LOCAL DEFINITIONS + + +

c     k      - Local loop variable
c     awcr   - Relative available soil wc, fraction (0-1.0)
c     wua    - Actual water use rate from soil layer (mm/day)
c     wup    - Potential water use rate from soil layer (mm/day)
c     twu    - Accumulated actual water use from the overlying
c              soil layers (mm)

c     + + + DATA INITIALIZATION + + +
c
c     + + + END SPECIFICATIONS + + +

      !handle special case of no roots yet
      if (prtd .le. 0.0) then
          awcr(1) = availwc(theta(1), thetaw(1), thetaf(1))
          twu = acplwu(awcr(1), ahzptp)
c         update soil water content in layer
          theta(1) = theta(1) - twu/bszlyt(1)
          goto 999  !finish up
      endif

      twu = 0.0
      wup_fac(0) = 0.0
      do 100 k= 1, layrsn
         !compute transpiration in root zone only
         if ((prtd*mtomm) .ge. bszlyd(k)) then
             depth = bszlyd(k)
         else if ((prtd*mtomm) .gt. (bszlyd(k)-bszlyt(k))) then
             depth = prtd*mtomm
         else
             goto 999  !we are done with layers having root mass
         endif
         wup_fac(k)=(1.0-exp(-wud*depth/(prtd*mtomm)))/(1.0-exp(-wud))
         wup(k) = ahzptp*(wup_fac(k)-(1.0-wuc)*wup_fac(k-1)) - wuc*twu
         awcr(k) = availwc(theta(k), thetaw(k), thetaf(k))
         wua(k) = acplwu(awcr(k), wup(k))
c        update soil water content in layer
c        prevent going beyond wilting point
         if( wua(k).gt.0.0 ) then
             temp = theta(k)-wua(k)/bszlyt(k)
             if( temp.lt.thetaw(k)) then
                 wua(k) = (theta(k)-thetaw(k))*bszlyt(k)
                 theta(k) = thetaw(k)
             else
                 theta(k) = temp
             end if
             twu = twu + wua(k)
          end if
100   continue

999   continue
      ahzpta = twu  !set WEPS global variable (accumulated actual water use)

      !Since ahzptp is set to zero when a crop is in dormancy, eg.,
      !winter wheat, we check that condition here.
      !However, this routine really shouldn't be called if a crop
      !is in dormancy because it isn't using any water - LEW 5/1/99
      if (ahzptp .eq. 0.0) then
         wsf = 1.0              !set it to not negatively influence growth
      else
         wsf = ahzpta / ahzptp  !set water stress factor
      endif

c      write(*,*) ' transp: ahzpta, ahzptp, wsf, prtd ',
c     *  ahzpta, ahzptp, wsf, prtd, bszlyd(1)

      return
      end