c__________________________________________________________________

      real*8 Function dFactor(iNfac)
      integer*4 iNfac, i 
c'
      dFactor = 1.
      do i = 1, iNfac
        dFactor = dFactor * i
      end do
c
      return
      end
c__________________________________________________________________

      real*8 Function dPrecipINC(dDELT, dK, iNfac)
      real*8 dDELT, dK, dDT
      integer*4 iNfac
c'
      dDT = dDELT / dK
      dPrecipINC = dFactor(iNfac - 1)
      dPrecipINC = Exp(-dDT) * dDT ** (iNfac - 1) / (dK * dPrecipINC)
c
      return
      end
c__________________________________________________________________

      real*8 Function dSUMPrecipINC(rTimeP, n1,n2, dK, iNfac, rDateNOW)
      include 'WLmodel.inc'
c      
      real*8 rTimeP(nMAX1D,2)
      integer*4 iNfac, n, i, n1, n2
      real*8 dWLadd, dDT, dPulse, dK, rDateNOW
c  
      dZERO = 1E-30
c'
      dSUMPrecipINC = 0.
      n = iMATCHrNset(rDateNOW, rTimeP, n1, n2)
      do i = n1, n
        dPulse = rTimeP(i, 2)
        If(dPulse ** 2 .gt. dZERO ) Then
          dPulse = dPulse * (rTimeP(i + 1, 1) - rTimeP(i, 1))
          dDT = rDateNOW - rTimeP(i, 1) + dZERO
          dWLadd = dPulse * dPrecipINC(dDT, dK, iNfac)
          dSUMPrecipINC = dSUMPrecipINC + dWLadd
        End If
      enddo
c
      return
      end
c__________________________________________________________________
