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

    MODULE d02nbfe_mod

!      D02NBF 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, 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 :: 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,p)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: d, h, t
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: p(neq,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(1,2) = -hxd*(beta*y(3))
          p(1,3) = -hxd*(beta*y(2))
          p(2,1) = -hxd*(alpha)
          p(2,2) = one - hxd*(-beta*y(3)-two*gamma*y(2))
          p(2,3) = -hxd*(-beta*y(2))
!         Do not need to set P(3,1) since Jacobian preset to zero
!         P(3,1) =       - HXD*(0.0E0)
          p(3,2) = -hxd*(two*gamma*y(2))
          p(3,3) = one - hxd*(0.0_nag_wp)
          RETURN
!dec$ attributes offload : mic :: jac
       END SUBROUTINE jac
    END MODULE d02nbfe_mod

    PROGRAM d02nbfe

!      D02NBF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02nbf, d02nby, d02nbz, d02nsf, d02nvf, d02nyf, &
                               x04abf
       USE d02nbfe_mod, ONLY : fcn, iset, itrace, jac, ldysav, nag_wp, neq,    &
                               nin, nout, nrw, nwkjac
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, h0, hmax, hmin, hu, t, tcrit, &
                                              tcur, tinit, tolsf, tout
       INTEGER                             :: i, icase, ifail, imxer, 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(:),      &
                                              wkjac(:), y(:), ydot(:),         &
                                              yinit(:), ysav(:,:)
       REAL (KIND=nag_wp)                  :: con(6)
       INTEGER                             :: inform(23)
       LOGICAL, ALLOCATABLE                :: algequ(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02NBF 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),wkjac(nwkjac),y(neq), &
          yinit(neq),ydot(neq),ysav(ldysav,sdysav),algequ(neq))
       outchn = nout
       CALL x04abf(iset,outchn)


       READ (nin,*) petzld
       READ (nin,*) hmin, hmax, h0
       READ (nin,*) tinit, tout
       READ (nin,*) itol
       READ (nin,*) yinit(1:neq)
       READ (nin,*) rtol(1), atol(1)

!      Two cases. In both cases:
!                 integrate to tout without passing tout;
!                 use B.D.F formulae with a Newton method;
!                 use default values for the array con; 
!                 use scalar tolerances;
!                 use NAG dummy routine D02NBY in place of MONITR subroutine.
       con(1:6) = 0.0_nag_wp
       tcrit = tout
       itask = 4

CASES: DO icase = 1, 2

          t = tinit
          y(1:neq) = yinit(1:neq)


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

          WRITE (nout,*)
          ifail = 0
          SELECT CASE (icase)
          CASE (1)
!            First case. The Jacobian is evaluated internally.
             CALL d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)

             WRITE (nout,*) '  Numerical Jacobian'
          CASE (2)
!            Second case. The Jacobian is evaluated by jac.
             CALL d02nsf(neq,neq,'Analytical',nwkjac,rwork,ifail)

             WRITE (nout,*) '  Analytic Jacobian'
          END SELECT

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

!         Soft fail and error messages only

          ifail = -1
          SELECT CASE (icase)
          CASE (1)
             CALL d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
                fcn,ysav,sdysav,d02nbz,wkjac,nwkjac,d02nby,itask,itrace,ifail)
          CASE (2)
             CALL d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
                fcn,ysav,sdysav,jac,wkjac,nwkjac,d02nby,itask,itrace,ifail)
          END SELECT

          IF (ifail==0) THEN
             WRITE (nout,99999) t, y(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 D02NBF with IFAIL = ', ifail, &
                '  and T = ', t
          END IF
       END DO CASES


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 d02nbfe
