!********************************************************************
!
!  This subroutine is called from WATBAL to calculate daily
!     drainage flux (m/d) if:
!       1) drainage system exists and,
!       2) water table  is above drainage tiles.
!
!*******************************************************************    
subroutine drain(soldep, nowcrp,    nsl,  drseq,  drainq,   &
         &   solthk, unsdep,    satdep,  drawat, st,  ssc, dg,  ddrain, sdrain, &
         &   drdiam,  drainc, por, thetfc,  coca, fc)

implicit none
include 'constants.inc'

!******************************************************************
!     soldep - soil profile depth
!     ddrain - depth from surface to drainage tiles (m)
!     drainc - drainage coefficient (m/day)
!     drdiam - drain tile diameter (m)
!     unsdep - unsaturated depth from surface to water table (m)
!     sdrain - drain spacing (m)
!     drainq - drainage flux (m/day)
!     satdep - saturated depth from soil surface (m)
!     satdep -
!     nowcrp - current crop
!     drseq  - drainage sequence based on which crop and overland flow
!              element in use
!
!******************************************************************
           
      integer, intent(in) :: nowcrp,nsl,drseq(mxnsl)
      real, intent(in) :: soldep, solthk(mxnsl)
      real, intent(in) :: ssc(mxnsl+1),dg(mxnsl)
      real, intent(in) :: ddrain,sdrain,drdiam,drainc,por(mxnsl)
      real, intent(in) :: thetfc(mxnsl),coca(mxnsl),fc(mxnsl)

      real, intent(inout) :: drainq, unsdep, satdep, drawat(mxnsl)
      real, intent(inout) :: st(mxnsl)

!
!     + + + LOCAL VARIABLES + + +
      real drfc(mxnsl), d, de, dranks, temp, totdg, totk, wattbl, watyld
      integer jj, mn
!
!     + + + LOCAL DEFINITIONS + + +
!     drfc   - layer's field capacity, corrected for entrapped air
!     dranks - hydraulic conductivity for drainage calculation (m/d)
!
!
!
!      Calculate average soil hydraulic conductivity in saturated zone
!      (below water table).
!
      totk = 0.0
      totdg = 0.0
      do 10 mn = 1, nsl
        if (solthk(mn).gt.unsdep) then
          totk = totk + (ssc(mn)*dg(mn))
          totdg = totdg + dg(mn)
        else
          totk = 0.0
          totdg = 0.0
        end if
   10 continue
!
!     Get average Ks of saturated layers and convert to m/day.
!     Assume horizontal flow KZ = KY, where KY = vertical K.
!
      dranks = (totk/totdg) * 86400.
!
!
!     Calculate drainage flux (m/day) for tile drains
!
      d = soldep - ddrain(drseq(nowcrp))
      temp = d / sdrain(drseq(nowcrp))
!
!c    Drain diameter is needed for DE calculations,
!c    0.1 m for Oregon; ie, rr=0.1
!
!     Calculate "equivalent depth (DE) using DRAINMOD equations
!     2-13 to 2/15.
!     (WEPP Equation 7b.2.6)
!
      if (temp.le..3) then
        de = d / (1.0+temp*((8.0/3.14)*log(d/drdiam(drseq(nowcrp)))-3.4))
      else
        de = (sdrain(drseq(nowcrp))*3.14) / (8.0*(log( &
     &      sdrain(drseq(nowcrp))/drdiam(drseq(nowcrp)))-1.15))
      end if
!
      if (unsdep.lt.ddrain(drseq(nowcrp))) then
        wattbl = ddrain(drseq(nowcrp)) - unsdep
!       sdrain=distance between tiles, m
!       ------ subsurface flow to drain tubes
!       (WEPP Equation 7b.2.5)
        drainq = (8.0*dranks*de*wattbl) + (4.0*dranks*(wattbl**2)) / (sdrain(drseq(nowcrp))**2)
      else
        drainq = 0.0
      end if
!
!     ---- Limit drain flux to the hydraulic capacity of the drainage system.
      if (drainq.gt.drainc(drseq(nowcrp))) drainq = drainc(drseq(nowcrp))
!
!c    (Need to work on drainage ditch, and where main tile comes to the
!c    surface.   reza june 1990.
!
!     Convert drainage flux to soil water depth using average
!     drainable porosity (water yield).
      watyld = por(1) - (thetfc(1)+(1.0-coca(1)))
      unsdep = unsdep + (drainq/watyld)
      if (unsdep.lt.0.0) unsdep = 0.0
      satdep = soldep - unsdep
!
!     *** L0 IF ***
!     If there is drainage from the tiles (DRAINQ > 0), update the
!     available water content (ST) in each soil layer to reflect
!     the drainage from each layer, starting at the top.
!
      if (drainq.gt.0.0) then
!       *** Begin L1 DO-LOOP ***
        jj = 0
   20   continue
        jj = jj + 1
!       -------- layer's field capacity, corrected for entrapped air
        drfc(jj) = fc(jj) + ((1-coca(1))*dg(jj))
!
!       *** L2 IF ***
!       -------- If layer is above field capacity, drain it.
        if (st(jj).gt.drfc(jj)) then
!         ---------- potential water which can be drained from this layer
          drawat(jj) = st(jj) - drfc(jj)
!
!         *** L3 IF ***
!         ---------- If water runs out the drain...
          if (drainq.gt.0.0) then
!
!           ------------ If water excess in this layer exceeds or equals what has
!           run from drain...
            if (drawat(jj).ge.drainq) then
!             -------------- subtract amount that has run out the drain fron this layer.
              st(jj) = st(jj) - drainq
              if (st(jj).lt.1e-10) st(jj) = 1e-10
              drainq = 0.0
!           ------------ If water excess in this layer is less than what has run
!           from drain...
            else
!             -------------- subtract from DRAINQ what this layer can contribute.
              drainq = drainq - drawat(jj)
!             -------------- adjust soil layer water content by same amount of water
              st(jj) = drfc(jj)
            end if
!
!         *** L3 ENDIF ***
          end if
!
!       *** L2 ENDIF ***
        end if
!
!       *** End L1 DO-LOOP ***
        if ((jj.lt.nsl).and.(drainq.gt.0.0)) &
        &    goto 20
!
!
!     *** L0 ENDIF ***
      end if
!
      return
      end subroutine
!      end module