!   D02LAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02lafe_mod

!      D02LAF 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       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: neq = 2, nin = 5, nout = 6
       INTEGER, PARAMETER                  :: lrwork = 16 + 20*neq
    CONTAINS
       SUBROUTINE fcn(neq,t,y,f)

!         Derivatives for two body problem in  y'' = f(t,y)  form

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: r
!         .. Intrinsic Functions ..
          INTRINSIC                              sqrt
!         .. Executable Statements ..
          r = sqrt(y(1)**2+y(2)**2)**3
          f(1) = -y(1)/r
          f(2) = -y(2)/r
          RETURN
       END SUBROUTINE fcn
    END MODULE d02lafe_mod

    PROGRAM d02lafe

!      D02LAF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02laf, d02lxf, d02lyf, d02lzf, nag_wp
       USE d02lafe_mod, ONLY : fcn, lrwork, neq, nin, nout, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, hnext, hstart, hused, t,      &
                                              tend, tinc, tnext, tol, tstart
       INTEGER                             :: i, ifail, itol, maxstp, natt,    &
                                              nfail, nsucc, nwant
       LOGICAL                             :: high, onestp, start
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: rwork(:), thres(:), thresp(:),   &
                                              y(:), ydp(:), yinit(:), yp(:),   &
                                              ypinit(:), ypwant(:), ywant(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02LAF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
!      neq: number of second-order ordinary differential equations
       READ (nin,*) nwant
       ALLOCATE (rwork(lrwork),thres(neq),thresp(neq),y(neq),ydp(neq), &
          yinit(neq),yp(neq),ypinit(neq),ypwant(nwant),ywant(nwant))
       READ (nin,*) high, onestp
       READ (nin,*) tinc

!      Initial conditions
       READ (nin,*) tstart, tend
       READ (nin,*) yinit(1:neq)
       READ (nin,*) ypinit(1:neq)

LOOP1: DO itol = 4, 5
          tol = 10.0_nag_wp**(-itol)
          WRITE (nout,*)

!         Call D02LXF with default THRES,THRESP,MAXSTP and H

          thres(1) = zero
          thresp(1) = zero
          h = zero
          maxstp = 0
          start = .TRUE.

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          CALL d02lxf(neq,h,tol,thres,thresp,maxstp,start,onestp,high,rwork, &
             lrwork,ifail)

          WRITE (nout,99999) 'Calculation with TOL = ', tol
          WRITE (nout,99995) (i,i=1,neq)

!         Set initial values

          y(1:neq) = yinit(1:neq)
          yp(1:neq) = ypinit(1:neq)
          t = tstart
          tnext = t + tinc
          WRITE (nout,99998) t, y(1:neq)

!         Loop point for onestep mode
LOOP2:    DO

             ifail = -1
             CALL d02laf(fcn,neq,t,tend,y,yp,ydp,rwork,lrwork,ifail)

             IF (ifail>0) THEN
                WRITE (nout,99997) ifail, t
                EXIT LOOP1
             END IF

!            Loop point for interpolation
             DO WHILE (tnext<=t)

                ifail = 0
                CALL d02lzf(neq,t,y,yp,neq,tnext,ywant,ypwant,rwork,lrwork, &
                   ifail)

                WRITE (nout,99998) tnext, ywant(1:neq)
                tnext = tnext + tinc
             END DO

             IF (t>=tend) EXIT LOOP2

          END DO LOOP2

          ifail = 0
          CALL d02lyf(neq,hnext,hused,hstart,nsucc,nfail,natt,thres,thresp, &
             rwork,lrwork,ifail)

          WRITE (nout,*)
          WRITE (nout,99996) ' Number of successful steps = ', nsucc
          WRITE (nout,99996) ' Number of   failed   steps = ', nfail
       END DO LOOP1

99999  FORMAT (1X,A,1P,E9.1)
99998  FORMAT (1X,F5.1,2(2X,F9.5))
99997  FORMAT (/1X,'D02LAF returned with IFAIL = ',I2,'  at T = ',1P,E10.3)
99996  FORMAT (1X,A,I5)
99995  FORMAT (/'   T ',2('     Y(',I1,')  '))
    END PROGRAM d02lafe