*     D03PWF Example Program Text
*     Mark 19 Revised. NAG Copyright 1999.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NPDE, NPTS, NCODE, NXI, NEQN, NIW, NWKRES,
     +                 LENODE, MLU, NW
      PARAMETER        (NPDE=3,NPTS=141,NCODE=0,NXI=0,
     +                 NEQN=NPDE*NPTS+NCODE,NIW=NEQN+24,
     +                 NWKRES=NPDE*(2*NPTS+3*NPDE+32)+7*NPTS+4,
     +                 LENODE=9*NEQN+50,MLU=3*NPDE-1,NW=(3*MLU+1)
     +                 *NEQN+NWKRES+LENODE)
*     .. Scalars in Common ..
      DOUBLE PRECISION EL0, ER0, GAMMA, RL0, RR0, UL0, UR0
*     .. Local Scalars ..
      DOUBLE PRECISION D, P, TOUT, TS, V
      INTEGER          I, IFAIL, IND, ITASK, ITOL, ITRACE, K
      CHARACTER        LAOPT, NORM
*     .. Local Arrays ..
      DOUBLE PRECISION ALGOPT(30), ATOL(1), RTOL(1), U(NPDE,NPTS),
     +                 UE(3,9), W(NW), X(NPTS), XI(1)
      INTEGER          IW(NIW)
*     .. External Subroutines ..
      EXTERNAL         BNDARY, D03PEK, D03PLF, D03PLP, NUMFLX
*     .. Common blocks ..
      COMMON           /INIT/EL0, ER0, RL0, RR0, UL0, UR0
      COMMON           /PARAMS/GAMMA
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D03PWF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
*
*     Problem parameters
*
      GAMMA = 1.4D0
      RL0 = 5.99924D0
      RR0 = 5.99242D0
      UL0 = 5.99924D0*19.5975D0
      UR0 = -5.99242D0*6.19633D0
      EL0 = 460.894D0/(GAMMA-1.0D0) + 0.5D0*RL0*19.5975D0**2
      ER0 = 46.095D0/(GAMMA-1.0D0) + 0.5D0*RR0*6.19633D0**2
*
*     Initialise mesh
*
      DO 20 I = 1, NPTS
         X(I) = 1.0D0*(I-1.0D0)/(NPTS-1.0D0)
   20 CONTINUE
*
*     Initial values
*
      DO 40 I = 1, NPTS
         IF (X(I).LT.0.5D0) THEN
            U(1,I) = RL0
            U(2,I) = UL0
            U(3,I) = EL0
         ELSE IF (X(I).EQ.0.5D0) THEN
            U(1,I) = 0.5D0*(RL0+RR0)
            U(2,I) = 0.5D0*(UL0+UR0)
            U(3,I) = 0.5D0*(EL0+ER0)
         ELSE
            U(1,I) = RR0
            U(2,I) = UR0
            U(3,I) = ER0
         END IF
   40 CONTINUE
*
      ITRACE = 0
      ITOL = 1
      NORM = '2'
      ATOL(1) = 0.5D-2
      RTOL(1) = 0.5D-3
      XI(1) = 0.0D0
      LAOPT = 'B'
      IND = 0
      ITASK = 1
      DO 60 I = 1, 30
         ALGOPT(I) = 0.0D0
   60 CONTINUE
*
*     Theta integration
*
      ALGOPT(1) = 2.0D0
      ALGOPT(6) = 2.0D0
      ALGOPT(7) = 2.0D0
*
*     Max. time step
*
      ALGOPT(13) = 0.5D-2
*
      TS = 0.0D0
      TOUT = 0.035D0
      IFAIL = 1
*
      CALL D03PLF(NPDE,TS,TOUT,D03PLP,NUMFLX,BNDARY,U,NPTS,X,NCODE,
     +            D03PEK,NXI,XI,NEQN,RTOL,ATOL,ITOL,NORM,LAOPT,ALGOPT,W,
     +            NW,IW,NIW,ITASK,ITRACE,IND,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,99998) TS
         WRITE (NOUT,99999)
*
*        Read exact data at output points
*
         DO 80 I = 1, 9
            READ (NIN,*) UE(1,I), UE(2,I), UE(3,I)
   80    CONTINUE
*
*        Calculate density, velocity and pressure
*
         K = 0
         DO 100 I = 15, NPTS - 14, 14
            D = U(1,I)
            V = U(2,I)/D
            P = D*(GAMMA-1.0D0)*(U(3,I)/D-0.5D0*V**2)
            K = K + 1
            WRITE (NOUT,99996) X(I), D, UE(1,K), V, UE(2,K), P, UE(3,K)
  100    CONTINUE
*
         WRITE (NOUT,99997) IW(1), IW(2), IW(3), IW(5)
      ELSE
         WRITE (NOUT,99995) IFAIL
      END IF
*
99999 FORMAT (4X,'X',6X,'APPROX D',3X,'EXACT D',4X,'APPROX V',3X,'EXAC',
     +       'T V',4X,'APPROX P',3X,'EXACT P')
99998 FORMAT (/' T = ',F6.3,/)
99997 FORMAT (/' Number of integration steps in time = ',I6,/' Number ',
     +       'of function evaluations = ',I6,/' Number of Jacobian ',
     +       'evaluations =',I6,/' Number of iterations = ',I6)
99996 FORMAT (1X,E8.2,6(1X,E10.4))
99995 FORMAT (1X,/1X,' ** D03PLF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE BNDARY(NPDE,NPTS,T,X,U,NCODE,V,VDOT,IBND,G,IRES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          IBND, IRES, NCODE, NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION G(NPDE), U(NPDE,NPTS), V(*), VDOT(*), X(NPTS)
*     .. Scalars in Common ..
      DOUBLE PRECISION EL0, ER0, RL0, RR0, UL0, UR0
*     .. Common blocks ..
      COMMON           /INIT/EL0, ER0, RL0, RR0, UL0, UR0
*     .. Executable Statements ..
      IF (IBND.EQ.0) THEN
         G(1) = U(1,1) - RL0
         G(2) = U(2,1) - UL0
         G(3) = U(3,1) - EL0
      ELSE
         G(1) = U(1,NPTS) - RR0
         G(2) = U(2,NPTS) - UR0
         G(3) = U(3,NPTS) - ER0
      END IF
      RETURN
      END
*
      SUBROUTINE NUMFLX(NPDE,T,X,NCODE,V,ULEFT,URIGHT,FLUX,IRES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T, X
      INTEGER          IRES, NCODE, NPDE
*     .. Array Arguments ..
      DOUBLE PRECISION FLUX(NPDE), ULEFT(NPDE), URIGHT(NPDE), V(*)
*     .. Scalars in Common ..
      DOUBLE PRECISION GAMMA
*     .. Local Scalars ..
      INTEGER          IFAIL
*     .. External Subroutines ..
      EXTERNAL         D03PWF
*     .. Common blocks ..
      COMMON           /PARAMS/GAMMA
*     .. Executable Statements ..
      IFAIL = 0
      CALL D03PWF(ULEFT,URIGHT,GAMMA,FLUX,IFAIL)
      RETURN
      END