Example description
!   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
      Read (nin,*)
      Read (nin,*) npts

      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

      Read (nin,*) itol
      Read (nin,*) norm
      Read (nin,*) atol(1), rtol(1)
      Read (nin,*) laopt

      ind = 0
      itask = 1
      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, &
        neqn,rtol,atol,itol,norm,laopt,algopt,w,nw,iw,niw,itask,itrace,ind,    &
        ifail)

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

!     Read exact data at output points

      Do i = 1, 9
        Read (nin,*) ue(1:3,i)
      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