!   D02LAF Example Program Text
!   Mark 29.3 Release. NAG Copyright 2023.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: neq = 2, nin = 5, nout = 6
      Integer, Parameter, Public       :: lrwork = 16 + 20*neq
    Contains
      Subroutine fcn(neq,t,y,f)

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

!       .. 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 Procedures ..
        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 d02lafe_mod, Only: fcn, lrwork, neq, nin, nout, zero
      Use nag_library, Only: d02laf, d02lxf, d02lyf, d02lzf, nag_wp
!     .. 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 one-step 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) Then
            Exit loop2
          End If

        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