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

    MODULE d02ncfe_mod

!      D02NCF 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
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
       INTEGER, PARAMETER                  :: iset = 1, itrace = 0, ml = 1,    &
                                              mu = 2, neq = 3, nin = 5
       INTEGER, PARAMETER                  :: njcpvt = neq
       INTEGER, PARAMETER                  :: nout = 6
       INTEGER, PARAMETER                  :: nrw = 50 + 4*neq
       INTEGER, PARAMETER                  :: nwkjac = neq*(2*ml+mu+1)
       INTEGER, PARAMETER                  :: ldysav = neq
    CONTAINS
!dec$ attributes offload : mic :: fcn
       SUBROUTINE fcn(neq,t,y,f,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)    :: f(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq)
!         .. Executable Statements ..
          f(1) = -alpha*y(1) + beta*y(2)*y(3)
          f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
          f(3) = gamma*y(2)*y(2)
          RETURN
!dec$ attributes offload : mic :: fcn
       END SUBROUTINE fcn

!dec$ attributes offload : mic :: jac
       SUBROUTINE jac(neq,t,y,h,d,ml,mu,p)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: d, h, t
          INTEGER, INTENT (IN)                :: ml, mu, neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: p(ml+mu+1,neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: hxd
!         .. Executable Statements ..
          hxd = h*d

          p(1,1) = one - hxd*(-alpha)
          p(2,1) = -hxd*(beta*y(3))
          p(3,1) = -hxd*(beta*y(2))

          p(1,2) = -hxd*(alpha)
          p(2,2) = one - hxd*(-beta*y(3)-two*gamma*y(2))
          p(3,2) = -hxd*(-beta*y(2))

          p(1,3) = -hxd*(two*gamma*y(2))
          p(2,3) = one - hxd*(0.0_nag_wp)
          RETURN
!dec$ attributes offload : mic :: jac
       END SUBROUTINE jac
    END MODULE d02ncfe_mod

    PROGRAM d02ncfe

!      D02NCF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02nby, d02ncf, d02ntf, d02nwf, d02nyf, d02nzf, &
                               nag_wp, x04abf
       USE d02ncfe_mod, ONLY : fcn, iset, itrace, jac, ldysav, ml, mu, neq,    &
                               nin, njcpvt, nout, nrw, nwkjac
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, h0, hmax, hmin, hu, t, tcrit, &
                                              tcur, tolsf, tout
       INTEGER                             :: i, ifail, imxer, itask, itol,    &
                                              maxord, maxstp, mxhnil, niter,   &
                                              nje, nq, nqu, nre, nst, outchn,  &
                                              sdysav
!      .. 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(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NCF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) maxord, maxstp, mxhnil
       sdysav = maxord + 3
       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
          ydot(neq),ysav(ldysav,sdysav),jacpvt(njcpvt),algequ(neq))
       outchn = nout
       CALL x04abf(iset,outchn)

!      Integrate to tout (itask=1 i.e. overshooting and internal interpolation)
!      using the blend method. Default values for the array con are used.
!      Employ scalar relative tolerance and vector absolute tolerance.
!      The Jacobian is evaluated by jac.
!      monitr subroutine replaced by NAG dummy routine D02NBY.

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

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

       ifail = 0
       CALL d02ntf(neq,neq,'Analytical',ml,mu,nwkjac,njcpvt,rwork,ifail)

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

!      Soft fail and error messages only

       ifail = -1
       CALL d02ncf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, &
          ysav,sdysav,jac,wkjac,nwkjac,jacpvt,njcpvt,d02nby,itask,itrace, &
          ifail)

       IF (ifail==0) THEN
          WRITE (nout,99999) t, (y(i),i=1,neq)
       ELSE
          WRITE (nout,*)
          WRITE (nout,99998) 'Exit D02NCF with IFAIL = ', ifail, '  and T = ', &
             t
       END IF
!      Reset tout and call D02NZF to override internal choice for stepsize.
!      No changes to other parameters.
       h = 0.7_nag_wp

       ifail = 0
       CALL d02nzf(neq,tcrit,h,hmin,hmax,maxstp,mxhnil,rwork,ifail)

       tout = 10.0_nag_wp

       ifail = -1
       CALL d02ncf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, &
          ysav,sdysav,jac,wkjac,nwkjac,jacpvt,njcpvt,d02nby,itask,itrace, &
          ifail)

       IF (ifail==0) THEN
          WRITE (nout,99999) t, (y(i),i=1,neq)

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

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