!   D03PWF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d03pwfe_mod

!      D03PWF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: alpha_l = 460.894_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: alpha_r = 46.095_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: beta_l = 19.5975_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: beta_r = 6.19633_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: half = 0.5_nag_wp
       INTEGER, PARAMETER                  :: itrace = 0, ncode = 0, nin = 5,  &
                                              nout = 6, npde = 3, nxi = 0
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: el0, er0, gamma, rl0, rr0, ul0,  &
                                              ur0
    CONTAINS
       SUBROUTINE bndary(npde,npts,t,x,u,ncode,v,vdot,ibnd,g,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 : d03pwf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 ..
          INTEGER                             :: ifail
!         .. Executable Statements ..
          ifail = 0
          CALL d03pwf(uleft,uright,gamma,flux,ifail)
          RETURN
       END SUBROUTINE numflx
    END MODULE d03pwfe_mod
    PROGRAM d03pwfe

!      D03PWF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pek, d03plf, d03plp
       USE d03pwfe_mod, ONLY : alpha_l, alpha_r, beta_l, beta_r, bndary, el0,  &
                               er0, gamma, half, itrace, nag_wp, 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 Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PWF 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 d03pwfe