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

    MODULE d03pdae_mod

!      D03PDA 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

!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,nptl,u,ux,p,q,r,ires,iuser,ruser)

!         .. 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 (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,nptl), ux(npde,nptl),  &
                                                 x(nptl)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. 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: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
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef

!dec$ attributes offload : mic :: bndary
       SUBROUTINE bndary(npde,t,u,ux,ibnd,beta,gamma,ires,iuser,ruser)

!         .. 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 (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. 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.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
!dec$ attributes offload : mic :: bndary
       END SUBROUTINE bndary

!dec$ attributes offload : mic :: uinit
       SUBROUTINE uinit(npde,npts,x,u,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: npde, npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npde,npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: piby2
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              sin
!         .. Executable Statements ..
          piby2 = ruser(1)
          DO i = 1, npts
             u(1,i) = -sin(piby2*x(i))
             u(2,i) = -piby2*piby2*u(1,i)
          END DO
          RETURN
!dec$ attributes offload : mic :: uinit
       END SUBROUTINE uinit
    END MODULE d03pdae_mod

    PROGRAM d03pdae

!      D03PDA Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pda, d03pyf, nag_wp, x01aaf
       USE d03pdae_mod, ONLY : bndary, nin, nout, npde, pdedef, uinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: acc, dx, pi, piby2, tout, ts
       INTEGER                             :: i, ifail, ind, intpts, it,       &
                                              itask, itrace, itype, lenode, m, &
                                              mu, nbkpts, nel, neqn, niw,      &
                                              npl1, npoly, npts, nw, nwkres
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: ruser(1), rwsav(1100)
       REAL (KIND=nag_wp), ALLOCATABLE     :: u(:,:), uout(:,:,:), w(:), x(:), &
                                              xbkpts(:), xout(:)
       INTEGER                             :: iuser(1), iwsav(505)
       INTEGER, ALLOCATABLE                :: iw(:)
       LOGICAL                             :: lwsav(100)
       CHARACTER (80)                      :: cwsav(10)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PDA 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
       niw = neqn + 24
       npl1 = npoly + 1
       nwkres = 3*npl1*npl1 + npl1*(npde*npde+6*npde+nbkpts+1) + 13*npde + 5
       lenode = (3*mu+1)*neqn
       nw = 11*neqn + 50 + nwkres + lenode

       ALLOCATE (u(npde,npts),uout(npde,intpts,itype),w(nw),x(npts), &
          xbkpts(nbkpts),xout(intpts),iw(niw))

       READ (nin,*) xout(1:intpts)
       READ (nin,*) acc
       READ (nin,*) m, itrace
       piby2 = 0.5_nag_wp*x01aaf(pi)
       ruser(1) = piby2

!      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 d03pda(npde,m,ts,tout,pdedef,bndary,u,nbkpts,xbkpts,npoly,npts, &
             x,uinit,acc,w,nw,iw,niw,itask,itrace,ind,iuser,ruser,cwsav,lwsav, &
             iwsav,rwsav,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, &
             w,nw,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) iw(1), iw(2), iw(3), iw(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 d03pdae
