!   D02HAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02hafe_mod

!     D02HAF 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        :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Integer, Parameter                   :: iset = 1, n = 3, nin = 5, nout = 6
      Integer, Parameter                   :: sdw = 3*n + 17 + max(11,n)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max
    Contains
      Subroutine fcn(x,y,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: f(*)
        Real (Kind=nag_wp), Intent (In)      :: y(*)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: cos, tan
!       .. Executable Statements ..
        f(1) = tan(y(3))
        f(2) = -0.032_nag_wp*tan(y(3))/y(2) - 0.02_nag_wp*y(2)/cos(y(3))
        f(3) = -0.032_nag_wp/y(2)**2
        Return
      End Subroutine fcn
    End Module d02hafe_mod

    Program d02hafe

!     D02HAF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02haf, nag_wp, x04abf
      Use d02hafe_mod, Only: fcn, iset, n, nin, nout, one, sdw, zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: a, b, dx, tol
      Integer                              :: i, ifail, l, m1, outchn
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: soln(:,:), x(:)
      Real (Kind=nag_wp)                   :: u(n,2), v(n,2), w(n,sdw)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02HAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
!     m1: solution is returned and printed for m1-1 grid points on [a, b].
      Read (nin,*) m1
      Allocate (soln(n,m1),x(m1))
!     a: left-hand boundary point, b: right-hand boundary point.
      Read (nin,*) a, b

!     Evaluate solution points x.
      x(1) = a
      dx = (b-a)/real(m1-1,kind=nag_wp)
      Do i = 2, m1 - 1
        x(i) = x(i-1) + dx
      End Do
      x(m1) = b

!     Set output channel for monitoring information.
      outchn = nout
      Call x04abf(iset,outchn)

!     Flag known (zero) and estimated (one) values in u 
      v(1:2,1:2) = zero
      v(2,2) = one
      v(3,1:2) = one
!     Set known values of u
      u(1,1:2) = zero
      u(2,1) = 0.5_nag_wp

loop: Do l = 4, 5
        tol = 5.0_nag_wp*10.0_nag_wp**(-l)
        Write (nout,*)
!       Set estimates of u
        u(2,2) = 0.46_nag_wp
        u(3,1) = 1.15_nag_wp
        u(3,2) = -1.2_nag_wp

!       ifail: behaviour on error exit
!              =1 for quiet-soft exit
!       * Set ifail to 111 to obtain monitoring information *
        ifail = 1
        Call d02haf(u,v,n,a,b,tol,fcn,soln,m1,w,sdw,ifail)

        If (ifail>=0) Then
          Write (nout,99999) 'Results with TOL = ', tol
          Write (nout,*)
          If (ifail==0) Then
            Write (nout,*) ' X-value and final solution'
            Do i = 1, m1
              If (l==4) Then
                Write (nout,99998) x(i), soln(1:n,i)
              Else
                Write (nout,99997) x(i), soln(1:n,i)
              End If
            End Do
          Else
            Write (nout,99996) ' IFAIL =', ifail
          End If
        Else
          Write (nout,99995) ifail
          Exit loop
        End If
      End Do loop

99999 Format (1X,A,E10.3)
99998 Format (1X,F4.1,3(1X,F9.3))
99997 Format (1X,F4.1,1X,3F10.4)
99996 Format (1X,A,I4)
99995 Format (1X/1X,' ** D02HAF returned with IFAIL = ',I5)
    End Program d02hafe
