*     D03RBF Example Program Text
*     Mark 19 Revised. NAG Copyright 1999.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
      INTEGER          MXLEV, NPDE, NPTS
      PARAMETER        (MXLEV=5,NPDE=2,NPTS=3000)
      INTEGER          LENIWK, LENRWK, LENLWK
      PARAMETER        (LENIWK=NPTS*(5*MXLEV+14)+2+7*MXLEV,
     +                 LENRWK=NPTS*NPDE*(5*MXLEV+9+18*NPDE)+2*NPTS,
     +                 LENLWK=2*NPTS)
*     .. Scalars in Common ..
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      DOUBLE PRECISION TOLS, TOLT, TOUT, TS
      INTEGER          I, IFAIL, IND, ITRACE, J, MAXLEV
*     .. Local Arrays ..
      DOUBLE PRECISION DT(3), OPTR(3,NPDE), RWK(LENRWK)
      INTEGER          IWK(LENIWK), OPTI(4)
      LOGICAL          LWK(LENLWK)
*     .. External Subroutines ..
      EXTERNAL         BNDARY, D03RBF, INIDOM, MONITR, PDEDEF, PDEIV
*     .. Common blocks ..
      COMMON           /OTIME/TWANT, IOUT
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D03RBF Example Program Results'
*
      IND = 0
      ITRACE = 0
      TS = 0.0D0
      TWANT(1) = 0.25D0
      TWANT(2) = 1.0D0
      DT(1) = 0.001D0
      DT(2) = 1.0D-7
      DT(3) = 0.0D0
      TOLS = 0.1D0
      TOLT = 0.05D0
      OPTI(1) = 5
      MAXLEV = OPTI(1)
      DO 20 I = 2, 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
*
*     Call main routine
*
      DO 120 IOUT = 1, 2
         IFAIL = 1
         TOUT = TWANT(IOUT)
         CALL D03RBF(NPDE,TS,TOUT,DT,TOLS,TOLT,INIDOM,PDEDEF,BNDARY,
     +               PDEIV,MONITR,OPTI,OPTR,RWK,LENRWK,IWK,LENIWK,LWK,
     +               LENLWK,ITRACE,IND,IFAIL)
         IF (IFAIL.EQ.0) THEN
*
*           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 = OPTI(1)
            DO 80 J = 1, MAXLEV
               IF (IWK(J+2).NE.0) THEN
                  WRITE (NOUT,'(I6,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,'(I6,2I14)') J, (IWK(J+2+I*MAXLEV),I=4,5)
               END IF
  100       CONTINUE
            WRITE (NOUT,*)
         ELSE
            WRITE (NOUT,99994) IFAIL
            GO TO 140
         END IF
*
  120 CONTINUE
*
  140 CONTINUE
*
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,' ** D03RBF returned with IFAIL = ',I5)
      END
*
      SUBROUTINE INIDOM(MAXPTS,XMIN,XMAX,YMIN,YMAX,NX,NY,NPTS,NROWS,
     +                  NBNDS,NBPTS,LROW,IROW,ICOL,LLBND,ILBND,LBND,
     +                  IERR)
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION XMAX, XMIN, YMAX, YMIN
      INTEGER          IERR, MAXPTS, NBNDS, NBPTS, NPTS, NROWS, NX, NY
*     .. Array Arguments ..
      INTEGER          ICOL(*), ILBND(*), IROW(*), LBND(*), LLBND(*),
     +                 LROW(*)
*     .. Local Scalars ..
      INTEGER          I, IFAIL, J, LENIWK
*     .. Local Arrays ..
      INTEGER          ICOLD(105), ILBNDD(28), IROWD(11), IWK(122),
     +                 LBNDD(72), LLBNDD(28), LROWD(11)
      CHARACTER*33     PGRID(11)
*     .. External Subroutines ..
      EXTERNAL         D03RYF
*     .. Data statements ..
      DATA             ICOLD/0, 1, 2, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
     +                 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4,
     +                 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 8, 9, 10, 0,
     +                 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5,
     +                 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
     +                 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6,
     +                 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8/
      DATA             ILBNDD/1, 2, 3, 4, 1, 4, 1, 2, 3, 4, 3, 4, 1, 2,
     +                 12, 23, 34, 41, 14, 41, 12, 23, 34, 41, 43, 14,
     +                 21, 32/
      DATA             IROWD/0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10/
      DATA             LBNDD/2, 4, 15, 26, 37, 46, 57, 68, 79, 88, 98,
     +                 99, 100, 101, 102, 103, 104, 96, 86, 85, 84, 83,
     +                 82, 70, 59, 48, 39, 28, 17, 6, 8, 9, 10, 11, 12,
     +                 13, 18, 29, 40, 49, 60, 72, 73, 74, 75, 76, 77,
     +                 67, 56, 45, 36, 25, 33, 32, 42, 52, 53, 43, 1,
     +                 97, 105, 87, 81, 3, 7, 71, 78, 14, 31, 51, 54,
     +                 34/
      DATA             LLBNDD/1, 2, 11, 18, 19, 24, 31, 37, 42, 48, 53,
     +                 55, 56, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
     +                 68, 69, 70, 71, 72/
      DATA             LROWD/1, 4, 15, 26, 37, 46, 57, 68, 79, 88, 97/
*     .. Executable Statements ..
      NX = 11
      NY = 11
*
*     Check MAXPTS against rough estimate of NPTS
*
      NPTS = NX*NY
      IF (MAXPTS.LT.NPTS) THEN
         IERR = -1
         RETURN
      END IF
*
      XMIN = 0.0D0
      YMIN = 0.0D0
      XMAX = 1.0D0
      YMAX = 1.0D0
*
      NROWS = 11
      NPTS = 105
      NBNDS = 28
      NBPTS = 72
*
      DO 20 I = 1, NROWS
         LROW(I) = LROWD(I)
         IROW(I) = IROWD(I)
   20 CONTINUE
*
      DO 40 I = 1, NBNDS
         LLBND(I) = LLBNDD(I)
         ILBND(I) = ILBNDD(I)
   40 CONTINUE
*
      DO 60 I = 1, NBPTS
         LBND(I) = LBNDD(I)
   60 CONTINUE
*
      DO 80 I = 1, NPTS
         ICOL(I) = ICOLD(I)
   80 CONTINUE
*
      WRITE (NOUT,*) 'Base grid:'
      WRITE (NOUT,*)
      LENIWK = 122
      IFAIL = -1
*
      CALL D03RYF(NX,NY,NPTS,NROWS,NBNDS,NBPTS,LROW,IROW,ICOL,LLBND,
     +            ILBND,LBND,IWK,LENIWK,PGRID,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,*) ' '
         DO 100 J = 1, NY
            WRITE (NOUT,*) PGRID(J)
            WRITE (NOUT,*) ' '
  100    CONTINUE
         WRITE (NOUT,*) ' '
      END IF
*
      RETURN
      END
*
      SUBROUTINE PDEIV(NPTS,NPDE,T,X,Y,U)
*     .. Parameters ..
      DOUBLE PRECISION EPS
      PARAMETER        (EPS=1D-3)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NPDE, NPTS
*     .. Array Arguments ..
      DOUBLE PRECISION U(NPTS,NPDE), X(NPTS), Y(NPTS)
*     .. Local Scalars ..
      DOUBLE PRECISION A
      INTEGER          I
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Executable Statements ..
      DO 20 I = 1, NPTS
         A = (-4.0D0*X(I)+4.0D0*Y(I)-T)/(32.0D0*EPS)
         IF (A.LE.0.0D0) THEN
            U(I,1) = 0.75D0 - 0.25D0/(1.0D0+EXP(A))
            U(I,2) = 0.75D0 + 0.25D0/(1.0D0+EXP(A))
         ELSE
            U(I,1) = 0.75D0 - 0.25D0*EXP(-A)/(EXP(-A)+1.0D0)
            U(I,2) = 0.75D0 + 0.25D0*EXP(-A)/(EXP(-A)+1.0D0)
         END IF
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE PDEDEF(NPTS,NPDE,T,X,Y,U,UT,UX,UY,UXX,UXY,UYY,RES)
*     .. Parameters ..
      DOUBLE PRECISION EPS
      PARAMETER        (EPS=1D-3)
*     .. 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)
*     .. Local Scalars ..
      INTEGER          I
*     .. Executable Statements ..
      DO 20 I = 1, NPTS
         RES(I,1) = UT(I,1) - (-U(I,1)*UX(I,1)-U(I,2)*UY(I,1)
     +              +EPS*(UXX(I,1)+UYY(I,1)))
         RES(I,2) = UT(I,2) - (-U(I,1)*UX(I,2)-U(I,2)*UY(I,2)
     +              +EPS*(UXX(I,2)+UYY(I,2)))
   20 CONTINUE
*
      RETURN
      END
      SUBROUTINE BNDARY(NPTS,NPDE,T,X,Y,U,UT,UX,UY,NBNDS,NBPTS,LLBND,
     +                  ILBND,LBND,RES)
*     .. Parameters ..
      DOUBLE PRECISION EPS
      PARAMETER        (EPS=1D-3)
*     .. Scalar Arguments ..
      DOUBLE PRECISION T
      INTEGER          NBNDS, 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          ILBND(NBNDS), LBND(NBPTS), LLBND(NBNDS)
*     .. Local Scalars ..
      DOUBLE PRECISION A
      INTEGER          I, K
*     .. Intrinsic Functions ..
      INTRINSIC        EXP
*     .. Executable Statements ..
      DO 20 K = LLBND(1), NBPTS
         I = LBND(K)
         A = (-4.0D0*X(I)+4.0D0*Y(I)-T)/(32.0D0*EPS)
         IF (A.LE.0.0D0) THEN
            RES(I,1) = U(I,1) - (0.75D0-0.25D0/(1.0D0+EXP(A)))
            RES(I,2) = U(I,2) - (0.75D0+0.25D0/(1.0D0+EXP(A)))
         ELSE
            RES(I,1) = U(I,1) - (0.75D0-0.25D0*EXP(-A)/(EXP(-A)+1.0D0))
            RES(I,2) = U(I,2) - (0.75D0+0.25D0*EXP(-A)/(EXP(-A)+1.0D0))
         END IF
   20 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE MONITR(NPDE,T,DT,DTNEW,TLAST,NLEV,XMIN,YMIN,DXB,DYB,
     +                  LGRID,ISTRUC,LSOL,SOL,IERR)
*     .. Parameters ..
      INTEGER          MAXPTS, NOUT
      PARAMETER        (MAXPTS=2500,NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION DT, DTNEW, DXB, DYB, T, XMIN, YMIN
      INTEGER          IERR, NLEV, NPDE
      LOGICAL          TLAST
*     .. Array Arguments ..
      DOUBLE PRECISION SOL(*)
      INTEGER          ISTRUC(*), LGRID(*), LSOL(NLEV)
*     .. Scalars in Common ..
      INTEGER          IOUT
*     .. Arrays in Common ..
      DOUBLE PRECISION TWANT(2)
*     .. Local Scalars ..
      INTEGER          IFAIL, IPSOL, IPT, LEVEL, NPTS
*     .. Local Arrays ..
      DOUBLE PRECISION UEX(105,2), X(MAXPTS), Y(MAXPTS)
*     .. External Subroutines ..
      EXTERNAL         D03RZF, PDEIV
*     .. Common blocks ..
      COMMON           /OTIME/TWANT, IOUT
*     .. Executable Statements ..
      IFAIL = -1
      IF (TLAST) THEN
         DO 40 LEVEL = 1, NLEV
            IPSOL = LSOL(LEVEL)
*
*           Get grid information
*
            CALL D03RZF(LEVEL,NLEV,XMIN,YMIN,DXB,DYB,LGRID,ISTRUC,NPTS,
     +                  X,Y,MAXPTS,IFAIL)
            IF (IFAIL.NE.0) THEN
               IERR = 1
               RETURN
            END IF
*
            IF (IOUT.EQ.2 .AND. LEVEL.EQ.1) THEN
*
*              Get exact solution
*
               CALL PDEIV(NPTS,NPDE,T,X,Y,UEX)
               WRITE (NOUT,*)
               WRITE (NOUT,
     +'('' Solution at every 2nd grid point '', ''in level 1 at time '',
     + F8.4,'':'')') T
               WRITE (NOUT,*)
               WRITE (NOUT,'(7X,A,10X,A,8X,A,5X,A,4X,A,4X,A)') 'x', 'y',
     +           'approx u', 'exact u', 'approx v', 'exact v'
               WRITE (NOUT,*)
               IPSOL = LSOL(LEVEL)
               DO 20 IPT = 1, NPTS, 2
                  WRITE (NOUT,'(6(1X,E11.4))') X(IPT), Y(IPT),
     +              SOL(IPSOL+IPT), UEX(IPT,1), SOL(IPSOL+NPTS+IPT),
     +              UEX(IPT,2)
   20          CONTINUE
               WRITE (NOUT,*)
            END IF
*
   40    CONTINUE
      END IF
*
      RETURN
      END