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

    MODULE d03pdfe_mod

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

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6, npde = 2
    CONTAINS
       SUBROUTINE uinit(npde,npts,x,u)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. 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 ..
          REAL (KIND=nag_wp)                  :: piby2
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              sin
!         .. Executable Statements ..
          piby2 = 0.5_nag_wp*x01aaf(piby2)
          DO i = 1, npts
             u(1,i) = -sin(piby2*x(i))
             u(2,i) = -piby2*piby2*u(1,i)
          END DO
          RETURN
       END SUBROUTINE uinit

       SUBROUTINE pdedef(npde,t,x,nptl,u,ux,p,q,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: npde, nptl
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: p(npde,npde,nptl),            &
                                                 q(npde,nptl), r(npde,nptl)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,nptl), ux(npde,nptl),  &
                                                 x(nptl)
!         .. Local Scalars ..
          INTEGER                             :: i
!         .. Executable Statements ..
          DO i = 1, nptl
             q(1,i) = u(2,i)
             q(2,i) = u(1,i)*ux(2,i) - ux(1,i)*u(2,i)
             r(1,i) = ux(1,i)
             r(2,i) = ux(2,i)
             p(1,1,i) = 0.0_nag_wp
             p(1,2,i) = 0.0_nag_wp
             p(2,1,i) = 0.0_nag_wp
             p(2,2,i) = 1.0_nag_wp
          END DO
          RETURN
       END SUBROUTINE pdedef

       SUBROUTINE bndary(npde,t,u,ux,ibnd,beta,gamma,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, npde
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: beta(npde), gamma(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde)
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             beta(1) = 1.0_nag_wp
             gamma(1) = 0.0_nag_wp
             beta(2) = 0.0_nag_wp
             gamma(2) = u(1) - 1.0_nag_wp
          ELSE
             beta(1) = 1.0E+0_nag_wp
             gamma(1) = 0.0_nag_wp
             beta(2) = 0.0_nag_wp
             gamma(2) = u(1) + 1.0_nag_wp
          END IF
          RETURN
       END SUBROUTINE bndary
    END MODULE d03pdfe_mod

    PROGRAM d03pdfe

!      D03PDF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pdf, d03pyf, nag_wp
       USE d03pdfe_mod, ONLY : bndary, nin, nout, npde, pdedef, uinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: acc, dx, tout, ts
       INTEGER                             :: i, ifail, ind, intpts, it,       &
                                              itask, itrace, itype, lenode,    &
                                              lisave, lrsave, m, mu, nbkpts,   &
                                              nel, neqn, npl1, npoly, npts,    &
                                              nwkres
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: rsave(:), u(:,:), uout(:,:,:),   &
                                              x(:), xbkpts(:), xout(:)
       INTEGER, ALLOCATABLE                :: isave(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PDF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) intpts, nbkpts, npoly, itype

       nel = nbkpts - 1
       npts = nel*npoly + 1
       mu = npde*(npoly+1) - 1
       neqn = npde*npts
       lisave = neqn + 24
       npl1 = npoly + 1
       nwkres = 3*npl1*npl1 + npl1*(npde*npde+6*npde+nbkpts+1) + 13*npde + 5
       lenode = (3*mu+1)*neqn
       lrsave = 11*neqn + 50 + nwkres + lenode

       ALLOCATE (u(npde,npts),uout(npde,intpts,itype),rsave(lrsave),x(npts), &
          xbkpts(nbkpts),xout(intpts),isave(lisave))

       READ (nin,*) xout(1:intpts)
       READ (nin,*) acc
       READ (nin,*) m, itrace

!      Set the break-points

       dx = 2.0_nag_wp/real(nbkpts-1,kind=nag_wp)
       xbkpts(1) = -1.0_nag_wp
       DO i = 2, nbkpts - 1
          xbkpts(i) = xbkpts(i-1) + dx
       END DO
       xbkpts(nbkpts) = 1.0_nag_wp

       ind = 0
       itask = 1
       READ (nin,*) ts, tout

!      Loop over output values of t

       DO it = 1, 5
          tout = 10.0_nag_wp*tout

!            ifail: behaviour on error exit   
!                   =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          CALL d03pdf(npde,m,ts,tout,pdedef,bndary,u,nbkpts,xbkpts,npoly,npts, &
             x,uinit,acc,rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

          IF (it==1) THEN
             WRITE (nout,99999) npoly, nel
             WRITE (nout,99998) acc, npts
             WRITE (nout,99997) xout(1:6)
          END IF

!         Interpolate at required spatial points

          ifail = 0
          CALL d03pyf(npde,u,nbkpts,xbkpts,npoly,npts,xout,intpts,itype,uout, &
             rsave,lrsave,ifail)

          WRITE (nout,99996) ts, uout(1,1:intpts,1)
          WRITE (nout,99995) uout(2,1:intpts,1)
       END DO

!      Print integration statistics

       WRITE (nout,99994) isave(1), isave(2), isave(3), isave(5)


99999  FORMAT (' Polynomial degree =',I4,'   No. of elements = ',I4)
99998  FORMAT (' Accuracy requirement =',E10.3,'  Number of points = ',I5/)
99997  FORMAT ('  T /   X    ',6F8.4/)
99996  FORMAT (1X,F7.4,' U(1)',6F8.4)
99995  FORMAT (9X,'U(2)',6F8.4/)
99994  FORMAT (' Number of integration steps in time                   ', &
          I4/' Number of residual evaluations of resulting ODE system', &
          I4/' Number of Jacobian evaluations                        ', &
          I4/' Number of iterations of nonlinear solver              ',I4)
    END PROGRAM d03pdfe