!   D03PSF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    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)

!       .. 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
!       .. 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 Procedures ..
        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)

!       .. 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

!       .. 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)

!       .. 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 Procedures ..
        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)

!       .. 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)

!       .. 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)

!       .. 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)

!       .. 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)

!       .. 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
!       .. 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 Procedures ..
        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, nag_wp
        Use d03psfe_mod, Only: bndry1, itrace, monit1, ncode, nin, nmflx1,     &
                               npde, nxfix, nxi, one, pdef1, uvin1, zero
!       .. 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 Procedures ..
        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, nag_wp
        Use d03psfe_mod, Only: bndry2, exact, itrace, monit2, ncode, nin,      &
                               nmflx2, npde, nxfix, nxi, one, pdef2, uvin2, zero
!       .. 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 Procedures ..
        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