```    Program d06cafe

!     D06CAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2017.

!     .. Use Statements ..
Use nag_library, Only: d06caf, g05kff, g05sqf, nag_wp, x01aaf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: delta, hx, hy, pi2, r, rad, sk,      &
theta, x1, x2, x3, y1, y2, y3
Integer                          :: genid, i, ifail, imax, ind, itrace,  &
j, jmax, k, liwork, lrwork, lseed,   &
lstate, me1, me2, me3, nedge, nelt,  &
nqint, nv, nvfix, reftk, subid
Character (1)                    :: pmesh
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: coor(:,:), rwork(:), variates(:)
Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:),      &
numfix(:), seed(:), state(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: cos, min, real, sin
!     .. Executable Statements ..
Write (nout,*) 'D06CAF Example Program Results'
Flush (nout)

!     Skip heading in data file

!     Read IMAX and JMAX, the number of vertices
!     in the x and y directions respectively.

nv = imax*jmax
nelt = 2*(imax-1)*(jmax-1)
nedge = 2*(imax-1) + 2*(jmax-1)
liwork = 8*nelt + 2*nv
lrwork = 2*nv + nelt

!     The array VARIATES will be used when distorting the mesh

Allocate (variates(2*nv),coor(2,nv),conn(3,nelt),edge(3,nedge),          &
iwork(liwork),rwork(lrwork))

!     of distortion neighbourhood so that cross-over
!     can only occur at 100% or greater.

hx = 1.0E0_nag_wp/real(imax-1,kind=nag_wp)
hy = 1.0E0_nag_wp/real(jmax-1,kind=nag_wp)
pi2 = 2.0E0_nag_wp*x01aaf(pi2)

!     GENID identifies the base generator

genid = 1
subid = 1

!     For GENID = 1 only one seed is required
!     The initializer is first called in query mode to get the value of
!     LSTATE for the chosen base generator

lseed = 1
lstate = -1
Allocate (seed(lseed),state(lstate))

!     Initialize the seed

seed(1:lseed) = (/1762541/)

ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

Deallocate (state)
Allocate (state(lstate))

!     Initialize the generator to a repeatable sequence

ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!     Generate two sets of uniform random variates

ifail = 0

ifail = 0
Call g05sqf(nv,0.0E0_nag_wp,pi2,state,variates(nv+1),ifail)

!     Generate a simple uniform mesh and then distort it
!     randomly within the distortion neighbourhood of each
!     node.

k = 0
ind = 0

Do j = 1, jmax

Do i = 1, imax
k = k + 1
r = variates(k)
theta = variates(nv+k)

If (i==1 .Or. i==imax .Or. j==1 .Or. j==jmax) Then
r = 0.E0_nag_wp
End If

coor(1,k) = real(i-1,kind=nag_wp)*hx + r*cos(theta)
coor(2,k) = real(j-1,kind=nag_wp)*hy + r*sin(theta)

If (i<imax .And. j<jmax) Then
ind = ind + 1
conn(1,ind) = k
conn(2,ind) = k + 1
conn(3,ind) = k + imax + 1
ind = ind + 1
conn(1,ind) = k
conn(2,ind) = k + imax + 1
conn(3,ind) = k + imax
End If

End Do

End Do

Write (nout,*)

Select Case (pmesh)
Case ('N')
Write (nout,*) 'The complete distorted mesh characteristics'
Write (nout,99999) 'Number of vertices =', nv
Write (nout,99999) 'Number of elements =', nelt
Case ('Y')

!       Output the mesh

Write (nout,99998) nv, nelt

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

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

reftk = 0

Do k = 1, nelt
me1 = conn(1,k)
me2 = conn(2,k)
me3 = conn(3,k)

x1 = coor(1,me1)
x2 = coor(1,me2)
x3 = coor(1,me3)
y1 = coor(2,me1)
y2 = coor(2,me2)
y3 = coor(2,me3)

sk = ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))/2.E0_nag_wp

If (sk<0.E0_nag_wp) Then
Write (nout,*) 'Error: the surface of the element is negative'
Write (nout,99999) 'element number  = ', k
Write (nout,99995) 'element surface = ', sk
Go To 100
End If

If (pmesh=='Y') Then
Write (nout,99996) conn(1,k), conn(2,k), conn(3,k), reftk
End If

End Do
Flush (nout)

!     Boundary edges

Do i = 1, imax - 1
edge(1,i) = i
edge(2,i) = i + 1
End Do

Do i = 1, jmax - 1
edge(1,imax-1+i) = i*imax
edge(2,imax-1+i) = (i+1)*imax
End Do

Do i = 1, imax - 1
edge(1,imax-1+jmax-1+i) = imax*jmax - i + 1
edge(2,imax-1+jmax-1+i) = imax*jmax - i
End Do

Do i = 1, jmax - 1
edge(1,2*(imax-1)+jmax-1+i) = (jmax-i)*imax + 1
edge(2,2*(imax-1)+jmax-1+i) = (jmax-i-1)*imax + 1
End Do

edge(3,1:nedge) = 0

nvfix = 0
Allocate (numfix(nvfix))

itrace = 1
nqint = 10

!     Call the smoothing routine

ifail = 0
Call d06caf(nv,nelt,nedge,coor,edge,conn,nvfix,numfix,itrace,nqint,      &
iwork,liwork,rwork,lrwork,ifail)

Select Case (pmesh)
Case ('N')
Write (nout,*)
Write (nout,*) 'The complete smoothed mesh characteristics'
Write (nout,99999) 'Number of vertices =', nv
Write (nout,99999) 'Number of elements =', nelt
Case ('Y')

!       Output the mesh

Write (nout,99998) nv, nelt

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

reftk = 0

Do k = 1, nelt
Write (nout,99996) conn(1,k), conn(2,k), conn(3,k), reftk
End Do

End Select

100   Continue

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