function uexa(x,t,kappa,l,io,task)

c General information:

c     Exact solution of the heat equation in a uniform bar of length l,
c               u_t = kappa u_xx
c     with a parabolic initial temperature profile at time t=0;
c               u(x,0) = 4 umax (1 - x/l) x/l   for  0 <= x <= l
c     (with umax a constant read in from file d_para.in),
c     and homogeneous Dirichlet boundary conditions at x=0 and x=l:
c               u(0,t) = u(l,t) = 0

c     Copyright 1996 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 12/16/96

c Usage information:

c     Function uexa looks for its parameters in a file 'd_para.in' of
c     the form:
c
c        <Value of umax>     <comment>
c
c     where umax is the initial temperature in the middle of the bar.
c        For better readability, comment lines may be inserted anywhere
c     in file d_para.in.  All lines are considered comments, unless the
c     first nonblank character is a digit, sign, or a decimal point.

c     Function uexa can be called in three different ways depending on
c     the chosen value of parameter task:

c     1. task = 0:
c        This is the standard call. In this case, uexa computes and
c     returns the exact temperature at position (x,t).
c        During such a call, all input parameters must have the
c     appropriate values.

c     2. task = -1:
c        This is the inquiry call.  Function uexa changes the value of
c     task to indicate its properties. To be precise, its sets task to
c     8+4=12 to indicate that the solution is antisymmetric around both
c     boundaries.

c     3. Any other value:
c        If task is positive, subroutine uexa merely initializes itself
c     and exits, otherwise it does nothing at all.

c Arguments:

c     Avoid typos:
      implicit none

c     Input: position x, (0 <= x <= l), at which the temperature is needed:
      double precision x

c     Input: time t, (0 <= t), at which the temperature is needed:
      double precision t

c     Input: conduction constant (0 <= kappa):
      double precision kappa

c     Input: length of the bar (0 < l):
      double precision l

c     Input: I/O unit of an already open output file or zero:
      integer io
c        Function uexa only writes to this unit if io is positive.

c     Input/Output: task to perform:
      integer task
c        See the usage information above for more information on task.

c     Output: The temperature at location x and time t (if task = 0):
      double precision uexa
c        If task is nonzero, the returned value uexa is always zero.

c External variables and info for compiling or changing function uexa:

c     The following utility routines from ../../../lib/util.f were used:

c        cwrite(module,text,io,show) writes line "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      external cwrite

c        fatal(module,text1,text2,text3) kills the program after a
c     fatal error, printing the lines "text1", "text2" and "text3".
      external fatal

c        fclose(module,filnam,io) closes the I/O unit "io" on which
c     file "filnam" is open.
      external fclose

c        fopen(module,filnam,io) opens the existing file "filnam" on
c     I/O unit "io". It selects "io" if initialized to zero.
      external fopen

c        rread(module,text,ioin,io,show) reads number "rread" from I/O
c     unit "ioin". It writes it and its description "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      double precision rread
      external rread

c        rwrite(module,text,val,io,show) writes number "val" and its
c     description "text" to I/O unit "io" if nonzero, and to the screen
c     if "show" is nonzero.
      external rwrite

c        warn(module,io,text1,text2,text3) writes the warning lines
c     "text1", "text2" and "text3" to the screen, and to I/O unit "io"
c     if nonzero.
      external warn

c Local variables:

c     An integer to keep track of whether uexa has been initialized:
      integer init

c     Integers that are returned during the inquiry call (see above):
      integer asyml,asymr
      parameter (asyml=4,asymr=8)

c     Maximum initial temperature:
      double precision umax
      save umax

c        The exact solution can be written as a Fourier series using the
c     standard separation-of-variable procedure.  Function uexa simply
c     sums that series. Used variables:

c     Index of the Fourier series:
      integer k
c     The same index as a floating point number (used since conversions
c     from integer to floating point may be very slow):
      double precision rk
c     Individual term in the Fourier series:
      double precision tk
c     Argument of the sine in the individual term:
      double precision args
c     Argument of the exponential factor in the individual term:
      double precision arge

c     A value for x above which exp(x) does not underflow (used since
c     underflow error messages can slow down execution greatly):
      double precision expufl
      parameter (expufl=-70.d0)

c     Maximum number of terms ever to sum in the Fourier series:
      integer klim
      parameter (klim=10000000)
c        Regardless of klim, the sum will terminate when converged.
c     Function uexa will crash if klim is exceeded.

c     Fortran I/O unit number for the input file:
      integer ioin

c     Constants defined to make changing precision easier:
      double precision zero,one,two,four,o2
      parameter (zero=0.d0,one=1.d0,two=2.d0,four=4.d0,o2=0.5d0)
      double precision pi
      save pi

c     Data statements:
      data init/0/

c Executable statements:

c     Default value of the temperature:
      uexa=zero

c For task = -1, return the properties of the solution:

      if(task.eq.-1)then
c        Both boundaries are anti-symmetric:
         task=asyml+asymr
         return
      endif

c Ignore other negative tasks:

      if(task.lt.0)return

c Initialization during the first time that out is called:

      if(init.eq.0)then

c        Show which function is being initialized:
         call cwrite('uexa',
     &   'Initialization of exact solution function uexa:',io,0)
         call cwrite('uexa',
     &   '- parabolic initial temperature distribution with',io,0)
         call cwrite('uexa',
     &   '  maximum temperature in the center of the bar;',io,0)
         call cwrite('uexa',
     &   '- zero temperature at the ends of the bar.',io,0)
         call cwrite('uexa',
     &   '- d_para/u_exa.f Copyright 1996 Leon van Dommelen.',io,0)

c        Value of pi:
         pi=four*atan(one)

c        Open the input file:
         ioin=0
         call fopen('uexa','d_para.in',ioin)

c        Read in the maximum temperature to use:
         umax=rread('uexa','The maximum initial temperature',ioin,io,1)

c        Accept a trivial solution:
         if(umax.eq.zero)call warn('uexa',io,
     &   'Zero initial temperature accepted.',' ',' ')

c        Close the input file:
         call fclose('uexa','d_para.in',ioin)

c        All done with the initialization:
         init=1

      endif

c Further ignore nonzero tasks:

      if(task.ne.0)return

c Check the arguments:

c     The Fourier series will not converge for negative kappa:
      if(kappa.lt.zero)then
         call rwrite('uexa',
     &   '    uexa: Heat conduction coefficient kappa',kappa,io,1)
         call fatal('uexa',
     &   'The heat conduction coefficient must be positive.',' ',' ')
      endif

c     Disallow negative length to simplify the code:
      if(l.le.zero)then
         call rwrite('uexa','    uexa: Bar length l',l,io,1)
         call fatal('uexa','The bar length must be positive.',' ',' ')
      endif

c     The Fourier series does not converge for negative times:
      if(t.lt.zero)then
         call rwrite('uexa','    uexa: Time t',t,io,1)
         call fatal('uexa','Negative times are not allowed.',' ',' ')
      endif

c     Disallow positions outside the bar to simplify the code:
      if(x.lt.zero .or. x.gt.l)then
         call rwrite('uexa','    uexa: Bar length l',l,io,1)
         call rwrite('uexa','    uexa:   Position x',x,io,1)
         call fatal('uexa',
     &   'Position x should be in the range from 0 to l.',' ',' ')
      endif

c Compute the temperature:

c     For umax = 0, return uexa as zero:
      if(umax.eq.zero)goto 900

c     For x = 0 or x = l, return the homogenuous boundary condition:
      if(x.eq.zero .or. x.eq.l)goto 900

c     For t = 0 or kappa = 0, return the parabolic initial profile:
      if(t.eq.zero .or. kappa.eq.zero)then
         uexa=four*umax*(one-x/l)*x/l
         goto 900
      endif

c     Else find the exact solution from summing the Fourier series:

c     Initialize the floating point value of the Fourier series index k:
      rk=zero

c     Sum up to klim terms:
      do 100 k=0,klim

c        Evaluate the argument of the sine, except for the factor x:
         args=(two*rk+one)*pi/l

c        Evaluate the argument of the exponential:
         arge=-args**2*kappa*t

c        Stop summing when underflow of the exponential is imminent:
         if(arge.lt.expufl)goto 900

c        The sin(args*x) term can make convergence hard to judge since
c        it can oscillate. So we will judge convergence based on its
c        amplitude alone, in the following way:

c        Evaluate the k-th term except for the sin(args*x) factor:
         tk=two*four**2*umax/(pi*(two*rk+one))**3*exp(arge)

c        Terminate when tk is too small to change the current value:
         if(tk+uexa.eq.uexa)goto 900

c        finish the k-th term:
         tk=tk*sin(args*x)

c        Add the term to the sum:
         uexa=uexa+tk

c        Increment the real value of k and go do the next term:
         rk=rk+one

 100  continue

c     We failed to get convergence since we did not jump out of the loop:
      call rwrite('uexa',
     &'    uexa: Last term in the Fourier series',tk,io,1)
      call rwrite('uexa',
     &'          Current value for the sum',uexa,io,1)
      call fatal('uexa','The Fourier series did not converge.',' ',' ')

c Exit:

c     Jump here when done:
 900  return
      end