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

    MODULE d03prfe_mod

!      D03PRF 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       :: half = 0.5_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
       INTEGER, PARAMETER                  :: itrace = 0, ncode = 0, nin = 5,  &
                                              nleft = 1, nout = 6, npde = 2,   &
                                              nxfix = 0, nxi = 0
    CONTAINS
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,nxi,x,xi,u,ncode,v)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ncode, npde, npts, nxi
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npde,npts), v(ncode)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts), xi(nxi)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: pi
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp, sin
!         .. Executable Statements ..
          pi = x01aaf(pi)
          DO i = 1, npts
             u(1,i) = exp(x(i))
             u(2,i) = x(i)**2 + sin(two*pi*x(i)**2)
          END DO
          RETURN
!dec$ attributes offload : mic :: uvinit
       END SUBROUTINE uvinit
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,udot,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), udot(npde),          &
                                                 ux(npde), v(ncode), vdot(ncode)
!         .. Executable Statements ..
          IF (ires==-1) THEN
             res(1) = udot(1)
             res(2) = udot(2)
          ELSE
             res(1) = udot(1) + ux(1) + ux(2)
             res(2) = udot(2) + 4.0_nag_wp*ux(1) + ux(2)
          END IF
          RETURN
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndary
       SUBROUTINE bndary(npde,t,ibnd,nobc,u,udot,ncode,v,vdot,res,ires)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. 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), udot(npde),          &
                                                 v(ncode), vdot(ncode)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: pp, ppt1, ppt3, t1, t3
!         .. Intrinsic Functions ..
          INTRINSIC                              exp, sin
!         .. Executable Statements ..
          IF (ires==-1) THEN
             res(1) = 0.0_nag_wp
          ELSE
             pp = two*x01aaf(pp)
             t1 = t
             t3 = -3.0_nag_wp*t
             IF (ibnd==0) THEN
                ppt3 = sin(pp*t3**2)
                ppt1 = sin(pp*t1**2)
                res(1) = u(1) - half*(exp(t3)+exp(t1)+half*(ppt3-ppt1))
                res(1) = res(1) - 2.0_nag_wp*t**2
             ELSE
                t3 = t3 + 1.0_nag_wp
                t1 = t1 + 1.0_nag_wp
                ppt3 = sin(pp*t3**2)
                ppt1 = sin(pp*t1**2)
                res(1) = u(2) - (exp(t3)-exp(t1)+half*(ppt3+ppt1))
                res(1) = res(1) - (1.0_nag_wp+5.0_nag_wp*t**2-2.0_nag_wp*t)
             END IF
          END IF
          RETURN
!dec$ attributes offload : mic :: bndary
       END SUBROUTINE bndary
!dec$ attributes offload : mic :: monitf
       SUBROUTINE monitf(t,npts,npde,x,u,fmon)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)     :: u(npde,npts), x(npts)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: d2x1, d2x2, h1, h2, h3
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              abs, max
!         .. Executable Statements ..
          DO i = 2, npts - 1
             h1 = x(i) - x(i-1)
             h2 = x(i+1) - x(i)
             h3 = half*(x(i+1)-x(i-1))
!            Second derivatives ..
             d2x1 = abs(((u(1,i+1)-u(1,i))/h2-(u(1,i)-u(1,i-1))/h1)/h3)
             d2x2 = abs(((u(2,i+1)-u(2,i))/h2-(u(2,i)-u(2,i-1))/h1)/h3)
             fmon(i) = max(d2x1,d2x2)
          END DO
          fmon(1) = fmon(2)
          fmon(npts) = fmon(npts-1)
          RETURN
!dec$ attributes offload : mic :: monitf
       END SUBROUTINE monitf
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(t,npde,npts,x,u)

!         Exact solution (for comparison purposes)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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)                  :: pp, ppt1, ppt3, x1, x3
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp, sin
!         .. Executable Statements ..
          pp = 2.0_nag_wp*x01aaf(pp)
          DO i = 1, npts
             x1 = x(i) + t
             x3 = x(i) - 3.0_nag_wp*t
             ppt3 = sin(pp*x3**2)
             ppt1 = sin(pp*x1**2)
             u(1,i) = half*(exp(x3)+exp(x1)+half*(ppt3-ppt1)) - two*x(i)*t + &
                two*t**2
             u(2,i) = (exp(x3)-exp(x1)+half*(ppt3+ppt1)) - two*x(i)*t + &
                x(i)**2 + 5.0_nag_wp*t**2
          END DO
          RETURN
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
    END MODULE d03prfe_mod
    PROGRAM d03prfe

!      D03PRF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pek, d03prf, d03pzf
       USE d03prfe_mod, ONLY : bndary, exact, itrace, monitf, nag_wp, ncode,   &
                               nin, nleft, nout, npde, nxfix, nxi, pdedef,     &
                               uvinit
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: con, dxmesh, tout, trmesh, ts,   &
                                              xratio
       INTEGER                             :: i, ifail, ind, intpts, ipminf,   &
                                              it, itask, itol, itype, lenode,  &
                                              lisave, lrsave, neqn, npts,      &
                                              nrmesh, nwkres
       LOGICAL                             :: remesh, theta
       CHARACTER (1)                       :: laopt, norm
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: algopt(30), atol(1), rtol(1),    &
                                              xfix(1), xi(1)
       REAL (KIND=nag_wp), ALLOCATABLE     :: rsave(:), u(:,:), ue(:,:),       &
                                              uout(:,:,:), x(:), xout(:)
       INTEGER, ALLOCATABLE                :: isave(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PRF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) npts, intpts, itype
       lisave = 25 + nxfix
       neqn = npde*npts + ncode
       nwkres = npde*(npts+21+3*npde) + 7*npts + nxfix + 3
       lenode = 11*neqn + 50
       lrsave = neqn*neqn + neqn + nwkres + lenode

       ALLOCATE (rsave(lrsave),u(npde,npts),ue(npde,npts), &
          uout(npde,intpts,itype),x(npts),xout(intpts),isave(lisave))
       READ (nin,*) itol
       READ (nin,*) atol(1), rtol(1)

!      Set remesh parameters
       remesh = .TRUE.
       nrmesh = 3
       dxmesh = 0.0_nag_wp
       con = 5.0_nag_wp/real(npts-1,kind=nag_wp)
       xratio = 1.2_nag_wp
       ipminf = 0

!      Initialise mesh
       DO i = 1, npts
          x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
       END DO
       xi(1) = 0.0_nag_wp

       READ (nin,*) xout(1:intpts)
       READ (nin,*) norm, laopt
       ind = 0
       itask = 1

!      Set theta to .TRUE. if the Theta integrator is required
       theta = .FALSE.
       algopt(1:30) = 0.0_nag_wp
       IF (theta) THEN
          algopt(1) = 2.0_nag_wp
          algopt(6) = 2.0_nag_wp
          algopt(7) = 1.0_nag_wp
       END IF

!      Loop over output value of t
       ts = 0.0_nag_wp
       tout = 0.0_nag_wp

       DO it = 1, 5
          tout = 0.05_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 d03prf(npde,ts,tout,pdedef,bndary,uvinit,u,npts,x,nleft,ncode, &
             d03pek,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix, &
             xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con,monitf,rsave,lrsave, &
             isave,lisave,itask,itrace,ind,ifail)

          IF (it==1) THEN
             WRITE (nout,99996) atol, npts
             WRITE (nout,99999) nrmesh
             WRITE (nout,99998) xout(1:intpts)
          END IF

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

!         Check against exact solution ..
          CALL exact(ts,npde,intpts,xout,ue)

          WRITE (nout,99997) ts
          WRITE (nout,99994) uout(1,1:intpts,1)
          WRITE (nout,99993) ue(1,1:intpts)
          WRITE (nout,99992) uout(2,1:intpts,1)
          WRITE (nout,99991) ue(2,1:intpts)

       END DO
       WRITE (nout,99995) isave(1), isave(2), isave(3), isave(5)

99999  FORMAT (' Remeshing every ',I3,' time steps'/)
99998  FORMAT (' X        ',5F10.4/)
99997  FORMAT (' T = ',F6.3)
99996  FORMAT (//'  Accuracy requirement =',E10.3,' Number of points = ',I3/)
99995  FORMAT (' Number of integration steps in time = ',I6/' Number o', &
          'f function evaluations = ',I6/' Number of Jacobian eval', &
          'uations =',I6/' Number of iterations = ',I6)
99994  FORMAT (' Approx U1',5F10.4)
99993  FORMAT (' Exact  U1',5F10.4)
99992  FORMAT (' Approx U2',5F10.4)
99991  FORMAT (' Exact  U2',5F10.4/)
    END PROGRAM d03prfe
