```!   D03PXF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.

Module d03pxfe_mod

!     D03PXF 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                           :: bndary, numflx
!     .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: alpha_l = 460.894_nag_wp
Real (Kind=nag_wp), Parameter, Public :: alpha_r = 46.095_nag_wp
Real (Kind=nag_wp), Parameter, Public :: beta_l = 19.5975_nag_wp
Real (Kind=nag_wp), Parameter, Public :: beta_r = 6.19633_nag_wp
Real (Kind=nag_wp), Parameter, Public :: half = 0.5_nag_wp
Integer, Parameter, Public       :: itrace = 0, nin = 5, nout = 6,       &
npde = 3, nv = 0, nxi = 0
!     .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: el0, er0, gamma, rl0, rr0, ul0, ur0
Contains
Subroutine bndary(npde,npts,t,x,u,nv,v,vdot,ibnd,g,ires)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (In)           :: ibnd, npde, npts, nv
Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g(npde)
Real (Kind=nag_wp), Intent (In) :: u(npde,npts), v(nv), vdot(nv),      &
x(npts)
!       .. Executable Statements ..
If (ibnd==0) Then
g(1) = u(1,1) - rl0
g(2) = u(2,1) - ul0
g(3) = u(3,1) - el0
Else
g(1) = u(1,npts) - rr0
g(2) = u(2,npts) - ur0
g(3) = u(3,npts) - er0
End If
Return
End Subroutine bndary
Subroutine numflx(npde,t,x,nv,v,uleft,uright,flux,ires)

!       .. Use Statements ..
Use nag_library, Only: d03pxf
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t, x
Integer, Intent (Inout)        :: ires
Integer, Intent (In)           :: npde, nv
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: flux(npde)
Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde), v(nv)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: tol
Integer                        :: ifail, niter
!       .. Executable Statements ..
tol = 0.0_nag_wp
niter = 0

ifail = 0
Call d03pxf(uleft,uright,gamma,tol,niter,flux,ifail)

Return
End Subroutine numflx
End Module d03pxfe_mod
Program d03pxfe

!     D03PXF Example Main Program

!     .. Use Statements ..
Use d03pxfe_mod, Only: alpha_l, alpha_r, beta_l, beta_r, bndary, el0,    &
er0, gamma, half, itrace, nin, nout, npde,        &
numflx, nv, nxi, rl0, rr0, ul0, ur0
Use nag_library, Only: d03pek, d03plf, d03plp, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: d, p, tout, ts, v
Integer                          :: i, ifail, ind, itask, itol, k,       &
lenode, mlu, neqn, niw, npts, nw,    &
nwkres
Character (1)                    :: laopt, norm
!     .. Local Arrays ..
Real (Kind=nag_wp)               :: algopt(30), atol(1), rtol(1),        &
ue(3,9), xi(1)
Real (Kind=nag_wp), Allocatable  :: u(:,:), w(:), x(:)
Integer, Allocatable             :: iw(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: real
!     .. Executable Statements ..
Write (nout,*) 'D03PXF Example Program Results'
!     Skip heading in data file

nwkres = npde*(2*npts+3*npde+32) + 7*npts + 4
mlu = 3*npde - 1
neqn = npde*npts + nv
niw = neqn + 24
lenode = 9*neqn + 50
nw = (3*mlu+1)*neqn + nwkres + lenode
Allocate (u(npde,npts),w(nw),x(npts),iw(niw))

Read (nin,*) gamma, rl0, rr0, ul0, ur0

el0 = alpha_l/(gamma-1.0_nag_wp) + half*rl0*beta_l**2
er0 = alpha_r/(gamma-1.0_nag_wp) + half*rr0*beta_r**2

!     Initialize mesh
Do i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
End Do
xi(1) = 0.0_nag_wp

!     Initial values
Do i = 1, npts
If (x(i)<half) Then
u(1,i) = rl0
u(2,i) = ul0
u(3,i) = el0
Else If (x(i)==half) Then
u(1,i) = half*(rl0+rr0)
u(2,i) = half*(ul0+ur0)
u(3,i) = half*(el0+er0)
Else
u(1,i) = rr0
u(2,i) = ur0
u(3,i) = er0
End If
End Do

ind = 0
algopt(1:30) = 0.0_nag_wp

!     Theta integration
algopt(1) = 2.0_nag_wp
algopt(6) = 2.0_nag_wp
algopt(7) = 2.0_nag_wp

!     Max. time step
algopt(13) = 0.5E-2_nag_wp

ts = 0.0_nag_wp
tout = 0.035_nag_wp

!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03plf(npde,ts,tout,d03plp,numflx,bndary,u,npts,x,nv,d03pek,nxi,xi, &
ifail)

Write (nout,99998) ts
Write (nout,99999)

!     Read exact data at output points

Do i = 1, 9
End Do

!     Calculate density, velocity and pressure

k = 0
Do i = 15, npts - 14, 14
d = u(1,i)
v = u(2,i)/d
p = d*(gamma-1.0_nag_wp)*(u(3,i)/d-half*v**2)
k = k + 1
Write (nout,99996) x(i), d, ue(1,k), v, ue(2,k), p, ue(3,k)
End Do

Write (nout,99997) iw(1), iw(2), iw(3), iw(5)

99999 Format (4X,'X',7X,'APPROX D',3X,'EXACT D',4X,'APPROX V',3X,'EXAC','T V', &
4X,'APPROX P',3X,'EXACT P')
99998 Format (/,' T = ',F6.3,/)
99997 Format (/,' Number of integration steps in time = ',I6,/,' Number ',     &
'of function evaluations = ',I6,/,' Number of Jacobian ',              &
'evaluations =',I6,/,' Number of iterations = ',I6)
99996 Format (1X,E9.2,6(E11.4))
End Program d03pxfe
```