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

    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
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(t,u,npde,x,npts)
!         Exact solution (for comparison and b.c. purposes)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,ux,ncode,v,vdot,p,c,d,s,ires)

!         .. 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)    :: 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
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndry1
       SUBROUTINE bndry1(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)
!         .. 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
!dec$ attributes offload : mic :: bndry1
       END SUBROUTINE bndry1
!dec$ attributes offload : mic :: nmflx1
       SUBROUTINE nmflx1(npde,t,x,ncode,v,uleft,uright,flux,ires)

!         .. 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 ..
          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
!dec$ attributes offload : mic :: nmflx1
       END SUBROUTINE nmflx1
!dec$ attributes offload : mic :: odedef
       SUBROUTINE odedef(npde,t,ncode,v,vdot,nxi,xi,ucp,ucpx,ucpt,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: odedef
       END SUBROUTINE odedef
!dec$ attributes offload : mic :: bndry2
       SUBROUTINE bndry2(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)
             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
!dec$ attributes offload : mic :: bndry2
       END SUBROUTINE bndry2
!dec$ attributes offload : mic :: nmflx2
       SUBROUTINE nmflx2(npde,t,x,ncode,v,uleft,uright,flux,ires)

!         .. Use Statements ..
          USE nag_library, ONLY : d03puf, d03pvf
!         .. 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
          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
!dec$ attributes offload : mic :: nmflx2
       END SUBROUTINE nmflx2
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,x,u)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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
!dec$ attributes offload : mic :: uvinit
       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
!dec$ attributes offload : mic :: ex1
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d03plf
          USE d03plfe_mod, ONLY : bndry1, exact, itrace, nag_wp, ncode1, nin,  &
                                  nmflx1, npde1, nxi1, odedef, pdedef
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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)
!dec$ attributes offload : mic :: ex1
       END SUBROUTINE ex1
!dec$ attributes offload : mic :: ex2
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d03pek, d03plf, d03plp
          USE d03plfe_mod, ONLY : bndry2, el0, er0, gamma, itrace, nag_wp,     &
                                  ncode2, nin, nmflx2, npde2, nxi2, rl0, rr0,  &
                                  uvinit
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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))
!dec$ attributes offload : mic :: ex2
       END SUBROUTINE ex2
    END PROGRAM d03plfe
