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

    Module d03ppae_mod

!     D03PPA 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        :: four = 4.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: half = 0.5_nag_wp
      Real (Kind=nag_wp), Parameter        :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: ptone = 0.1_nag_wp
      Real (Kind=nag_wp), Parameter        :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: zero = 0.0_nag_wp
      Integer, Parameter                   :: itrace = 0, m = 0, ncode = 0,    &
                                              nin = 5, nout = 6, npde = 1,     &
                                              nxfix = 0, nxi = 0
    Contains
      Subroutine uvinit(npde,npts,nxi,x,xi,u,ncode,v,iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: ncode, npde, npts, nxi
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (Out)     :: u(npde,npts), v(ncode)
        Real (Kind=nag_wp), Intent (In)      :: x(npts), xi(nxi)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: a, b, c, e, t
        Integer                              :: i
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        e = ruser(1)
        t = zero
        Do i = 1, npts
          a = (x(i)-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x(i)-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x(i)-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(1,i) = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x(i)+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(1,i) = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            u(1,i) = (one+half*a+ptone*b)/(one+a+b)
          End If
        End Do
        Return
      End Subroutine uvinit
      Subroutine pdedef(npde,t,x,u,ux,ncode,v,vdot,p,q,r,ires,iuser,ruser)

!       .. 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)     :: p(npde,npde), q(npde), r(npde)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: u(npde), ux(npde), v(ncode),   &
                                                vdot(ncode)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: e
!       .. Executable Statements ..
        e = ruser(1)
        p(1,1) = one
        r(1) = e*ux(1)
        q(1) = u(1)*ux(1)
        Return
      End Subroutine pdedef
      Subroutine bndary(npde,t,u,ux,ncode,v,vdot,ibnd,beta,gamma,ires,iuser, &
        ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: ibnd, ncode, npde
        Integer, Intent (Inout)              :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: beta(npde), gamma(npde)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: u(npde), ux(npde), v(ncode),   &
                                                vdot(ncode)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: a, b, c, e, ue, x
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        e = ruser(1)
        beta(1) = zero
        If (ibnd==0) Then
          x = zero
          a = (x-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            ue = (one+half*a+ptone*b)/(one+a+b)
          End If
        Else
          x = one
          a = (x-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            ue = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            ue = (one+half*a+ptone*b)/(one+a+b)
          End If
        End If
        gamma(1) = u(1) - ue
        Return
      End Subroutine bndary
      Subroutine monitf(t,npts,npde,x,u,r,fmon,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npde, npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: fmon(npts)
        Real (Kind=nag_wp), Intent (In)      :: r(npde,npts), u(npde,npts),    &
                                                x(npts)
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: drdx, h
        Integer                              :: i, k, l
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs, max, min
!       .. Executable Statements ..
        Do i = 1, npts - 1
          k = max(1,i-1)
          l = min(npts,i+1)
          h = (x(l)-x(k))*half
!         Second derivative ..
          drdx = (r(1,i+1)-r(1,i))/h
          fmon(i) = abs(drdx)
        End Do
        fmon(npts) = fmon(npts-1)
        Return
      End Subroutine monitf
      Subroutine exact(t,x,npts,u,iuser,ruser)
!       Exact solution (for comparison purposes)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: npts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (Out)     :: u(npts)
        Real (Kind=nag_wp), Intent (In)      :: x(npts)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: a, b, c, e
        Integer                              :: i
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        e = ruser(1)
        Do i = 1, npts
          a = (x(i)-0.25_nag_wp-0.75_nag_wp*t)/(four*e)
          b = (0.9_nag_wp*x(i)-0.325_nag_wp-0.495_nag_wp*t)/(two*e)
          If (a>zero .And. a>b) Then
            a = exp(-a)
            c = (0.8_nag_wp*x(i)-0.4_nag_wp-0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(i) = (half+ptone*c+a)/(one+c+a)
          Else If (b>zero .And. b>=a) Then
            b = exp(-b)
            c = (-0.8_nag_wp*x(i)+0.4_nag_wp+0.24_nag_wp*t)/(four*e)
            c = exp(c)
            u(i) = (ptone+half*c+b)/(one+c+b)
          Else
            a = exp(a)
            b = exp(b)
            u(i) = (one+half*a+ptone*b)/(one+a+b)
          End If
        End Do
        Return
      End Subroutine exact
    End Module d03ppae_mod
    Program d03ppae

!     D03PPA Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d03ppa, d03pzf, d53pck, nag_wp
      Use d03ppae_mod, Only: bndary, exact, half, itrace, m, monitf, ncode,    &
                             nin, nout, npde, nxfix, nxi, pdedef, two, uvinit, &
                             zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: con, dx, dxmesh, e, tout,        &
                                              trmesh, ts, x0, xmid, xratio
      Integer                              :: i, ifail, ind, intpts, ipminf,   &
                                              it, itask, itol, itype, lenode,  &
                                              neqn, niw, npts, nrmesh, nw,     &
                                              nwkres
      Logical                              :: remesh, theta
      Character (1)                        :: laopt, norm
!     .. Local Arrays ..
      Real (Kind=nag_wp)                   :: algopt(30), atol(1), rtol(1),    &
                                              ruser(1), rwsav(1100), xfix(1),  &
                                              xi(1)
      Real (Kind=nag_wp), Allocatable      :: u(:), ue(:), uout(:,:,:), w(:),  &
                                              x(:), xout(:)
      Integer                              :: iuser(1), iwsav(505)
      Integer, Allocatable                 :: iw(:)
      Logical                              :: lwsav(100)
      Character (80)                       :: cwsav(10)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: min, real
!     .. Executable Statements ..
      Write (nout,*) 'D03PPA Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) npts, intpts, itype
      neqn = npde*npts + ncode
      niw = 25 + nxfix
      nwkres = npde*(npts+3*npde+21) + 7*npts + nxfix + 3
      lenode = 11*neqn + 50
      nw = neqn*neqn + neqn + nwkres + lenode

      Allocate (u(neqn),ue(intpts),uout(npde,intpts,itype),w(nw),x(npts), &
        xout(intpts),iw(niw))
      Read (nin,*) itol
      Read (nin,*) atol(1), rtol(1)
      Read (nin,*) e
      ruser(1) = e

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

!     Set remesh parameters
      remesh = .True.
      nrmesh = 3
      dxmesh = half
      con = two/real(npts-1,kind=nag_wp)
      xratio = 1.5_nag_wp
      ipminf = 0

      xi(1) = zero
      norm = 'A'
      laopt = 'F'
      ind = 0
      itask = 1

!     Set theta to .TRUE. if the Theta integrator is required
      theta = .False.
      algopt(1:30) = zero
      If (theta) Then
        algopt(1) = two
      End If

!     Loop over output value of t
      ts = zero
      tout = zero
      Do it = 1, 5
        xmid = half + half*tout
        tout = 0.2_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 d03ppa(npde,m,ts,tout,pdedef,bndary,uvinit,u,npts,x,ncode,d53pck, &
          nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix,xfix, &
          nrmesh,dxmesh,trmesh,ipminf,xratio,con,monitf,w,nw,iw,niw,itask, &
          itrace,ind,iuser,ruser,cwsav,lwsav,iwsav,rwsav,ifail)

        If (it==1) Then
          Write (nout,99998) atol, npts
          Write (nout,99993) nrmesh
          Write (nout,99992) e
          Write (nout,*)
        End If

!       Set output points
        dx = 0.1_nag_wp
        If (tout>half) dx = 0.05_nag_wp

        x0 = xmid - half*real(intpts-1,kind=nag_wp)*dx
        Do i = 1, intpts
          xout(i) = x0
          x0 = x0 + dx
        End Do
        xout(intpts) = min(xout(intpts),x(npts))

        Write (nout,99999) ts
        Write (nout,99996) xout(1:intpts)

!       Interpolate at output points
        ifail = 0
        Call d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

!       Check against exact solution

        Call exact(ts,xout,intpts,ue,iuser,ruser)

        Write (nout,99995) uout(1,1:intpts,1)
        Write (nout,99994) ue(1:intpts)

      End Do
      Write (nout,99997) iw(1), iw(2), iw(3), iw(5)

99999 Format (' T = ',F6.3)
99998 Format (//'  Accuracy requirement =',E10.3,' Number of points = ',I3/)
99997 Format (' Number of integration steps in time = ',I6/' Number o', &
        'f function evaluations = ',I6/' Number of Jacobian eval','uations =', &
        I6/' Number of iterations = ',I6)
99996 Format (1X,'X           ',5F9.4)
99995 Format (1X,'Approx sol. ',5F9.4)
99994 Format (1X,'Exact  sol. ',5F9.4/)
99993 Format (2X,'Remeshing every',I3,' time steps'/)
99992 Format (2X,'E =',F8.3)
    End Program d03ppae
