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

    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
!      .. 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                  :: iset = 1, m = 4, n = 3, nin = 5, &
                                              nout = 6
!      .. Local Scalars ..
       INTEGER, SAVE                       :: icap
    CONTAINS
!dec$ attributes offload : mic :: eqn
       SUBROUTINE eqn(e,q,p,m)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: eqn
       END SUBROUTINE eqn
!dec$ attributes offload : mic :: fcn
       SUBROUTINE fcn(x,y,f,n,p,m,i)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: fcn
       END SUBROUTINE fcn
!dec$ attributes offload : mic :: bc
       SUBROUTINE bc(g1,g2,p,m,n)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: bc
       END SUBROUTINE bc
!dec$ attributes offload : mic :: range
       SUBROUTINE range(x,npoint,p,m)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: range
       END SUBROUTINE range
!dec$ attributes offload : mic :: prsol
       SUBROUTINE prsol(z,y,n)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          INTRINSIC                              abs
!         .. Executable Statements ..
          IF (icap/=1) THEN
             icap = 1
             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) z = xend
          RETURN

99999     FORMAT (1X,F9.3,3F10.4)
!dec$ attributes offload : mic :: prsol
       END SUBROUTINE prsol
!dec$ attributes offload : mic :: constr
       FUNCTION constr(p,m)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          LOGICAL                             :: constr
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: m
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: p(m)
!         .. Intrinsic Functions ..
          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
!dec$ attributes offload : mic :: constr
       END FUNCTION constr
    END MODULE d02safe_mod
    PROGRAM d02safe

!      D02SAF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02saf, d02sas, x04abf
       USE d02safe_mod, ONLY : bc, constr, eqn, fcn, icap, iset, m, n, nag_wp, &
                               nin, nout, prsol, range
!      .. 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 Functions ..
       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)
       icap = 0

!      * 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
