c__________________________________________________________________

      real*8 Function rSTheis( Q, Tran, S, r, Time)
      real*8 Q, Tran, S, r, Time, Pi, Wu, t
      Pi = 3.141592654
c
      t = Time * Tran / S
c  W_1=Abs(h)*T/Q     r'=r/L    t'=tT/(SL**2)
      Wu = rW1(r, t)
      rSTheis = Q * Wu / Tran
c
      return
      end
c__________________________________________________________________
      real*8 Function rSUMTheis(rTimeQ, n1, n2, Tran,S,r, DateNOW, Qcon)
      include 'WLmodel.inc'
c      
      real*8 rTimeQ(nMAX1D,2)
      real*8 Tran,S,r, DateNOW
      integer n1, n2, n, i 
      real*8 DQ, Sadd, DT, Qprev, Qcon
c
      rSUMTheis = 0.000000000000000000000000
      If(rTimeQ(n1, 1) .lt. DateNOW )Then
        n = iMATCHrNset(DateNOW, rTimeQ, n1, n2)
        Qprev = 0.000000000000000000000000
        do i = n1, n
          DT = DateNOW - rTimeQ(i, 1) + 0.000000001
          DQ = (rTimeQ(i, 2) - Qprev) * Qcon
          Sadd = rSTheis(DQ, Tran, S, r, DT)
          rSUMTheis = rSUMTheis + Sadd
          Qprev = rTimeQ(i, 2)
        enddo 
      End If
c
      return
      end
c__________________________________________________________________
c  This calculates drawdowns for flow to a well from the Theis solution.
c  All variables are dimensionless and are defined as follows:
c        W_1=Abs(h)*T/Q     r'=r/L    t'=tT/(SL^2)
c  where Abs(h)=drawdown, Q=well flow, r=radial distance from the pumped well,
c  t=time, T=transmissivity, S=storage coefficient and L=length that may be
c  chosen arbitrarily since it cancels out in the Theis solution. For example,
c  choose L=one metre if units of metre are used in the problem. All input and
c  output variables are dimensionless, and the prime superscript has been omitted
c  for notational convenience.
      real*8 Function rW1(r, t)
      real*8 A0,A1,A2,A3,A4,A5, B0,B1,B2,B3,B4, C0,c1,c2,c3,C4
      real*8 P1, P2
      real*8 r, t, Pi, x
      Pi = 3.141592654
      A0 = -0.57721566
      A1 = 0.99999193
      A2 = -0.24991055
      A3 = 0.05519968 
      A4 = -0.00976004
      A5 = 0.00107857
      B0 = 0.2677737343
      B1 = 8.6347608925
      B2 = 18.059016973
      B3 = 8.5733287401
      B4 = 1
      C0 = 3.9584969228
      c1 = 21.0996530827
      c2 = 25.6329561486
      c3 = 9.5733223454
      C4 = 1
c
      If( t .le. 0 )Then
        rW1 = 0.000000000000000
      Else
        x = r**2 / (4. * t)
        If( x.le.1.0000 )Then
          x = -dLog(x) + A0 + x*(A1+x*(A2+x*(A3+x*(A4+x*A5))))
        Else
          P1 = B0 + x * (B1 + x * (B2 + x * (B3 + x * B4)))
          P2 = C0 + x * (c1 + x * (c2 + x * (c3 + x * C4)))
          x = (p1 / P2) * dExp(-x) / x
        End If
        rW1 = x / (4. * Pi)
c        rW_1 = rExp1(r**2 / (4. * t)) / (4. * Pi)
      End If
c
      return
      end      
