!   D02ZAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02zafe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      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
      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) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3)
        f(2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - &
          3.0E7_nag_wp*y(2)*y(2)
        f(3) = 3.0E7_nag_wp*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) = 1.0E0_nag_wp - hxd*(-0.04E0_nag_wp)
        p(1,2) = -hxd*(1.0E4_nag_wp*y(3))
        p(1,3) = -hxd*(1.0E4_nag_wp*y(2))
        p(2,1) = -hxd*(0.04E0_nag_wp)
        p(2,2) = 1.0E0_nag_wp - hxd*(-1.0E4_nag_wp*y(3)-6.0E7_nag_wp*y(2))
        p(2,3) = -hxd*(-1.0E4_nag_wp*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*(6.0E7_nag_wp*y(2))
        p(3,3) = 1.0E0_nag_wp - hxd*(0.0E0_nag_wp)
        Return
      End Subroutine jac
      Subroutine monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon,inln, &
        hmin,hmax,nqu)

!       .. Use Statements ..
        Use nag_library, Only: d02zaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: hlast, t
        Real (Kind=nag_wp), Intent (Inout)   :: hmax, hmin, hnext
        Integer, Intent (Inout)              :: imon
        Integer, Intent (Out)                :: inln
        Integer, Intent (In)                 :: ldysav, neq, nqu
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: acor(neq,2), r(neq),           &
                                                ydot(neq), ysav(ldysav,*)
        Real (Kind=nag_wp), Intent (Inout)   :: y(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: errloc
        Integer                              :: i, ifail
!       .. Executable Statements ..
        inln = 3
        If (imon==1) Then

          ifail = -1
          errloc = d02zaf(neq,acor(1,2),acor(1,1),ifail)

          If (ifail/=0) Then
            imon = -2
          Else If (errloc>5.0E0_nag_wp) Then
            Write (nout,99999) t, (y(i),i=1,neq), errloc
          Else
            Write (nout,99998) t, (y(i),i=1,neq)
          End If
        End If

        Return

99999   Format (1X,F10.6,3(F13.7,2X)/1X,' ** WARNING scaled local error = ', &
          F13.5)
99998   Format (1X,F10.6,3(F13.7,2X))
      End Subroutine monitr
    End Module d02zafe_mod
    Program d02zafe

!     D02ZAF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02nbf, d02nsf, d02nvf, d02nyf, nag_wp, x04abf
      Use d02zafe_mod, Only: fcn, iset, itrace, jac, ldysav, monitr, neq, nin, &
                             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
      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)
      Logical, Allocatable                 :: algequ(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02ZAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
!     neq: number of differential equations
      Read (nin,*) maxord, maxstp, mxhnil
      sdysav = maxord + 1
      Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), &
        ysav(ldysav,sdysav),algequ(neq))
      outchn = nout
      Call x04abf(iset,outchn)

!     Set algorithmic and problem parameters

      Read (nin,*) hmin, hmax, h0, t, tout
      Read (nin,*) petzld

!     Initialisation

!     Integrate to tout without passing tout.
      tcrit = tout
      itask = 4

!     Use default values for the array con.
      con(1:6) = 0.0_nag_wp

!     Use BDF formulae with modified Newton method.
!     Use averaged L2 norm for local error control.
      ifail = 0
      Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, &
        maxstp,mxhnil,'Average-L2',rwork,ifail)

!     Setup for using analytic Jacobian
      ifail = 0
      Call d02nsf(neq,neq,'Analytical',nwkjac,rwork,ifail)

      Write (nout,*)
      Write (nout,*) '  Analytic Jacobian'
      Write (nout,*)

!     Set tolerances.
      Read (nin,*) itol
      Read (nin,*) rtol(1), atol(1)

!     Initial values for Y.
      Read (nin,*) y(1:neq)

      Write (nout,*) '    X          Y(1)           Y(2)           Y(3)'
      Write (nout,99999) t, (y(i),i=1,neq)

!     Solve the problem using MONITR to output intermediate results.
      ifail = -1
      Call d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, &
        ysav,sdysav,jac,wkjac,nwkjac,monitr,itask,itrace,ifail)

      If (ifail==0) Then

!       Get integration statistics.
        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

99999 Format (1X,F10.6,3(F13.7,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 d02zafe
