```!   D06BAF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.
Module d06bafe_mod

!     D06BAF 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                           :: fbnd
!     .. Parameters ..
Integer, Parameter, Public       :: nin = 5, nout = 6
Contains
Function fbnd(i,x,y,ruser,iuser)

!       .. Function Return Value ..
Real (Kind=nag_wp)             :: fbnd
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x, y
Integer, Intent (In)           :: i
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: radius2, x0, xa, xb, y0
!       .. Executable Statements ..
xa = ruser(1)
xb = ruser(2)
x0 = ruser(3)
y0 = ruser(4)

fbnd = 0.0_nag_wp

Select Case (i)
Case (1)

!         line 1,2,3, and 4: ellipse centred in (X0,Y0) with
!         XA and XB as coefficients

fbnd = ((x-x0)/xa)**2 + ((y-y0)/xb)**2 - 1.0_nag_wp
Case (2)

!         line 7, lower arc on letter n, is a circle centred in (X0,Y0)

x0 = 0.5_nag_wp
y0 = 6.25_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (3)

!         line 11, upper arc on letter n, is a circle centred in (X0,Y0)

x0 = 1.0_nag_wp
y0 = 4.0_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (4)

!         line 15, upper arc on letter a, is a circle centred in (X0,Y0)

x0 = 8.5_nag_wp
y0 = 2.75_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (5)

!         line 25, lower arc on hat of 'a', is a circle centred in (X0,Y0)

x0 = 8.5_nag_wp
y0 = 4.0_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (6)

!         lines 20, 21 and 22, belly of letter a, is an ellipse centered
!         in (X0, Y0) with semi-axes 3.5 and 2.75.

x0 = 8.5_nag_wp
y0 = 5.75_nag_wp
fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
Case (7)

!         lines 43, 44 and 45, outer curve on bottom of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 3.5 and 2.5.

x0 = 17.5_nag_wp
y0 = 2.5_nag_wp
fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.5_nag_wp)**2 - 1.0_nag_wp
Case (8)

!         lines 28, 29 and 30, inner curve on bottom of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 2.0 and 1.5.

x0 = 17.5_nag_wp
y0 = 2.5_nag_wp
fbnd = ((x-x0)/2.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
Case (9)

!         line 42, inner curve on lower middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.5 and 0.5.

x0 = 17.5_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/0.5_nag_wp)**2 - 1.0_nag_wp
Case (10)

!         line 31, outer curve on lower middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 2.0 and 1.5.

x0 = 17.5_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
Case (11)

!         line 41, inner curve on upper middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.0 and 1.0.

x0 = 17.0_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.0_nag_wp)**2 + ((y-y0)/1.0_nag_wp)**2 - 1.0_nag_wp
Case (12)

!         line 32, outer curve on upper middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.5 and 1.1573.

x0 = 16.0_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/1.1573_nag_wp)**2 -          &
1.0_nag_wp
Case (13)

!         lines 33, 33, 34, 39 and 40, upper portion of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 3.0 and 2.75.
x0 = 17.0_nag_wp
y0 = 9.25_nag_wp
fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
End Select

Return

End Function fbnd
End Module d06bafe_mod
Program d06bafe

!     D06BAF Example Main Program

!     .. Use Statements ..
Use d06bafe_mod, Only: fbnd, nin, nout
Use nag_library, Only: d06abf, d06acf, d06baf, f16dnf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: x0, xa, xb, xmax, xmin, y0, ymax,    &
ymin
Integer                          :: i, ifail, itrace, j, k, liwork,      &
lrwork, maxind, maxval, ncomp,       &
nedge, nedmx, nelt, nlines, npropa,  &
nv, nvb, nvint, nvmax, reftk, sdcrus
Character (1)                    :: pmesh
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: coor(:,:), coorch(:,:), crus(:,:),   &
rate(:), rwork(:), weight(:)
Real (Kind=nag_wp)               :: ruser(4)
Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:),      &
lcomp(:), lined(:,:), nlcomp(:)
Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
Intrinsic                        :: abs
!     .. Executable Statements ..
Write (nout,*) 'D06BAF Example Program Results'
Flush (nout)

!     Skip heading in data file

!     Initialize boundary mesh inputs:
!     the number of line and of the characteristic points of
!     the boundary mesh

Allocate (coor(2,nvmax),coorch(2,nlines),rate(nlines),edge(3,nedmx),     &
lcomp(nlines),lined(4,nlines))

!     The Lines of the boundary mesh

sdcrus = 0

Do i = 1, nlines

If (lined(4,i)<0) Then
sdcrus = sdcrus + lined(1,i) - 2
End If

End Do

liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus

!     Get max(LINED(1,:)) for computing LRWORK

Call f16dnf(nlines,lined,4,maxind,maxval)

lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines

Allocate (crus(2,sdcrus),iwork(liwork),rwork(lrwork))

!     The ellipse boundary which envelops the NAG Logo
!     the N, the A and the G

!     The number of connected components to the boundary
!     and their information

Allocate (nlcomp(ncomp))

j = 1

Do i = 1, ncomp
k = j + abs(nlcomp(i)) - 1
j = k + 1
End Do

!     Data passed to the user-supplied function

xmin = coorch(1,4)
xmax = coorch(1,2)
ymin = coorch(2,1)
ymax = coorch(2,3)

xa = (xmax-xmin)/2.0_nag_wp
xb = (ymax-ymin)/2.0_nag_wp

x0 = (xmin+xmax)/2.0_nag_wp
y0 = (ymin+ymax)/2.0_nag_wp

ruser(1:4) = (/xa,xb,x0,y0/)
iuser(1) = 0

itrace = -1

Write (nout,*)
Flush (nout)

!     Call to the boundary mesh generator

ifail = 0
Call d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp,      &
lcomp,nvmax,nedmx,nvb,coor,nedge,edge,itrace,ruser,iuser,rwork,lrwork, &
iwork,liwork,ifail)

Select Case (pmesh)
Case ('N')
Write (nout,*) 'Boundary mesh characteristics'
Write (nout,99999) 'NVB   =', nvb
Write (nout,99999) 'NEDGE =', nedge
Case ('Y')

!       Output the mesh

Write (nout,99998) nvb, nedge

Do i = 1, nvb
Write (nout,99997) i, coor(1:2,i)
End Do

Do i = 1, nedge
Write (nout,99996) i, edge(1:3,i)
End Do

Flush (nout)

Case Default
Write (nout,*) 'Problem with the printing option Y or N'
Go To 100
End Select

Deallocate (rwork,iwork)

!     Initialize mesh control parameters

itrace = 0
npropa = 1
nvint = 0
lrwork = 12*nvmax + 15
liwork = 6*nedge + 32*nvmax + 2*nvb + 78
Allocate (weight(nvint),rwork(lrwork),iwork(liwork),conn(3,2*nvmax+5))

!     Call to the 2D Delaunay-Voronoi mesh generator

ifail = 0
Call d06abf(nvb,nvint,nvmax,nedge,edge,nv,nelt,coor,conn,weight,npropa,  &
itrace,rwork,lrwork,iwork,liwork,ifail)

Select Case (pmesh)
Case ('N')
Write (nout,*) 'Complete mesh (via the 2D Delaunay-Voronoi'
Write (nout,*) 'mesh generator) characteristics'
Write (nout,99999) 'NV (rounded to nearest 10)   =', 10*((nv+5)/10)
Write (nout,99999) 'NELT (rounded to nearest 10) =', 10*((nelt+5)/10)
Case ('Y')

!       Output the mesh

Write (nout,99998) nv, nelt

Do i = 1, nv
Write (nout,99995) coor(1:2,i)
End Do

reftk = 0

Do k = 1, nelt
Write (nout,99994) conn(1:3,k), reftk
End Do

Flush (nout)

End Select

Deallocate (rwork,iwork)
lrwork = 12*nvmax + 30015
liwork = 8*nedge + 53*nvmax + 2*nvb + 10078
Allocate (rwork(lrwork),iwork(liwork))

!     Call to the 2D Advancing front mesh generator

ifail = 0
Call d06acf(nvb,nvint,nvmax,nedge,edge,nv,nelt,coor,conn,weight,itrace,  &
rwork,lrwork,iwork,liwork,ifail)

Select Case (pmesh)
Case ('N')
Write (nout,*) 'Complete mesh (via the 2D Advancing front mesh'
Write (nout,*) 'generator) characteristics'
Write (nout,99999) 'NV (rounded to nearest 10)   =', 10*((nv+5)/10)
Write (nout,99999) 'NELT (rounded to nearest 10) =', 10*((nelt+5)/10)
Case ('Y')

!       Output the mesh

Write (nout,99998) nv, nelt

Do i = 1, nv
Write (nout,99995) coor(1:2,i)
End Do

reftk = 0

Do k = 1, nelt
Write (nout,99994) conn(1:3,k), reftk
End Do

End Select

100   Continue

99999 Format (1X,A,I6)
99998 Format (1X,2I10)
99997 Format (2X,I4,2(2X,E13.6))
99996 Format (1X,4I4)
99995 Format (2(2X,E13.6))
99994 Format (1X,4I10)
End Program d06bafe
```