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

    Module d03pffe_mod

!     D03PFF 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                               :: bndry1, bndry2, exact, nmflx1,   &
                                              nmflx2, pdedef
!     .. Parameters ..
      Integer, Parameter, Public           :: nin = 5, nout = 6, npde1 = 2,    &
                                              npde2 = 1
    Contains
      Subroutine exact(t,u,npde,x,npts)
!       Exact solution (for comparison and b.c. purposes)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
        Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
!       .. 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(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: pi, px1, px2, x1, x2
        Integer                              :: i
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp, sin
!       .. Executable Statements ..
        pi = x01aaf(pi)

        Do i = 1, npts
          x1 = x(i) + t
          x2 = x(i) - 3.0_nag_wp*t
          px1 = half*sin(two*pi*x1**2)
          px2 = half*sin(two*pi*x2**2)
          u(1,i) = half*(exp(x2)+exp(x1)+px2-px1) - t*(x1+x2)
          u(2,i) = (exp(x2)-exp(x1)+px2+px1) + x1*x2 + 8.0_nag_wp*t**2
        End Do
        Return
      End Subroutine exact
      Subroutine bndry1(npde,npts,t,x,u,ibnd,g,ires)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
        Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: ibnd, npde, npts
        Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: g(npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npde,3), x(npts)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: c, exu1, exu2
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: ue(2,1)
!       .. Executable Statements ..
        If (ibnd==0) Then
          Call exact(t,ue,npde,x(1),1)
          c = (x(2)-x(1))/(x(3)-x(2))
          exu1 = (one+c)*u(1,2) - c*u(1,3)
          exu2 = (one+c)*u(2,2) - c*u(2,3)
          g(1) = two*u(1,1) + u(2,1) - two*ue(1,1) - ue(2,1)
          g(2) = two*u(1,1) - u(2,1) - two*exu1 + exu2
        Else
          Call exact(t,ue,npde,x(npts),1)
          c = (x(npts)-x(npts-1))/(x(npts-1)-x(npts-2))
          exu1 = (one+c)*u(1,2) - c*u(1,3)
          exu2 = (one+c)*u(2,2) - c*u(2,3)
          g(1) = two*u(1,1) - u(2,1) - two*ue(1,1) + ue(2,1)
          g(2) = two*u(1,1) + u(2,1) - two*exu1 - exu2
        End If
        Return
      End Subroutine bndry1
      Subroutine nmflx1(npde,t,x,uleft,uright,flux,ires)

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t, x
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: npde
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: flux(npde)
        Real (Kind=nag_wp), Intent (In)      :: uleft(npde), uright(npde)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: ltmp, rtmp
!       .. Executable Statements ..
        ltmp = 3.0_nag_wp*uleft(1) + 1.5_nag_wp*uleft(2)
        rtmp = uright(1) - half*uright(2)
        flux(1) = half*(ltmp-rtmp)
        flux(2) = ltmp + rtmp
        Return
      End Subroutine nmflx1
      Subroutine pdedef(npde,t,x,u,ux,p,c,d,s,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t, x
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: 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)
!       .. Executable Statements ..
        p(1,1) = 1.0_nag_wp
        c(1) = 0.01_nag_wp
        d(1) = ux(1)
        s(1) = u(1)
        Return
      End Subroutine pdedef
      Subroutine bndry2(npde,npts,t,x,u,ibnd,g,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: ibnd, npde, npts
        Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: g(npde)
        Real (Kind=nag_wp), Intent (In)      :: u(npde,3), x(npts)
!       .. Executable Statements ..
        If (ibnd==0) Then
          g(1) = u(1,1) - 3.0_nag_wp
        Else
          g(1) = u(1,1) - 5.0_nag_wp
        End If
        Return
      End Subroutine bndry2
      Subroutine nmflx2(npde,t,x,uleft,uright,flux,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t, x
        Integer, Intent (Inout)              :: ires
        Integer, Intent (In)                 :: npde
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: flux(npde)
        Real (Kind=nag_wp), Intent (In)      :: uleft(npde), uright(npde)
!       .. Executable Statements ..
        If (x>=0.0E0_nag_wp) Then
          flux(1) = x*uleft(1)
        Else
          flux(1) = x*uright(1)
        End If
        Return
      End Subroutine nmflx2
    End Module d03pffe_mod
    Program d03pffe

!     D03PFF Example Main Program

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

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use nag_library, Only: d03pff, d03pfp, nag_wp
        Use d03pffe_mod, Only: bndry1, exact, nin, nmflx1, npde1
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: dx, tout, ts, tsmax
        Integer                              :: i, ifail, inc, ind, it, itask, &
                                                itrace, lisave, lrsave,        &
                                                nfuncs, niters, njacs, nop,    &
                                                npde, npts, nsteps, outpts
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: acc(2)
        Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:,:), ue(:,:),     &
                                                x(:), 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
        npde = npde1
        Read (nin,*)
        Read (nin,*) npts, inc, outpts
        lisave = 24 + npde*npts
        lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

        Allocate (rsave(lrsave),u(npde,npts),ue(npde,outpts),x(npts), &
          xout(outpts),isave(lisave))
        Read (nin,*) acc(1:2)
        Read (nin,*) itrace
        Read (nin,*) tsmax

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

!       Set initial values
        Read (nin,*) ts
        Call exact(ts,u,npde,x,npts)

        ind = 0
        itask = 1

        Write (nout,99992) npts
        Write (nout,99991) acc(1)
        Write (nout,99990) acc(2)
        Write (nout,99999)

        Do it = 1, 2
          tout = 0.1E0_nag_wp*real(it,kind=nag_wp)

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          Call d03pff(npde,ts,tout,d03pfp,nmflx1,bndry1,u,npts,x,acc,tsmax, &
            rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

!         Set output points

          nop = 0
          Do i = 1, npts, inc
            nop = nop + 1
            xout(nop) = x(i)
          End Do

!         Check against exact solution

          Call exact(tout,ue,npde,xout,nop)
          nop = 1
          i = 1
          Write (nout,99998) ts, xout(nop), u(1,i), ue(1,nop), u(2,i), &
            ue(2,nop)
          Do i = inc + 1, npts, inc
            nop = nop + 1
            Write (nout,99997) xout(nop), u(1,i), ue(1,nop), u(2,i), ue(2,nop)
          End Do
        End Do

!       Print integration statistics (reasonably rounded)
        nsteps = 5*((isave(1)+2)/5)
        nfuncs = 50*((isave(2)+25)/50)
        njacs = 10*((isave(3)+5)/10)
        niters = 50*((isave(5)+25)/50)
        Write (nout,99996) nsteps
        Write (nout,99995) nfuncs
        Write (nout,99994) njacs
        Write (nout,99993) niters

        Return

99999   Format (/5X,'t',9X,'x',8X,'Approx u',4X,'Exact u',5X,'Approx v',4X, &
          'Exact v')
99998   Format (/1X,F6.3,5(3X,F9.4))
99997   Format (7X,5(3X,F9.4))
99996   Format (/' Number of time steps           (nearest 5)  = ',I6)
99995   Format (' Number of function evaluations (nearest 50) = ',I6)
99994   Format (' Number of Jacobian evaluations (neaerst 10) = ',I6)
99993   Format (' Number of iterations           (nearest 50) = ',I6)
99992   Format (/' Number of mesh points used = ',I4)
99991   Format (' Relative tolerance used    = ',E10.3)
99990   Format (' Absolute tolerance used    = ',E10.3)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use nag_library, Only: d03pff, nag_wp
        Use d03pffe_mod, Only: bndry2, nin, nmflx2, npde2, pdedef
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: dx, tout, ts, tsmax
        Integer                              :: i, ifail, ind, it, itask,      &
                                                itrace, lisave, lrsave,        &
                                                nfuncs, niters, njacs, npde,   &
                                                npts, nsteps, outpts
!       .. Local Arrays ..
        Real (Kind=nag_wp)                   :: acc(2)
        Real (Kind=nag_wp), Allocatable      :: rsave(:), u(:,:), x(:), xout(:)
        Integer, Allocatable                 :: iout(:), isave(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: real
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 2'
        Write (nout,*)
        npde = npde2
        Read (nin,*)
        Read (nin,*) npts, outpts
        lisave = 24 + npde*npts
        lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

        Allocate (rsave(lrsave),u(npde,npts),x(npts),xout(outpts), &
          isave(lisave),iout(outpts))
        Read (nin,*) iout(1:outpts)
        Read (nin,*) acc(1:2)
        Read (nin,*) itrace
        Read (nin,*) tsmax

!       Initialise mesh
        dx = 2.0_nag_wp/real(npts-1,kind=nag_wp)
        Do i = 1, npts
          x(i) = -1.0_nag_wp + real(i-1,kind=nag_wp)*dx
        End Do

!       Set initial values
        u(1,1:npts) = x(1:npts) + 4.0_nag_wp

        ind = 0
        itask = 1

!       Set output points from inout indices
        Do i = 1, outpts
          xout(i) = x(iout(i))
        End Do

!       Two output value of t: tout read from file and 10.0

        Read (nin,*) ts, tout

        Write (nout,99995) npts
        Write (nout,99994) acc(1)
        Write (nout,99993) acc(2)
        Write (nout,99992) xout(1:outpts)

        Do it = 1, 2

          ifail = 0
          Call d03pff(npde,ts,tout,pdedef,nmflx2,bndry2,u,npts,x,acc,tsmax, &
            rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

          Write (nout,99991) ts, (u(1,iout(i)),i=1,outpts)
          tout = 10.0_nag_wp

        End Do

!       Print integration statistics (reasonably rounded)
        nsteps = 5*((isave(1)+2)/5)
        nfuncs = 50*((isave(2)+25)/50)
        njacs = 10*((isave(3)+5)/10)
        niters = 50*((isave(5)+25)/50)
        Write (nout,99999) nsteps
        Write (nout,99998) nfuncs
        Write (nout,99997) njacs
        Write (nout,99996) niters

        Return

99999   Format (/' Number of time steps           (nearest 5)  = ',I6)
99998   Format (' Number of function evaluations (nearest 50) = ',I6)
99997   Format (' Number of Jacobian evaluations (neaerst 10) = ',I6)
99996   Format (' Number of iterations           (nearest 50) = ',I6)
99995   Format (/' Number of mesh points used = ',I4)
99994   Format (' Relative tolerance used    = ',E10.3)
99993   Format (' Absolute tolerance used    = ',E10.3)
99992   Format (/1X,'x',10X,7F9.4)
99991   Format (1X,'u(t=',F6.3,')',7F9.4)
      End Subroutine ex2
    End Program d03pffe