Program d06abfe

!     D06ABF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. 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
      Read (nin,*)

!     Reading of the geometry:
!          number of points on RAE aerofoil data;
!          number of points on circular boundary;
!          maximum number of vertices.

      Read (nin,*) nrae, nvint, nvmax
      Read (nin,*) pmesh

      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
        Read (nin,*) i1, coor(1,nvint+i), coor(2,nvint+i)
      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