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

    Module d02qffe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: neqf = 2, neqg = 2, nin = 5,     &
                                              nout = 6
      Integer, Parameter                   :: latol = neqf
      Integer, Parameter                   :: liwork = 21 + 4*neqg
      Integer, Parameter                   :: lrtol = neqf
      Integer, Parameter                   :: lrwork = 23 + 23*neqf + 14*neqg
    Contains
      Subroutine fcn(neqf,x,y,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neqf
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: f(neqf)
        Real (Kind=nag_wp), Intent (In)      :: y(neqf)
!       .. Executable Statements ..
        f(1) = y(2)
        f(2) = -y(1)
        Return
      End Subroutine fcn

      Function g(neqf,x,y,yp,k)

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: g
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: k, neqf
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: y(neqf), yp(neqf)
!       .. Executable Statements ..
        If (k==1) Then
          g = yp(1)
        Else
          g = y(1)
        End If
        Return
      End Function g
    End Module d02qffe_mod

    Program d02qffe

!     D02QFF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02qff, d02qwf, d02qxf, d02qyf, nag_wp
      Use d02qffe_mod, Only: fcn, g, latol, liwork, lrtol, lrwork, neqf, neqg, &
                             nin, nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: hlast, hmax, hnext, t, tcrit,    &
                                              tcurr, tolfac, tout, tstart
      Integer                              :: badcmp, i, ifail, index, maxstp, &
                                              nfail, nsucc, odlast, odnext, type
      Logical                              :: alterg, crit, onestp, root,      &
                                              sophst, vectol
      Character (1)                        :: statef
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: atol(:), resids(:), rtol(:),     &
                                              rwork(:), y(:), yp(:)
      Integer, Allocatable                 :: events(:), iwork(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02QFF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Allocate (atol(latol),resids(neqg),rtol(lrtol),rwork(lrwork),y(neqf), &
        yp(neqf),events(neqg),iwork(liwork))
      Read (nin,*) hmax, tstart, tcrit
      Read (nin,*) statef
      Read (nin,*) vectol, onestp, crit, sophst
      Read (nin,*) maxstp
      Read (nin,*) rtol(1:neqf)
      Read (nin,*) atol(1:neqf)

!     Initialize
      ifail = 0
      Call d02qwf(statef,neqf,vectol,atol,latol,rtol,lrtol,onestp,crit,tcrit, &
        hmax,maxstp,neqg,alterg,sophst,rwork,lrwork,iwork,liwork,ifail)

      t = tstart
      tout = tcrit
      Read (nin,*) y(1:neqf)

!      Cycle through roots and print info when encountered.
findr: Do
        ifail = -1
        Call d02qff(fcn,neqf,t,y,tout,g,neqg,root,rwork,lrwork,iwork,liwork, &
          ifail)

        If (ifail/=0) Exit findr

        ifail = 0
        Call d02qxf(neqf,yp,tcurr,hlast,hnext,odlast,odnext,nsucc,nfail, &
          tolfac,badcmp,rwork,lrwork,iwork,liwork,ifail)

        If (.Not. root) Exit findr

        ifail = 0
        Call d02qyf(neqg,index,type,events,resids,rwork,lrwork,iwork,liwork, &
          ifail)

        Write (nout,99999) t
        Write (nout,99998) index, type, resids(index)
        Write (nout,99997) y(1), yp(1)

        Do i = 1, neqg
          If (i/=index) Then
            If (events(i)/=0) Then
              Write (nout,99996) i, events(i), resids(i)
            End If
          End If
        End Do

        If (tcurr>=tout) Exit findr

      End Do findr

99999 Format (/1X,'Root at ',1P,E13.5)
99998 Format (1X,'for event equation ',I2,' with type',I3,' and residual ',1P, &
        E13.5)
99997 Format (1X,' Y(1) = ',1P,E13.5,'   Y''(1) = ',1P,E13.5)
99996 Format (1X,'and also for event equation ',I2,' with type',I3, &
        ' and residual ',1P,E13.5)
    End Program d02qffe
