*     D06CAF Example Program Text
*     Mark 20 Release. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NBEDMX, NVMAX, NELTMAX, NVFIXMX, LNUME, LIWORK,
     +                 LRWORK
      PARAMETER        (NBEDMX=100,NVMAX=400,NELTMAX=2*NVMAX-1,
     +                 NVFIXMX=20,LNUME=3*NELTMAX,
     +                 LIWORK=2*NVMAX+5*NELTMAX+LNUME,
     +                 LRWORK=2*NVMAX+NELTMAX)
*     .. Local Scalars ..
      DOUBLE PRECISION DELTA, HX, HY, PI, R, RAD, SK, THETA, X1, X2, X3,
     +                 Y1, Y2, Y3
      INTEGER          I, IFAIL, IMAX, IND, ITRACE, J, JMAX, K, ME1,
     +                 ME2, ME3, NEDGE, NELT, NQINT, NV, NVFIX, REFTK
      CHARACTER        PMESH
*     .. Local Arrays ..
      DOUBLE PRECISION COOR(2,NVMAX), RWORK(LRWORK)
      INTEGER          CONN(3,NELTMAX), EDGE(3,NBEDMX), IWORK(LIWORK),
     +                 NUMFIX(NVFIXMX)
*     .. External Functions ..
      DOUBLE PRECISION G05DAF
      EXTERNAL         G05DAF
*     .. External Subroutines ..
      EXTERNAL         D06CAF, G05CBF
*     .. Intrinsic Functions ..
*
      INTRINSIC        ATAN, COS, DBLE, MIN, SIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D06CAF Example Program Results'
      WRITE (NOUT,*)
*
*     Skip heading in data file
*
      READ (NIN,*)
*
*     Read IMAX and JMAX, the number of vertices
*     in the x and y directions respectively.
*
      READ (NIN,*) IMAX, JMAX
*
*     Read distortion percentage and calculate radius
*     of distortion neighbourhood so that cross-over
*     can only occur at 100% or greater.
*
      READ (NIN,*) DELTA
*
      NV = IMAX*JMAX
      IF (NV.GT.NVMAX) THEN
         WRITE (NOUT,99999) 'Dimension problem NV MAX ', NV, NVMAX
         STOP
      END IF
*
      READ (NIN,*) PMESH
*
      HX = 1.D0/DBLE(IMAX-1)
      HY = 1.D0/DBLE(JMAX-1)
      RAD = 0.01D0*DELTA*MIN(HX,HY)/2.D0
      PI = 4.D0*ATAN(1.D0)
      CALL G05CBF(0)
      IND = 0
*
*     Generate a simple uniform mesh and then distort it
*     randomly within the distortion neighbourhood of each
*     node.
*
      DO 40 J = 1, JMAX
         DO 20 I = 1, IMAX

            R = G05DAF(0.D0,RAD)
            THETA = G05DAF(0.D0,2*PI)
            IF (I.EQ.1 .OR. I.EQ.IMAX .OR. J.EQ.1 .OR. J.EQ.JMAX)
     +          R = 0.D0

            K = (J-1)*IMAX + I
            COOR(1,K) = DBLE(I-1)*HX + R*COS(THETA)
            COOR(2,K) = DBLE(J-1)*HY + R*SIN(THETA)

            IF (I.LT.IMAX .AND. J.LT.JMAX) THEN
               IND = IND + 1
               CONN(1,IND) = K
               CONN(2,IND) = K + 1
               CONN(3,IND) = K + IMAX + 1
               IND = IND + 1
               CONN(1,IND) = K
               CONN(2,IND) = K + IMAX + 1
               CONN(3,IND) = K + IMAX
            END IF
   20    CONTINUE
   40 CONTINUE
*
      NELT = IND
*
      IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'The complete distorted mesh characteristics'
         WRITE (NOUT,99998) 'NV   =', NV
         WRITE (NOUT,99998) 'NELT =', NELT
      ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99997) NV, NELT
         DO 60 I = 1, NV
            WRITE (NOUT,99996) COOR(1,I), COOR(2,I)
   60    CONTINUE
      ELSE
         WRITE (NOUT,*) 'Problem with the printing option Y or N'
         STOP
      END IF
*
      REFTK = 0
      DO 80 K = 1, NELT
         ME1 = CONN(1,K)
         ME2 = CONN(2,K)
         ME3 = CONN(3,K)
*
         X1 = COOR(1,ME1)
         X2 = COOR(1,ME2)
         X3 = COOR(1,ME3)
         Y1 = COOR(2,ME1)
         Y2 = COOR(2,ME2)
         Y3 = COOR(2,ME3)
*
         SK = ((X2-X1)*(Y3-Y1)-(Y2-Y1)*(X3-X1))/2.D0
         IF (SK.LT.0.D0) THEN
            WRITE (NOUT,*)
     +        'Error the surface of the element is negative'
            WRITE (NOUT,99998) 'K = ', K
            WRITE (NOUT,99994) 'SK = ', SK
            STOP
         END IF
         IF (PMESH.EQ.'Y') WRITE (NOUT,99995) CONN(1,K), CONN(2,K),
     +       CONN(3,K), REFTK
   80 CONTINUE
*
*     Boundary edges
*
      NEDGE = 0
      DO 100 I = 1, IMAX - 1
         NEDGE = NEDGE + 1
         EDGE(1,NEDGE) = I
         EDGE(2,NEDGE) = I + 1
         EDGE(3,NEDGE) = 0
  100 CONTINUE
*
      DO 120 I = 1, JMAX - 1
         NEDGE = NEDGE + 1
         EDGE(1,NEDGE) = I*IMAX
         EDGE(2,NEDGE) = (I+1)*IMAX
         EDGE(3,NEDGE) = 0
  120 CONTINUE
*
      DO 140 I = 1, IMAX - 1
         NEDGE = NEDGE + 1
         EDGE(1,NEDGE) = IMAX*JMAX - I + 1
         EDGE(2,NEDGE) = IMAX*JMAX - I
         EDGE(3,NEDGE) = 0
  140 CONTINUE
*
      DO 160 I = 1, JMAX - 1
         NEDGE = NEDGE + 1
         EDGE(1,NEDGE) = (JMAX-I)*IMAX + 1
         EDGE(2,NEDGE) = (JMAX-I-1)*IMAX + 1
         EDGE(3,NEDGE) = 0
  160 CONTINUE
*
      NVFIX = 0
      NUMFIX(1) = 0
      ITRACE = 1
      NQINT = 10
      IFAIL = 0
*
*     Call the smoothing routine
*
      CALL D06CAF(NV,NELT,NEDGE,COOR,EDGE,CONN,NVFIX,NUMFIX,ITRACE,
     +            NQINT,IWORK,LIWORK,RWORK,LRWORK,IFAIL)
*
      IF (PMESH.EQ.'N') THEN
         WRITE (NOUT,*) 'The complete smoothed mesh characteristics'
         WRITE (NOUT,99998) 'NV   =', NV
         WRITE (NOUT,99998) 'NELT =', NELT
      ELSE IF (PMESH.EQ.'Y') THEN
*
*     Output the mesh to view it using the NAG Graphics Library
*
         WRITE (NOUT,99997) NV, NELT
         DO 180 I = 1, NV
            WRITE (NOUT,99996) COOR(1,I), COOR(2,I)
  180    CONTINUE
*
         REFTK = 0
         DO 200 K = 1, NELT
            WRITE (NOUT,99995) CONN(1,K), CONN(2,K), CONN(3,K), REFTK
  200    CONTINUE
      END IF
*
      STOP
*
99999 FORMAT (1X,A,2I6)
99998 FORMAT (1X,A,I6)
99997 FORMAT (1X,2I10)
99996 FORMAT (2(2X,E12.6))
99995 FORMAT (1X,4I10)
99994 FORMAT (1X,A,E12.6)
      END