*     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
      STOP
      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
*     .. Save statement ..
      SAVE             /OTIME1/, /PARAM1/
*     .. 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)
*
*        Print statistics
*
         WRITE (NOUT,'('' Statistics:'')')
         WRITE (NOUT,'('' Time = '',F8.4)') TS
         WRITE (NOUT,
     +     '('' Total number of accepted timesteps ='', I5)') IWK(1)
         WRITE (NOUT,
     +     '('' Total number of rejected timesteps ='', I5)') IWK(2)
         WRITE (NOUT,*)
         WRITE (NOUT,
     +     '(''             T o t a l  n u m b e r  o f    '')')
         WRITE (NOUT,
     + '(''          Residual  Jacobian    Newton ''   , ''  Lin sys'')'
     +     )
         WRITE (NOUT,
     + '(''             evals     evals     iters ''   , ''    iters'')'
     +     )
         WRITE (NOUT,'('' At level '')')
         MAXLEV = 3
         DO 80 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) WRITE (NOUT,'(I8,4I10)') J, IWK(J+2),
     +          IWK(J+2+MAXLEV), IWK(J+2+2*MAXLEV), IWK(J+2+3*MAXLEV)
*
   80    CONTINUE
         WRITE (NOUT,*)
         WRITE (NOUT,
     +     '(''             M a x i m u m  n u m b e r '', '' o f'')')
         WRITE (NOUT,
     +     '(''              Newton iters    Lin sys iters '')')
         WRITE (NOUT,'('' At level '')')
         DO 100 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) WRITE (NOUT,'(I8,2I14)') J,
     +          IWK(J+2+4*MAXLEV), IWK(J+2+5*MAXLEV)
  100    CONTINUE
         WRITE (NOUT,*)
*
  120 CONTINUE
*
      RETURN
      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
*     .. Save statement ..
      SAVE             /PARAM1/
*     .. 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
*     .. Save statement ..
      SAVE              /OTIME1/
*     .. Executable Statements ..
*
      IF (TLAST) THEN
*
*        Print solution
*
         IF (IOUT.EQ.2) THEN
            WRITE (NOUT,
     +'('' Solution at every 4th grid point '', ''in level 1 at time '',
     + F8.4,'':'')') T
            WRITE (NOUT,*)
            WRITE (NOUT,'(7X,''x'',10X,''y'',8X,''approx u'')')
            WRITE (NOUT,*)
            LEVEL = 1
            NPTS = NGPTS(LEVEL)
            IPSOL = LSOL(LEVEL)
            IPT = 1
            DO 20 I = 1, NPTS, 4
               WRITE (NOUT,'(3(1X,D11.4))') XPTS(IPT+I-1),
     +           YPTS(IPT+I-1), SOL(IPSOL+I)
   20       CONTINUE
            WRITE (NOUT,*)
         END IF
      END IF
*
      RETURN
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          MXLEV, NPDE, NPTS
      PARAMETER        (MXLEV=4,NPDE=2,NPTS=8000)
      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, 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
*     .. Save statement ..
      SAVE             /OTIME2/, /PARAM2/
*     .. 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)
*
*        Print statistics
*
         MAXLEV = OPTI(1)
         WRITE (NOUT,'('' Statistics:'')')
         WRITE (NOUT,'('' Time = '',F8.4)') TS
         WRITE (NOUT,
     +     '('' Total number of accepted timesteps ='', I5)') IWK(1)
         WRITE (NOUT,
     +     '('' Total number of rejected timesteps ='', I5)') IWK(2)
         WRITE (NOUT,*)
         WRITE (NOUT,
     +     '(''             T o t a l  n u m b e r  o f    '')')
         WRITE (NOUT,
     + '(''          Residual  Jacobian    Newton ''   , ''  Lin sys'')'
     +     )
         WRITE (NOUT,
     + '(''             evals     evals     iters ''   , ''    iters'')'
     +     )
         WRITE (NOUT,'('' At level '')')
         MAXLEV = OPTI(1)
         DO 80 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) WRITE (NOUT,'(I6,4I10)') J, IWK(J+2),
     +          IWK(J+2+MAXLEV), IWK(J+2+2*MAXLEV), IWK(J+2+3*MAXLEV)
*
   80    CONTINUE
         WRITE (NOUT,*)
         WRITE (NOUT,
     +     '(''             M a x i m u m  n u m b e r '', '' o f'')')
         WRITE (NOUT,
     +     '(''              Newton iters    Lin sys iters '')')
         WRITE (NOUT,'('' At level '')')
         DO 100 J = 1, MAXLEV
            IF (IWK(J+2).NE.0) WRITE (NOUT,'(I6,2I14)') J,
     +          IWK(J+2+4*MAXLEV), IWK(J+2+5*MAXLEV)
  100    CONTINUE
         WRITE (NOUT,*)
*
  120 CONTINUE
*
      RETURN
      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
*     .. Save statement ..
      SAVE              /PARAM2/
*     .. 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
*     .. Save statement ..
      SAVE             /PARAM2/
*     .. 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
*     .. Save statement ..
      SAVE              /OTIME2/
*     .. Executable Statements ..
*
      IF (TLAST) THEN
*
*        Print solution
*
         IF (IOUT.EQ.2) THEN
            WRITE (NOUT,
     +'('' Solution at every 2nd grid point '', ''in level 1 at time '',
     + F8.4,'':'')') T
            WRITE (NOUT,*)
            WRITE (NOUT,
     +        '(7X,''x'',10X,''y'',9X,''approx c1'',3X,''approx c2'')')
            WRITE (NOUT,*)
            LEVEL = 1
            NPTS = NGPTS(LEVEL)
            IPSOL = LSOL(LEVEL)
            IPT = 1
            DO 20 I = 1, NPTS, 2
               WRITE (NOUT,'(2(1X,D11.4),2X,D11.4,2X,D11.4)')
     +           XPTS(IPT+I-1), YPTS(IPT+I-1), SOL(IPSOL+I),
     +           SOL(IPSOL+NPTS+I)
   20       CONTINUE
            WRITE (NOUT,*)
         END IF
      END IF
*
      RETURN
      END