Example description
!   D06BAF Example Program Text
!   Mark 26.2 Release. NAG Copyright 2017.
    Module d06bafe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fbnd
!     .. Parameters ..
      Integer, Parameter, Public       :: 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)             :: radius2, x0, xa, xb, y0
!       .. Executable Statements ..
        xa = ruser(1)
        xb = ruser(2)
        x0 = ruser(3)
        y0 = ruser(4)

        fbnd = 0.0_nag_wp

        Select Case (i)
        Case (1)

!         line 1,2,3, and 4: ellipse centred in (X0,Y0) with
!         XA and XB as coefficients

          fbnd = ((x-x0)/xa)**2 + ((y-y0)/xb)**2 - 1.0_nag_wp
        Case (2)

!         line 7, lower arc on letter n, is a circle centred in (X0,Y0)
!         with radius SQRT(RADIUS2)

          x0 = 0.5_nag_wp
          y0 = 6.25_nag_wp
          radius2 = 20.3125_nag_wp
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        Case (3)

!         line 11, upper arc on letter n, is a circle centred in (X0,Y0)
!         with radius SQRT(RADIUS2)

          x0 = 1.0_nag_wp
          y0 = 4.0_nag_wp
          radius2 = 9.0_nag_wp + (11.0_nag_wp-y0)**2
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        Case (4)

!         line 15, upper arc on letter a, is a circle centred in (X0,Y0)
!         with radius SQRT(RADIUS2) touching point (5,11).

          x0 = 8.5_nag_wp
          y0 = 2.75_nag_wp
          radius2 = (x0-5.0_nag_wp)**2 + (11.0_nag_wp-y0)**2
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        Case (5)

!         line 25, lower arc on hat of 'a', is a circle centred in (X0,Y0)
!         with radius SQRT(RADIUS2) touching point (11,10).

          x0 = 8.5_nag_wp
          y0 = 4.0_nag_wp
          radius2 = 2.5_nag_wp**2 + (10.0_nag_wp-y0)**2
          fbnd = (x-x0)**2 + (y-y0)**2 - radius2
        Case (6)

!         lines 20, 21 and 22, belly of letter a, is an ellipse centered
!         in (X0, Y0) with semi-axes 3.5 and 2.75.

          x0 = 8.5_nag_wp
          y0 = 5.75_nag_wp
          fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
        Case (7)

!         lines 43, 44 and 45, outer curve on bottom of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 3.5 and 2.5.

          x0 = 17.5_nag_wp
          y0 = 2.5_nag_wp
          fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.5_nag_wp)**2 - 1.0_nag_wp
        Case (8)

!         lines 28, 29 and 30, inner curve on bottom of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 2.0 and 1.5.

          x0 = 17.5_nag_wp
          y0 = 2.5_nag_wp
          fbnd = ((x-x0)/2.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
        Case (9)

!         line 42, inner curve on lower middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.5 and 0.5.

          x0 = 17.5_nag_wp
          y0 = 5.5_nag_wp
          fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/0.5_nag_wp)**2 - 1.0_nag_wp
        Case (10)

!         line 31, outer curve on lower middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 2.0 and 1.5.

          x0 = 17.5_nag_wp
          y0 = 5.5_nag_wp
          fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
        Case (11)

!         line 41, inner curve on upper middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.0 and 1.0.

          x0 = 17.0_nag_wp
          y0 = 5.5_nag_wp
          fbnd = ((x-x0)/1.0_nag_wp)**2 + ((y-y0)/1.0_nag_wp)**2 - 1.0_nag_wp
        Case (12)

!         line 32, outer curve on upper middle of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 1.5 and 1.1573.

          x0 = 16.0_nag_wp
          y0 = 5.5_nag_wp
          fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/1.1573_nag_wp)**2 -          &
            1.0_nag_wp
        Case (13)

!         lines 33, 33, 34, 39 and 40, upper portion of 'g', is an ellipse
!         centered in (X0, Y0) with semi-axes 3.0 and 2.75.
          x0 = 17.0_nag_wp
          y0 = 9.25_nag_wp
          fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
        End Select

        Return

      End Function fbnd
    End Module d06bafe_mod
    Program d06bafe

!     D06BAF Example Main Program

!     .. Use Statements ..
      Use d06bafe_mod, Only: fbnd, nin, nout
      Use nag_library, Only: d06abf, d06acf, d06baf, f16dnf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: x0, xa, xb, xmax, xmin, y0, ymax,    &
                                          ymin
      Integer                          :: i, ifail, itrace, j, k, liwork,      &
                                          lrwork, maxind, maxval, ncomp,       &
                                          nedge, nedmx, nelt, nlines, npropa,  &
                                          nv, nvb, nvint, nvmax, reftk, sdcrus
      Character (1)                    :: pmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: coor(:,:), coorch(:,:), crus(:,:),   &
                                          rate(:), rwork(:), weight(:)
      Real (Kind=nag_wp)               :: ruser(4)
      Integer, Allocatable             :: conn(:,:), edge(:,:), iwork(:),      &
                                          lcomp(:), lined(:,:), nlcomp(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs
!     .. Executable Statements ..
      Write (nout,*) 'D06BAF Example Program Results'
      Flush (nout)

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

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

!     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

      Allocate (crus(2,sdcrus),iwork(liwork),rwork(lrwork))

!     The ellipse boundary which envelops the NAG Logo
!     the N, the A and the G

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

      Read (nin,*) crus(1,1:sdcrus)
      Read (nin,*) crus(2,1:sdcrus)

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

      Read (nin,*) ncomp
      Allocate (nlcomp(ncomp))

      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

      xmin = coorch(1,4)
      xmax = coorch(1,2)
      ymin = coorch(2,1)
      ymax = coorch(2,3)

      xa = (xmax-xmin)/2.0_nag_wp
      xb = (ymax-ymin)/2.0_nag_wp

      x0 = (xmin+xmax)/2.0_nag_wp
      y0 = (ymin+ymax)/2.0_nag_wp

      ruser(1:4) = (/xa,xb,x0,y0/)
      iuser(1) = 0

      itrace = -1

      Write (nout,*)
      Flush (nout)

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

      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

        Flush (nout)

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

      Deallocate (rwork,iwork)

!     Initialize mesh control parameters

      itrace = 0
      npropa = 1
      nvint = 0
      lrwork = 12*nvmax + 15
      liwork = 6*nedge + 32*nvmax + 2*nvb + 78
      Allocate (weight(nvint),rwork(lrwork),iwork(liwork),conn(3,2*nvmax+5))

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

      Select Case (pmesh)
      Case ('N')
        Write (nout,*) 'Complete mesh (via the 2D Delaunay-Voronoi'
        Write (nout,*) 'mesh generator) characteristics'
        Write (nout,99999) 'NV (rounded to nearest 10)   =', 10*((nv+5)/10)
        Write (nout,99999) 'NELT (rounded to nearest 10) =', 10*((nelt+5)/10)
      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

        Flush (nout)

      End Select

      Deallocate (rwork,iwork)
      lrwork = 12*nvmax + 30015
      liwork = 8*nedge + 53*nvmax + 2*nvb + 10078
      Allocate (rwork(lrwork),iwork(liwork))

!     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 (via the 2D Advancing front mesh'
        Write (nout,*) 'generator) characteristics'
        Write (nout,99999) 'NV (rounded to nearest 10)   =', 10*((nv+5)/10)
        Write (nout,99999) 'NELT (rounded to nearest 10) =', 10*((nelt+5)/10)
      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 d06bafe