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

    Module d02mzfe_mod

!     D02MZF 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        :: tstep = 0.02_nag_wp
      Integer, Parameter                   :: iset = 1, neq = 3, nin = 5,      &
                                              nout = 6
      Integer, Parameter                   :: nrw = 50 + 4*neq
      Integer, Parameter                   :: nwkjac = neq*(neq+1)
      Integer, Parameter                   :: sdysav = 6
      Integer, Parameter                   :: ldysav = neq
    Contains
      Subroutine resid(neq,t,y,ydot,r,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)     :: r(neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq), ydot(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: alpha = 0.04_nag_wp
        Real (Kind=nag_wp)                   :: beta = 1.0E4_nag_wp
        Real (Kind=nag_wp)                   :: gamma = 3.0E7_nag_wp
!       .. Executable Statements ..
        r(1:3) = -ydot(1:3)
        If (ires==1) Then
          r(1) = -alpha*y(1) + beta*y(2)*y(3) + r(1)
          r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) + r(2)
          r(3) = gamma*y(2)*y(2) + r(3)
        End If
        Return
      End Subroutine resid
    End Module d02mzfe_mod

    Program d02mzfe

!     D02MZF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02mzf, d02nby, d02ngf, d02ngz, d02nsf, d02nvf,   &
                             d02nyf, nag_wp, x04abf
      Use d02mzfe_mod, Only: iset, ldysav, neq, nin, nout, nrw, nwkjac, resid, &
                             sdysav, tstep
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: h, h0, hmax, hmin, hu, t, tcrit, &
                                              tcur, tolsf, tout, xout
      Integer                              :: i, ifail, imxer, iout, itask,    &
                                              itol, itrace, maxord, maxstp,    &
                                              mxhnil, niter, nje, nq, nqu,     &
                                              nre, nst, outchn, pstat
      Logical                              :: petzld
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: atol(:), rtol(:), rwork(:),      &
                                              sol(:), wkjac(:), y(:), ydot(:), &
                                              ysav(:,:)
      Real (Kind=nag_wp)                   :: con(6)
      Integer                              :: inform(23)
      Logical, Allocatable                 :: algequ(:)
      Logical                              :: lderiv(2)
!     .. Executable Statements ..
      Write (nout,*) 'D02MZF Example Program Results'
!     Skip heading in data file
      Read (nin,*)

!     Allocations based on number of differential equations (neq)
      Allocate (atol(neq),rtol(neq),rwork(nrw),sol(neq),wkjac(nwkjac),y(neq), &
        ydot(neq),ysav(ldysav,sdysav),algequ(neq))

!     Read algorithmic parameters
      Read (nin,*) maxord, maxstp, mxhnil, pstat
      Read (nin,*) petzld
      Read (nin,*) hmin, hmax, h0, tcrit
      Read (nin,*) t, tout
      Read (nin,*) itol
      Read (nin,*) y(1:neq)
      Read (nin,*) lderiv(1:2)
      Read (nin,*) rtol(1:neq)
      Read (nin,*) atol(1:neq)

      outchn = nout
      Call x04abf(iset,outchn)
      con(1:6) = 0.0_nag_wp
      itask = 2
      itrace = 0

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

!     Linear algebra setup.
      ifail = 0
      Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail)


      Write (nout,99994)(i,i=1,neq)
      Write (nout,99999) t, y(1:neq)

!     First value of t to interpolate and print solution at.
      iout = 1
      xout = tstep

!      Integrate one step at a time and overshoot tout (itask=2).
steps: Do
        ifail = 0
        Call d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, &
          resid,ysav,sdysav,d02ngz,wkjac,nwkjac,d02nby,lderiv,itask,itrace, &
          ifail)

!       Get Current value of t (tcur) and last step used (hu)
        Call d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, &
          imxer,algequ,inform,ifail)

interp: Do
          If (tcur-hu<xout .And. xout<=tcur) Then
!           xout is in interval of last step, so interpolate.
            ifail = 0
            Call d02mzf(xout,sol,neq,ldysav,neq,ysav,sdysav,rwork,ifail)

            Write (nout,99999) xout, sol(1:neq)
            iout = iout + 1
            If (iout>=6) Then
!             Final solution point printed. Print stats if required.
              If (pstat==1) Then
                Write (nout,*)
                Write (nout,99998) hu, h, tcur
                Write (nout,99997) nst, nre, nje
                Write (nout,99996) nqu, nq, niter
                Write (nout,99995) ' Max err comp = ', imxer
              End If
              Exit steps
            End If
!           Set next xout. This might also be in last step, so keep
!           looping for interpolation.
            xout = xout + tstep
            Cycle interp
          Else
!           Take another step.
            Cycle steps
          End If
        End Do interp
      End Do steps

99999 Format (1X,F8.3,3(F13.5,2X))
99998 Format (1X,' HUSED = ',E12.5,'  HNEXT = ',E12.5,'  TCUR = ',E12.5)
99997 Format (1X,' NST = ',I6,'    NRE = ',I6,'    NJE = ',I6)
99996 Format (1X,' NQU = ',I6,'    NQ  = ',I6,'  NITER = ',I6)
99995 Format (1X,A,I4)
99994 Format (/1X,'    X ',3('         Y(',I1,')  '))
    End Program d02mzfe