!   D02PGF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module d02pgfe_mod

!     D02PGF Example Program Module:
!            Parameters

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: tol = 1.0E-5_nag_wp
      Integer, Parameter, Public       :: liwsav = 130, n = 2, nin = 5,        &
                                          nout = 6
      Integer, Parameter, Public       :: lrwsav = 350 + 32*n
      Integer, Parameter, Public       :: lwcomm = 6*n
    End Module d02pgfe_mod

    Program d02pgfe

!     D02PGF Example Main Program

!     .. Use Statements ..
      Use d02pgfe_mod, Only: liwsav, lrwsav, lwcomm, n, nin, nout, tol
      Use nag_library, Only: c05azf, d02pgf, d02phf, d02pjf, d02pqf, d02ptf,   &
                             nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: hnext, hstart, t, t1, t2, tend,      &
                                          tnow, tout, tprev, waste
      Integer                          :: i, icheck, ifail, ind, irevcm, j,    &
                                          method, nchange, stpcst, stpsok,     &
                                          totf
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: c(17)
      Real (Kind=nag_wp), Allocatable  :: rwsav(:), thres(:), troot(:),        &
                                          wcomm(:), y(:), ynow(:), yout(:),    &
                                          yp(:), ypnow(:), yprev(:)
      Integer, Allocatable             :: iroot(:), iwsav(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02PGF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) method
      Allocate (thres(n),iwsav(liwsav),rwsav(lrwsav),ynow(n),ypnow(n),         &
        yprev(n),wcomm(lwcomm),yout(n),iroot(n),y(n),yp(n),troot(n))

!     Set initial conditions and input for D02PQF

      Read (nin,*) t, tend
      Read (nin,*) ynow(1:n)
      Read (nin,*) hstart
      Read (nin,*) thres(1:n)

!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d02pqf(n,t,tend,ynow,tol,thres,method,hstart,iwsav,rwsav,ifail)

      Write (nout,99999) tol
      Write (nout,99998)
      Write (nout,99997) t, ynow(1:n)

      tout = 0.1_nag_wp
      tnow = t
      Do While (tnow<tend)
!       Integrate by one time step using reverse communication
        tprev = tnow
        yprev(1:n) = ynow(1:n)
        ifail = 1
        irevcm = 0
        Do While (irevcm>=0)
          Call d02pgf(irevcm,n,tnow,ynow,ypnow,iwsav,rwsav,ifail)
          If (irevcm>0) Then
            ypnow(1) = ynow(2)
            ypnow(2) = -ynow(1)
          End If
        End Do
        If (ifail==1 .Or. ifail>4) Then
          Write (6,99996) ifail, tnow
          Go To 100
        End If

!       Detect sign changes in last step
        iroot(1:n) = 0
        nchange = 0
        Do i = 1, n
          If (ynow(i)*yprev(i)<0.0_nag_wp) Then
            nchange = nchange + 1
            iroot(nchange) = i
          End If
        End Do

!       Interpolate for values of t at increments of 0.1 within step and
!       find roots of components of y.
        If (tnow>=tout .Or. nchange>0) Then
!         Set up interpolation by reverse communication
          ifail = 0
          irevcm = 0
          Do While (irevcm>=0)
            ifail = 0
            Call d02phf(irevcm,n,n,t,y,yp,wcomm,lwcomm,iwsav,rwsav,ifail)
            If (irevcm>0) Then
              yp(1) = y(2)
              yp(2) = -y(1)
            End If
          End Do

          icheck = 1
!         If there are sign changes: find roots
          Do i = 1, nchange
            j = iroot(i)
            t1 = tprev
            t2 = tnow
            ind = 1
            ifail = 0
            Do While (ind/=0)
              Call c05azf(t1,t2,y(j),tol,1,c,ind,ifail)
              If (ind>1) Then
                ifail = 0
                Call d02pjf(icheck,n,n,t1,0,y,wcomm,lwcomm,iwsav,rwsav,ifail)
                icheck = 0
              End If
            End Do
            troot(i) = t1
          End Do
!         Evaluate Interpolant at increments and print any contained roots
          ifail = 0
          Do While (tnow>=tout)
            Do i = 1, nchange
              If (troot(i)<tout .And. iroot(i)>0) Then
                Write (nout,99995) iroot(i), troot(i)
                iroot(i) = -iroot(i)
              End If
            End Do
            Call d02pjf(icheck,n,n,tout,0,yout,wcomm,lwcomm,iwsav,rwsav,ifail)
            icheck = 0
            Write (nout,99997) tout, yout(1:n)
            tout = tout + 0.1_nag_wp
          End Do
!         Print any remaining roots up to current time
          Do i = 1, nchange
            If (iroot(i)>0) Then
              Write (nout,99995) iroot(i), troot(i)
            End If
          End Do
        End If
      End Do

!     Print some integration statistics
      ifail = 0
      Call d02ptf(totf,stpcst,waste,stpsok,hnext,iwsav,rwsav,ifail)
      Write (nout,99994) totf
100   Continue

99999 Format (/,' Calculation with TOL = ',E8.1)
99998 Format (/,'    t         y1        y2',/)
99997 Format (1X,F6.3,2(3X,F8.4))
99996 Format (/,'D02PGF returned with ifail = ',I3,' at t = ',E12.5,'.')
99995 Format ('Component ',I2,' has a root at t = ',F7.4)
99994 Format (/,' Cost of the integration in evaluations of F is',I6)
    End Program d02pgfe