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

    MODULE d03ppae_mod

!      D03PPA 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
    CONTAINS
!dec$ attributes offload : mic :: uvinit
       SUBROUTINE uvinit(npde,npts,nxi,x,xi,u,ncode,v,iuser,ruser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: ncode, npde, npts, nxi
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npde,npts), v(ncode)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts), xi(nxi)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: a, b, c, e, t
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          e = ruser(1)
          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,iuser,ruser)

!         .. 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 (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: e
!         .. Executable Statements ..
          e = ruser(1)
          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,iuser, &
             ruser)

!         .. 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 (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde), v(ncode),  &
                                                 vdot(ncode)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: a, b, c, e, ue, x
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          e = ruser(1)
          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,iuser,ruser)

!         .. 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)
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(*)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: drdx, h
          INTEGER                             :: i, k, l
!         .. Intrinsic Functions ..
          INTRINSIC                              abs, max, min
!         .. Executable Statements ..
          DO i = 1, npts - 1
             k = max(1,i-1)
             l = min(npts,i+1)
             h = (x(l)-x(k))*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,iuser,ruser)
!         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 (INOUT)  :: ruser(*)
          REAL (KIND=nag_wp), INTENT (OUT)    :: u(npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(npts)
          INTEGER, INTENT (INOUT)             :: iuser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: a, b, c, e
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp
!         .. Executable Statements ..
          e = ruser(1)
          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 d03ppae_mod
    PROGRAM d03ppae

!      D03PPA Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d03ppa, d03pzf, d53pck
       USE d03ppae_mod, ONLY : bndary, 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, e, tout,        &
                                              trmesh, ts, x0, xmid, xratio
       INTEGER                             :: i, ifail, ind, intpts, ipminf,   &
                                              it, itask, itol, itype, lenode,  &
                                              neqn, niw, npts, nrmesh, nw,     &
                                              nwkres
       LOGICAL                             :: remesh, theta
       CHARACTER (1)                       :: laopt, norm
!      .. Local Arrays ..
       REAL (KIND=nag_wp)                  :: algopt(30), atol(1), rtol(1),    &
                                              ruser(1), rwsav(1100), xfix(1),  &
                                              xi(1)
       REAL (KIND=nag_wp), ALLOCATABLE     :: u(:), ue(:), uout(:,:,:), w(:),  &
                                              x(:), xout(:)
       INTEGER                             :: iuser(1), iwsav(505)
       INTEGER, ALLOCATABLE                :: iw(:)
       LOGICAL                             :: lwsav(100)
       CHARACTER (80)                      :: cwsav(10)
!      .. Intrinsic Functions ..
       INTRINSIC                              min, real
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PPA Example Program Results'
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) npts, intpts, itype
       neqn = npde*npts + ncode
       niw = 25 + nxfix
       nwkres = npde*(npts+3*npde+21) + 7*npts + nxfix + 3
       lenode = 11*neqn + 50
       nw = neqn*neqn + neqn + nwkres + lenode

       ALLOCATE (u(neqn),ue(intpts),uout(npde,intpts,itype),w(nw),x(npts), &
          xout(intpts),iw(niw))
       READ (nin,*) itol
       READ (nin,*) atol(1), rtol(1)
       READ (nin,*) e
       ruser(1) = 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
       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 d03ppa(npde,m,ts,tout,pdedef,bndary,uvinit,u,npts,x,ncode, &
             d53pck,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt,remesh,nxfix, &
             xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con,monitf,w,nw,iw,niw, &
             itask,itrace,ind,iuser,ruser,cwsav,lwsav,iwsav,rwsav,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,iuser,ruser)

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

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