!   D03PLF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d03plfe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: itrace = 0, ncode1 = 2,          &
                                              ncode2 = 0, nin = 5, nout = 6,   &
                                              npde1 = 2, npde2 = 3, nxi1 = 2,  &
                                              nxi2 = 0
!     .. Local Scalars ..
      Real (Kind=nag_wp), Save             :: el0, er0, gamma, rl0, rr0
    Contains
      Subroutine exact(t,u,npde,x,npts)
!       Exact solution (for comparison and b.c. purposes)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: u(npde,npts)
        Real (Kind=nag_wp), Intent (In)      :: x(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: f, g, pi, pi2, x1, x3
        Integer                              :: i
!       .. Intrinsic Procedures ..
        Intrinsic                            :: cos, exp, sin
!       .. Executable Statements ..
        f = 0.0_nag_wp
        pi = x01aaf(f)
        pi2 = 2.0_nag_wp*pi
        Do i = 1, npts
          x1 = x(i) + t
          x3 = x(i) - 3.0_nag_wp*t
          f = exp(pi*x3)*sin(pi2*x3)
          g = exp(-pi2*x1)*cos(pi2*x1)
          u(1,i) = f + g
          u(2,i) = f - g
        End Do
        Return
      End Subroutine exact
      Subroutine pdedef(npde,t,x,u,ux,ncode,v,vdot,p,c,d,s,ires)

!       .. 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)     :: c(npde), d(npde),              &
                                                p(npde,npde), s(npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npde), ux(npde), v(ncode),   &
                                                vdot(ncode)
!       .. Local Scalars ..
        Integer                              :: i
!       .. Executable Statements ..
        c(1:npde) = 1.0_nag_wp
        d(1:npde) = 0.0_nag_wp
        s(1:npde) = 0.0_nag_wp
        p(1:npde,1:npde) = 0.0_nag_wp
        Do i = 1, npde
          p(i,i) = 1.0_nag_wp
        End Do
        Return
      End Subroutine pdedef
      Subroutine bndry1(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)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: dudx
        Integer                              :: i
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: ue(2,1)
!       .. Executable Statements ..
        If (ibnd==0) Then
          i = 1
          Call exact(t,ue,npde,x(i),1)
          g(1) = u(1,i) + u(2,i) - ue(1,1) - ue(2,1)
          dudx = (u(1,i+1)-u(2,i+1)-u(1,i)+u(2,i))/(x(i+1)-x(i))
          g(2) = vdot(1) - dudx
        Else
          i = npts
          Call exact(t,ue,npde,x(i),1)
          g(1) = u(1,i) - u(2,i) - ue(1,1) + ue(2,1)
          dudx = (u(1,i)+u(2,i)-u(1,i-1)-u(2,i-1))/(x(i)-x(i-1))
          g(2) = vdot(2) + 3.0_nag_wp*dudx
        End If
        Return
      End Subroutine bndry1
      Subroutine nmflx1(npde,t,x,ncode,v,uleft,uright,flux,ires)

!       .. 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)                   :: tmpl, tmpr
!       .. Executable Statements ..
        tmpl = 3.0_nag_wp*(uleft(1)+uleft(2))
        tmpr = uright(1) - uright(2)
        flux(1) = 0.5_nag_wp*(tmpl-tmpr)
        flux(2) = 0.5_nag_wp*(tmpl+tmpr)
        Return
      End Subroutine nmflx1
      Subroutine odedef(npde,t,ncode,v,vdot,nxi,xi,ucp,ucpx,ucpt,r,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: ncode, npde, nxi
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: r(ncode)
        Real (Kind=nag_wp), Intent (In)      :: ucp(npde,*), ucpt(npde,*),     &
                                                ucpx(npde,*), v(ncode),        &
                                                vdot(ncode), xi(nxi)
!       .. Executable Statements ..
        If (ires==-1) Then
          r(1) = 0.0_nag_wp
          r(2) = 0.0_nag_wp
        Else
          r(1) = v(1) - ucp(1,1) + ucp(2,1)
          r(2) = v(2) - ucp(1,2) - ucp(2,2)
        End If
        Return
      End Subroutine odedef
      Subroutine bndry2(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)
          g(3) = u(3,1) - el0
        Else
          g(1) = u(1,npts) - rr0
          g(2) = u(2,npts)
          g(3) = u(3,npts) - er0
        End If
        Return
      End Subroutine bndry2
      Subroutine nmflx2(npde,t,x,ncode,v,uleft,uright,flux,ires)

!       .. Use Statements ..
        Use nag_library, Only: d03puf, d03pvf
!       .. 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
        Character (1)                        :: path, solver
!       .. Executable Statements ..
        ifail = 0
        solver = 'R'
        If (solver=='R') Then
!         ROE scheme ..
          Call d03puf(uleft,uright,gamma,flux,ifail)
        Else
!         OSHER scheme ..
          path = 'P'
          Call d03pvf(uleft,uright,gamma,path,flux,ifail)
        End If
        Return
      End Subroutine nmflx2
      Subroutine uvinit(npde,npts,x,u)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: u(npde,npts)
        Real (Kind=nag_wp), Intent (In)      :: x(npts)
!       .. Local Scalars ..
        Integer                              :: i
!       .. Executable Statements ..
        Do i = 1, npts
          If (x(i)<0.5_nag_wp) Then
            u(1,i) = rl0
            u(2,i) = 0.0_nag_wp
            u(3,i) = el0
          Else If (x(i)==0.5_nag_wp) Then
            u(1,i) = 0.5_nag_wp*(rl0+rr0)
            u(2,i) = 0.0_nag_wp
            u(3,i) = 0.5_nag_wp*(el0+er0)
          Else
            u(1,i) = rr0
            u(2,i) = 0.0_nag_wp
            u(3,i) = er0
          End If
        End Do
        Return
      End Subroutine uvinit
    End Module d03plfe_mod
    Program d03plfe

!     D03PLF Example Main Program

!     .. Use Statements ..
      Use d03plfe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'D03PLF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use nag_library, Only: d03plf, nag_wp
        Use d03plfe_mod, Only: bndry1, exact, itrace, ncode1, nin, nmflx1,     &
                               npde1, nxi1, odedef, pdedef
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: tout, ts
        Integer                              :: i, ifail, ind, itask, itol,    &
                                                lenode, lisave, lrsave, ncode, &
                                                neqn, nop, npde, npts, nwkres, &
                                                nxi, outpts
        Character (1)                        :: laopt, norm
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: algopt(30), atol(1), rtol(1)
        Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:), ue(:,:),       &
                                                uout(:,:), x(:), xi(:), xout(:)
        Integer, Allocatable                 :: isave(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: real
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 1'
        Write (nout,*)
!       Skip heading in data file
        Read (nin,*)
        Read (nin,*) npts
        npde = npde1
        ncode = ncode1
        nxi = nxi1
        neqn = npde*npts + ncode
        lisave = 25*neqn + 24
        nwkres = npde*(2*npts+6*nxi+3*npde+26) + nxi + ncode + 7*npts + 2
        lenode = 11*neqn + 50
        lrsave = 4*neqn + 11*neqn/2 + 1 + nwkres + lenode
        lisave = lisave*4
        lrsave = lrsave*4
        outpts = npts/20 + 1
        Allocate (rsave(lrsave),u(neqn),ue(npde,outpts),uout(npde,outpts), &
          x(npts),xi(nxi),xout(outpts),isave(lisave))

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

!       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
        xi(2) = 1.0_nag_wp

!       Set initial values ..
        ts = 0.0_nag_wp
        Call exact(ts,u,npde,x,npts)
        u(neqn-1) = u(1) - u(2)
        u(neqn) = u(neqn-2) + u(neqn-3)

        laopt = 'S'
        ind = 0
        itask = 1

        algopt(1:30) = 0.0_nag_wp
!       Theta integration
        algopt(1) = 1.0_nag_wp
!       Sparse matrix algebra parameters
        algopt(29) = 0.1_nag_wp
        algopt(30) = 1.1_nag_wp

        tout = 0.5_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,pdedef,nmflx1,bndry1,u,npts,x,ncode,odedef, &
          nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
          lisave,itask,itrace,ind,ifail)

        Write (nout,99995) npts, atol, rtol

!       Set output points ..
        nop = 0
        Do i = 1, npts, 20
          nop = nop + 1
          xout(nop) = x(i)
        End Do

        Write (nout,99996) ts
        Write (nout,99999)

        nop = 0
        Do i = 1, npts, 20
          nop = nop + 1
          uout(1,nop) = u(npde*(i-1)+1)
          uout(2,nop) = u(npde*(i-1)+2)
        End Do

!       Check against exact solution ..
        Call exact(tout,ue,npde,xout,nop)
        Write (nout,99998)(xout(i),uout(1,i),ue(1,i),uout(2,i),ue(2,i),i=1, &
          nop)
        Write (nout,99997)
        Write (nout,99994) isave(1), isave(2), isave(3), isave(5)

        Return

99999   Format (8X,'X',8X,'Approx U1',3X,'Exact U1',4X,'Approx U2',3X, &
          'Exact U2'/)
99998   Format (5F12.4)
99997   Format (E11.4,4E14.4)
99996   Format (' T = ',F6.3)
99995   Format (/' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3/)
99994   Format (' Number of integration steps in time = ',I6/' Number ', &
          'of function evaluations = ',I6/' Number of Jacobian ', &
          'evaluations =',I6/' Number of iterations = ',I6)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use nag_library, Only: d03pek, d03plf, d03plp, nag_wp
        Use d03plfe_mod, Only: bndry2, el0, er0, gamma, itrace, ncode2, nin,   &
                               nmflx2, npde2, nxi2, rl0, rr0, uvinit
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: d, p, tout, ts, v
        Integer                              :: i, ifail, ind, it, itask,      &
                                                itol, k, lenode, lisave,       &
                                                lrsave, mlu, ncode, neqn,      &
                                                npde, npts, nskip, nwkres,     &
                                                nxi, outpts
        Character (1)                        :: laopt, norm
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: algopt(30), atol(1), rtol(1),  &
                                                xi(1)
        Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:,:), ue(:,:), x(:)
        Integer, Allocatable                 :: isave(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: real
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 2'
        Write (nout,*)
        Read (nin,*)
        Read (nin,*) npts, nskip
        npde = npde2
        ncode = ncode2
        nxi = nxi2
        Read (nin,*) el0, er0, gamma, rl0, rr0
        nwkres = npde*(2*npts+3*npde+32) + 7*npts + 4
        mlu = 3*npde - 1
        neqn = npde*npts + ncode
        lenode = 9*neqn + 50
        lisave = neqn + 24
        lrsave = (3*mlu+1)*neqn + nwkres + lenode
        outpts = (npts-2*nskip-1)/nskip
        Allocate (rsave(lrsave),u(npde,npts),x(npts),isave(lisave), &
          ue(npde,outpts))

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

        Write (nout,99995) gamma, el0, er0, rl0, rr0
        Write (nout,99997) npts, atol, rtol

!       Initialise mesh
        Do i = 1, npts
          x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
        End Do

!       Initial values of variables
        Call uvinit(npde,npts,x,u)

        xi(1) = 0.0_nag_wp
        laopt = 'B'
        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
        Write (nout,99999)
        Do it = 1, 2
          tout = real(it,kind=nag_wp)*0.1_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,nmflx2,bndry2,u,npts,x,ncode,d03pek, &
            nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
            lisave,itask,itrace,ind,ifail)

          Write (nout,99998) ts

!         Read exact data (to 4 d.p.) at output points ..
          Read (nin,*)
          Do i = 1, outpts
            Read (nin,*) ue(1,i), ue(2,i), ue(3,i)
          End Do

!         Calculate density, velocity and pressure ..
          k = 0
          Do i = 2*nskip + 1, npts - nskip, nskip
            d = u(1,i)
            v = u(2,i)/d
            p = d*(gamma-1.0_nag_wp)*(u(3,i)/d-0.5_nag_wp*v**2)
            k = k + 1
            Write (nout,99994) x(i), d, ue(1,k), v, ue(2,k), p, ue(3,k)
          End Do
        End Do

        Write (nout,99996) isave(1), isave(2), isave(3), isave(5)

        Return

99999   Format (4X,'X',5X,'APPROX D',1X,'EXACT D',2X,'APPROX V',1X,'EXAC', &
          'T V',2X,'APPROX P',1X,'EXACT P')
99998   Format (/' T = ',F6.3/)
99997   Format (/' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',E10.3/)
99996   Format (/' Number of integration steps in time = ',I6/' Number ', &
          'of function evaluations = ',I6/' Number of Jacobian ', &
          'evaluations =',I6/' Number of iterations = ',I6)
99995   Format (/' GAMMA =',F6.3,'  EL0 =',F6.3,'  ER0 =',F6.3,'  RL0 =',F6.3, &
          '  RR0 =',F6.3)
99994   Format (1X,F7.4,6(2X,F7.4))
      End Subroutine ex2
    End Program d03plfe
