!   D06BAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
       FUNCTION fbnd(i,x,y,ruser,iuser)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 nag_library, ONLY : d06abf, d06acf, d06baf, f16dnf, nag_wp
       USE d06bafe_mod, ONLY : fbnd, nin, nout
!      .. 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 Functions ..
       INTRINSIC                              abs
!      .. Executable Statements ..
       WRITE (nout,*) 'D06BAF Example Program Results'
       FLUSH (nout)

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

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

       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 20
       END SELECT

       DEALLOCATE (rwork,iwork)

!      Initialise 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   =', 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

          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   =', 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

20     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