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

    MODULE d03ppfe_mod

!      D03PPF 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
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: e
    CONTAINS
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,nxi,x,xi,u,ncode,v)

!         .. 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)                  :: a, b, c, t
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          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
!dec$ attributes offload : mic :: uvinit
       END SUBROUTINE uvinit
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,ux,ncode,v,vdot,p,q,r,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)    :: p(npde,npde), q(npde), r(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
!         .. Executable Statements ..
          p(1,1) = one
          r(1) = e*ux(1)
          q(1) = u(1)*ux(1)
          RETURN
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndary
       SUBROUTINE bndary(npde,t,u,ux,ncode,v,vdot,ibnd,beta,gamma,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: a, b, c, ue, x
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          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
!dec$ attributes offload : mic :: bndary
       END SUBROUTINE bndary
!dec$ attributes offload : mic :: monitf
       SUBROUTINE monitf(t,npts,npde,x,u,r,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)     :: r(npde,npts), u(npde,npts),   &
                                                 x(npts)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: drdx, h
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              abs
!         .. Executable Statements ..
          fmon(1) = abs((r(1,2)-r(1,1))/((x(2)-x(1))*half))
          DO i = 2, npts - 1
             h = (x(i+1)-x(i-1))*half
!            Second derivative ..
             drdx = (r(1,i+1)-r(1,i))/h
             fmon(i) = abs(drdx)
          END DO
          fmon(npts) = fmon(npts-1)
          RETURN
!dec$ attributes offload : mic :: monitf
       END SUBROUTINE monitf
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(t,x,npts,u)

!         Exact solution (for comparison purposes)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: npts
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: a, b, c
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          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
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
    END MODULE d03ppfe_mod
    PROGRAM d03ppfe

!      D03PPF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03pck, d03ppf, d03pzf
       USE d03ppfe_mod, ONLY : bndary, e, exact, half, itrace, m, monitf,      &
                               nag_wp, ncode, nin, nout, npde, nxfix, nxi,     &
                               pdedef, two, uvinit, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: con, dx, dxmesh, tout, trmesh,   &
                                              ts, x0, xmid, 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                              min, real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PPF 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+3*npde+21) + 7*npts + nxfix + 3
       lenode = 11*neqn + 50
       lrsave = neqn*neqn + neqn + nwkres + lenode

       ALLOCATE (u(neqn),ue(intpts),uout(npde,intpts,itype),rsave(lrsave), &
          x(npts),xout(intpts),isave(lisave))
       READ (nin,*) itol
       READ (nin,*) atol(1), rtol(1)
       READ (nin,*) 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
       ELSE
          algopt(1) = zero
       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 d03ppf(npde,m,ts,tout,pdedef,bndary,uvinit,u,npts,x,ncode, &
             d03pck,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,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)

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

       END DO
       WRITE (nout,99997) isave(1), isave(2), isave(3), isave(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 d03ppfe
