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

    MODULE d03psfe_mod

!      D03PSF 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       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: itrace = 0, ncode = 0, nin = 5,  &
                                              nout = 6, npde = 1, nxfix = 0,   &
                                              nxi = 0
    CONTAINS
       SUBROUTINE exact(t,u,x,npts)
!         Exact solution (for comparison and b.c. 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(1,npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: x(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: del, psi, rm, rn, s
          INTEGER                             :: i
!         .. Executable Statements ..
          s = 0.1_nag_wp
          del = 0.01_nag_wp
          rm = -one/del
          rn = one + s/del
          DO i = 1, npts
             psi = x(i) - t
             IF (psi<s) THEN
                u(1,i) = one
             ELSE IF (psi>(del+s)) THEN
                u(1,i) = zero
             ELSE
                u(1,i) = rm*psi + rn
             END IF
          END DO
          RETURN
       END SUBROUTINE exact
       SUBROUTINE uvin1(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, tmp
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              sin
!         .. Executable Statements ..
          tmp = zero
          pi = x01aaf(tmp)
          DO i = 1, npts
             IF (x(i)>0.2_nag_wp .AND. x(i)<=0.4_nag_wp) THEN
                tmp = pi*(5.0_nag_wp*x(i)-one)
                u(1,i) = sin(tmp)
             ELSE
                u(1,i) = zero
             END IF
          END DO
          RETURN
       END SUBROUTINE uvin1
       SUBROUTINE pdef1(npde,t,x,u,ux,ncode,v,vdot,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)                :: ncode, 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), v(ncode),  &
                                                 vdot(ncode)
!         .. Executable Statements ..
          p(1,1) = one
          c(1) = 0.002_nag_wp
          d(1) = ux(1)
          s(1) = zero
          RETURN
       END SUBROUTINE pdef1
       SUBROUTINE bndry1(npde,npts,t,x,u,ncode,v,vdot,ibnd,g,ires)

!         Zero solution at both boundaries

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, ncode, npde, npts
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,npts), v(ncode),       &
                                                 vdot(ncode), x(npts)
!         .. Executable Statements ..
          IF (ibnd==0) THEN
             g(1) = u(1,1)
          ELSE
             g(1) = u(1,npts)
          END IF
          RETURN
       END SUBROUTINE bndry1
       SUBROUTINE monit1(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)                  :: h1, h2, h3
          INTEGER                             :: i
!         .. Intrinsic Functions ..
          INTRINSIC                              abs
!         .. 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 ..
             fmon(i) = abs(((u(1,i+1)-u(1,i))/h2-(u(1,i)-u(1,i-1))/h1)/h3)
          END DO
          fmon(1) = fmon(2)
          fmon(npts) = fmon(npts-1)
          RETURN
       END SUBROUTINE monit1
       SUBROUTINE nmflx1(npde,t,x,ncode,v,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)                :: ncode, npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: flux(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: uleft(npde), uright(npde),    &
                                                 v(ncode)
!         .. Executable Statements ..
          flux(1) = uleft(1)
          RETURN
       END SUBROUTINE nmflx1
       SUBROUTINE uvin2(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)                  :: t
!         .. Executable Statements ..
          t = zero
          CALL exact(t,u,x,npts)
          RETURN
       END SUBROUTINE uvin2
       SUBROUTINE pdef2(npde,t,x,u,ux,ncode,v,vdot,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)                :: ncode, 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), v(ncode),  &
                                                 vdot(ncode)
!         .. Executable Statements ..
          p(1,1) = one
          c(1) = zero
          d(1) = zero
          s(1) = -100.0_nag_wp*u(1)*(u(1)-one)*(u(1)-half)
          RETURN
       END SUBROUTINE pdef2
       SUBROUTINE bndry2(npde,npts,t,x,u,ncode,v,vdot,ibnd,g,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (IN)                :: ibnd, ncode, npde, npts
          INTEGER, INTENT (INOUT)             :: ires
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: g(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,npts), v(ncode),       &
                                                 vdot(ncode), x(npts)
!         .. Local Arrays ..
          REAL (KIND=nag_wp)                  :: ue(1,1)
!         .. Executable Statements ..
!         Solution known to be constant at both boundaries
          IF (ibnd==0) THEN
             CALL exact(t,ue,x(1),1)
             g(1) = ue(1,1) - u(1,1)
          ELSE
             CALL exact(t,ue,x(npts),1)
             g(1) = ue(1,1) - u(1,npts)
          END IF
          RETURN
       END SUBROUTINE bndry2
       SUBROUTINE nmflx2(npde,t,x,ncode,v,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)                :: ncode, npde
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: flux(npde)
          REAL (KIND=nag_wp), INTENT (IN)     :: uleft(npde), uright(npde),    &
                                                 v(ncode)
!         .. Executable Statements ..
          flux(1) = uleft(1)
          RETURN
       END SUBROUTINE nmflx2
       SUBROUTINE monit2(t,npts,npde,x,u,fmon)

!         .. 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)    :: fmon(npts)
          REAL (KIND=nag_wp), INTENT (IN)     :: u(npde,npts), x(npts)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: h1, pi, ux, uxmax, xl, xleft, &
                                                 xmax, xr, xright, xx
          REAL (KIND=nag_wp), SAVE            :: xa = zero
          INTEGER                             :: i
          INTEGER, SAVE                       :: icount = 0
!         .. Intrinsic Functions ..
          INTRINSIC                              abs, cos
!         .. Executable Statements ..
          xx = zero
          pi = x01aaf(xx)
!         Locate shock ..
          uxmax = zero
          xmax = zero
          DO i = 2, npts - 1
             h1 = x(i) - x(i-1)
             ux = abs((u(1,i)-u(1,i-1))/h1)
             IF (ux>uxmax) THEN
                uxmax = ux
                xmax = x(i)
             END IF
          END DO
!         Assign width (on first call only) ..
          IF (icount==0) THEN
             icount = 1
             xleft = xmax - x(1)
             xright = x(npts) - xmax
             IF (xleft>xright) THEN
                xa = xright
             ELSE
                xa = xleft
             END IF
          END IF
          xl = xmax - xa
          xr = xmax + xa
!         Assign monitor function ..
          DO i = 1, npts
             IF (x(i)>xl .AND. x(i)<xr) THEN
                fmon(i) = one + cos(pi*(x(i)-xmax)/xa)
             ELSE
                fmon(i) = zero
             END IF
          END DO
          RETURN
       END SUBROUTINE monit2
    END MODULE d03psfe_mod
    PROGRAM d03psfe

!      D03PSF Example Main Program

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

       CALL ex1

       CALL ex2

    CONTAINS
       SUBROUTINE ex1

!         .. Use Statements ..
          USE nag_library, ONLY : d03pek, d03psf, d03pzf
          USE d03psfe_mod, ONLY : bndry1, itrace, monit1, nag_wp, ncode, nin,  &
                                  nmflx1, npde, nxfix, nxi, one, pdef1, uvin1, &
                                  zero
!         .. 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, m, mlu, neqn, npts,   &
                                                 nrmesh, nwkres
          LOGICAL                             :: remesh
          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(:,:),             &
                                                 uout(:,:,:), x(:), xout(:)
          INTEGER, ALLOCATABLE                :: isave(:)
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
          WRITE (nout,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 1'
!         Skip heading in data file
          READ (nin,*)
          READ (nin,*) npts, intpts, itype
          nwkres = npde*(3*npts+3*npde+32) + 7*npts + 3
          mlu = 3*npde - 1
          neqn = npde*npts + ncode
          lenode = 11*neqn + 50
          lisave = 25 + nxfix + neqn
          lrsave = (3*mlu+1)*neqn + nwkres + lenode

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

!         Initialise mesh
          DO i = 1, npts
             x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
          END DO
          xfix(1) = zero

!         Set remesh parameters
          remesh = .TRUE.
          nrmesh = 3
          dxmesh = zero
          trmesh = zero
          con = 2.0_nag_wp/real(npts-1,kind=nag_wp)
          xratio = 1.5_nag_wp
          ipminf = 0

          xi(1) = zero
          laopt = 'B'
          ind = 0
          itask = 1

          algopt(1:30) = zero
!         b.d.f. integration
          algopt(1) = one
          algopt(13) = 0.005_nag_wp

!         Loop over output value of t

          ts = zero
          tout = zero
          DO it = 1, 3
             tout = real(it,kind=nag_wp)*0.1_nag_wp

!            ifail: behaviour on error exit   
!                   =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft  
             ifail = 0
             CALL d03psf(npde,ts,tout,pdef1,nmflx1,bndry1,uvin1,u,npts,x, &
                ncode,d03pek,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt, &
                remesh,nxfix,xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con, &
                monit1,rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

             IF (it==1) THEN
                WRITE (nout,*)
                WRITE (nout,99998) npts, atol, rtol
             END IF

             WRITE (nout,99999) ts
             WRITE (nout,99996) xout(1:intpts)
!            Interpolate at output points ..
             m = 0

             ifail = 0
             CALL d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

             WRITE (nout,99995) uout(1,1:intpts,1)
          END DO

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

          RETURN

99999     FORMAT (' T = ',F6.3)
99998     FORMAT (/'  NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',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,'Approx U ',7F9.4/)
       END SUBROUTINE ex1
       SUBROUTINE ex2

!         .. Use Statements ..
          USE nag_library, ONLY : d03pek, d03psf, d03pzf
          USE d03psfe_mod, ONLY : bndry2, exact, itrace, monit2, nag_wp,       &
                                  ncode, nin, nmflx2, npde, nxfix, nxi, one,   &
                                  pdef2, uvin2, zero
!         .. 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, m, mlu, neqn, npts,   &
                                                 nrmesh, nwkres
          LOGICAL                             :: remesh
          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,*)
          WRITE (nout,*)
          WRITE (nout,*) 'Example 2'
!         Skip heading in data file
          READ (nin,*)
          READ (nin,*) npts, intpts, itype
          nwkres = npde*(3*npts+3*npde+32) + 7*npts + 3
          mlu = 3*npde - 1
          neqn = npde*npts + ncode
          lenode = 11*neqn + 50
          lisave = 25 + nxfix + neqn
          lrsave = (3*mlu+1)*neqn + nwkres + lenode

          ALLOCATE (rsave(lrsave),u(neqn),ue(1,intpts),uout(1,intpts,itype), &
             x(npts),xout(intpts),isave(lisave))

          READ (nin,*) xout(1:intpts)
          READ (nin,*) itol
          READ (nin,*) norm
          READ (nin,*) atol(1), rtol(1)

!         Initialise mesh
          DO i = 1, npts
             x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
          END DO
          xfix(1) = zero

!         Set remesh parameters
          remesh = .TRUE.
          nrmesh = 5
          dxmesh = zero
          con = one/real(npts-1,kind=nag_wp)
          xratio = 1.5_nag_wp
          ipminf = 0

          xi(1) = zero
          laopt = 'B'
          ind = 0
          itask = 1

          algopt(1:30) = zero
!         Theta integration ..
          algopt(1) = 2.0_nag_wp
          algopt(6) = 2.0_nag_wp
          algopt(7) = 2.0_nag_wp
!         Max. time step ..
          algopt(13) = 2.5E-3_nag_wp

          ts = zero
          tout = zero
          DO it = 1, 2
             tout = real(it,kind=nag_wp)*0.2_nag_wp

             ifail = 0
             CALL d03psf(npde,ts,tout,pdef2,nmflx2,bndry2,uvin2,u,npts,x, &
                ncode,d03pek,nxi,xi,neqn,rtol,atol,itol,norm,laopt,algopt, &
                remesh,nxfix,xfix,nrmesh,dxmesh,trmesh,ipminf,xratio,con, &
                monit2,rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)

             IF (it==1) THEN
                WRITE (nout,*)
                WRITE (nout,99998) npts, atol, rtol
             END IF

             WRITE (nout,99999) ts
             WRITE (nout,99996)
!            Interpolate at output points ..
             m = 0

             ifail = 0
             CALL d03pzf(npde,m,u,npts,x,xout,intpts,itype,uout,ifail)

!            Check against exact solution ..
             CALL exact(tout,ue,xout,intpts)
             DO i = 1, intpts
                WRITE (nout,99995) xout(i), uout(1,i,1), ue(1,i)
             END DO
          END DO

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

          RETURN

99999     FORMAT (' T = ',F6.3)
99998     FORMAT (/' NPTS = ',I4,' ATOL = ',E10.3,' RTOL = ',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 (8X,'X',8X,'Approx U',4X,'Exact U'/)
99995     FORMAT (3(3X,F9.4))
       END SUBROUTINE ex2
    END PROGRAM d03psfe