```    Program d06aafe

!     D06AAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2017.

!     .. Use Statements ..
Use nag_library, Only: d06aaf, nag_wp, x01aaf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: meshout = 7, nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: coef, pi2, power, r, theta, theta_i, &
x0, y0
Integer                          :: i, ifail, itrace, liwork, lrwork,    &
nedge, nelt, nv, nvb, nvb1, nvb2,    &
nvb3, nvmax
Logical                          :: smooth
Character (1)                    :: pmesh
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: bspace(:), coor(:,:), rwork(:)
Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: cos, max, real, sin
!     .. Executable Statements ..
Write (nout,*) 'D06AAF Example Program Results'

!     Skip heading in data file

!     Coordinates of the boundary mesh vertices and
!     edges references.

Read (nin,*) nvb1, nvb2, nvb3, nvmax
nvb = nvb1 + nvb2 + nvb3
nedge = nvb

lrwork = nvmax
liwork = 16*nvmax + 2*nedge + max(4*nvmax+2,nedge-14)
Allocate (bspace(nvb),coor(2,nvmax),rwork(lrwork),conn(3,2*(nvmax-       &
1)),edge(3,nedge),iwork(liwork))

!     Outer circle
pi2 = 2.0_nag_wp*x01aaf(theta)
theta = pi2/real(nvb1,kind=nag_wp)
r = 1.0_nag_wp
x0 = 0.0_nag_wp
y0 = 0.0_nag_wp
Do i = 1, nvb1
theta_i = theta*real(i,kind=nag_wp)
coor(1,i) = x0 + r*cos(theta_i)
coor(2,i) = y0 + r*sin(theta_i)
End Do
!     Larger inner circle
theta = pi2/real(nvb2,kind=nag_wp)
r = 0.49_nag_wp
x0 = -0.5_nag_wp
y0 = 0.0_nag_wp
Do i = 1, nvb2
theta_i = theta*real(i,kind=nag_wp)
coor(1,nvb1+i) = x0 + r*cos(theta_i)
coor(2,nvb1+i) = y0 + r*sin(theta_i)
End Do
!     Smaller inner circle
theta = pi2/real(nvb3,kind=nag_wp)
r = 0.15_nag_wp
x0 = -0.5_nag_wp
y0 = 0.65_nag_wp
Do i = 1, nvb3
theta_i = theta*real(i,kind=nag_wp)
coor(1,nvb1+nvb2+i) = x0 + r*cos(theta_i)
coor(2,nvb1+nvb2+i) = y0 + r*sin(theta_i)
End Do

!     Boundary edges

Do i = 1, nedge
edge(1,i) = i
edge(2,i) = i + 1
edge(3,i) = 1
End Do
edge(2,nvb1) = 1
edge(2,nvb1+nvb2) = nvb1 + 1
edge(2,nvb) = nvb1 + nvb2 + 1

!     Initialize mesh control parameters

bspace(1:nvb) = 0.05E0_nag_wp
smooth = .True.
itrace = 0
coef = 0.75E0_nag_wp
power = 0.25E0_nag_wp

!     Call to the mesh generator

ifail = 0
Call d06aaf(nvb,nvmax,nedge,edge,nv,nelt,coor,conn,bspace,smooth,coef,   &
power,itrace,rwork,lrwork,iwork,liwork,ifail)

Write (nout,*)

Select Case (pmesh)
Case ('N')
Write (nout,99999) 'NV   =', nv
Write (nout,99999) 'NELT =', nelt
Case ('Y')

!       Output the mesh in a form suitable for printing

Write (meshout,*) '# D06ABF Example Program Mesh results'
Do i = 1, nelt
Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,99998) coor(1,conn(2,i)), coor(2,conn(2,i))
Write (meshout,99998) coor(1,conn(3,i)), coor(2,conn(3,i))
Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
Case Default
Write (nout,*) 'Problem with the printing option Y or N'
End Select

99999 Format (1X,A,I6)
99998 Format (2(2X,E13.6))
End Program d06aafe
```