!   D06ACF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.
    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
!dec$ attributes offload : mic :: fbnd
       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)                  :: c, radius, x0, x1, y0, y1
!         .. Intrinsic Functions ..
          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

!dec$ attributes offload : mic :: fbnd
       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 Functions ..
       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 20
       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

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 d06acfe
