!   D03PXF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

    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, ncode = 0, nin = 5,  &
                                              nout = 6, npde = 3, 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,ncode,v,vdot,ibnd,g,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: ibnd, ncode, npde, npts
        Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: g(npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npde,npts), v(ncode),        &
                                                vdot(ncode), 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,ncode,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)                 :: ncode, npde
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: flux(npde)
        Real (Kind=nag_wp), Intent (In)      :: uleft(npde), uright(npde),     &
                                                v(ncode)
!       .. 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 nag_library, Only: d03pek, d03plf, d03plp, nag_wp
      Use d03pxfe_mod, Only: alpha_l, alpha_r, beta_l, beta_r, bndary, el0,    &
                             er0, gamma, half, itrace, ncode, nin, nout, npde, &
                             numflx, nxi, rl0, rr0, ul0, ur0
!     .. 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 + ncode
      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

!     Initialise 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,ncode,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