!   D02NGF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02ngfe_mod

!      D02NGF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: alpha = 0.04_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: beta = 1.0E4_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: gamma = 3.0E7_nag_wp
       INTEGER, PARAMETER                  :: iset = 1, itrace = 0, neq = 3,   &
                                              nin = 5, nout = 6
       INTEGER, PARAMETER                  :: nrw = 50 + 4*neq
       INTEGER, PARAMETER                  :: nwkjac = neq*(neq+1)
       INTEGER, PARAMETER                  :: ldysav = neq
    CONTAINS
!dec$ attributes offload : mic :: resid
       SUBROUTINE resid(neq,t,y,ydot,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: r(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq), ydot(neq)
!         .. Executable Statements ..
          r(1) = -ydot(1)
          r(2) = -ydot(2)
          r(3) = -ydot(3)

          IF (ires==1) THEN
             r(1) = -alpha*y(1) + beta*y(2)*y(3) + r(1)
             r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2)
             r(3) = gamma*y(2)*y(2) + r(3)
          END IF
          RETURN
!dec$ attributes offload : mic :: resid
       END SUBROUTINE resid
    END MODULE d02ngfe_mod

    PROGRAM d02ngfe

!      D02NGF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02nby, d02ngf, d02ngz, d02nsf, d02nvf, d02nyf, &
                               d02xjf, nag_wp, x04abf
       USE d02ngfe_mod, ONLY : iset, itrace, ldysav, neq, nin, nout, nrw,      &
                               nwkjac, resid
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, h0, hmax, hmin, hu, t, tcrit, &
                                              tcur, tolsf, tout, xout
       INTEGER                             :: i, ifail, imxer, iout, itask,    &
                                              itol, maxord, maxstp, mxhnil,    &
                                              niter, nje, nq, nqu, nre, nst,   &
                                              outchn, sdysav
       LOGICAL                             :: petzld
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), rtol(:), rwork(:),      &
                                              sol(:), wkjac(:), y(:), ydot(:), &
                                              ysav(:,:)
       REAL (KIND=nag_wp)                  :: con(6)
       INTEGER                             :: inform(23)
       LOGICAL, ALLOCATABLE                :: algequ(:)
       LOGICAL                             :: lderiv(2)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NGF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) maxord, maxstp, mxhnil
       sdysav = maxord + 1
       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),sol(neq),wkjac(nwkjac),y(neq), &
          ydot(neq),ysav(ldysav,sdysav),algequ(neq))
       outchn = nout
       CALL x04abf(iset,outchn)

!      Integrate to tout by overshooting tout in one step mode (itask=2)
!      using B.D.F formulae with a functional iteration method.
!      Default values for the array con are used. Employ vector
!      tolerances and the Jacobian is evaluated internally, if necessary.
!      monitr subroutine replaced by NAG dummy routine D02NBY.
!      Interpolation outside D02NGF using D02XJF.

       itask = 2
       READ (nin,*) petzld
       READ (nin,*) hmin, hmax, h0, tcrit
       READ (nin,*) t, tout
       READ (nin,*) y(1:neq)
       READ (nin,*) lderiv(1:2)
       READ (nin,*) itol
       READ (nin,*) rtol(1:neq)
       READ (nin,*) atol(1:neq)
       con(1:6) = 0.0_nag_wp

!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL d02nvf(neq,sdysav,maxord,'Functional-iteration',petzld,con,tcrit, &
          hmin,hmax,h0,maxstp,mxhnil,'Average-L2',rwork,ifail)

!      Linear algebra setup required (in case functional iteration
!      encounters any difficulty).
       ifail = 0
       CALL d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)

       xout = 0.02E0_nag_wp
       iout = 1

       WRITE (nout,99993) (i,i=1,neq)
       WRITE (nout,99999) t, (y(i),i=1,neq)

!      Soft fail and error messages only

LOOP1: DO
          ifail = -1
          CALL d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
             resid,ysav,sdysav,d02ngz,wkjac,nwkjac,d02nby,lderiv,itask,itrace, &
             ifail)

          IF (ifail==0) THEN

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

LOOP2:       DO
                IF (tcur-hu<xout .AND. xout<=tcur) THEN

!                  C0 interpolation
                   ifail = 0
                   CALL d02xjf(xout,sol,neq,ysav,ldysav,sdysav,neq,tcur,nqu, &
                      hu,h,ifail)

                   WRITE (nout,99999) xout, (sol(i),i=1,neq)
                   iout = iout + 1
                   xout = real(iout,kind=nag_wp)*0.02E0_nag_wp
                   IF (iout>=6) THEN
                      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
                      EXIT LOOP1
                   END IF
                ELSE
                   EXIT LOOP2
                END IF
             END DO LOOP2
          ELSE
             WRITE (nout,*)
             WRITE (nout,99998) 'Exit D02NGF with IFAIL = ', ifail, &
                '  and T = ', t
          END IF
       END DO LOOP1

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)
99993  FORMAT (/1X,'    X ',3('         Y(',I1,')  '))
    END PROGRAM d02ngfe
