Program d02qgfe

!     D02QGF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
      Use nag_library, Only: d02qgf, d02qwf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: alpha = -0.032_nag_wp
      Real (Kind=nag_wp), Parameter    :: beta = -0.02_nag_wp
      Integer, Parameter               :: neqf = 3, neqg = 1, nin = 5,         &
                                          nout = 6
      Integer, Parameter               :: liwork = 21 + 4*neqg
      Integer, Parameter               :: lrwork = 23 + 23*neqf + 14*neqg
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: grvcm, hmax, t, tcrit, tinc, tout,   &
                                          trvcm, tstart
      Integer                          :: i, ifail, irevcm, j, kgrvcm, latol,  &
                                          lrtol, maxstp, yprvcm, yrvcm
      Logical                          :: alterg, crit, onestp, root, sophst,  &
                                          vectol
      Character (1)                    :: statef
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: atol(:), rtol(:), rwork(:), y(:)
      Integer, Allocatable             :: iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cos, tan
!     .. Executable Statements ..
      Write (nout,*) 'D02QGF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) latol, lrtol

      Allocate (atol(latol),rtol(lrtol),rwork(lrwork),y(neqf),iwork(liwork))

      Read (nin,*) hmax, tstart, tcrit, tinc
      Read (nin,*) statef
      Read (nin,*) vectol, onestp, crit, sophst
      Read (nin,*) maxstp
      Read (nin,*) rtol(1:lrtol), atol(1:latol)
      Read (nin,*) y(1:neqf)

      t = tstart

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d02qwf(statef,neqf,vectol,atol,latol,rtol,lrtol,onestp,crit,tcrit,  &
        hmax,maxstp,neqg,alterg,sophst,rwork,lrwork,iwork,liwork,ifail)

      Write (nout,*)
      Write (nout,*) '  T         Y(1)     Y(2)     Y(3)'
      Write (nout,99999) t, (y(i),i=1,neqf)
      j = 1
      tout = tinc
      irevcm = 0

revcm: Do
        ifail = -1
        Call d02qgf(neqf,t,y,tout,neqg,root,irevcm,trvcm,yrvcm,yprvcm,grvcm,   &
          kgrvcm,rwork,lrwork,iwork,liwork,ifail)

        Select Case (irevcm)
        Case (0)
          If (ifail==0) Then
!           Print solution at current t
            Write (nout,99999) t, (y(i),i=1,neqf)
            If (t==tout .And. j<5) Then
!             Increment tout and cycle to find solution at this new time.
              j = j + 1
              tout = tout + tinc
              Cycle revcm
            End If
          End If
          Exit revcm
        Case (1:7)
          If (yrvcm==0) Then
            rwork(yprvcm) = tan(y(3))
            rwork(yprvcm+1) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3))
            rwork(yprvcm+2) = alpha/y(2)**2
          Else
            rwork(yprvcm) = tan(rwork(yrvcm+2))
            rwork(yprvcm+1) = alpha*tan(rwork(yrvcm+2))/rwork(yrvcm+1) +       &
              beta*rwork(yrvcm+1)/cos(rwork(yrvcm+2))
            rwork(yprvcm+2) = alpha/rwork(yrvcm+1)**2
          End If
        Case (9:)
          grvcm = y(1)
        End Select
      End Do revcm

99999 Format (1X,F7.4,2X,3(F7.4,2X))
    End Program d02qgfe