!   D02NBF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn, jac
!     .. 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, Public       :: iset = 1, itrace = 0, neq = 3,       &
                                          nin = 5, nout = 6
      Integer, Parameter, Public       :: nrw = 50 + 4*neq
      Integer, Parameter, Public       :: nwkjac = neq*(neq+1)
      Integer, Parameter, Public       :: ldysav = neq
    Contains
      Subroutine fcn(neq,t,y,f,ires)

!       .. 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
      End Subroutine fcn

      Subroutine jac(neq,t,y,h,d,p)

!       .. 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
      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,   &
                             nag_wp, x04abf
      Use d02nbfe_mod, Only: fcn, iset, itrace, jac, ldysav, 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