    PROGRAM d02nmfe

!      D02NMF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : d02nmf, d02nsf, d02nvf, d02nyf, d02xkf, nag_wp, &
                               x04abf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER   :: tstep = 2.0E0_nag_wp
       INTEGER, PARAMETER              :: iset = 1, neq = 3, nin = 5,          &
                                          njcpvt = 1, nout = 6
       INTEGER, PARAMETER              :: nrw = 50 + 4*neq
       INTEGER, PARAMETER              :: nwkjac = neq*(neq+1)
       INTEGER, PARAMETER              :: sdysav = 6
       INTEGER, PARAMETER              :: ldysav = neq
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: h, h0, hlast, hmax, hmin, hnext, hu, &
                                          t, tc, tcrit, tcur, tolsf, tout, xout
       INTEGER                         :: i, ifail, iflag, imon, imxer, inln,  &
                                          ires, irevcm, itask, itol, itrace,   &
                                          lacor1, lacor2, lacor3, lacorb,      &
                                          lsavr1, lsavr2, lsavr3, lsavrb,      &
                                          maxord, maxstp, mxhnil, niter, nje,  &
                                          nq, nqu, nre, nst, outchn
       LOGICAL                         :: petzld
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: atol(:), rtol(:), rwork(:),          &
                                          wkjac(:), y(:), ydot(:), ysav(:,:)
       REAL (KIND=nag_wp)              :: con(6)
       INTEGER                         :: inform(23)
       INTEGER, ALLOCATABLE            :: jacpvt(:)
       LOGICAL, ALLOCATABLE            :: algequ(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          int
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NMF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
!      neq: number of differential equations
       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
          ydot(neq),ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq))


!      Integrate to tout by overshooting tout (itask=1) using B.D.F.
!      formulae with a Newton method. Default values for the array con
!      are used. Employ scalar tolerances and the Jacobian is evaluated
!      internally. On the reverse communication call equivalent to the
!      monitr call in forward communication routines carry out
!      interpolation using D02XKF.

       READ (nin,*) maxord, maxstp, mxhnil
       READ (nin,*) hmin, hmax, h0, tcrit
       READ (nin,*) petzld
       READ (nin,*) t, tout
       READ (nin,*) itol
       READ (nin,*) y(1:neq)
       READ (nin,*) rtol(1), atol(1)


       outchn = nout
       CALL x04abf(iset,outchn)
       itask = 1
       xout = tstep
       con(1:6) = 0.0E0_nag_wp

       ifail = 0
       CALL d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, &
          maxstp,mxhnil,'Average-L2',rwork,ifail)

       ifail = 0
       CALL d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)

       lacorb = 50 + neq
       lacor1 = lacorb + 1
       lacor2 = lacorb + 2
       lacor3 = lacorb + 3
       lsavrb = lacorb + neq
       lsavr1 = lsavrb + 1
       lsavr2 = lsavrb + 2
       lsavr3 = lsavrb + 3
       WRITE (nout,*) '    X          Y(1)           Y(2)           Y(3)'
       WRITE (nout,99999) t, (y(i),i=1,neq)

!      Soft fail and error messages only
       irevcm = 0
       itrace = 0

STEPS: DO
          ifail = 1
          CALL d02nmf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
             ysav,sdysav,wkjac,nwkjac,jacpvt,njcpvt,imon,inln,ires,irevcm, &
             itask,itrace,ifail)

          SELECT CASE (irevcm)
          CASE (0)
             EXIT STEPS
          CASE (1,3)
!            Equivalent to fcn evaluation in forward communication
!            routines
             rwork(lsavr1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
             rwork(lsavr2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - &
                3.0E7_nag_wp*y(2)*y(2)
             rwork(lsavr3) = 3.0E7_nag_wp*y(2)*y(2)
          CASE (4)
!            Equivalent to fcn evaluation in forward communication
!            routines
             rwork(lacor1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
             rwork(lacor2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - &
                3.0E7_nag_wp*y(2)*y(2)
             rwork(lacor3) = 3.0E7_nag_wp*y(2)*y(2)
          CASE (5)
!            Equivalent to fcn evaluation in forward communication
!            routines
             ydot(1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
             ydot(2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - &
                3.0E7_nag_wp*y(2)*y(2)
             ydot(3) = 3.0E7_nag_wp*y(2)*y(2)
          CASE (9)
!            Equivalent to monitr call in forward communication routines
             IF (imon==1) THEN
                tc = rwork(19)
                hlast = rwork(15)
                hnext = rwork(16)
                nqu = int(rwork(10))
INTERP:         DO
                   IF (tc-hlast<xout .AND. xout<=tc) THEN
                      iflag = 1

                      CALL d02xkf(xout,rwork(lsavr1),neq,ysav,ldysav,sdysav, &
                         rwork(lacor1),neq,tc,nqu,hlast,hnext,iflag)

                      IF (iflag/=0) THEN
                         imon = -2
                         EXIT INTERP
                      ELSE
                         WRITE (nout,99999) xout, (rwork(lsavrb+i),i=1,neq)
                         xout = xout + tstep
                         IF (xout>tout) EXIT INTERP
                      END IF
                   ELSE
                      EXIT INTERP
                   END IF
                END DO INTERP
             END IF
          CASE (2,6:8)
             WRITE (nout,*)
             WRITE (nout,99994) 'Illegal value of IREVCM = ', irevcm
             EXIT STEPS
          END SELECT
       END DO STEPS
       IF (ifail==0) THEN

          iflag = 0
          CALL d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, &
             imxer,algequ,inform,iflag)

          WRITE (nout,*)
          WRITE (nout,99997) hu, h, tcur
          WRITE (nout,99996) nst, nre, nje
          WRITE (nout,99995) nqu, nq, niter
          WRITE (nout,99994) ' Max err comp = ', imxer
          WRITE (nout,*)
       ELSE
          WRITE (nout,*)
          WRITE (nout,99998) 'Exit D02NMF with IFAIL = ', ifail, '  and T = ', &
             t
       END IF

99999  FORMAT (1X,F8.3,3(F13.5,2X))
99998  FORMAT (1X,A,I2,A,E12.5)
99997  FORMAT (1X,' HUSED = ',E12.5,'  HNEXT = ',E12.5,'  TCUR = ',E12.5)
99996  FORMAT (1X,' NST = ',I6,'    NRE = ',I6,'    NJE = ',I6)
99995  FORMAT (1X,' NQU = ',I6,'    NQ  = ',I6,'  NITER = ',I6)
99994  FORMAT (1X,A,I4)
    END PROGRAM d02nmfe
