!   D06ACF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.
    Module d06acfe_mod

!     D06ACF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: 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)                   :: c, radius, x0, x1, y0, y1
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sqrt
!       .. Executable Statements ..
        fbnd = 0.0_nag_wp

        Select Case (i)
        Case (1)

!         upper NACA0012 wing beginning at the origin

          c = 1.008930411365_nag_wp
          fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*x)-0.126_nag_wp*c*x- &
            0.3516_nag_wp*(c*x)**2+0.2843_nag_wp*(c*x)**3- &
            0.1015_nag_wp*(c*x)**4) - c*y
        Case (2)

!         lower NACA0012 wing beginning at the origin

          c = 1.008930411365_nag_wp
          fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*x)-0.126_nag_wp*c*x- &
            0.3516_nag_wp*(c*x)**2+0.2843_nag_wp*(c*x)**3- &
            0.1015_nag_wp*(c*x)**4) + c*y
        Case (3)
          x0 = ruser(1)
          y0 = ruser(2)
          radius = ruser(3)
          fbnd = (x-x0)**2 + (y-y0)**2 - radius**2
        Case (4)

!         upper NACA0012 wing beginning at (X1;Y1)

          c = 1.008930411365_nag_wp
          x1 = ruser(4)
          y1 = ruser(5)
          fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*(x- &
            x1))-0.126_nag_wp*c*(x-x1)-0.3516_nag_wp*(c*(x- &
            x1))**2+0.2843_nag_wp*(c*(x-x1))**3-0.1015_nag_wp*(c*(x-x1))**4) - &
            c*(y-y1)
        Case (5)

!         lower NACA0012 wing beginning at (X1;Y1)

          c = 1.008930411365_nag_wp
          x1 = ruser(4)
          y1 = ruser(5)
          fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*(x- &
            x1))-0.126_nag_wp*c*(x-x1)-0.3516_nag_wp*(c*(x- &
            x1))**2+0.2843_nag_wp*(c*(x-x1))**3-0.1015_nag_wp*(c*(x-x1))**4) + &
            c*(y-y1)
        End Select

        Return

      End Function fbnd
    End Module d06acfe_mod
    Program d06acfe

!     D06ACF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d06acf, d06baf, f16dnf, nag_wp
      Use d06acfe_mod, Only: fbnd, nin, nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: dnvint, radius, x0, x1, y0, y1
      Integer                              :: i, ifail, itrace, j, k, liwork,  &
                                              lrwork, maxind, maxval, ncomp,   &
                                              nedge, nedmx, nelt, nlines, nv,  &
                                              nvb, nvint, nvint2, nvmax,       &
                                              reftk, sdcrus
      Character (1)                        :: pmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: coor(:,:), coorch(:,:),          &
                                              crus(:,:), rate(:), rwork(:),    &
                                              weight(:)
      Real (Kind=nag_wp)                   :: ruser(5)
      Integer, Allocatable                 :: conn(:,:), edge(:,:), iwork(:),  &
                                              lcomp(:), lined(:,:), nlcomp(:)
      Integer                              :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: abs, real
!     .. Executable Statements ..
      Write (nout,*) 'D06ACF Example Program Results'

!     Skip heading in data file
      Read (nin,*)

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

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

      Read (nin,*) coorch(1,1:nlines)
      Read (nin,*) coorch(2,1:nlines)

!     The Lines of the boundary mesh

      Read (nin,*)(lined(1:4,j),rate(j),j=1,nlines)

      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

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

      Read (nin,*) ncomp
      Allocate (crus(2,sdcrus),nlcomp(ncomp),iwork(liwork),rwork(lrwork))

      j = 1

      Do i = 1, ncomp
        Read (nin,*) nlcomp(i)
        k = j + abs(nlcomp(i)) - 1
        Read (nin,*) lcomp(j:k)
        j = k + 1
      End Do

!     Data passed to the user-supplied function

      x0 = 1.5_nag_wp
      y0 = 0.0_nag_wp
      radius = 4.5_nag_wp
      x1 = 0.8_nag_wp
      y1 = -0.3_nag_wp

      ruser(1:5) = (/x0,y0,radius,x1,y1/)
      iuser(1) = 0

      itrace = 0

!     Call to the 2D 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)

      Write (nout,*)
      Read (nin,*) pmesh

      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
      Case Default
        Write (nout,*) 'Problem with the printing option Y or N'
        Go To 100
      End Select

      Deallocate (rwork,iwork)

!     Initialise mesh control parameters

      itrace = 0

!     Generation of interior vertices
!     for the wake of the first NACA

      nvint = 40
      lrwork = 12*nvmax + 30015
      liwork = 8*nedge + 53*nvmax + 2*nvb + 10078
      Allocate (weight(nvint),rwork(lrwork),conn(3,2*nvmax+5),iwork(liwork))

      nvint2 = 20
      dnvint = 5.0_nag_wp/real(nvint2+1,kind=nag_wp)

      Do i = 1, nvint2
        reftk = nvb + i
        coor(1,reftk) = 1.0_nag_wp + real(i,kind=nag_wp)*dnvint
        coor(2,reftk) = 0.0_nag_wp
      End Do

      weight(1:nvint2) = 0.05_nag_wp

!     for the wake of the second one

      dnvint = 4.19_nag_wp/real(nvint2+1,kind=nag_wp)

      Do i = nvint2 + 1, nvint
        reftk = nvb + i
        coor(1,reftk) = 1.8_nag_wp + real(i-nvint2,kind=nag_wp)*dnvint
        coor(2,reftk) = -0.3_nag_wp
      End Do

      weight((nvint2+1):nvint) = 0.05_nag_wp

!     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 characteristics'
        Write (nout,99999) 'NV   =', nv
        Write (nout,99999) 'NELT =', nelt
      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 d06acfe
