Example description
!   D03PSF Example Program Text
!   Mark 27.0 Release. NAG Copyright 2019.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: bndry1, bndry2, exact, monit1,       &
                                          monit2, nmflx1, nmflx2, pdef1,       &
                                          pdef2, uvin1, uvin2
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: half = 0.5_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: itrace = 0, nin = 5, nout = 6,       &
                                          npde = 1, nv = 0, 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,nv,v)

!       .. Use Statements ..
        Use nag_library, Only: x01aaf
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: npde, npts, nv, nxi
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: u(npde,npts), v(nv)
        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,nv,v,vdot,p,c,d,s,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. 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(nv), vdot(nv)
!       .. 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,nv,v,vdot,ibnd,g,ires)

!       Zero solution at both boundaries

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: ibnd, npde, npts, nv
        Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: g(npde)
        Real (Kind=nag_wp), Intent (In) :: u(npde,npts), v(nv), vdot(nv),      &
                                          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,nv,v,uleft,uright,flux,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: flux(npde)
        Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde), v(nv)
!       .. Executable Statements ..
        flux(1) = uleft(1)
        Return
      End Subroutine nmflx1
      Subroutine uvin2(npde,npts,nxi,x,xi,u,nv,v)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: npde, npts, nv, nxi
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: u(npde,npts), v(nv)
        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,nv,v,vdot,p,c,d,s,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. 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(nv), vdot(nv)
!       .. 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,nv,v,vdot,ibnd,g,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t
        Integer, Intent (In)           :: ibnd, npde, npts, nv
        Integer, Intent (Inout)        :: ires
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: g(npde)
        Real (Kind=nag_wp), Intent (In) :: u(npde,npts), v(nv), vdot(nv),      &
                                          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,nv,v,uleft,uright,flux,ires)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: t, x
        Integer, Intent (Inout)        :: ires
        Integer, Intent (In)           :: npde, nv
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: flux(npde)
        Real (Kind=nag_wp), Intent (In) :: uleft(npde), uright(npde), v(nv)
!       .. 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, xa, xl, xmax, xr, &
                                          xx
        Integer                        :: i
!       .. 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

!       Desired width of nonzero region of monitor function
        xa = 7.0_nag_wp/60.0_nag_wp

        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 d03psfe_mod, Only: bndry1, itrace, monit1, nin, nmflx1, npde, nv,  &
                               nxfix, nxi, one, pdef1, uvin1, zero
        Use nag_library, Only: d03pek, d03psf, d03pzf, nag_wp
!       .. 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(nxfix), xi(nxi)
        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 + nv
        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)

!       Initialize 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 = zero
        trmesh = zero
        con = 2.0_nag_wp/real(npts-1,kind=nag_wp)
        xratio = 1.5_nag_wp
        ipminf = 0

        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,nv,      &
            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 d03psfe_mod, Only: bndry2, exact, itrace, monit2, nin, nmflx2,     &
                               npde, nv, nxfix, nxi, one, pdef2, uvin2, zero
        Use nag_library, Only: d03pek, d03psf, d03pzf, nag_wp
!       .. 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(nxfix), xi(nxi)
        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 + nv
        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)

!       Initialize 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 = 5
        dxmesh = zero
        con = one/real(npts-1,kind=nag_wp)
        xratio = 1.5_nag_wp
        ipminf = 0

        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,nv,      &
            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