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

    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
!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) = -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
!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) = 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
!dec$ attributes offload : mic :: jac
       END SUBROUTINE jac
!dec$ attributes offload : mic :: monitr
       SUBROUTINE monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon, &
             inln,hmin,hmax,nqu)

!         .. Use Statements ..
          USE nag_library, ONLY : d02zaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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))
!dec$ attributes offload : mic :: monitr
       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
