*     D06ACF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NEDMX, NVMAX, NLINESX, NUS, NCOMPX, MAXCAN,
     +                 NVIMX, LRWORK, LIWORK
      PARAMETER        (NEDMX=200,NVMAX=2000,NLINESX=50,NUS=100,
     +                 NCOMPX=5,MAXCAN=10000,NVIMX=40,
     +                 LRWORK=12*NVMAX+3*MAXCAN+15,
     +                 LIWORK=8*NEDMX+55*NVMAX+MAXCAN+78)
*     .. Local Scalars ..
      DOUBLE PRECISION DNVINT, RADIUS, X0, X1, Y0, Y1
      INTEGER          I, IFAIL, ITRACE, J, K, NCOMP, NEDGE, NELT,
     +                 NLINES, NV, NVB, NVINT, NVINT2, REFTK
      CHARACTER        PMESH
*     .. Local Arrays ..
      DOUBLE PRECISION COOR(2,NVMAX), COORCH(2,NLINESX), COORUS(2,NUS),
     +                 RATE(NLINESX), RUSER(5), RWORK(LRWORK),
     +                 WEIGHT(NVIMX)
      INTEGER          CONN(3,2*NVMAX+5), EDGE(3,NEDMX), IUSER(1),
     +                 IWORK(LIWORK), LCOMP(NLINESX), LINE(4,NLINESX),
     +                 NLCOMP(NCOMPX)
*     .. External Functions ..
      DOUBLE PRECISION FBND
      EXTERNAL         FBND
*     .. External Subroutines ..
      EXTERNAL         D06ACF, D06BAF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS, DBLE
*     .. Executable Statements ..
*
      WRITE (NOUT,*) 'D06ACF Example Program Results'
      WRITE (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
*
      READ (NIN,*) (COORCH(1,J),J=1,NLINES)
*
      READ (NIN,*) (COORCH(2,J),J=1,NLINES)
*
*     The Lines of the boundary mesh
*
      READ (NIN,*) ((LINE(I,J),I=1,4),RATE(J),J=1,NLINES)
*
*     The number of connected components to the boundary
*     and their informations
*
      READ (NIN,*) NCOMP
      J = 1
      DO 20 I = 1, NCOMP
         READ (NIN,*) NLCOMP(I)
*
         READ (NIN,*) (LCOMP(K),K=J,J+ABS(NLCOMP(I))-1)
         J = J + ABS(NLCOMP(I))
   20 CONTINUE
*
      READ (NIN,*) PMESH
*
*     Data passed to the user-supplied function
*
      X0 = 1.5D0
      Y0 = 0.D0
      RADIUS = 4.5D0
      X1 = 0.8D0
      Y1 = -0.3D0
*
      RUSER(1) = X0
      RUSER(2) = Y0
      RUSER(3) = RADIUS
      RUSER(4) = X1
      RUSER(5) = Y1
      IUSER(1) = 0
*
      ITRACE = 0
*
*     Call to the 2D boundary mesh generator
*
      IFAIL = 0
*
      CALL D06BAF(NLINES,COORCH,LINE,FBND,COORUS,NUS,RATE,NCOMP,NLCOMP,
     +            LCOMP,NVMAX,NEDMX,NVB,COOR,NEDGE,EDGE,ITRACE,RUSER,
     +            IUSER,RWORK,LRWORK,IWORK,LIWORK,IFAIL)
*
      IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'Boundary mesh characteristics'
         WRITE (NOUT,99999) 'NVB   =', NVB
         WRITE (NOUT,99999) 'NEDGE =', NEDGE
      ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NVB, NEDGE
*
         DO 40 I = 1, NVB
            WRITE (NOUT,99997) I, COOR(1,I), COOR(2,I)
   40    CONTINUE
*
         DO 60 I = 1, NEDGE
            WRITE (NOUT,99996) I, EDGE(1,I), EDGE(2,I), EDGE(3,I)
   60    CONTINUE
      ELSE
         WRITE (NOUT,*) 'Problem with the printing option Y or N'
         STOP
      END IF
*
*     Initialise mesh control parameters
*
      ITRACE = 0
*
*     Generation of interior vertices
*     for the wake of the first NACA
*
      NVINT = 40
      NVINT2 = 20
      DNVINT = 5.D0/DBLE(NVINT2+1)
      DO 80 I = 1, NVINT2
         REFTK = NVB + I
         COOR(1,REFTK) = 1.D0 + DBLE(I)*DNVINT
         COOR(2,REFTK) = 0.D0
         WEIGHT(I) = 0.05D0
   80 CONTINUE
*
*     for the wake of the second one
*
      DNVINT = 4.19D0/DBLE(NVINT2+1)
      DO 100 I = NVINT2 + 1, NVINT
         REFTK = NVB + I
         COOR(1,REFTK) = 1.8D0 + DBLE(I-NVINT2)*DNVINT
         COOR(2,REFTK) = -0.3D0
         WEIGHT(I) = 0.05D0
  100 CONTINUE
*
*     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)
*
      IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'Complete mesh characteristics'
         WRITE (NOUT,99999) 'NV   =', NV
         WRITE (NOUT,99999) 'NELT =', NELT
      ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99998) NV, NELT
         DO 120 I = 1, NV
            WRITE (NOUT,99995) COOR(1,I), COOR(2,I)
  120    CONTINUE
*
         REFTK = 0
         DO 140 K = 1, NELT
            WRITE (NOUT,99994) CONN(1,K), CONN(2,K), CONN(3,K),
     +        REFTK
  140    CONTINUE
      END IF
*
      STOP
*
99999 FORMAT (1X,A,I6)
99998 FORMAT (1X,2I10)
99997 FORMAT (2X,I4,2(2X,E12.6))
99996 FORMAT (1X,4I4)
99995 FORMAT (2(2X,E12.6))
99994 FORMAT (1X,4I10)
      END
*
      DOUBLE PRECISION FUNCTION FBND(I,X,Y,RUSER,IUSER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION               X, Y
      INTEGER                        I
*     .. Array Arguments ..
      DOUBLE PRECISION               RUSER(*)
      INTEGER                        IUSER(*)
*     .. Local Scalars ..
      DOUBLE PRECISION               C, RADIUS, X0, X1, Y0, Y1
*     .. Intrinsic Functions ..
      INTRINSIC                      SQRT
*     .. Executable Statements ..
      FBND = 0.D0
      IF (I.EQ.1) THEN
*
*     upper NACA0012 wing beginning at the origin
*
         C = 1.008930411365D0
         FBND = 0.6D0*(0.2969D0*SQRT(C*X)-0.126D0*C*X-0.3516D0*(C*X)
     +          **2+0.2843D0*(C*X)**3-0.1015D0*(C*X)**4) - C*Y
      ELSE IF (I.EQ.2) THEN
*
*     lower NACA0012 wing beginning at the origin
*
         C = 1.008930411365D0
         FBND = 0.6D0*(0.2969D0*SQRT(C*X)-0.126D0*C*X-0.3516D0*(C*X)
     +          **2+0.2843D0*(C*X)**3-0.1015D0*(C*X)**4) + C*Y
      ELSE IF (I.EQ.3) THEN
         X0 = RUSER(1)
         Y0 = RUSER(2)
         RADIUS = RUSER(3)
         FBND = (X-X0)**2 + (Y-Y0)**2 - RADIUS**2
      ELSE IF (I.EQ.4) THEN
*
*     upper NACA0012 wing beginning at (X1;Y1)
*
         C = 1.008930411365D0
         X1 = RUSER(4)
         Y1 = RUSER(5)
         FBND = 0.6D0*(0.2969D0*SQRT(C*(X-X1))-0.126D0*C*(X-X1)
     +          -0.3516D0*(C*(X-X1))**2+0.2843D0*(C*(X-X1))
     +          **3-0.1015D0*(C*(X-X1))**4) - C*(Y-Y1)
      ELSE IF (I.EQ.5) THEN
*
*     lower NACA0012 wing beginning at (X1;Y1)
*
         C = 1.008930411365D0
         X1 = RUSER(4)
         Y1 = RUSER(5)
         FBND = 0.6D0*(0.2969D0*SQRT(C*(X-X1))-0.126D0*C*(X-X1)
     +          -0.3516D0*(C*(X-X1))**2+0.2843D0*(C*(X-X1))
     +          **3-0.1015D0*(C*(X-X1))**4) + C*(Y-Y1)
      END IF
*
      RETURN
      END