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

    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
    CONTAINS
       SUBROUTINE fcn(x,y,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              cos, tan
!         .. Executable Statements ..
          f(1) = tan(y(3))
          f(2) = -0.032E0_nag_wp*tan(y(3))/y(2) - 0.02E0_nag_wp*y(2)/cos(y(3))
          f(3) = -0.032E0_nag_wp/y(2)**2
          RETURN
       END SUBROUTINE fcn
    END MODULE d02hafe_mod

    PROGRAM d02hafe

!      D02HAF Example Main Program

!      N.B the definition of SDW must be changed for N.GT.11

!      .. Use Statements ..
       USE nag_library, ONLY : d02haf, nag_wp, x04abf
       USE d02hafe_mod, ONLY : fcn, iset, n, nin, nout, one, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: a, b, tol
       INTEGER                             :: i, ifail, j, l, m1, outchn, sdw
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: soln(:,:), u(:,:), v(:,:), w(:,:)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02HAF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
!      m1: controls output.
       READ (nin,*) m1
       sdw = 3*n + 17 + 11
       ALLOCATE (soln(n,m1),u(n,2),v(n,2),w(n,sdw))
!      a: left-hand boundary point, b: right-hand boundary point.
       READ (nin,*) a, b
       outchn = nout
       CALL x04abf(iset,outchn)
LOOP:  DO l = 3, 4
          tol = 5.0_nag_wp*10.0_nag_wp**(-l)
          WRITE (nout,*)
          u(1,1:2) = zero
          u(2,1) = 0.5_nag_wp
          u(2,2) = 0.46_nag_wp
          u(3,1) = 1.15_nag_wp
          u(3,2) = -1.2_nag_wp
          v(1:2,1:2) = zero
          v(2,2) = one
          v(3,1:2) = one

!         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
                   WRITE (nout,99998) i - 1, (soln(j,i),j=1,n)
                END DO
             ELSE
                WRITE (nout,99997) ' IFAIL =', ifail
             END IF
          ELSE
             WRITE (nout,99996) ifail
             EXIT LOOP
          END IF
       END DO LOOP

99999  FORMAT (1X,A,E10.3)
99998  FORMAT (1X,I3,3F10.4)
99997  FORMAT (1X,A,I4)
99996  FORMAT (1X/1X,' ** D02HAF returned with IFAIL = ',I5)
    END PROGRAM d02hafe