NAG Library Manual, Mark 29
```!   D02SAF Example Program Text
!   Mark 29.0 Release. NAG Copyright 2023.

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
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
Call x04abf(iset,outchn)
swp(1:npoint-1,1:3) = 0.0_nag_wp

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