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

    MODULE d03pffe_mod

!      D03PFF 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, npde1 = 2,    &
                                              npde2 = 1
    CONTAINS
!dec$ attributes offload : mic :: exact
       SUBROUTINE exact(t,u,npde,x,npts)
!         Exact solution (for comparison and b.c. purposes)

!         .. Use Statements ..
          USE nag_library, ONLY : x01aaf
!         .. 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
!         .. 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)                  :: pi, px1, px2, x1, x2
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              exp, sin
!         .. Executable Statements ..
          pi = x01aaf(pi)

          DO i = 1, npts
             x1 = x(i) + t
             x2 = x(i) - 3.0_nag_wp*t
             px1 = half*sin(two*pi*x1**2)
             px2 = half*sin(two*pi*x2**2)
             u(1,i) = half*(exp(x2)+exp(x1)+px2-px1) - t*(x1+x2)
             u(2,i) = (exp(x2)-exp(x1)+px2+px1) + x1*x2 + 8.0_nag_wp*t**2
          END DO
          RETURN
!dec$ attributes offload : mic :: exact
       END SUBROUTINE exact
!dec$ attributes offload : mic :: bndry1
       SUBROUTINE bndry1(npde,npts,t,x,u,ibnd,g,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
          REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, npde, npts
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,3), x(npts)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: c, exu1, exu2
!         .. Local Arrays ..
          REAL (KIND=nag_wp)                  :: ue(2,1)
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             CALL exact(t,ue,npde,x(1),1)
             c = (x(2)-x(1))/(x(3)-x(2))
             exu1 = (one+c)*u(1,2) - c*u(1,3)
             exu2 = (one+c)*u(2,2) - c*u(2,3)
             g(1) = two*u(1,1) + u(2,1) - two*ue(1,1) - ue(2,1)
             g(2) = two*u(1,1) - u(2,1) - two*exu1 + exu2
          ELSE
             CALL exact(t,ue,npde,x(npts),1)
             c = (x(npts)-x(npts-1))/(x(npts-1)-x(npts-2))
             exu1 = (one+c)*u(1,2) - c*u(1,3)
             exu2 = (one+c)*u(2,2) - c*u(2,3)
             g(1) = two*u(1,1) - u(2,1) - two*ue(1,1) + ue(2,1)
             g(2) = two*u(1,1) + u(2,1) - two*exu1 - exu2
          END IF
          RETURN
!dec$ attributes offload : mic :: bndry1
       END SUBROUTINE bndry1
!dec$ attributes offload : mic :: nmflx1
       SUBROUTINE nmflx1(npde,t,x,uleft,uright,flux,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER       :: half = 0.5_nag_wp
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t, x
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: flux(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: uleft(npde), uright(npde)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: ltmp, rtmp
!         .. Executable Statements ..
          ltmp = 3.0_nag_wp*uleft(1) + 1.5_nag_wp*uleft(2)
          rtmp = uright(1) - half*uright(2)
          flux(1) = half*(ltmp-rtmp)
          flux(2) = ltmp + rtmp
          RETURN
!dec$ attributes offload : mic :: nmflx1
       END SUBROUTINE nmflx1
!dec$ attributes offload : mic :: pdedef
       SUBROUTINE pdedef(npde,t,x,u,ux,p,c,d,s,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t, x
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: c(npde), d(npde),             &
                                                 p(npde,npde), s(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde), ux(npde)
!         .. Executable Statements ..
          p(1,1) = 1.0_nag_wp
          c(1) = 0.01_nag_wp
          d(1) = ux(1)
          s(1) = u(1)
          RETURN
!dec$ attributes offload : mic :: pdedef
       END SUBROUTINE pdedef
!dec$ attributes offload : mic :: bndry2
       SUBROUTINE bndry2(npde,npts,t,x,u,ibnd,g,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, npde, npts
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,3), x(npts)
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             g(1) = u(1,1) - 3.0_nag_wp
          ELSE
             g(1) = u(1,1) - 5.0_nag_wp
          END IF
          RETURN
!dec$ attributes offload : mic :: bndry2
       END SUBROUTINE bndry2
!dec$ attributes offload : mic :: nmflx2
       SUBROUTINE nmflx2(npde,t,x,uleft,uright,flux,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t, x
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: flux(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: uleft(npde), uright(npde)
!         .. Executable Statements ..
          IF (x>=0.0E0_nag_wp) THEN
             flux(1) = x*uleft(1)
          ELSE
             flux(1) = x*uright(1)
          END IF
          RETURN
!dec$ attributes offload : mic :: nmflx2
       END SUBROUTINE nmflx2
    END MODULE d03pffe_mod
    PROGRAM d03pffe

!      D03PFF Example Main Program

!      .. Use Statements ..
       USE d03pffe_mod, ONLY : nout
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Executable Statements ..
       WRITE (nout,*) 'D03PFF Example Program Results'

       CALL ex1

       CALL ex2

    CONTAINS
!dec$ attributes offload : mic :: ex1
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d03pff, d03pfp, nag_wp
          USE d03pffe_mod, ONLY : bndry1, exact, nin, nmflx1, npde1
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: dx, tout, ts, tsmax
          INTEGER                             :: i, ifail, inc, ind, it,       &
                                                 itask, itrace, lisave,        &
                                                 lrsave, nop, npde, npts, outpts
!         .. Local Arrays ..
          REAL (KIND=nag_wp)                  :: acc(2)
          REAL (KIND=nag_wp), ALLOCATABLE     :: rsave(:), u(:,:), ue(:,:),    &
                                                 x(:), xout(:)
          INTEGER, ALLOCATABLE                :: isave(:)
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 1'
          WRITE (nout,*)
!         Skip heading in data file
          npde = npde1
          READ (nin,*)
          READ (nin,*) npts, inc, outpts
          lisave = 24 + npde*npts
          lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

          ALLOCATE (rsave(lrsave),u(npde,npts),ue(npde,outpts),x(npts), &
             xout(outpts),isave(lisave))
          READ (nin,*) acc(1:2)
          READ (nin,*) itrace
          READ (nin,*) tsmax

!         Initialise mesh
          dx = 1.0_nag_wp/real(npts-1,kind=nag_wp)
          DO i = 1, npts
             x(i) = real(i-1,kind=nag_wp)*dx
          END DO

!         Set initial values
          READ (nin,*) ts
          CALL exact(ts,u,npde,x,npts)

          ind = 0
          itask = 1

          DO it = 1, 2
             tout = 0.1E0_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 d03pff(npde,ts,tout,d03pfp,nmflx1,bndry1,u,npts,x,acc,tsmax, &
                rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

             IF (it==1) THEN
                WRITE (nout,99996) npts, acc(1), acc(2)
                WRITE (nout,99999)
             END IF

!            Set output points

             nop = 0
             DO i = 1, npts, inc
                nop = nop + 1
                xout(nop) = x(i)
             END DO

             WRITE (nout,99995) ts

!            Check against exact solution

             CALL exact(tout,ue,npde,xout,nop)
             nop = 0
             DO i = 1, npts, inc
                nop = nop + 1
                WRITE (nout,99998) xout(nop), u(1,i), ue(1,nop), u(2,i), &
                   ue(2,nop)
             END DO
          END DO

          WRITE (nout,99997) isave(1), isave(2), isave(3), isave(5)

          RETURN

99999     FORMAT (8X,'X',8X,'Approx U',4X,'Exact U',5X,'Approx V',4X,'Exac', &
             't V')
99998     FORMAT (5(3X,F9.4))
99997     FORMAT (/' Number of integration steps in time = ',I6/' Number ', &
             'of function evaluations = ',I6/' Number of Jacobian ', &
             'evaluations =',I6/' Number of iterations = ',I6)
99996     FORMAT (/' NPTS = ',I4,' ACC(1) = ',E10.3,' ACC(2) = ',E10.3/)
99995     FORMAT (/' T = ',F6.3/)
!dec$ attributes offload : mic :: ex1
       END SUBROUTINE ex1
!dec$ attributes offload : mic :: ex2
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d03pff, nag_wp
          USE d03pffe_mod, ONLY : bndry2, nin, nmflx2, npde2, pdedef
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: dx, tout, ts, tsmax
          INTEGER                             :: i, ifail, ind, it, itask,     &
                                                 itrace, lisave, lrsave, npde, &
                                                 npts, outpts
!         .. Local Arrays ..
          REAL (KIND=nag_wp)                  :: acc(2)
          REAL (KIND=nag_wp), ALLOCATABLE     :: rsave(:), u(:,:), x(:), xout(:)
          INTEGER, ALLOCATABLE                :: iout(:), isave(:)
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 2'
          WRITE (nout,*)
          npde = npde2
          READ (nin,*)
          READ (nin,*) npts, outpts
          lisave = 24 + npde*npts
          lrsave = (11+9*npde)*npde*npts + (32+3*npde)*npde + 7*npts + 54

          ALLOCATE (rsave(lrsave),u(npde,npts),x(npts),xout(outpts), &
             isave(lisave),iout(outpts))
          READ (nin,*) iout(1:outpts)
          READ (nin,*) acc(1:2)
          READ (nin,*) itrace
          READ (nin,*) tsmax

!         Initialise mesh
          dx = 2.0_nag_wp/real(npts-1,kind=nag_wp)
          DO i = 1, npts
             x(i) = -1.0_nag_wp + real(i-1,kind=nag_wp)*dx
          END DO

!         Set initial values
          u(1,1:npts) = x(1:npts) + 4.0_nag_wp

          ind = 0
          itask = 1

!         Set output points from inout indices
          DO i = 1, outpts
             xout(i) = x(iout(i))
          END DO

!         Loop over output value of t

          READ (nin,*) ts, tout
          DO it = 1, 2

             ifail = 0
             CALL d03pff(npde,ts,tout,pdedef,nmflx2,bndry2,u,npts,x,acc,tsmax, &
                rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

             IF (it==1) THEN
                WRITE (nout,99998) npts, acc(1), acc(2)
                WRITE (nout,99996) xout(1:outpts)
                tout = 10.0_nag_wp
             END IF

             WRITE (nout,99999) ts
             WRITE (nout,99995) (u(1,iout(i)),i=1,outpts)
          END DO

          WRITE (nout,99997) isave(1), isave(2), isave(3), isave(5)

          RETURN

99999     FORMAT (' T = ',F6.3)
99998     FORMAT (/' NPTS = ',I4,'  ACC(1) = ',E10.3,' ACC(2) = ',E10.3/)
99997     FORMAT (' Number of integration steps in time = ',I6/' Number ', &
             'of function evaluations = ',I6/' Number of Jacobian ', &
             'evaluations =',I6/' Number of iterations = ',I6)
99996     FORMAT (1X,'X    ',7F9.4/)
99995     FORMAT (1X,'U    ',7F9.4/)
!dec$ attributes offload : mic :: ex2
       END SUBROUTINE ex2
    END PROGRAM d03pffe
