*     D03RAF Example Program Text
*     Mark 19 Revised. NAG Copyright 1999.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D03RAF Example Program Results'
      CALL EX1
      CALL EX2
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          MXLEV, NPDE, NPTS
      PARAMETER        (MXLEV=3,NPDE=1,NPTS=2000)
      INTEGER          LENIWK, LENRWK, LENLWK
      PARAMETER        (LENIWK=NPTS*(5*MXLEV+14)+2+7*MXLEV,
     +                 LENRWK=NPTS*NPDE*(5*MXLEV+9+18*NPDE)+NPTS*2,
     +                 LENLWK=NPTS+1)
*     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA, D, DELTA, DIFF, REAC
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      DOUBLE PRECISION TOLS, TOLT, TOUT, TS, XMAX, XMIN, YMAX, YMIN
      INTEGER          I, IFAIL, IND, ITRACE, J, MAXLEV, NX, NY
*     .. Local Arrays ..
      DOUBLE PRECISION DT(3), OPTR(3,NPDE), RWK(LENRWK)
      INTEGER          IWK(LENIWK), OPTI(4)
      LOGICAL          LWK(LENLWK)
*     .. External Subroutines ..
      EXTERNAL         BNDRY1, D03RAF, MONIT1, PDEF1, PDEIV1
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Common blocks ..
      COMMON           /OTIME1/TWANT, IOUT
      COMMON           /PARAM1/ALPHA, DELTA, REAC, DIFF, D
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 1'
      WRITE (NOUT,*)
*
*     Problem Parameters
*
      ALPHA = 1.0D0
      DELTA = 20.0D0
      REAC = 5.0D0
      DIFF = 0.1D0
      D = REAC*EXP(DELTA)/(ALPHA*DELTA)
*
      IND = 0
      ITRACE = 0
      TS = 0.0D0
      DT(1) = 0.1D-2
      DT(2) = 0.0D0
      DT(3) = 0.0D0
      TOUT = 0.24D0
      TWANT(1) = 0.24D0
      TWANT(2) = 0.25D0
      XMIN = 0.0D0
      XMAX = 1.0D0
      YMIN = 0.0D0
      YMAX = 1.0D0
      NX = 21
      NY = 21
      TOLS = 0.5D0
      TOLT = 0.01D0
      DO 20 I = 1, 4
         OPTI(I) = 0
   20 CONTINUE
      DO 60 J = 1, NPDE
         DO 40 I = 1, 3
            OPTR(I,J) = 1.0D0
   40    CONTINUE
   60 CONTINUE
*
      DO 120 IOUT = 1, 2
         IFAIL = 1
         TOUT = TWANT(IOUT)
         CALL D03RAF(NPDE,TS,TOUT,DT,XMIN,XMAX,YMIN,YMAX,NX,NY,TOLS,
     +               TOLT,PDEF1,BNDRY1,PDEIV1,MONIT1,OPTI,OPTR,RWK,
     +               LENRWK,IWK,LENIWK,LWK,LENLWK,ITRACE,IND,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99994) IFAIL
            GO TO 140
         END IF
*
*        Print statistics
*
         WRITE (NOUT,99999) 'Statistics:'
         WRITE (NOUT,99998) 'Time = ', TS
         WRITE (NOUT,99997) 'Total number of accepted timesteps =',
     +     IWK(1)
         WRITE (NOUT,99997) 'Total number of rejected timesteps =',
     +     IWK(2)
         WRITE (NOUT,99996)
         MAXLEV = 3
         DO 80 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) THEN
               WRITE (NOUT,'(I8,4I10)') J, (IWK(J+2+I*MAXLEV),I=0,3)
            END IF
   80    CONTINUE
         WRITE (NOUT,99995)
         DO 100 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) THEN
               WRITE (NOUT,'(I8,2I14)') J, (IWK(J+2+I*MAXLEV),I=4,5)
            END IF
  100    CONTINUE
         WRITE (NOUT,*)
*
  120 CONTINUE
*
  140 CONTINUE
*
      RETURN
*
99999 FORMAT (1X,A)
99998 FORMAT (1X,A,F8.4)
99997 FORMAT (1X,A,I5)
99996 FORMAT (1X,/'             T o t a l  n u m b e r  o f    ',/'   ',
     +       '       Residual  Jacobian    Newton   Lin sys',/'       ',
     +       '      evals     evals     iters     iters',/' At level ')
99995 FORMAT (1X,/'             M a x i m u m  n u m b e r  o f',/'   ',
     +       '           Newton iters    Lin sys iters ',/' At level ')
99994 FORMAT (1X,/1X,' ** D03RAF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE PDEIV1(NPTS,NPDE,T,X,Y,U)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION U(NPTS,NPDE), X(NPTS), Y(NPTS)
*     .. Local Scalars ..
      INTEGER          I
*     .. Executable Statements ..
      DO 20 I = 1, NPTS
         U(I,1) = 1.0D0
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE PDEF1(NPTS,NPDE,T,X,Y,U,UT,UX,UY,UXX,UXY,UYY,RES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION RES(NPTS,NPDE), U(NPTS,NPDE), UT(NPTS,NPDE),
     +                 UX(NPTS,NPDE), UXX(NPTS,NPDE), UXY(NPTS,NPDE),
     +                 UY(NPTS,NPDE), UYY(NPTS,NPDE), X(NPTS), Y(NPTS)
*     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA, D, DELTA, DIFF, REAC
*     .. Local Scalars ..
      INTEGER          I
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Common blocks ..
      COMMON           /PARAM1/ALPHA, DELTA, REAC, DIFF, D
*     .. Executable Statements ..
      DO 20 I = 1, NPTS
         RES(I,1) = UT(I,1) - DIFF*(UXX(I,1)+UYY(I,1)) -
     +              D*(1.0D0+ALPHA-U(I,1))*EXP(-DELTA/U(I,1))
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE BNDRY1(NPTS,NPDE,T,X,Y,U,UT,UX,UY,NBPTS,LBND,RES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NBPTS, NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION RES(NPTS,NPDE), U(NPTS,NPDE), UT(NPTS,NPDE),
     +                 UX(NPTS,NPDE), UY(NPTS,NPDE), X(NPTS), Y(NPTS)
      INTEGER          LBND(NBPTS)
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, J
*     .. External Functions ..
      DOUBLE PRECISION X02AJF
      EXTERNAL         X02AJF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS
*     .. Executable Statements ..
      TOL = 10.D0*X02AJF()
*
      DO 20 I = 1, NBPTS
         J = LBND(I)
         IF (ABS(X(J)).LE.TOL) THEN
            RES(J,1) = UX(J,1)
         ELSE IF (ABS(X(J)-1.0D0).LE.TOL) THEN
            RES(J,1) = U(J,1) - 1.0D0
         ELSE IF (ABS(Y(J)).LE.TOL) THEN
            RES(J,1) = UY(J,1)
         ELSE IF (ABS(Y(J)-1.0D0).LE.TOL) THEN
            RES(J,1) = U(J,1) - 1.0D0
         END IF
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE MONIT1(NPDE,T,DT,DTNEW,TLAST,NLEV,NGPTS,XPTS,YPTS,LSOL,
     +                  SOL,IERR)
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION DT, DTNEW, T
      INTEGER          IERR, NLEV, NPDE
      LOGICAL          TLAST
*     .. Array Arguments ..
      DOUBLE PRECISION SOL(*), XPTS(*), YPTS(*)
      INTEGER          LSOL(NLEV), NGPTS(NLEV)
*     .. Scalars in Common ..
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      INTEGER          I, IPSOL, IPT, LEVEL, NPTS
*     .. Common blocks ..
      COMMON           /OTIME1/TWANT, IOUT
*     .. Executable Statements ..
      IF (TLAST) THEN
*
*        Print solution
*
         IF (IOUT.EQ.2) THEN
            WRITE (NOUT,99999) T
            WRITE (NOUT,99998)
            LEVEL = 1
            NPTS = NGPTS(LEVEL)
            IPSOL = LSOL(LEVEL)
            IPT = 1
            DO 20 I = 1, NPTS, 4
               WRITE (NOUT,99997) XPTS(IPT+I-1), YPTS(IPT+I-1),
     +           SOL(IPSOL+I)
   20       CONTINUE
            WRITE (NOUT,*)
         END IF
      END IF
*
      RETURN
*
99999 FORMAT (1X,'Solution at every 4th grid point in level 1 at time ',
     +       F8.4,':')
99998 FORMAT (1X,/7X,'x',10X,'y',8X,'approx u',/)
99997 FORMAT (3(1X,E11.4))
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          MXLEV, NPDE, NPTS
      PARAMETER        (MXLEV=4,NPDE=2,NPTS=2000)
      INTEGER          LENIWK, LENRWK, LENLWK
      PARAMETER        (LENIWK=2*NPTS*(5*MXLEV+14)+2+7*MXLEV,
     +                 LENRWK=2*NPTS*NPDE*(5*MXLEV+9+18*NPDE)+NPTS*2,
     +                 LENLWK=2*NPTS+400)
*     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA, BETA, PI
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      DOUBLE PRECISION TOLS, TOLT, TOUT, TS, XMAX, XMIN, XX, YMAX, YMIN
      INTEGER          I, IFAIL, IND, ITRACE, J, MAXLEV, NX, NY
*     .. Local Arrays ..
      DOUBLE PRECISION DT(3), OPTR(3,NPDE), RWK(LENRWK)
      INTEGER          IWK(LENIWK), OPTI(4)
      LOGICAL          LWK(LENLWK)
*     .. External Functions ..
      DOUBLE PRECISION X01AAF
      EXTERNAL         X01AAF
*     .. External Subroutines ..
      EXTERNAL         BNDRY2, D03RAF, MONIT2, PDEF2, PDEIV2
*     .. Common blocks ..
      COMMON           /OTIME2/TWANT, IOUT
      COMMON           /PARAM2/ALPHA, BETA, PI
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 2'
      WRITE (NOUT,*)
*
      XX = 0.0D0
      PI = X01AAF(XX)
      ALPHA = 50.0D0
      BETA = 300.0D0
*
      IND = 0
      ITRACE = 0
      TS = 0.0D0
      TWANT(1) = 0.01D0
      TWANT(2) = 0.025D0
      DT(1) = 0.5D-3
      DT(2) = 1.0D-6
      DT(3) = 0.0D0
      XMIN = 0.0D0
      XMAX = 1.0D0
      YMIN = 0.0D0
      YMAX = 1.0D0
      TOLS = 0.075D0
      TOLT = 0.1D0
      NX = 11
      NY = 11
      OPTI(1) = 4
      DO 20 I = 2, 4
         OPTI(I) = 0
   20 CONTINUE
      OPTR(1,1) = 250.0D0
      OPTR(1,2) = 1.5D6
      DO 60 J = 1, NPDE
         DO 40 I = 2, 3
            OPTR(I,J) = 1.0D0
   40    CONTINUE
   60 CONTINUE
*
      DO 120 IOUT = 1, 2
         IFAIL = 1
         TOUT = TWANT(IOUT)
         CALL D03RAF(NPDE,TS,TOUT,DT,XMIN,XMAX,YMIN,YMAX,NX,NY,TOLS,
     +               TOLT,PDEF2,BNDRY2,PDEIV2,MONIT2,OPTI,OPTR,RWK,
     +               LENRWK,IWK,LENIWK,LWK,LENLWK,ITRACE,IND,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,99994) IFAIL
            GO TO 140
         END IF
*
*        Print statistics
*
         MAXLEV = OPTI(1)
         WRITE (NOUT,99999) 'Statistics:'
         WRITE (NOUT,99998) 'Time = ', TS
         WRITE (NOUT,99997) 'Total number of accepted timesteps =',
     +     IWK(1)
         WRITE (NOUT,99997) 'Total number of rejected timesteps =',
     +     IWK(2)
         WRITE (NOUT,99996)
         MAXLEV = OPTI(1)
         DO 80 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) THEN
               WRITE (NOUT,'(I8,4I10)') J, (IWK(J+2+I*MAXLEV),I=0,3)
            END IF
   80    CONTINUE
         WRITE (NOUT,99995)
         DO 100 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) THEN
               WRITE (NOUT,'(I8,2I14)') J, (IWK(J+2+I*MAXLEV),I=4,5)
            END IF
  100    CONTINUE
         WRITE (NOUT,*)
*
  120 CONTINUE
*
  140 CONTINUE
*
      RETURN
*
99999 FORMAT (1X,A)
99998 FORMAT (1X,A,F8.4)
99997 FORMAT (1X,A,I5)
99996 FORMAT (1X,/'             T o t a l  n u m b e r  o f    ',/'   ',
     +       '       Residual  Jacobian    Newton   Lin sys',/'       ',
     +       '      evals     evals     iters     iters',/' At level ')
99995 FORMAT (1X,/'             M a x i m u m  n u m b e r  o f',/'   ',
     +       '           Newton iters    Lin sys iters ',/' At level ')
99994 FORMAT (1X,/1X,' ** D03RAF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE PDEIV2(NPTS,NPDE,T,X,Y,U)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION U(NPTS,NPDE), X(NPTS), Y(NPTS)
*     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA, BETA, PI
*     .. Local Scalars ..
      DOUBLE PRECISION B2, FP
      INTEGER          I
*     .. Intrinsic Functions ..
      INTRINSIC        SIN
*     .. Common blocks ..
      COMMON           /PARAM2/ALPHA, BETA, PI
*     .. Executable Statements ..
      FP = 4.0D0*PI
*
      DO 20 I = 1, NPTS
         B2 = -1.0D0 - ALPHA*X(I)*Y(I) - BETA*SIN(FP*X(I))*SIN(FP*Y(I))
         U(I,1) = 1.0D1 + (16.0D0*X(I)*(1.0D0-X(I))*Y(I)*(1.0D0-Y(I)))
     +            **2
         U(I,2) = B2 + 1.0D4*U(I,1)
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE PDEF2(NPTS,NPDE,T,X,Y,U,UT,UX,UY,UXX,UXY,UYY,RES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION RES(NPTS,NPDE), U(NPTS,NPDE), UT(NPTS,NPDE),
     +                 UX(NPTS,NPDE), UXX(NPTS,NPDE), UXY(NPTS,NPDE),
     +                 UY(NPTS,NPDE), UYY(NPTS,NPDE), X(NPTS), Y(NPTS)
*     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA, BETA, PI
*     .. Local Scalars ..
      DOUBLE PRECISION B1, B2, FP
      INTEGER          I
*     .. Intrinsic Functions ..
      INTRINSIC        SIN
*     .. Common blocks ..
      COMMON           /PARAM2/ALPHA, BETA, PI
*     .. Executable Statements ..
      FP = 4.0D0*PI
*
      DO 20 I = 1, NPTS
         B1 = 1.0D0 + ALPHA*X(I)*Y(I) + BETA*SIN(FP*X(I))*SIN(FP*Y(I))
         B2 = -B1
         RES(I,1) = UT(I,1) - (UXX(I,1)+UYY(I,1)) - U(I,1)*(B1-U(I,1)
     +              -0.5D-6*U(I,2))
         RES(I,2) = -0.05D0*(UXX(I,2)+UYY(I,2)) - U(I,2)
     +              *(B2+1.0D4*U(I,1)-U(I,2))
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE BNDRY2(NPTS,NPDE,T,X,Y,U,UT,UX,UY,NBPTS,LBND,RES)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NBPTS, NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION RES(NPTS,NPDE), U(NPTS,NPDE), UT(NPTS,NPDE),
     +                 UX(NPTS,NPDE), UY(NPTS,NPDE), X(NPTS), Y(NPTS)
      INTEGER          LBND(NBPTS)
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, J
*     .. External Functions ..
      DOUBLE PRECISION X02AJF
      EXTERNAL         X02AJF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS
*     .. Executable Statements ..
      TOL = 10.D0*X02AJF()
*
      DO 20 I = 1, NBPTS
         J = LBND(I)
         IF (ABS(X(J)).LE.TOL .OR. ABS(X(J)-1.0D0).LE.TOL) THEN
            RES(J,1) = UX(J,1)
            RES(J,2) = UX(J,2)
         ELSE IF (ABS(Y(J)).LE.TOL .OR. ABS(Y(J)-1.0D0).LE.TOL) THEN
            RES(J,1) = UY(J,1)
            RES(J,2) = UY(J,2)
         END IF
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE MONIT2(NPDE,T,DT,DTNEW,TLAST,NLEV,NGPTS,XPTS,YPTS,LSOL,
     +                  SOL,IERR)
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION DT, DTNEW, T
      INTEGER          IERR, NLEV, NPDE
      LOGICAL          TLAST
*     .. Array Arguments ..
      DOUBLE PRECISION SOL(*), XPTS(*), YPTS(*)
      INTEGER          LSOL(NLEV), NGPTS(NLEV)
*     .. Scalars in Common ..
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      INTEGER          I, IPSOL, IPT, LEVEL, NPTS
*     .. Common blocks ..
      COMMON           /OTIME2/TWANT, IOUT
*     .. Executable Statements ..
      IF (TLAST) THEN
*
*        Print solution
*
         IF (IOUT.EQ.2) THEN
            WRITE (NOUT,99999) T
            WRITE (NOUT,99998)
            LEVEL = 1
            NPTS = NGPTS(LEVEL)
            IPSOL = LSOL(LEVEL)
            IPT = 1
            DO 20 I = 1, NPTS, 2
               WRITE (NOUT,99997) XPTS(IPT+I-1), YPTS(IPT+I-1),
     +           SOL(IPSOL+I), SOL(IPSOL+NPTS+I)
   20       CONTINUE
            WRITE (NOUT,*)
         END IF
      END IF
*
      RETURN
*
99999 FORMAT (1X,'Solution at every 2nd grid point in level 1 at time ',
     +       F8.4,':')
99998 FORMAT (1X,/7X,'x',10X,'y',9X,'approx c1',3X,'approx c2',/)
99997 FORMAT (2(1X,E11.4),2X,E11.4,2X,E11.4)
      END