!   D02SAF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module d02safe_mod

!     D02SAF 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                           :: bc, constr, eqn, fcn, prsol, range
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: alpha = 0.032_nag_wp
      Real (Kind=nag_wp), Parameter    :: beta = 0.02_nag_wp
      Real (Kind=nag_wp), Parameter    :: xend = 5.0_nag_wp
      Integer, Parameter, Public       :: iset = 1, m = 4, n = 3, nin = 5,     &
                                          nout = 6
    Contains
      Subroutine eqn(e,q,p,m)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m, q
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: e(q)
        Real (Kind=nag_wp), Intent (In) :: p(m)
!       .. Executable Statements ..
        e(1) = 0.02_nag_wp - p(4) - 1.0E-5_nag_wp*p(3)
        Return
      End Subroutine eqn
      Subroutine fcn(x,y,f,n,p,m,i)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: i, m, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(n)
        Real (Kind=nag_wp), Intent (In) :: p(m), y(n)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, tan
!       .. Executable Statements ..
        f(1) = tan(y(3))
        If (i==1) Then
          f(2) = -alpha*tan(y(3))/y(2) - beta*y(2)/cos(y(3))
          f(3) = -alpha/y(2)**2
        Else
          f(2) = -p(2)*tan(y(3))/y(2) - p(4)*y(2)/cos(y(3))
          f(3) = -p(2)/y(2)**2
        End If
        Return
      End Subroutine fcn
      Subroutine bc(g1,g2,p,m,n)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: g1(n), g2(n)
        Real (Kind=nag_wp), Intent (In) :: p(m)
!       .. Executable Statements ..
        g1(1) = 0.0_nag_wp
        g1(2) = 0.5_nag_wp
        g1(3) = p(1)
        g2(1) = 0.0_nag_wp
        g2(2) = 0.45_nag_wp
        g2(3) = -1.2_nag_wp
        Return
      End Subroutine bc
      Subroutine range(x,npoint,p,m)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m, npoint
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: p(m)
        Real (Kind=nag_wp), Intent (Out) :: x(npoint)
!       .. Executable Statements ..
        x(1) = 0.0_nag_wp
        x(2) = p(3)
        x(3) = xend
        Return
      End Subroutine range
      Subroutine prsol(z,y,n)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: z
        Integer, Intent (In)           :: n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: y(n)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. Executable Statements ..
        If (z==0.0E0_nag_wp) Then
          Write (nout,*)
          Write (nout,*) '      Z       Y(1)      Y(2)      Y(3)'
        End If
        Write (nout,99999) z, (y(i),i=1,n)
        z = z + 0.5_nag_wp
        If (abs(z-xend)<0.25_nag_wp) Then
          z = xend
        End If
        Return

99999   Format (1X,F9.3,3F10.4)
      End Subroutine prsol
      Function constr(p,m)

!       .. Function Return Value ..
        Logical                        :: constr
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: p(m)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: any
!       .. Executable Statements ..
        If (any(p(1:m)<0.0_nag_wp) .Or. p(3)>5.0_nag_wp) Then
          constr = .False.
        Else
          constr = .True.
        End If
        Return
      End Function constr
    End Module d02safe_mod
    Program d02safe

!     D02SAF Example Main Program

!     .. Use Statements ..
      Use d02safe_mod, Only: bc, constr, eqn, fcn, iset, m, n, nin, nout,      &
                             prsol, range
      Use nag_library, Only: d02saf, d02sas, nag_wp, x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: ymax
      Integer                          :: i, icount, ifail, ldswp, ldw, n1,    &
                                          npoint, outchn, sdw
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dp(:), e(:), p(:), pe(:), pf(:),     &
                                          swp(:,:), w(:,:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'D02SAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) npoint
      n1 = n
      sdw = 3*m + 23
      ldswp = npoint
      ldw = max(m,n)
      Allocate (dp(m),e(n),p(m),pe(m),pf(m),swp(ldswp,6),w(ldw,sdw))

      outchn = nout
      Read (nin,*) icount
      Read (nin,*) ymax
      Read (nin,*) pe(1:m)
      Read (nin,*) pf(1:m)
      Read (nin,*) dp(1:m)
      Read (nin,*) e(1:n)
      Call x04abf(iset,outchn)
      swp(1:npoint-1,1:3) = 0.0_nag_wp
      Read (nin,*) p(1:m)

!     * To obtain monitoring information, replace the name d02sas by d02hbx
!     in the next statement and USE nag_library : d02hbx

!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 1
      Call d02saf(p,m,n,n1,pe,pf,e,dp,npoint,swp,ldswp,icount,range,bc,fcn,    &
        eqn,constr,ymax,d02sas,prsol,w,ldw,sdw,ifail)

      If (ifail/=0) Then
        Write (nout,99999) ifail
      End If
      If (ifail>=4 .And. ifail<=12) Then
        Write (nout,99998) 'SWP(NPOINT,1) = ', swp(npoint,1)
        If (ifail<=6) Then
          Write (nout,99998) 'SWP(1,5) = ', swp(1,5)
          Write (nout,*) '  i     W(i,1) '
          Write (nout,99997)(i,w(i,1),i=1,n)
        End If
      End If

99999 Format (1X,/,1X,' ** D02SAF returned with IFAIL = ',I5)
99998 Format (1X,A,F10.4)
99997 Format (1X,I4,1X,E10.3)
    End Program d02safe