```!   D01ESF Example Program Text
!   Mark 27.0 Release. NAG Copyright 2019.

Module d01esfe_mod

!     D01ESF Example Program Module:
!            User-defined Routines

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                           :: f
Contains
Subroutine f(ni,ndim,nx,xtr,nntr,icolzp,irowix,xs,qs,fm,iflag,iuser,     &
ruser)

!       .. Use Statements ..
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: xtr
Integer, Intent (Inout)        :: iflag
Integer, Intent (In)           :: ndim, ni, nntr, nx
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: fm(ni,nx)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: xs(nntr)
Integer, Intent (In)           :: icolzp(nx+1), irowix(nntr), qs(nntr)
Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: s_ntr, s_tr
Integer                        :: i, j, logs_hi, logs_lo, s_hi, s_lo,  &
tid
!       .. Intrinsic Procedures ..
Intrinsic                      :: log, real, sin, sum
!       .. Executable Statements ..

!       For each evaluation point x_i, i = 1, ..., nx, return in fm the
!       computed values of the ni integrands f_j, j = 1, ..., ni defined by

!         fm(j,i) = f_j(x_i)
!                                                         ndim
!                 = sin(j + S(i))*log(S(i)), where S(i) =  Sum  k*x_i(k).
!                                                          k=1

!       Split the S expression into two components, one involving only the
!       'trivial' value xtr:

!                ndim             ndim
!         S(i) =  Sum  (k*xtr)  +  Sum  (k*(x_i(k)-xtr))
!                 k=1              k=1

!                      ndim*(ndim+1)   ndim
!              = xtr * ------------- +  Sum  (k*(x_i(k)-xtr))
!                           2           k=1

!             := s_tr                + s_ntr(i)

!       By definition the summands in the s_ntr(i) term on the right-hand side
!       are zero for those k outside the range of indices defined in irowix.

!       As a demonstration of safely operating with the user arrays iuser and
!       ruser when running in parallel, 'partition' these based on the current
!       thread number. Store some of the s_tr and s_ntr computations in these
!       array sections.

!       The thread number, converted to 1-based numbering.

s_lo = iuser(tid)
s_hi = s_lo + nx - 1
logs_lo = s_hi + 1
logs_hi = logs_lo + nx - 1

If (iflag==0) Then
!         First call: nx=1, no non-trivial dimensions.
!         The constant s_tr can be reused by all subsequent calculations.
s_tr = 0.5E0_nag_wp*xtr*real(ndim*(ndim+1),kind=nag_wp)
ruser(1) = s_tr
ruser(s_lo) = s_tr
ruser(logs_lo) = log(s_tr)
Else
!         Calculate S(i) = s_tr + s_ntr(i).
s_tr = ruser(1)
Do i = 1, nx
s_ntr = sum(real(irowix(icolzp(i):icolzp(i+1)-                     &
1),kind=nag_wp)*(xs(icolzp(i):icolzp(i+1)-1)-xtr))
ruser(s_lo+i-1) = s_ntr + s_tr
ruser(logs_lo+i-1) = log(s_ntr+s_tr)
End Do
End If

!       Finally we obtain fm(j,:) = sin(j+S(:))*log(S(:))
Do j = 1, ni
fm(j,:) = sin(real(j,kind=nag_wp)+ruser(s_lo:s_hi))*                 &
ruser(logs_lo:logs_hi)
End Do

Return
End Subroutine f
End Module d01esfe_mod
Program d01esfe

!     .. Use Statements ..
Use d01esfe_mod, Only: f
Use nag_library, Only: d01esf, d01zkf, d01zlf, nag_wp, x06acf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: rvalue
Integer                          :: ifail, j, liuser, lruser, maxnx,     &
ndim, ni, optype, smpthd
Character (16)                   :: cvalue
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: dinest(:), errest(:), opts(:),       &
ruser(:)
Integer, Allocatable             :: iopts(:), iuser(:), ivalid(:),       &
maxdlv(:)
!     .. Executable Statements ..
Write (nout,*) 'D01ESF Example Program Results'
Write (nout,*)
ni = 10
ndim = 4

Allocate (iopts(100),opts(100),ivalid(ni),dinest(ni),errest(ni),         &
maxdlv(ndim))

!     Initialize option arrays.
ifail = 0
Call d01zkf('Initialize = D01ESF',iopts,100,opts,100,ifail)

!     Set any required options.
Call d01zkf('Absolute Tolerance = 0.0',iopts,100,opts,100,ifail)
Call d01zkf('Relative Tolerance = 1.0e-3',iopts,100,opts,100,ifail)
Call d01zkf('Maximum Level = 6',iopts,100,opts,100,ifail)
Call d01zkf('Index Level = 5',iopts,100,opts,100,ifail)

!     Set any required maximum dimension levels.
maxdlv(:) = 0

!     As a demonstration of safely operating with the user arrays iuser and
!     ruser when running in parallel, we will 'partition' these in the user-
!     supplied function f based on the current thread number.
!     The size of these arrays is a function of Maximum Nx and the maximum
!     allowed number of OpenMP threads.

ifail = 0
Call d01zlf('Maximum Nx',maxnx,rvalue,cvalue,optype,iopts,opts,ifail)

smpthd = x06acf()

lruser = 1 + 2*maxnx*smpthd
liuser = smpthd
Allocate (iuser(liuser),ruser(lruser))

!     iuser stores the partition indices for ruser:
iuser(1) = 2
Do j = 2, smpthd
iuser(j) = iuser(j-1) + 2*maxnx
End Do

!     Approximate the integrals.
ifail = -1
Call d01esf(ni,ndim,f,maxdlv,dinest,errest,ivalid,iopts,opts,iuser,      &
ruser,ifail)
Select Case (ifail)
Case (0,1,2,-1)
!       0: The result returned satisfies the requested accuracy requirements.
!       1, 2: The result returned is inaccurate for at least one integral.
!       -1: Exit was requested by setting iflag negative in f.
!           A result will be returned if at least one call to f was
!           successful.
Write (nout,99999)
Do j = 1, ni
Write (nout,99998) j, dinest(j), errest(j), ivalid(j)
End Do
Case Default
!       If internal memory allocation failed consider reducing the options
!       'Maximum Nx' and 'Index Level', or run with fewer threads.
Write (nout,99997) ifail
End Select

99999 Format (1X,'Integral # | Estimated value | Error estimate | ',           &
'Final state of integral')
99998 Format (1X,I11,'|',Es17.5,'|',Es16.5,'|',I8)
99997 Format (1X,'D01ESF exited with IFAIL = ',I8)
End Program d01esfe
```