NAG Library Manual, Mark 28.4
!   D03PSF Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

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
!       There are no coupled ODEs in this problem (nv = 0):
v(:) = 0._nag_wp
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)
!       There are no coupled ODEs in this problem (nv = 0):
v(:) = 0._nag_wp
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
Use nag_library, Only: nag_wp
!     .. 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
!       .. 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
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))

!       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

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

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

!       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

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

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