Program d06abfe

!     D06ABF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2017.

!     .. Use Statements ..
Use nag_library, Only: d06abf, f06epf, nag_wp, x01aaf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: meshout = 7, nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: c, d, dnvint, r, s, theta, theta_i
Integer                          :: i, i1, ic, ifail, itrace, j, liwork, &
lrwork, nearest, nedge, nelt,        &
nelt_near, npropa, nrae, nv, nvb,    &
nvint, nvmax, nv_near
Character (1)                    :: pmesh
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: coor(:,:), rwork(:), weight(:)
Real (Kind=nag_wp)               :: t(2)
Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: cos, real, sin, tan
!     .. Executable Statements ..
Write (nout,*) 'D06ABF Example Program Results'

!     Skip heading in data file

!          number of points on RAE aerofoil data;
!          number of points on circular boundary;
!          maximum number of vertices.

nvb = nvint + 2*nrae
nedge = nvb

lrwork = 12*nvmax + 15
liwork = 6*nedge + 32*nvmax + 2*nvb + 78
Allocate (coor(2,nvmax),rwork(lrwork),weight(nvint),conn(3,2*nvmax+5),   &
edge(3,nedge),iwork(liwork))

!     Circular outer boundary, radius 3 and centre (1,0)
theta = 2.0_nag_wp*x01aaf(r)/real(nvint,kind=nag_wp)
r = 3.0_nag_wp
t(1) = 1.0_nag_wp
Do i = 1, nvint
theta_i = theta*real(i-1,kind=nag_wp)
coor(1,i) = r*cos(theta_i) + t(1)
coor(2,i) = r*sin(theta_i)
End Do

!     Read data for aerofoil RAE 2822
Do i = 1, nrae
End Do

!     Transform RAE 2822 for secondary foil
theta = x01aaf(theta)/12.0_nag_wp
c = cos(theta)
s = sin(theta)
ic = nvint + nrae
!     Copy and rotate coordinates by theta = pi/12
coor(1:2,ic+1:ic+nrae) = coor(1:2,nvint+1:ic)
Call f06epf(nrae,coor(1,ic+1),2,coor(2,ic+1),2,c,s)
!     Reduce by 0.4 and translate to distance 0.25 from intercept at (0.75,0)
d = 0.4_nag_wp
t(1) = 0.75_nag_wp + 0.25_nag_wp*c
t(2) = -0.25_nag_wp*s
Do i = 1, nrae
coor(1:2,ic+i) = d*coor(1:2,ic+i) + t(1:2)
End Do

!     Boundary edges
Do i = 1, nedge
edge(1,i) = i
edge(2,i) = i + 1
edge(3,i) = 0
End Do
!     Tie up end of three boundary edges
edge(2,nvint) = 1
edge(2,nvint+nrae) = nvint + 1
edge(2,nedge) = nvint + nrae + 1

!     Initialize mesh control parameters

itrace = 0

!     Generation of interior vertices on the
!     RAE airfoil's wake

dnvint = 2.5E0_nag_wp/real(nvint+1,kind=nag_wp)

Do i = 1, nvint
i1 = nvb + i
coor(1,i1) = 1.38E0_nag_wp + real(i,kind=nag_wp)*dnvint
coor(2,i1) = -tan(theta)*(coor(1,i1)-0.75_nag_wp)
End Do

weight(1:nvint) = 0.01E0_nag_wp

Write (nout,*)

!      Loop on the propagation coef

pcoef: Do j = 1, 4

nearest = 250
Select Case (j)
Case (1)
npropa = -5
Case (2)
npropa = -1
Case (3)
npropa = 1
Case Default
npropa = 5
End Select

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

Write (nout,99999) 'Mesh characteristics with NPROPA =', npropa
nv_near = ((nv+nearest/2)/nearest)*nearest
nelt_near = ((nelt+nearest/2)/nearest)*nearest
Write (nout,99998) 'NV  ', nv_near, nearest
Write (nout,99998) 'NELT', nelt_near, nearest

If (pmesh=='Y') Then
!         Output the mesh in a form suitable for printing

If (j==1) Then
Write (meshout,*) '# D06ABF Example Program Mesh results'
End If
Write (meshout,99999) '# Mesh line segments for NPROPA =', npropa
Do i = 1, nelt
Write (meshout,99997) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,99997) coor(1,conn(2,i)), coor(2,conn(2,i))
Write (meshout,99997) coor(1,conn(3,i)), coor(2,conn(3,i))
Write (meshout,99997) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
End If

End Do pcoef

99999 Format (1X,A,I6)
99998 Format (1X,A5,' = ',I10,' to the nearest ',I3)
99997 Format (2(2X,E13.6))
End Program d06abfe