      real*8 Function RORAq( t, h0, Tran, Sy, a )
      real*8  t, h0, Tran, Sy, a 
      real*8 c1,c2, aD, tD, Pi, dZero
      integer m 
      Pi = 3.141592654
      dZero = 1.0E-34
c'
      c2 = -Pi ** 2 * Tran * t / 4 / a ** 2 / Sy
      RORAq = 0.000000
      m = 1
      aD = 1.
      do while( Abs(aD).gt.dZero .AND. m.lt.333)
        c1 = c2 * m ** 2
        aD = dExp(c1)
        RORAq = RORAq + aD
        m = m+2
      end do
c'
      RORAq = RORAq * 2 * Tran * h0 / a
c'
      return
      end
C ______________________________________________________________________________________________
c
      real*8 Function RORAhd( t, x, h0, Tran, Sy, a )
      real*8  t, x, h0, Tran, Sy, a 
      real*8 c1,c2, aD, tD, Pi, dZero
      integer m 
      Pi = 3.141592654
      dZero = 1.0E-34
c'
c'c    Solves for a 1-D, dimensionless solution of head as a function of time and
c'c  distance with no-flow boundary at x=0 and a constant head boundary at x=1.0.
c'c  The initial conditionis h = 1.0 everywhere.
c'
      xD = 1.000000 - x / a
      tD = t * Tran / Sy / a ** 2
c'c
      If( tD.gt.dZero ) Then
        RORAhd = 0.000000
        m = 0
        aD = 1.
        do while( Abs(aD).gt.dZero .AND. m.lt.333)
          m = m + 1
          e = (m - 0.5) * Pi
          aD = (-1) ** (m + 1) / e * Exp(-e ** 2 * tD) * Cos(e * xD)
          RORAhd = RORAhd + aD
        end do
        RORAhd = h0 * RORAhd * 2
      Else
        RORAhd = h0
      End If
c'
      return
      end
c______________________________________________________________________________________________
c
      real*8 Function SpcHD1D(t, x, h0, Tran, Sy, a)
      real*8  t, x, h0, Tran, Sy, a 
      real*8 c1, aD, e, tD, Pi, dZero
      integer m, mMAX 
      Pi = 3.141592654
      dZero = 1.0E-34
      mMAX = 1000
c'
c'c    Solves for a 1-D, dimensionless solution of head as a function of time and
c'c  distance with constant head boundaries and initial conditions of h0
c'c  everywhere.
c'
      xD = 1.00000000000000 - x / a
      tD = t * Tran / Sy / a ** 2
c'c
      If( tD.gt.dZero )Then
        SpcHD1D = 0.000000000000000000000
        m = 0
        aD = 1.
        do while( Abs(aD).gt.dZero .AND. m.lt.mMAX ) 
          m = m + 1
          e = m * Pi
          aD = (-1) ** (m + 1) / e * dExp(-e ** 2 * tD) * dSin(e * xD)
          SpcHD1D = SpcHD1D + aD
        end do 
        SpcHD1D = 1.000000 - xD + SpcHD1D * 2
        SpcHD1D = h0 * SpcHD1D
      Else
        SpcHD1D = h0
      End If
c'
      return
      end
C ______________________________________________________________________________________________
c
      real*8 Function SpcHD1Dq( t, h0, Tran, Sy, a)
      real*8  t, h0, Tran, Sy, a 
      real*8 c1, aD, tD, Pi, dZero
      integer m 
      Pi = 3.141592654
      dZero = 1.0E-34
c'
      tD = t * Tran / Sy / a ** 2
c'
      SpcHD1Dq = 0.000000
      m = 1
      aD = 1.
      do while( Abs(aD).gt.dZero .AND. m.lt.100)
        e = (m * Pi) ** 2
        aD = dExp(-e * tD)
        SpcHD1Dq = SpcHD1Dq + aD
        m = m+2
      end do
c'
      SpcHD1Dq = SpcHD1Dq * 2 * Tran * h0 / a
c'
      return
      end
C ______________________________________________________________________________________________
c
      real*8 Function rSUMSpcHD1D(rTimeQ, n1, n2, Tran,Sy,x,a, DateNOW)
      include 'WLmodel.inc'
c      
      real*8 rTimeQ(nMAX1D,2)
      real*8 Tran,Sy,x,a, DateNOW
      integer n1, n2, n, i 
      real*8 Dh0, Sadd, DT, H0prev
c
      rSUMSpcHD1D = 0.000000000000000000000000
      If(rTimeQ(n1, 1) .lt. DateNOW )Then
        n = iMATCHrNset(DateNOW, rTimeQ, n1, n2)
        H0prev = 0.000000000000000000000000
        do i = n1, n
          DT = DateNOW - rTimeQ(i, 1) + 0.000000001
          Dh0 = rTimeQ(i, 2) - H0prev
          Sadd = Dh0 - SpcHD1D(DT, x, Dh0, Tran, Sy, a)
          rSUMSpcHD1D = rSUMSpcHD1D + Sadd
          H0prev = rTimeQ(i, 2)
        enddo 
      End If
c
      return
      end
C ______________________________________________________________________________________________
c
      real*8 Function rSUMroraHD(rTimeQ, n1, n2, Tran,Sy,x,a, DateNOW)
      include 'WLmodel.inc'
c      
      real*8 rTimeQ(nMAX1D,2)
      real*8 Tran,Sy,x,a, DateNOW
      integer n1, n2, n, i 
      real*8 Dh0, Sadd, DT, H0prev
c
      rSUMroraHD = 0.000000000000000000000000
      If(rTimeQ(n1, 1) .lt. DateNOW )Then
        n = iMATCHrNset(DateNOW, rTimeQ, n1, n2)
        H0prev = rTimeQ(n1, 2)
        do i = n1, n
          DT = DateNOW - rTimeQ(i, 1) + 0.000000001
          Dh0 = rTimeQ(i, 2) - H0prev
          Sadd = Dh0 - RORAhd(DT, x, Dh0, Tran, Sy, a)
          rSUMroraHD = rSUMroraHD + Sadd
          H0prev = rTimeQ(i, 2)
        enddo 
      End If
c
      return
      end
C ______________________________________________________________________________________________
c
      real*8 Function rSUMroraQ(rTimeQ, n1, n2, Tran,Sy,a, DateNOW)
      include 'WLmodel.inc'
c      
      real*8 rTimeQ(nMAX1D,2)
      real*8 Tran,Sy, a, DateNOW
      integer n1, n2, n, i 
      real*8 Dh0, Qadd, DT, H0prev
c
      rSUMroraQ = 0.000000000000000000000000
      If(rTimeQ(n1, 1) .lt. DateNOW )Then
        n = iMATCHrNset(DateNOW, rTimeQ, n1, n2)
        H0prev = 0.000000000000000000000000
        do i = n1, n
          DT = DateNOW - rTimeQ(i, 1) + 0.000000001
          Dh0 = rTimeQ(i, 2) - H0prev
          Qadd = RORAq(DT, Dh0, Tran, Sy, a)
          rSUMroraQ = rSUMroraQ + Sadd
          H0prev = rTimeQ(i, 2)
        enddo 
      End If
c
      return
      end
