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

    MODULE d03pkfe_mod

!      D03PKF 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       :: one = 1.0_nag_wp
       INTEGER, PARAMETER                  :: itrace = 0, ncode = 1, nin = 5,  &
                                              nleft = 1, nout = 6, npde = 2,   &
                                              nxi = 1
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: ts
    CONTAINS
!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) = vdot(1)
          ELSE
             r(1) = vdot(1) - v(1)*ucp(1,1) - ucp(2,1) - one - t
          END IF
          RETURN
!dec$ attributes offload : mic :: odedef
       END SUBROUTINE odedef
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,ut,ux,ncode,v,vdot,res,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)    :: res(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ut(npde), ux(npde),  &
                                                 v(ncode), vdot(ncode)
!         .. Executable Statements ..
          IF (ires==-1) THEN
             res(1) = v(1)*v(1)*ut(1) - x*u(2)*v(1)*vdot(1)
             res(2) = 0.0_nag_wp
          ELSE
             res(1) = v(1)*v(1)*ut(1) - x*u(2)*v(1)*vdot(1) - ux(2)
             res(2) = u(2) - ux(1)
          END IF
          RETURN
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndary
       SUBROUTINE bndary(npde,t,ibnd,nobc,u,ut,ncode,v,vdot,res,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, ncode, nobc, npde
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: res(nobc)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ut(npde), v(ncode),  &
                                                 vdot(ncode)
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             IF (ires==-1) THEN
                res(1) = 0.0_nag_wp
             ELSE
                res(1) = u(2) + v(1)*exp(t)
             END IF
          ELSE
             IF (ires==-1) THEN
                res(1) = v(1)*vdot(1)
             ELSE
                res(1) = u(2) + v(1)*vdot(1)
             END IF
          END IF
          RETURN
!dec$ attributes offload : mic :: bndary
       END SUBROUTINE bndary
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,x,u,ncode,neqn)

!         Routine for PDE initial values

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ncode, neqn, npde, npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(neqn)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          INTEGER                             :: i, k
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          k = 1
          DO i = 1, npts
             u(k) = exp(ts*(one-x(i))) - one
             u(k+1) = -ts*exp(ts*(one-x(i)))
             k = k + 2
          END DO
          u(neqn) = ts
          RETURN
!dec$ attributes offload : mic :: uvinit
       END SUBROUTINE uvinit
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(time,neqn,npts,x,u)
!         Exact solution (for comparison purposes)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: time
          INTEGER, INTENT (IN)                :: neqn, npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(neqn)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          INTEGER                             :: i, k
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          k = 1
          DO i = 1, npts
             u(k) = exp(time*(one-x(i))) - one
             k = k + 2
          END DO
          RETURN
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
    END MODULE d03pkfe_mod
    PROGRAM d03pkfe

!      D03PKF Example Program Text

!      .. Use Statements ..
       USE nag_library, ONLY : d03pkf
       USE d03pkfe_mod, ONLY : bndary, exact, itrace, nag_wp, ncode, nin,      &
                               nleft, nout, npde, nxi, odedef, one, pdedef,    &
                               ts, uvinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: tout
       INTEGER                             :: i, ifail, ind, it, itask, itol,  &
                                              latol, lenode, lisave, lrsave,   &
                                              lrtol, neqn, npts, nwkres
       CHARACTER (1)                       :: laopt, norm
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: algopt(30), xi(nxi)
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), exy(:), rsave(:),       &
                                              rtol(:), u(:), x(:)
       INTEGER, ALLOCATABLE                :: isave(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              mod, real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PKF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) npts
       neqn = npde*npts + ncode
       nwkres = npde*(npts+6*nxi+3*npde+15) + ncode + nxi + 7*npts + 2
       lenode = 11*neqn + 50
       lisave = 25*neqn + 24
       lrsave = neqn*neqn + neqn + nwkres + lenode
       ALLOCATE (exy(neqn),u(neqn),rsave(lrsave),x(npts),isave(lisave))

       READ (nin,*) itol
       latol = 1
       lrtol = 1
       IF (itol>2) latol = neqn
       IF (mod(itol,2)==0) lrtol = neqn
       ALLOCATE (atol(latol),rtol(lrtol))
       READ (nin,*) atol(1:latol), rtol(1:lrtol)
       READ (nin,*) ts

!      Set spatial-mesh points
       DO i = 1, npts
          x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
       END DO

       READ (nin,*) xi(1:nxi)
       READ (nin,*) norm, laopt
       ind = 0
       itask = 1

       algopt(1:30) = 0.0_nag_wp
       algopt(1) = one
       algopt(13) = 0.005_nag_wp

!      Loop over output value of t

       CALL uvinit(npde,npts,x,u,ncode,neqn)

       tout = 0.2_nag_wp
       DO it = 1, 5
!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
          ifail = 0
          CALL d03pkf(npde,ts,tout,pdedef,bndary,u,npts,x,nleft,ncode,odedef, &
             nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,rsave,lrsave,isave, &
             lisave,itask,itrace,ind,ifail)

          IF (it==1) THEN
             WRITE (nout,99997) atol, npts
             WRITE (nout,99999) (x(i),i=1,13,4), x(npts)
          END IF

!         Check against the exact solution

          CALL exact(tout,neqn,npts,x,exy)

          WRITE (nout,99998) ts
          WRITE (nout,99995) (u(i),i=1,25,4*npde), u(neqn-2), u(neqn)
          WRITE (nout,99994) (exy(i),i=1,25,4*npde), exy(neqn-2), ts
          tout = 2.0_nag_wp*tout
       END DO
       WRITE (nout,99996) isave(1), isave(2), isave(3), isave(5)

99999  FORMAT ('  X        ',5F9.3/)
99998  FORMAT (' T = ',F6.3)
99997  FORMAT (//'  Accuracy requirement =',E10.3,' Number of points = ',I3/)
99996  FORMAT (' Number of integration steps in time = ',I6/' Number o', &
          'f function evaluations = ',I6/' Number of Jacobian eval', &
          'uations =',I6/' Number of iterations = ',I6)
99995  FORMAT (1X,'App.  sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3)
99994  FORMAT (1X,'Exact sol.  ',F7.3,4F9.3,'  ODE sol. =',F8.3/)
    END PROGRAM d03pkfe
