PROGRAM d02qgfe

!      D02QGF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       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