real function exact(x,t,nu,l,V_0,V_1,V_2,tol)

c     The exact solution to the unidirectional flow Navier-Stokes equation:

c                            v_t = nu v_xx

c     with boundary conditions:

c                v(0,t) = V_1            v(l,t) = V_2

c     and initial condition:

c              v(x,0) = V_0  4 (x/l) (1 - (x/l)) + V_2

c     the exact solution is normally found by summing the Separation
c     of Variables series to the desired accuracy tol.

c     avoid typos and duplicate names
      implicit none

c input arguments

c     value of x and t at which the velocity is desired
      real x,t

c     coefficient of kinematic viscosity
      real nu

c     length of the flow region
      real l

c     constants in the given boundary conditions and initial condition
      real V_0,V_1,V_2

c     allowable error.  The simplest is to specify this as zero; then
c     the solution will computed as accurately as possible.
      real tol

c local variables

c     Fourier mode number
      integer k
c     the value of k as a real
      real rk

c     the amplitude of the k-th Fourier mode
      real amplit

c     the total k-th term of the Fourier series
      real term

c     pi  (the correct value 3.14... of pi will be set later)
      real pi
      data pi/0./
      save pi

c executable statements

c     compute pi if it has not been done yet
      if(pi.eq.0.)pi=4.*atan(1.)

c     on the boundaries, the solution is given
      if(x.le.0)then
         exact=V_1
         return
      endif
      if(x.ge.l)then
         exact=V_2
         return
      endif

c     else, at t=0, the solution is also given
      if(t.le.0.)then
         exact=4.*V_0*(x/l)*(1.-(x/l))+V_2
         return
      endif

c     initialize the exact solution to zero
      exact=0.

c     we will be computing the exact solution in 3 parts; the first part
c     is the SOV series due to the V_0 part of the initial condition
c     alone:

c     initial value of k as a real number (why do we need it?)
      rk=1.

c     sum up to 1,000,000 Fourier modes together
      do 100 k=1,2000000,2

c        the k-th term of the sum, except for the sin(k pi x/l) factor
         amplit=(32.*V_0/pi**3)/rk**3*exp(-t*nu*(pi/l)**2*rk**2)

c        the total k-th term to add (why is pi*x/l between parentheses??)
         term=amplit*sin(rk*(pi*x/l))

c        done if we have have the requested accuracy
         if(abs(amplit*rk).le.tol)goto 110
c        Why look at the amplitude alone, not the complete 'term' below??
c        Why is the *rk there???

c        done in any case if maximum possible accuracy has been achieved
         if(amplit+exact.eq.exact)goto 110
c        Does this make sense?????????????????????????  Explain!

c        add the k-th term to the Fourier sum
         exact=exact+term

c        increment the real value of k (is this correct?)
         rk=rk+2.

 100  continue

      print*,'*** exact: First sum not converged after 1,000,000 terms'
      stop

 110  continue

c     the second part of the solution is that due to the V_2 initial
c     condition, V_2 boundary condition at x=l, and also assuming a
c     V_2 boundary condition at x=0.  The solution is simply V_2.

      exact=exact+V_2

c     the last part is the solution due to the missing V_1-V_2 boundary
c     condition at x=0, with zero initial conditions and zero boundary
c     condition at x=l.  This solution itself consists of two parts;
c     the first is the large time solution satisfying the inhomogeneous
c     boundary condition:

      exact=exact+(V_1-V_2)*(1.-(x/l))

c     the other, final part is the SOV solution to the problem with
c     homogeneous boundary conditions and a -(V_1-V_2)*(1.-(x/l))
c     initial condition

c     initial value of k as a real number
      rk=1.

c     sum up to 1,000,000 Fourier modes together
      do 200 k=1,1000000

c        the k-th term of the sum, except for the sin(k pi x/l) factor
         amplit=-(2.*(V_1-V_2)/pi)/rk*exp(-t*nu*(pi/l)**2*rk**2)

c        the total k-th term to add
         term=amplit*sin(rk*(pi*x/l))

c        done if we have have the requested accuracy
         if(abs(amplit*rk).le.tol)goto 210

c        done in any case if maximum possible accuracy has been achieved
         if(amplit+exact.eq.exact)goto 210

c        add the k-th term to the Fourier sum
         exact=exact+term

c        increment the real value of k
         rk=rk+1.

 200  continue

      print*,'*** exact: 2nd sum not converged after 1,000,000 terms'
      stop

 210  continue

      return
      end