      SUBROUTINE TRICP (XD,YD,ZD,ND,C,NC,WK,IWK,MODE)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     USER INTERFACE   (VERSION 0.1)
C
C     COPYRIGHT(C): A. PREUSSER, 1984
C                   SEN. F. WISS. U. FORSCH.
C                   PROJEKTGRUPPE PARALLELRECHNEN
C                   HEILBRONNER STR. 10
C                   D-1000 BERLIN 31
C
      DIMENSION XD(ND),YD(ND),ZD(ND),C(NC),WK(*),IWK(*)
C
C     XD,YD,ZD   COORDINATES OF DATA POINTS
C                XD AND YD MUST BE GIVEN IN CM.
C                FOR USING INCH,  SET INSTALLATION PARAMETER
C                CMSCAL= 2.54 .
C     ND         NUMBER OF DATA POINTS
C     C          CONTOUR LEVELS
C     NC         NUMBER OF CONTOUR LEVELS
C     WK         REAL WORK ARRAY, DIMENSION 5*ND
C
C                ON EXIT, WK CONTAINS THE PARTIAL DERIVATIVES
C                (ZX(I),ZY(I),ZXX(I),ZXY(I),ZYY(I),I=1,ND)
C
C     IWK        INTEGER WORK ARRAY. FOR MODE=0, THE DIMENSION IS
C                                  31*ND  .
C                FOR OTHER MODES, SEE BELOW.
C
C                ON EXIT, IWK CONTAINS
C                A) IWK(1)= NT= NUMBER OF TRIANGLES
C                B) IWK(2...3*NT+1) POINT NUMBERS FOR THE TRIANGLES.
C                   THE FIRST 3 NUMBERS DETERMINE THE VERTICES OF THE
C                   FIRST TRIANGLE, THE NEXT 3 FOR THE SECOND AND SO ON.
C                   THE NUMBERS COORESPOND TO THE INDICES OF THE
C                   XD,YD,ZD ARRAYS. THEY ARE ARRANGED COUNTER-
C                   CLOCKWISE WITHIN A TRIANGLE.
C
C                C) IWK(3*NT+2  ...  3*NT+NCP*ND+1)
C                   A LIST OF NCP*ND POINT NUMBERS, REPRESENTING
C                   NCP POINTS FOR EACH DATA POINT THAT ARE
C                   THE CLOSEST TO THIS POINT. THE FIRST NCP
C                   NUMBERS ARE FOR THE FIRST DATA POINT, THE NEXT
C                   NCP FOR THE SECOND AND SO ON. THESE NUMBERS
C                   WERE USED FOR THE COMPUTATION OF THE PARTIAL
C                   DERIVATIVES.
C                   NCP IS AN INSTALLATION PARAMETER AND WILL BE SET
C                   TO 4.
C
C     MODE       MODE OF USAGE
C                0, NORMAL MODE, SEE ABOVE.
C                1, NO TRIANGULATION REQUESTED.
C                   IWK MUST CONTAIN THE INFORMATION ABOUT
C                   THE TRIANGLES AS DESCRIBED UNDER A) AND B)
C                   ON ENTRY.
C                   THESE LOCATIONS OF IWK WILL NOT BE CHANGED
C                   AND IN ADDITION, THE INFORMATION DESCRIBED
C                   UNDER C) WILL BE AVAILABLE ON EXIT.
C                2, NO TRIANGULATION AND NO DETERMINATION OF THE
C                   NCP CLOSEST POINTS FOR EACH DATA POINT.
C                   IWK MUST CONTAIN THE INFORMATION AS DESCRIBED
C                   UNDER A), B), AND C) AS INPUT INFORMATION.
C                   THE CONTENTS OF IWK WILL NOT BE CHANGED.
C                3, NO TRIANGULATION AND NO COMPUTATION OF THE
C                   PARTIAL DERIVATIVES.
C                   IWK MUST CONTAIN THE INFORMATION A) AND B)
C                   AND
C                   WK MUST CONTAIN THE PARTIAL DERIVATIVES
C                   (ZX(I),ZY(I),ZXX(I),ZXY(I),ZYY(I),I=1,ND)
C                   ON ENTRY.
C                   THE CONTENTS OF WK AND IWK WILL NOT BE CHANGED.
C                   THIS MODE IS ESPECIALLY USEFUL WHEN TRICP IS
C                   CALLED AGAIN AFTER A PREVIOUS CALL. FOR INSTANCE,
C                   ONLY THE CONTOUR LEVELS MAY HAVE CHANGED.
C                   WHEN DESIGNING A SURFACE, IT MAY BE APPROPRIATE
C                   TO CHANGE THE XD,YD,ZD PARAMETERS AND THE PARTIAL
C                   DERIVATIVES INTERACTIVELY AND TO CALL TRICP AGAIN
C                   USING THIS MODE.
C
C                   FOR MODES 1,2, AND 3 THE DIMENSION OF IWK CAN
C                   BE REDUCED TO
C                                 3*NT + NCP*ND + 1.
C
C
      COMMON /TRPCOC/ DSPMIN,CMSCAL,NPMAX
C
C *** DSPMIN     RESOLUTION OF PLOTTING DEVICE IN USE
C                (MINIMAL DISTANCE OF TWO POINTS TO BE PLOTTED)
C *** CMSCAL     VARIABLE FOR SWITCHING BETWEEN CM AND INCH
C     NPMAX      MAXIMUM NUMBER OF POINTS PER LINE FOR A TRIANGLE
C
C
C     SET INSTALLATION PARAMETERS
      NCP= 4
      CMSCAL= 1.d0
C *** SET  CMSCAL= 2.54  FOR INCH CALIBRATED PLOTTERS
      DSPMIN= 0.02d0/CMSCAL
      NPMAX= 500
C
C     SET STARTING ADDRESSES FOR IWK WORK ARRAY
      N1= 2
      N2= N1 + 6*ND - 15
      IF (MODE.GE.1) N2= N1 + 3*IWK(1)
      N3= N2 + 6*ND
      N4= N3 + 18*ND
C
      IF (MODE.EQ.1) GOTO 50
      IF (MODE.EQ.2) GOTO 60
      IF (MODE.EQ.3) GOTO 70
C
C     TRIANGULATION
      CALL IDTANG (ND,XD,YD,IWK,IWK(N1),NL,IWK(N2),IWK(N3),IWK(N4),WK)
C     IDTANG IS PART OF ACM ALG. 526 BY H. AKIMA
      N2= N1+3*IWK(1)
C
C     FIND NCP POINTS CLOSEST TO EACH DATA POINT
   50 CONTINUE
      CALL IDCLDP (ND,XD,YD,NCP,IWK(N2))
C     IDCLDP IS PART OF ACM ALG. 526 BY H. AKIMA
C
C     COMPUTE PARTIAL DERIVATIVES
   60 CONTINUE
      CALL IDPDRV (ND,XD,YD,ZD,NCP,IWK(N2),WK)
C     IDPDRV IS PART OF ACM ALG. 526 BY H. AKIMA
C
C     COMPUTE CONTOURS BY SUCCESSIVE SOLUTION OF
C     QUINTIC POLYNOMIAL EQUATIONS
   70 CONTINUE
      CALL TRP00 (XD,YD,ZD,C,NC,IWK(N1),WK,IWK(1))
C
      RETURN
      END



      SUBROUTINE TRP00 (XD,YD,ZD,CN,NC,IPT,PDD,NT)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     MAIN ROUTINE
C
C     COPYRIGHT(C): A. PREUSSER, 1984
C                   SEN. F. WISS. U. FORSCH.
C                   PROJEKTGRUPPE PARALLELRECHNEN
C                   HEILBRONNER STR. 10
C                   D-1000 BERLIN 31
C
      DIMENSION XD(*),YD(*),ZD(*),CN(NC),IPT(*),PDD(*)
C
C     XD,YD,ZD  COORDINATES OF DATA POINTS
C     CN        CONTOUR LEVELS
C     NC        NUMBER OF CONTOUR LEVELS
C     IPT       POINT NUMBERS STORED AS 3*I-2, 3*I-1, 3*I TH
C               ELEMENT FOR TRIANGLE I
C     PDD       PARTIAL DERIVATIVES
C               (ZX,ZY,ZXX,ZXY,ZYY)
C
C
C     FOR PLOTTING THE CALCOMP ROUTINE PLOT(X,Y,IPEN) IS USED
C     IPEN= 3 MEANS     MOVE TO X,Y WITH PEN IN UPWARD POSITION
C     IPEN= 2 MEANS     MOVE TO X,Y WITH PEN DOWN
C
C
      COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
     1,               SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
     2,               AP,BP,CP,DP
C
C     X,Y,Z      COORDINATES AT THE THREE VERTICES OF A TRIANGLE
C     DX,DY      DIFFERENCES OF X,Y
C     DX2,DY2    (DIFFERENCES OF X,Y)**2
C     SL         SIDE LENGTHS
C     SL2        (SIDE LENGTHS)**2
C     ZT(IV,K)   FIRST PARTIAL DERIVATIVE WITH RESPECT TO VARIABLE IV
C                (U,V,W FOR IV=1,2,3) AT POINT K
C     ZTT        SECOND PARTIAL DERIVATIVES
C     ZUV        MIXED PARTIAL DERIVATIVES
C     AP,BP,CP,DP  CONSTANTS DERIVED FROM DX,DY
C                NAMES DUE TO ALG. 526 CACM
C
      COMMON /TRPCOC/ DSPMIN,CMSCAL,NPMAX
C     FOR EXPLANATION OF VARIABLES IN /TRPCOC/ SEE SUBROUTINE TRICP
C
      COMMON /TRPCOF/ KK,KSE,XX4F,YY4F,SIR,COR,CL
C
C     /TRPCOF/ CONTAINS VARIABLES WHICH ARE PASSED TO FUNCTION
C              TRICPF AS PARAMETERS
C
C     KK          NUMBER OF FUNCTION TO BE EVALUATED BY TRICPF
C     KSE         ACTUAL SIDE NUMBER
C     XX4F,YY4F   COORDINATES FOR POINT P4F (PRELEMINARY POSITION
C                 OF POINT P4)
C     SIR,COR     COSINUS OF DIRECTION NORMAL TO CURVE DIRECTION
C     CL          ACTUAL SCALED CONTOUR LEVEL
C
      DIMENSION TS1(6,3),Z1(6,3),ZMAX(3),ZMIN(3),TS2(6),IN(3)
     1,         TZR(5,3),X0(5,3),Y0(5,3),NZ(3),TPER(3),DT(3),THETAS(3)
     2,         SI(3),CO(3)
C
C     TS1          U,V,W-VALUES AT ENDPOINTS OF INTERVALS
C     Z1           Z-VALUES AT ENDPOINTS OF INTERVALS
C     ZMAX,ZMIN    MAXIMUM AND MINIMUM Z VALUE ALONG A SIDE
C     TS2          WORKING ARRAY FOR COMPUTING TS1
C     IN           NUMBER OF INTERVALS
C     TZR          U,V,W VALUES FOR ZEROS FOUND ON SIDES
C     X0,Y0        X,Y COORDINATES FOR ZEROS FOUND
C     NZ           NUMBER OF ZEROS FOUND
C     TPER         PERMITTED POSITION ERROR FOR ZEROS ON SIDES
C     DT           DIFFERENCE IN U,V,W FOR ESTIMATING STARTING
C                  DIRECTION OF A CONTOUR LINE
C     THETAS       ANGLE OF DIRECTION FOR SIDES
C     SI,CO        COSINUS OF DIRECTION FOR SIDES
C
      PI= 4.d0*ATAN2(1.d0,1.d0)
C
C     START MAIN LOOP OVER TRIANGLES
      DO 5000  IT=1,NT
        icb_found = 0
C
C     LOAD COORDINATES AT VERTICES
      JI= 3*(IT-1)
      DO 50 J=1,3
        JI= JI + 1
        ID= IPT(JI)
        X(J)= XD(ID)
        Y(J)= YD(ID)
        Z(J)= ZD(ID)
   50 CONTINUE
      ZS= (Z(1)+Z(2)+Z(3))/3.d0
C
C     SOME BASIC GEOMETRY FOR THE TRIANGLE
      DO 60 J=1,3
        Z(J)= Z(J) - ZS
        NP1= 3
        NP2= 2
        IF (J.EQ.3) NP1= 1
        IF (J.EQ.2) NP2= 1
        DX(J)= X(NP2) - X(NP1)
        DY(J)= Y(NP2) - Y(NP1)
        DX2(J)= DX(J)*DX(J)
        DY2(J)= DY(J)*DY(J)
        SL2(J)= DX2(J) + DY2(J)
        SL(J)= SQRT(SL2(J))
        CO(J)= DX(J)/SL(J)
        SI(J)= DY(J)/SL(J)
   60 CONTINUE
      CO(1)= -CO(1)
      SI(1)= -SI(1)
      SLMAX= MAX(SL(1),SL(2),SL(3))
      DI1= ABS(DY(3)*CO(1) - DX(3)*SI(1))
      DI2= ABS(DY(1)*CO(2) - DX(1)*SI(2))
      DI3= ABS(DY(2)*CO(3) - DX(2)*SI(3))
      HMIN= MIN(DI1,DI2,DI3)
C         = SHORTEST DISTANCE BETWEEN A VERTEX AND
C           ITS OPPOSITE SIDE (MINIMUM HEIGHT)
C
C *** ISSUE A WARNING MESSAGE HERE, IF HMIN IS LESS THAN,
C     SAY 0.01/CMSCAL
C
C     DEFINE CONSTANTS DEPENDING ON HMIN
C
      RMAX  = MIN(0.02d0/CMSCAL,HMIN*0.01d0)
C           = DISTANCE NORMAL TO CURVE DIRECTION WITHIN WHICH
C             A ZERO MUST BE FOUND
      RMIN  = RMAX*0.2d0
C             IF ZERO HAS BEEN FOUND WITHIN A DISTANCE SMALLER
C             THAN RMIN, THEN STEPSIZE DS IS MULTIPLIED BY 2.
      POSERR= MIN(1.d-3/CMSCAL,1.d-3*HMIN*HMIN/SLMAX)
C           = PERMITTED POSITION ERROR
      DSMAX = HMIN*0.2d0
C           = MAXIMUM STEP SIZE
      FSTEP = RMAX*10.d0
C           = STARTING STEP SIZE
      EPS   = HMIN*0.01d0
C           = DIFFERENCE FOR ESTIMATING STARTING DIRECTION
      RMAX2 = RMAX*2.d0
C           = DISTANCE NORMAL TO CURVE DIRECTION. A ZERO MUST BE
C             FOUND WITHIN WHEN CROSSING TRIANGLE BORDER.
C
C     COMPUTE COEFFICIENTS FOR POLYNOMIALS ALONG SIDES
      CALL TRP001 (IT,IPT,PDD)
C
      CL= 0.d0
C
C     LOOP OVER SIDES
      DO 200 JSE=1,3
        KSE= JSE
        TPER(KSE)= POSERR/SL(KSE)
C
C       COMPUTE ENDPOINTS OF INTERVALS
        TS1(1,KSE)= 0.d0
        TS1(2,KSE)= 1.d0
        I= 2
        DO 150 K=1,4
          KK= K
          TS2(1)= 0.d0
          II= 2
          TB= 0.d0
          F2= TRICPF(TB)
          DO 100 J=2,I
            TA= TB
            F1= F2
            TB= TS1(J,KSE)
            F2= TRICPF(TB)
            IF (F1*F2.GT.0.) GOTO 100
            IF (F1.EQ.0. .AND. F2.EQ.0.d0) GOTO 100
            TS2(II)= TRICPZ(TA,TB,F1,F2,TPER(KSE))
            II= II + 1
  100     CONTINUE
          TS2(II)= 1.d0
          DO 120 J=1,II
  120       TS1(J,KSE)= TS2(J)
          I= II
  150   CONTINUE
C       IN(KSE)= NUMBER OF INTERVALS
        I= I-1
        IN(KSE)= I
C       (E.G. IF IN(KSE)=1, THERE IS NO POINT FOR WHICH 1ST DER.=0)
C
C       COMPUTE MAXIMA AND MINIMA FOR EACH SIDE
        NP1= 3
        NP2= 2
        IF (KSE.EQ.3) NP1=1
        IF (KSE.EQ.2) NP2=1
        ZMAX(KSE)= MAX(Z(NP1),Z(NP2))
        ZMIN(KSE)= MIN(Z(NP1),Z(NP2))
        Z1(1,KSE)= Z(NP1)
        Z1(I+1,KSE)= Z(NP2)
        IF (I.EQ.1) GOTO 170
        KK= 5
        DO 160 J=2,I
          Z1(J,KSE)= TRICPF(TS1(J,KSE))
          IF (Z1(J,KSE).GT.ZMAX(KSE)) ZMAX(KSE)= Z1(J,KSE)
          IF (Z1(J,KSE).LT.ZMIN(KSE)) ZMIN(KSE)= Z1(J,KSE)
  160   CONTINUE
  170   CONTINUE
  200 CONTINUE
C
      NPTRI= 0
      KA= 1
      KE= 1
      KC= 0
C
C     START LOOP OVER CONTOUR LEVELS
      DO 4000 K=1,NC
C
        CL= CN(K) - ZS
C
C       COMPUTE ZEROS (IF ANY) ON THE THREE SIDES FOR LEVEL CL
        KK= 5
        DO 1200 JSE=1,3
          NZ(JSE)= 0
          IF (CL.LT.ZMIN(JSE) .OR. CL.GT.ZMAX(JSE)) GOTO 1200
          KSE= JSE
          JN=0
          NI= IN(KSE)
C         CHECK INTERVALS FOR ZEROS
          DO 1100 J= 1,NI
            F1= Z1(J,KSE) - CL
            F2= Z1(J+1,KSE) - CL
            F1F2= F1*F2
            IF (F1F2.GT.0.) GOTO 1100
            IF (F1F2.LT.0.) GOTO 1090
C
C           SPECIAL SITUATIONS
            IF (NI.EQ.1 .AND. F1.EQ.0. .AND. F2.EQ.0.) GOTO 1070
            IF (J.EQ.1  .AND. F1.EQ.0. .OR.
     A          J.EQ.NI .AND. F2.EQ.0.) GOTO 1080
            GOTO 1090
 1070       CONTINUE
C           CONTOURLINE = SIDE KSE
            NP1= 3
            NP2= 2
            IF (KSE.EQ.3) NP1= 1
            IF (KSE.EQ.2) NP2= 1
            CALL PLOT_CALLBACK (X(NP1),Y(NP1),3 ,icb_found ,cn(k))
            CALL PLOT_CALLBACK (X(NP2),Y(NP2),2 ,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
            GOTO 1100
 1080       CONTINUE
C           LINE PASSES THROUGH DATA POINT
C           ONLY ONE ZERO FOR THIS TRIANGLE, SKIP EITHER SIDE
            IF (KSE.EQ.3 .OR. (F1.EQ.0. .AND. KSE.EQ.2)) GOTO 1100
            JN= JN+1
            IF (F2.EQ.0.) TZR(JN,KSE)= 1.-TPER(KSE)*0.5d0
            IF (F1.EQ.0.) TZR(JN,KSE)= TPER(KSE)*0.5d0
            GOTO 1100
C
 1090       JN= JN+1
            TZR(JN,KSE)= TRICPZ(TS1(J,KSE),TS1(J+1,KSE),F1,F2,TPER(KSE))
 1100     CONTINUE
          NZ(KSE)= JN
C                = NUMBER OF ZEROS FOR SIDE KSE AND CONTOUR LEVEL CL
 1200   CONTINUE
C
        IF (NZ(1)+NZ(2)+NZ(3).LT.2) GOTO 4000
C       IF (.TRUE.) THEN STOP PROCESSING THIS CONTOUR LEVEL
C
C       COMPUTE X0,Y0 FOR EACH ZERO (RELATIVE TO X(3),Y(3))
        DO 1300 JSE=1,3
          JN= NZ(JSE)
          IF (JN.EQ.0) GOTO 1300
          DO 1280 JZR=1,JN
            T= TZR(JZR,JSE)
            IF (JSE.EQ.3) GOTO 1230
            X0(JZR,JSE)= DX(JSE)*T
            Y0(JZR,JSE)= DY(JSE)*T
            GOTO 1280
 1230       X0(JZR,JSE)= DX(2) + DX(3)*T
            Y0(JZR,JSE)= DY(2) + DY(3)*T
 1280     CONTINUE
 1300   CONTINUE
C
        KC= KC+1
        IF (KC.NE.1) GOTO 1350
C       COMPUTE DIRECTION OF SIDES
        THETAS(1)= ATAN2(-DY(1),-DX(1))
        THETAS(2)= ATAN2( DY(2), DX(2))
        THETAS(3)= ATAN2( DY(3), DX(3))
C       COMPUTE DIFFERENCES FOR ESTIMATING START DIRECTION
        DT(1)= -EPS/SL(1)
        DT(2)=  EPS/SL(2)
        DT(3)=  EPS/SL(3)
C       COMPUTE COEFFICIENTS FOR POLYNOMIAL INSIDE TRIANGLE
        CALL TRP002
 1350   CONTINUE
C
C       FOLLOW CONTOURS
C
        RMA= RMAX
C       LOOP OVER SIDES, KSE= SIDE NUMBER
        DO 3000 JSE=1,3
          KSE= KA
          JN= NZ(KSE)
          IF (JN.EQ.0) GOTO 2010
          RMA= SIGN (RMAX,-RMA)
C         LOOP OVER ZEROS
          DO 2000 JZR= 1,JN
            IF (TZR(JZR,KSE).GT.5.) GOTO 2000
C           IF ()  THEN THIS ZERO HAS ALREADY BEEN PROCESSED
C
C           ESTIMATE STARTING DIRECTION
C
C           COMPUTE F1 ON SIDE USING FIRST DERIVATIVE FD1
            KK=4
            FD1= TRICPF(TZR(JZR,KSE))
            F1= FD1 * DT(KSE)
C
C           COMPUTE F2 NORMAL TO SIDE
            KK= 6
            XX4F= X0(JZR,KSE)
            YY4F= Y0(JZR,KSE)
            SIR= CO(KSE)
            COR= - SI(KSE)
            F2= TRICPF(EPS)
C
C           COMPUTE ANGLE
            IF (F2.EQ.0.d0) GOTO 1470
            THETSC= ATAN2 (-F1,F2)
            IF (THETSC.LE.0.d0) THETSC= THETSC + PI
            GOTO 1480
 1470       CONTINUE
            IF (FD1.EQ.0.d0) GOTO 2000
            THETSC= PI*0.5d0
 1480       THETAC= THETAS(KSE) + THETSC
C
C           MOVE PEN TO START
            CALL PLOT_CALLBACK (X0(JZR,KSE)+X(3),Y0(JZR,KSE)+Y(3),3
     $           ,icb_found ,cn(k))
C
C           INITIALIZE TRACING
            DS= FSTEP
            DX12= DS*COS(THETAC)
            DY12= DS*SIN(THETAC)
            XX3= XX4F
            YY3= YY4F
            XX2= XX3 - DX12
            YY2= YY3 - DY12
            XX1= XX2 - DX12
            YY1= YY2 - DY12
            DS23= DS
            DS12= DS
C           POINTS P1,P2,P3,P4 WITH COORDINATES XX1...XX4, YY1...YY4
C           ARE REFERRED TO AS *QUEUE*. DS12...DS34 ARE THE DISTANCES
C           BETWEEN POINTS IN THE QUEUE.
C           EVERY POINT NORMALLY MUST PASS THROUGH THE QUEUE
C           BEFORE IT IS PLOTTED.
C           FIRST (OLDEST) POINT IN THE QUEUE IS P1, WHICH IS NORMALLY
C           (NPQ = 4) THE FIRST CANDIDATE TO BE PLOTTED.
C           P4 IS THE NEXT POINT TO BE COMPUTED.
            NPQ= 0
C           NPQ= NUMBER OF POINTS IN THE QUEUE COUNTED FROM THE END
C                WHICH CAN BE PLOTTED.
C                E.G. NPQ=1, P4 IS STILL TO BE PLOTTED
            NP=0
C             = NUMBER OF POINTS COMPUTED FOR THIS CONTOUR LINE
            DSP=0.d0
C              = DISTANCE OF POINT TO BE PLOTTED TO LAST POINT PLOTTED
            NSTOP= 0
            NFOUND= 0
            NOST= 0
C
C           COMPUTE NEW POINT FOR CONTOUR LINE
C
C           COMPUTE CURVE DIRECTION
 1500       CONTINUE
            DS13= DS23 + DS12
            PL0=  DS23/(DS12*DS13)
            PL1= -DS13/(DS12*DS23)
            PL2= (DS13+DS23)/(DS13*DS23)
            DXDS= PL0*XX1 + PL1*XX2 + PL2*XX3
            DYDS= PL0*YY1 + PL1*YY2 + PL2*YY3
            SQ= SQRT(DXDS*DXDS+DYDS*DYDS)
            DXDS= DXDS/SQ
            DYDS= DYDS/SQ
            COR= -DYDS
            SIR=  DXDS
C
C           SEARCH FOR TWO POINTS WITH OPPOSITE SIGN
            RMA= SIGN (RMAX,RMA)
 1550       CONTINUE
            XX4F= XX3 + DXDS*DS
            YY4F= YY3 + DYDS*DS
            F1= TRICPF(0.d0)
            F2= TRICPF(RMA)
            IF (F1*F2 .LE. 0.d0) GOTO 1600
            IF (ABS(F2) .LT. ABS(F1)) GOTO 1560
            RMA= -RMA
            F2= TRICPF(RMA)
            IF (F1*F2 .LE. 0.d0) GOTO 1600
C
C           DIVIDE STEPSIZE IN CURVE DIRECTION BY 2.
 1560       CONTINUE
            DS= DS*0.5d0
            IF (DS.GT.POSERR) GOTO 1550
C
C           DIVIDE STEPSIZE NORMAL TO CURVE BY 2.
            NOST= NOST + 1
            DS= DS*2.d0
            RMA= RMA*0.5d0
 1570       F2= TRICPF(RMA)
            IF (F1*F2.LE.0.d0) GOTO 1600
            RMA= -RMA
            F2= TRICPF(RMA)
            IF (F1*F2.LE.0.d0) GOTO 1600
            RMA= RMA*0.5d0
            IF (ABS(RMA).GT.POSERR) GOTO 1570
            NSTOP= 1
            GOTO 1900
C
C
C           FIND ZERO FOR NEW POINT
 1600       NPQ= NPQ + 1
            NP= NP + 1
            IF (NP.GT.NPMAX) GOTO 1990
            R= TRICPZ(0.d0,RMA,F1,F2,POSERR)
            DS34= SQRT(DS*DS + R*R)
            XX4= XX4F + COR*R
            YY4= YY4F + SIR*R
            U= AP*XX4  + BP*YY4
            V= CP*XX4  + DP*YY4
C
C           CHECK IF POINT IS OUTSIDE THE TRIANGLE
            KE=1
            IF (U.LT.0.d0) GOTO 1700
            KE=2
            IF (V.LT.0.d0) GOTO 1700
            KE= 3
            IF (V .GT. 1.d0-U) GOTO 1700
C
C           POINT IS INSIDE
C           PLOT LAST POINT OF QUEUE AND UPDATE QUEUE
            IF (NPQ.NE.4) GOTO 1610
C ***       INSERT CALL FOR LABELING ROUTINE HERE
C ***       USING XX1...XX4, YY1...YY4, DS12...DS34 AS PARAMETERS
C ***       AND REDUCE NPQ BY THE NUMBER OF POINTS WHICH HAVE BEEN
C ***       USED FROM THE QUEUE FOR LABELING
            NPQ= NPQ -1
            DSP= DSP + DS01
            IF (DSP.LT.DSPMIN) GOTO 1610
            CALL PLOT_CALLBACK (XX1+X(3),YY1+Y(3),2 ,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
            DSP= 0.d0
 1610       DS01= DS12
            DS12= DS23
            DS23= DS34
            XX1= XX2
            YY1= YY2
            XX2= XX3
            YY2= YY3
            XX3= XX4
            YY3= YY4
C           SET NEW STEP SIZE
            IF (ABS(R).LT.RMIN) DS= MIN(DSMAX,DS*2.d0)
            GOTO 1500
C
C           POINT IS OUTSIDE
C
 1700       CONTINUE
C           NO SEARCH IF FIRST POINT WAS COMPUTED AND
C           START WAS AT DATA POINT
            IF (NP.NE.1) GOTO 1710
            IF (TZR(JZR,KSE).GT.1.d0-TPER(KSE) .OR.
     A          TZR(JZR,KSE).LT.TPER(KSE)) GOTO 2000
 1710       CONTINUE
C
C           FLAG ZERO WHERE LINE STARTED
            TZR(JZR,KSE)= TZR(JZR,KSE) + 10.d0
C
C           SEARCH FOR CORRESPONDING ZERO ON SIDE KE
            CHECK= 5.d0
            DIST= 99999.d0
            DO 1790 N=1,2
              DO 1770 JESE=1,3
                JJ= NZ(KE)
                IF (JJ.EQ.0) GOTO 1760
                DO 1750 J=1,JJ
                  IF (TZR(J,KE).GT.CHECK) GOTO 1750
                  DI1= ABS((Y0(J,KE)-YY3)*DXDS - (X0(J,KE)-XX3)*DYDS)
                  IF (DI1.GE.DIST) GOTO 1750
                  J1= J
                  DIST= DI1
 1750           CONTINUE
 1760           IF (DIST.LT.RMAX2) GOTO 1800
                KE= MOD(KE,3) + 1
 1770         CONTINUE
C             FOR N=2, DISREGARD FLAGS
              CHECK= 15.d0
 1790       CONTINUE
C           NO ACCEPTABLE ZERO ON ALL THREE SIDES
            NFOUND= 1
            GOTO 1900
C
C           FLAG ZERO FOUND
 1800       TZR(J1,KE)= TZR(J1,KE) + 10.d0
C
C           CLEAR QUEUE AND FINISH CONTOUR LINE
 1900       CONTINUE
            IF (NPQ.GT.3) GOTO 1910
            IF (NPQ.GT.2) GOTO 1920
            IF (NPQ.GT.1) GOTO 1930
            GOTO               1940
 1910       DSP= DSP + DS01
            IF (DSP.LT.DSPMIN) GOTO 1920
            CALL PLOT_CALLBACK (XX1+X(3),YY1+Y(3),2,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
            DSP= 0.d0
 1920       DSP= DSP + DS12
            IF (DSP.LT.DSPMIN) GOTO 1930
            CALL PLOT_CALLBACK (XX2+X(3),YY2+Y(3),2,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
            DSP= 0.d0
 1930       DSP= DSP + DS23
            IF (DSP.LT.DSPMIN) GOTO 1940
            CALL PLOT_CALLBACK (XX3+X(3),YY3+Y(3),2,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
 1940       IF (NFOUND+ NSTOP .EQ.0)
     A      CALL PLOT_CALLBACK (X0(J1,KE)+X(3),Y0(J1,KE)+Y(3),2
     $           ,icb_found ,cn(k))
            IF (icb_found .EQ. 1) RETURN 
C
 1990       CONTINUE
            NPTRI= NPTRI + NP
 2000     CONTINUE
C         END OF LOOP OVER ZEROS
 2010     CONTINUE
          KA= MOD(KA,3) + 1
 3000   CONTINUE
C       END OF LOOP OVER SIDES
        KA= KE
 4000 CONTINUE
C     END OF LOOP OVER CONTOUR LEVELS
 5000 CONTINUE
C     END OF LOOP OVER TRIANGLES
      RETURN
      END



      SUBROUTINE TRP001 (IT0,IPT,PDD)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     COMPUTATION OF COEFFICIENTS FOR POLYNOMIALS ALONG SIDES
C
C     COPYRIGHT(C): A. PREUSSER, 1984
C                   SEN. F. WISS. U. FORSCH.
C                   PROJEKTGRUPPE PARALLELRECHNEN
C                   HEILBRONNER STR. 10
C                   D-1000 BERLIN 31
C
      DIMENSION IPT(*),PDD(*)
C
C     IT0     NUMBER OF TRIANGLE IN USE
C     IPT     POINT NUMBERS STORED AS 3*I-2, 3*I-1, 3*I TH
C             ELEMENT FOR TRIANGLE I
C     PDD     PARTIAL DERIVATIVES AT DATA POINTS
C             (ZX,ZY,ZXX,ZXY,ZYY)
C
      COMMON /TRPCOP/ P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
     1,     Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
     2,     R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
     3,     P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C     P0...P5     COEFFICIENTS OF POLYNOMIALS ALONG THE SIDES
C     Q0...Q4     COEFFICIENTS OF 1ST DERIVATIVES ALONG THE SIDES
C     R0...R3     COEFFICIENTS OF 2ND DERIVATIVES ALONG THE SIDES
C     S0...S2     COEFFICIENTS OF 3RD DERIVATIVES ALONG THE SIDES
C     T0...T1     COEFFICIENTS OF 4TH DERIVATIVES ALONG THE SIDES
C     P11...P41   COEFFICIENTS FOR BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
      COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
     1,               SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
     2,               AP,BP,CP,DP
C     FOR EXPLANATION OF VARIABLES IN /TRPCOT/ SEE SUBROUTINE TRP00
C
      DIMENSION PD(5,3)
C     PARTIAL DERIVATIVES IN X-Y DIRECTIONS AT THE THREE VERTICES
C
C     LOADS PARTIAL DERIVATIVES AT THE VERTICES
      JIPT=3*(IT0-1)
      DO 230 K=1,3
        JIPT=JIPT+1
        IDP=IPT(JIPT)
        JPDD= 5*(IDP-1)
        DO 220 KPD=1,5
          JPDD=JPDD+1
          PD(KPD,K)=PDD(JPDD)
  220   CONTINUE
  230 CONTINUE
C     DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
C     TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
      AD=  DX(2)*DY(1)
      BC=  DX(1)*DY(2)
      DLT= AD-BC
      AP=  DY(1)/DLT
      BP= -DX(1)/DLT
      CP= -DY(2)/DLT
      DP=  DX(2)/DLT
C     CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE
C     TRIANGLE FOR THE U-V COORDINATE SYSTEM.
      AB=  DX(2)*DX(1)
      ADBC= AD+BC
      CD=  DY(1)*DY(2)
      DXDY1= 2.0d0*DX(1)*DY(1)
      DXDY2= 2.0d0*DX(2)*DY(2)
      DXDY3= 2.0d0*DX(3)*DY(3)
      DO 260 K=1,3
        ZT(1,K) = DX(2)*PD(1,K)  + DY(2)*PD(2,K)
        ZT(2,K) = DX(1)*PD(1,K)  + DY(1)*PD(2,K)
        ZTT(1,K)= DX2(2)*PD(3,K) + DXDY2*PD(4,K) + DY2(2)*PD(5,K)
        ZUV(K)  =   AB  *PD(3,K) +  ADBC*PD(4,K) +   CD  *PD(5,K)
        ZTT(2,K)= DX2(1)*PD(3,K) + DXDY1*PD(4,K) + DY2(1)*PD(5,K)
  260 CONTINUE
      DO 270 K=1,2
        ZT(3,K) = DX(3)*PD(1,K)  + DY(3)*PD(2,K)
        ZTT(3,K)= DX2(3)*PD(3,K) + DXDY3*PD(4,K) + DY2(3)*PD(5,K)
  270 CONTINUE
C     CALCULATES THE COEFFICIENTS OF THE POLYNOMIALS ALONG
C     THE THREE SIDES OF THE TRIANGLE
      DO 280 JSE=1,3
        NP1= 3
        NP2= 2
        IF (JSE.EQ.3) NP1= 1
        IF (JSE.EQ.2) NP2= 1
        IV= NP2
        IF (JSE.EQ.3) IV= 3
        P0(JSE)= Z(NP1)
        P1(JSE)= ZT(IV,NP1)
        P2(JSE)= 0.5d0*ZTT(IV,NP1)
        H1= Z(NP2) - P0(JSE) - P1(JSE) - P2(JSE)
        H2= ZT(IV,NP2) - P1(JSE) - ZTT(IV,NP1)
        H3= ZTT(IV,NP2) - ZTT(IV,NP1)
        P3(JSE)= 10.0d0*H1-4.0d0*H2+0.5d0*H3
        P4(JSE)=-15.0d0*H1+7.0d0*H2    -H3
        P5(JSE)=  6.0d0*H1-3.0d0*H2+0.5d0*H3
C       CALCULATES COEFFICIENTS FOR DERIVATIVES ALONG SIDES
        Q0(JSE)=     P1(JSE)
        Q1(JSE)= 2.d0 *P2(JSE)
        Q2(JSE)= 3.d0 *P3(JSE)
        Q3(JSE)= 4.d0 *P4(JSE)
        Q4(JSE)= 5.d0 *P5(JSE)
        R0(JSE)=     Q1(JSE)
        R1(JSE)= 2.d0 *Q2(JSE)
        R2(JSE)= 3.d0 *Q3(JSE)
        R3(JSE)= 4.d0 *Q4(JSE)
        S0(JSE)=     R1(JSE)
        S1(JSE)= 2.d0 *R2(JSE)
        S2(JSE)= 3.d0 *R3(JSE)
        T0(JSE)=     S1(JSE)
        T1(JSE)= 2.d0 *S2(JSE)
  280 CONTINUE
      RETURN
      END



      SUBROUTINE TRP002
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     COMPUTATION OF POLYNOMIAL COEFFICIENTS P11...P41 FOR
C     BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
C     COPYRIGHT(C): A. PREUSSER, 1984
C                   SEN. F. WISS. U. FORSCH.
C                   PROJEKTGRUPPE PARALLELRECHNEN
C                   HEILBRONNER STR. 10
C                   D-1000 BERLIN 31
C
      COMMON /TRPCOP/  P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
     1,     Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
     2,     R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
     3,     P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C     FOR EXPLANATION OF VARIABLES IN /TRPCOP/ SEE  SUBROUTINE TRP001
C
      COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
     1,               SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
     2,               AP,BP,CP,DP
C     FOR EXPLANATION OF VARIABLES IN /TRPCOT/ SEE  SUBROUTINE TRP00
C
C
      P14= P5(1) * (2.5d0*(SL2(2)-SL2(3))/SL2(1) + 2.5d0)
      P41= P5(2) * (2.5d0*(SL2(1)-SL2(3))/SL2(2) + 2.5d0)
      P11= ZUV(3)
      H1= ZT(2,1)-P1(1)-P11-P41
      H2= ZUV(1)-P11-4.0d0*P41
      P21= 3.0d0*H1-H2
      P31=-2.0d0*H1+H2
      H1= ZT(1,2)-P1(2)-P11-P14
      H2= ZUV(2)-P11-4.0d0*P14
      P12= 3.0d0*H1-H2
      P13=-2.0d0*H1+H2
      H1= 0.5d0*ZTT(2,1)-P2(1)-P12
      H2= 0.5d0*ZTT(1,2)-P2(2)-P21
      E1= 2.5d0*(SL2(2)-SL2(1))/SL2(3) + 2.5d0
      G1= 3.d0 - E1
      G2=-2.d0 + E1
      G3= E1*(P5(1) - P5(2) + P41 -P14)
     1    + P14 - 4.d0*P41 + 5.d0*P5(2)
      P22= G1*H1 + G2*H2 + G3
      P32= H1-P22
      P23= H2-P22
      RETURN
      END



      FUNCTION TRICPZ (TA,TB,F1,F2,ER)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     COMPUTE ZERO BETWEEN TA AND TB
C
C     F1= FUNCTION VALUE AT TA
C     F2= FUNCTION VALUE AT TB
C         F1 AND F2 MUST HAVE OPPOSITE SIGN
C         THIS MUST BE CHECKED BEFORE ENTRY
C     ER= PERMITTED ERROR FOR SOLUTION TRICPZ
C     NAME OF FUNCTION = TRICPF
C
C     THE METHOD IS A COMBINATION OF THE REGULA FALSI
C     AND THE MIDPOINT METHOD
C
C     IT IS A MODIFIED VERSION OF THE VIM- (CONTROL DATA
C     USER GROUP) ROUTINE WITH CATALOG IDENTIFICATION
C                C2BKYZERO
C     WRITTEN BY LOREN P. MEISSNER, 1965
C
      A=TA
      B=TB
      FA=F1
      FB=F2
      C=A
      FC=FA
      S=C
      FS=FC
C
   10 CONTINUE
      H=0.5d0*(B+C)
      IF(ABS(H-B) .LE.ER) GO TO 110
      IF (ABS(FB) .LE. ABS(FC)) GO TO 15
      Y=B
      FY=FB
      G=B
      FG=FB
      S=C
      FS=FC
      GO TO 20
   15 Y=S
      FY=FS
      G=C
      FG=FC
      S=B
      FS=FB
   20 CONTINUE
      IF (FY .NE. FS) GO TO 21
      B=H
      GO TO 29
   21 CONTINUE
      E=(S*FY-Y*FS)/(FY-FS)
      IF (ABS(E-S) .LE.ER) E=S+SIGN(ER,G-S)
      IF ((E-H)*(S-E) .LT. 0.0d0) GO TO 28
      B=E
      GO TO 29
   28 B=H
C
C *** FUNCTION CALL
   29 FB=TRICPF(B)
C
      IF (FG*FB .LT. 0.0d0) GO TO 35
      C=S
      FC=FS
      GO TO 10
   35 CONTINUE
      C=G
      FC=FG
      GO TO 10
C
  110 TRICPZ= H
C
      RETURN
      END



      FUNCTION TRICPF(T)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C
C     T R I ANGLE    C  ONTOUR    P  LOTTING
C
C     FUNCTION EVALUATION
C
C     COPYRIGHT(C): A. PREUSSER, 1984
C                   SEN. F. WISS. U. FORSCH.
C                   PROJEKTGRUPPE PARALLELRECHNEN
C                   HEILBRONNER STR. 10
C                   D-1000 BERLIN 31
C
      COMMON /TRPCOF/ KK,KSE,XX4F,YY4F,SIR,COR,CL
C
C     VARIABELS IN /TRPCOF/ ARE USED AS ARGUMENTS
C     FOR AN EXPLANATION SEE SUBROUTINE TRP00
C
C     KK      NUMBER OF FUNCTION TO BE EVALUATED
C     KK=1    4TH DERIVATIVE ALONG SIDE KSE
C     KK=2    3RD DERIVATIVE ALONG SIDE KSE
C     KK=3    2ND DERIVATIVE ALONG SIDE KSE
C     KK=4    1ST DERIVATIVE ALONG SIDE KSE
C     KK=5    ORIGINAL POLYNOMIAL ALONG SIDE KSE
C     KK=6    BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
      COMMON /TRPCOP/  P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
     1,     Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
     2,     R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
     3,     P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C     FOR EXPLANATION OF VARIABLES IN /TRPCOP/ SEE SUBROUTINE TRP001
C
      COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
     1,               SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
     2,               AP,BP,CP,DP
C     FOR AN EXPLANATION OF VARIABLES IN /TRPCOT/ SEE SUBROUTINE TRP00
C
      IF (KK.EQ.6) GOTO 60
      GOTO (10,20,30,40,50),KK
   10 TRICPF= T0(KSE) + T*T1(KSE)
      RETURN
   20 TRICPF= S0(KSE) + T*(S1(KSE)+T*S2(KSE))
      RETURN
   30 TRICPF= R0(KSE) + T*(R1(KSE)+T*(R2(KSE)+T*R3(KSE)))
      RETURN
   40 TRICPF= Q0(KSE)+T*(Q1(KSE)+T*(Q2(KSE)+T*(Q3(KSE)+T*Q4(KSE))))
      RETURN
   50 TRICPF= P0(KSE)+T*(P1(KSE)+T*(P2(KSE)+T*(P3(KSE)+T*(P4(KSE)+
     1        T*P5(KSE))))) - CL
      RETURN
   60 XX4= XX4F + COR*T
      YY4= YY4F + SIR*T
      U= AP*XX4 + BP*YY4
      V= CP*XX4 + DP*YY4
      H0=P0(1)+V*(P1(1)+V*(P2(1)+V*(P3(1)+V*(P4(1)+V*P5(1)))))
      H1=P1(2)+V*(P11+V*(P12+V*(P13+V*P14)))
      H2=P2(2)+V*(P21+V*(P22+V*P23))
      H3=P3(2)+V*(P31+V*P32)
      H4=P4(2)+V*P41
      TRICPF=H0+U*(H1+U*(H2+U*(H3+U*(H4+U*P5(2)))))  - CL
      RETURN
      END



      SUBROUTINE  IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C THIS SUBROUTINE PERFORMS TRIANGULATION.  IT DIVIDES THE X-Y
C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA
C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE
C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS
C CORRESPONDING TO THE BORDER LINE SEGMENTS.
C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE.  POINT NUMBERS OF THE END POINTS
C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE,
C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.
C THIS SUBROUTINE CALLS THE IDXCHG FUNCTION.
C THE INPUT PARAMETERS ARE
C     NDP = NUMBER OF DATA POINTS,
C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE
C           X COORDINATES OF THE DATA POINTS,
C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE
C           Y COORDINATES OF THE DATA POINTS.
C THE OUTPUT PARAMETERS ARE
C     NT  = NUMBER OF TRIANGLES,
C     IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE
C           POINT NUMBERS OF THE VERTEXES OF THE (IT)TH
C           TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND,
C           (3*IT-1)ST, AND (3*IT)TH ELEMENTS,
C           IT=1,2,...,NT,
C     NL  = NUMBER OF BORDER LINE SEGMENTS,
C     IPL = INTEGER ARRAY OF DIMENSION 6*NDP, WHERE THE
C           POINT NUMBERS OF THE END POINTS OF THE (IL)TH
C           BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE
C           NUMBER ARE TO BE STORED AS THE (3*IL-2)ND,
C           (3*IL-1)ST, AND (3*IL)TH ELEMENTS,
C           IL=1,2,..., NL.
C THE OTHER PARAMETERS ARE
C     IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
C           INTERNALLY AS A WORK AREA,
C     IWP = INTEGER ARRAY OF DIMENSION NDP USED
C           INTERNALLY AS A WORK AREA,
C     WK  = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
C           WORK AREA.
C DECLARATION STATEMENTS
      DIMENSION   XD(100),YD(100),IPT(585),IPL(600),
     1            IWL(1800),IWP(100),WK(100)
      DIMENSION   ITF(2)
      DATA  RATIO/1.0d-6/, NREP/100/
C STATEMENT FUNCTIONS
      DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
      SIDE(U1,V1,U2,V2,U3,V3)=(V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
C PRELIMINARY PROCESSING
   10 NDP0=NDP
      NDPM1=NDP0-1
      IF(NDP0.LT.4)       GO TO 90
C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
   20 DSQMN=DSQF(XD(1),YD(1),XD(2),YD(2))
      IPMN1=1
      IPMN2=2
      DO 22  IP1=1,NDPM1
        X1=XD(IP1)
        Y1=YD(IP1)
        IP1P1=IP1+1
        DO 21  IP2=IP1P1,NDP0
          DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
          IF(DSQI.EQ.0.0d0)      GO TO 91
          IF(DSQI.GE.DSQMN)    GO TO 21
          DSQMN=DSQI
          IPMN1=IP1
          IPMN2=IP2
   21   CONTINUE
   22 CONTINUE
      DSQ12=DSQMN
      XDMP=(XD(IPMN1)+XD(IPMN2))/2.0d0
      YDMP=(YD(IPMN1)+YD(IPMN2))/2.0d0
C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
C NUMBERS IN THE IWP ARRAY.
   30 JP1=2
      DO 31  IP1=1,NDP0
        IF(IP1.EQ.IPMN1.OR.IP1.EQ.IPMN2)      GO TO 31
        JP1=JP1+1
        IWP(JP1)=IP1
        WK(JP1)=DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
   31 CONTINUE
      DO 33  JP1=3,NDPM1
        DSQMN=WK(JP1)
        JPMN=JP1
        DO 32  JP2=JP1,NDP0
          IF(WK(JP2).GE.DSQMN)      GO TO 32
          DSQMN=WK(JP2)
          JPMN=JP2
   32   CONTINUE
        ITS=IWP(JP1)
        IWP(JP1)=IWP(JPMN)
        IWP(JPMN)=ITS
        WK(JPMN)=WK(JP1)
   33 CONTINUE
C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
   35 AR=DSQ12*RATIO
      X1=XD(IPMN1)
      Y1=YD(IPMN1)
      DX21=XD(IPMN2)-X1
      DY21=YD(IPMN2)-Y1
      DO 36  JP=3,NDP0
        IP=IWP(JP)
        IF(ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21).GT.AR)
     1               GO TO 37
   36 CONTINUE
      GO TO 92
   37 IF(JP.EQ.3)    GO TO 40
      JPMX=JP
      JP=JPMX+1
      DO 38  JPC=4,JPMX
        JP=JP-1
        IWP(JP)=IWP(JP-1)
   38 CONTINUE
      IWP(3)=IP
C FORMS THE FIRST TRIANGLE.  STORES POINT NUMBERS OF THE VER-
C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
C THE IPL ARRAY.
   40 IP1=IPMN1
      IP2=IPMN2
      IP3=IWP(3)
      IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
     1     .GE.0.0d0)       GO TO 41
      IP1=IPMN2
      IP2=IPMN1
   41 NT0=1
      NTT3=3
      IPT(1)=IP1
      IPT(2)=IP2
      IPT(3)=IP3
      NL0=3
      NLT3=9
      IPL(1)=IP1
      IPL(2)=IP2
      IPL(3)=1
      IPL(4)=IP2
      IPL(5)=IP3
      IPL(6)=1
      IPL(7)=IP3
      IPL(8)=IP1
      IPL(9)=1
C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
   50 DO 79  JP1=4,NDP0
        IP1=IWP(JP1)
        X1=XD(IP1)
        Y1=YD(IP1)
C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
        IP2=IPL(1)
        JPMN=1
        DXMN=XD(IP2)-X1
        DYMN=YD(IP2)-Y1
        DSQMN=DXMN**2+DYMN**2
        ARMN=DSQMN*RATIO
        JPMX=1
        DXMX=DXMN
        DYMX=DYMN
        DSQMX=DSQMN
        ARMX=ARMN
        DO 52  JP2=2,NL0
          IP2=IPL(3*JP2-2)
          DX=XD(IP2)-X1
          DY=YD(IP2)-Y1
          AR=DY*DXMN-DX*DYMN
          IF(AR.GT.ARMN)       GO TO 51
          DSQI=DX**2+DY**2
          IF(AR.GE.(-ARMN).AND.DSQI.GE.DSQMN)      GO TO 51
          JPMN=JP2
          DXMN=DX
          DYMN=DY
          DSQMN=DSQI
          ARMN=DSQMN*RATIO
   51     AR=DY*DXMX-DX*DYMX
          IF(AR.LT.(-ARMX))    GO TO 52
          DSQI=DX**2+DY**2
          IF(AR.LE.ARMX.AND.DSQI.GE.DSQMX)    GO TO 52
          JPMX=JP2
          DXMX=DX
          DYMX=DY
          DSQMX=DSQI
          ARMX=DSQMX*RATIO
   52   CONTINUE
        IF(JPMX.LT.JPMN)  JPMX=JPMX+NL0
        NSH=JPMN-1
        IF(NSH.LE.0)      GO TO 60
C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
        NSHT3=NSH*3
        DO 53  JP2T3=3,NSHT3,3
          JP3T3=JP2T3+NLT3
          IPL(JP3T3-2)=IPL(JP2T3-2)
          IPL(JP3T3-1)=IPL(JP2T3-1)
          IPL(JP3T3)  =IPL(JP2T3)
   53   CONTINUE
        DO 54  JP2T3=3,NLT3,3
          JP3T3=JP2T3+NSHT3
          IPL(JP2T3-2)=IPL(JP3T3-2)
          IPL(JP2T3-1)=IPL(JP3T3-1)
          IPL(JP2T3)  =IPL(JP3T3)
   54   CONTINUE
        JPMX=JPMX-NSH
C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
   60   JWL=0
        DO 64  JP2=JPMX,NL0
          JP2T3=JP2*3
          IPL1=IPL(JP2T3-2)
          IPL2=IPL(JP2T3-1)
          IT  =IPL(JP2T3)
C - - ADDS A TRIANGLE TO THE IPT ARRAY.
          NT0=NT0+1
          NTT3=NTT3+3
          IPT(NTT3-2)=IPL2
          IPT(NTT3-1)=IPL1
          IPT(NTT3)  =IP1
C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
          IF(JP2.NE.JPMX)      GO TO 61
          IPL(JP2T3-1)=IP1
          IPL(JP2T3)  =NT0
   61     IF(JP2.NE.NL0)       GO TO 62
          NLN=JPMX+1
          NLNT3=NLN*3
          IPL(NLNT3-2)=IP1
          IPL(NLNT3-1)=IPL(1)
          IPL(NLNT3)  =NT0
C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
C - - LINE SEGMENTS.
   62     ITT3=IT*3
          IPTI=IPT(ITT3-2)
          IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2)   GO TO 63
          IPTI=IPT(ITT3-1)
          IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2)   GO TO 63
          IPTI=IPT(ITT3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
   63     IF(IDXCHG(XD,YD,IP1,IPTI,IPL1,IPL2).EQ.0)     GO TO 64
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
          IPT(ITT3-2)=IPTI
          IPT(ITT3-1)=IPL1
          IPT(ITT3)  =IP1
          IPT(NTT3-1)=IPTI
          IF(JP2.EQ.JPMX)      IPL(JP2T3)=IT
          IF(JP2.EQ.NL0.AND.IPL(3).EQ.IT)     IPL(3)=NT0
C - - SETS FLAGS IN THE IWL ARRAY.
          JWL=JWL+4
          IWL(JWL-3)=IPL1
          IWL(JWL-2)=IPTI
          IWL(JWL-1)=IPTI
          IWL(JWL)  =IPL2
   64   CONTINUE
        NL0=NLN
        NLT3=NLNT3
        NLF=JWL/2
        IF(NLF.EQ.0)      GO TO 79
C - IMPROVES TRIANGULATION.
   70   NTT3P3=NTT3+3
        DO 78  IREP=1,NREP
          DO 76  ILF=1,NLF
            ILFT2=ILF*2
            IPL1=IWL(ILFT2-1)
            IPL2=IWL(ILFT2)
C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
C - - THE FLAGGED LINE SEGMENT.
            NTF=0
            DO 71  ITT3R=3,NTT3,3
              ITT3=NTT3P3-ITT3R
              IPT1=IPT(ITT3-2)
              IPT2=IPT(ITT3-1)
              IPT3=IPT(ITT3)
              IF(IPL1.NE.IPT1.AND.IPL1.NE.IPT2.AND.
     1           IPL1.NE.IPT3)      GO TO 71
              IF(IPL2.NE.IPT1.AND.IPL2.NE.IPT2.AND.
     1           IPL2.NE.IPT3)      GO TO 71
              NTF=NTF+1
              ITF(NTF)=ITT3/3
              IF(NTF.EQ.2)     GO TO 72
   71       CONTINUE
            IF(NTF.LT.2)       GO TO 76
C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
C - - ON THE LINE SEGMENT.
   72       IT1T3=ITF(1)*3
            IPTI1=IPT(IT1T3-2)
            IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2)    GO TO 73
            IPTI1=IPT(IT1T3-1)
            IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2)    GO TO 73
            IPTI1=IPT(IT1T3)
   73       IT2T3=ITF(2)*3
            IPTI2=IPT(IT2T3-2)
            IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2)    GO TO 74
            IPTI2=IPT(IT2T3-1)
            IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2)    GO TO 74
            IPTI2=IPT(IT2T3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
   74       IF(IDXCHG(XD,YD,IPTI1,IPTI2,IPL1,IPL2).EQ.0)
     1         GO TO 76
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
            IPT(IT1T3-2)=IPTI1
            IPT(IT1T3-1)=IPTI2
            IPT(IT1T3)  =IPL1
            IPT(IT2T3-2)=IPTI2
            IPT(IT2T3-1)=IPTI1
            IPT(IT2T3)  =IPL2
C - - SETS NEW FLAGS.
            JWL=JWL+8
            IWL(JWL-7)=IPL1
            IWL(JWL-6)=IPTI1
            IWL(JWL-5)=IPTI1
            IWL(JWL-4)=IPL2
            IWL(JWL-3)=IPL2
            IWL(JWL-2)=IPTI2
            IWL(JWL-1)=IPTI2
            IWL(JWL)  =IPL1
            DO 75  JLT3=3,NLT3,3
              IPLJ1=IPL(JLT3-2)
              IPLJ2=IPL(JLT3-1)
              IF((IPLJ1.EQ.IPL1.AND.IPLJ2.EQ.IPTI2).OR.
     1           (IPLJ2.EQ.IPL1.AND.IPLJ1.EQ.IPTI2))
     2                         IPL(JLT3)=ITF(1)
              IF((IPLJ1.EQ.IPL2.AND.IPLJ2.EQ.IPTI1).OR.
     1           (IPLJ2.EQ.IPL2.AND.IPLJ1.EQ.IPTI1))
     2                         IPL(JLT3)=ITF(2)
   75       CONTINUE
   76     CONTINUE
          NLFC=NLF
          NLF=JWL/2
          IF(NLF.EQ.NLFC)      GO TO 79
C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
          JWL=0
          JWL1MN=(NLFC+1)*2
          NLFT2=NLF*2
          DO 77  JWL1=JWL1MN,NLFT2,2
            JWL=JWL+2
            IWL(JWL-1)=IWL(JWL1-1)
            IWL(JWL)  =IWL(JWL1)
   77     CONTINUE
          NLF=JWL/2
   78   CONTINUE
   79 CONTINUE
C REARRANGES THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE.
   80 DO 81  ITT3=3,NTT3,3
        IP1=IPT(ITT3-2)
        IP2=IPT(ITT3-1)
        IP3=IPT(ITT3)
        IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
     1       .GE.0.0d0)     GO TO 81
        IPT(ITT3-2)=IP2
        IPT(ITT3-1)=IP1
   81 CONTINUE
      NT=NT0
      NL=NL0
      RETURN
C ERROR EXIT
   90 CALL ERROR_CALLBACK ("Less than 4 data points.")
      GO TO 93
   91 CALL ERROR_CALLBACK
     1("Two data points have identical X and Y values.")
      GO TO 93
   92 CALL ERROR_CALLBACK ("Data points are all colinear.")
   93 NT=0
      RETURN
      END



      FUNCTION  IDXCHG(X,Y,I1,I2,I3,I4)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO
C TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION
C BY C. L. LAWSON.
C THE INPUT PARAMETERS ARE
C     X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA
C           POINTS,
C     I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1, P2,
C           P3, AND P4 THAT FORM A QUADRILATERAL WITH P3
C           AND P4 CONNECTED DIAGONALLY.
C THIS FUNCTION RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EX-
C CHANGE IS NECESSARY, AND 0 (ZERO) OTHERWISE.
C DECLARATION STATEMENTS
      DIMENSION   X(100),Y(100)
      EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ),
     1            (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ)
C PRELIMINARY PROCESSING
   10 X1=X(I1)
      Y1=Y(I1)
      X2=X(I2)
      Y2=Y(I2)
      X3=X(I3)
      Y3=Y(I3)
      X4=X(I4)
      Y4=Y(I4)
C CALCULATION
   20 IDX=0
      U3=(Y2-Y3)*(X1-X3)-(X2-X3)*(Y1-Y3)
      U4=(Y1-Y4)*(X2-X4)-(X1-X4)*(Y2-Y4)
      IF(U3*U4.LE.0.0d0)    GO TO 30
      U1=(Y3-Y1)*(X4-X1)-(X3-X1)*(Y4-Y1)
      U2=(Y4-Y2)*(X3-X2)-(X4-X2)*(Y3-Y2)
      A1SQ=(X1-X3)**2+(Y1-Y3)**2
      B1SQ=(X4-X1)**2+(Y4-Y1)**2
      C1SQ=(X3-X4)**2+(Y3-Y4)**2
      A2SQ=(X2-X4)**2+(Y2-Y4)**2
      B2SQ=(X3-X2)**2+(Y3-Y2)**2
      C3SQ=(X2-X1)**2+(Y2-Y1)**2
      S1SQ=U1*U1/(C1SQ*MAX(A1SQ,B1SQ))
      S2SQ=U2*U2/(C2SQ*MAX(A2SQ,B2SQ))
      S3SQ=U3*U3/(C3SQ*MAX(A3SQ,B3SQ))
      S4SQ=U4*U4/(C4SQ*MAX(A4SQ,B4SQ))
      IF(MIN(S1SQ,S2SQ).LT.MIN(S3SQ,S4SQ))     IDX=1
   30 IDXCHG=IDX
      RETURN
      END



      SUBROUTINE  IDCLDP(NDP,XD,YD,NCP,IPC)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST
C TO EACH OF THE DATA POINT.
C THE INPUT PARAMETERS ARE
C     NDP = NUMBER OF DATA POINTS,
C     XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C           COORDINATES OF THE DATA POINTS,
C     NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA
C           POINTS.
C THE OUTPUT PARAMETER IS
C     IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE
C           POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C           EACH OF THE NDP DATA POINTS ARE TO BE STORED.
C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST
C NOT EXCEED 25.
C DECLARATION STATEMENTS
      DIMENSION   XD(100),YD(100),IPC(400)
      DIMENSION   DSQ0(25),IPC0(25)
      DATA  NCPMX/25/
C STATEMENT FUNCTION
      DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
C PRELIMINARY PROCESSING
   10 NDP0=NDP
      NCP0=NCP
      IF(NDP0.LT.2)  GO TO 90
      IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0)    GO TO 90
C CALCULATION
   20 DO 59  IP1=1,NDP0
C - SELECTS NCP POINTS.
        X1=XD(IP1)
        Y1=YD(IP1)
        J1=0
        DSQMX=0.0d0
        DO 22  IP2=1,NDP0
          IF(IP2.EQ.IP1)  GO TO 22
          DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
          J1=J1+1
          DSQ0(J1)=DSQI
          IPC0(J1)=IP2
          IF(DSQI.LE.DSQMX)    GO TO 21
          DSQMX=DSQI
          JMX=J1
   21     IF(J1.GE.NCP0)  GO TO 23
   22   CONTINUE
   23   IP2MN=IP2+1
        IF(IP2MN.GT.NDP0)      GO TO 30
        DO 25  IP2=IP2MN,NDP0
          IF(IP2.EQ.IP1)  GO TO 25
          DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
          IF(DSQI.GE.DSQMX)    GO TO 25
          DSQ0(JMX)=DSQI
          IPC0(JMX)=IP2
          DSQMX=0.0d0
          DO 24  J1=1,NCP0
            IF(DSQ0(J1).LE.DSQMX)   GO TO 24
            DSQMX=DSQ0(J1)
            JMX=J1
   24     CONTINUE
   25   CONTINUE
C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR.
   30   IP2=IPC0(1)
        DX12=XD(IP2)-X1
        DY12=YD(IP2)-Y1
        DO 31  J3=2,NCP0
          IP3=IPC0(J3)
          DX13=XD(IP3)-X1
          DY13=YD(IP3)-Y1
          IF((DY13*DX12-DX13*DY12)**2/(DSQ0(1)*DSQ0(J3))
     A       .GT.0.06698d0)       GOTO 50
C ----- 0.06698 CORRESPONDS TO AN ANGLE OF ABOUT 15. DEGREES(=SIN**2)
   31   CONTINUE
C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT.
   40   NCLPT=0
        DO 43  IP3=1,NDP0
          IF(IP3.EQ.IP1)       GO TO 43
          DO 41  J4=1,NCP0
            IF(IP3.EQ.IPC0(J4))     GO TO 43
   41     CONTINUE
          DX13=XD(IP3)-X1
          DY13=YD(IP3)-Y1
          DSQI=DSQF(X1,Y1,XD(IP3),YD(IP3))
          IF((DY13*DX12-DX13*DY12)**2/(DSQ0(1)*DSQI)
     A        .LE.0.06698d0)     GO TO 43
          IF(NCLPT.EQ.0)       GO TO 42
          IF(DSQI.GE.DSQMN)    GO TO 43
   42     NCLPT=1
          DSQMN=DSQI
          IP3MN=IP3
   43   CONTINUE
        IF(NCLPT.EQ.0)    GO TO 91
        DSQMX=DSQMN
        IPC0(JMX)=IP3MN
C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY.
   50   J1=(IP1-1)*NCP0
        DO 51  J2=1,NCP0
          J1=J1+1
          IPC(J1)=IPC0(J2)
   51   CONTINUE
   59 CONTINUE
      RETURN
C ERROR EXIT
   90 CALL ERROR_CALLBACK ("Improper input parameter value(s).")
      GO TO 92
   91 CALL ERROR_CALLBACK ("Data points are all colinear.")
   92 IPC(1)=0
      RETURN
      END



      SUBROUTINE  IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z), INTEGER(I-N)
C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND
C SECOND ORDER AT THE DATA POINTS.
C THE INPUT PARAMETERS ARE
C     NDP = NUMBER OF DATA POINTS,
C     XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
C           Y, AND Z COORDINATES OF THE DATA POINTS,
C     NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C           MATING PARTIAL DERIVATIVES AT EACH DATA POINT,
C     IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING
C           THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C           EACH OF THE NDP DATA POINTS.
C THE OUTPUT PARAMETER IS
C     PD  = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED
C           ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA
C           POINTS ARE TO BE STORED.
C DECLARATION STATEMENTS
      DIMENSION   XD(100),YD(100),ZD(100),IPC(400),PD(500)
      REAL        NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY
C PRELIMINARY PROCESSING
   10 NDP0=NDP
      NCP0=NCP
      NCPM1=NCP0-1
C ESTIMATION OF ZX AND ZY
  20  DO 24  IP0=1,NDP0
        X0=XD(IP0)
        Y0=YD(IP0)
        Z0=ZD(IP0)
        NMX=0.0d0
        NMY=0.0d0
        NMZ=0.0d0
        JIPC0=NCP0*(IP0-1)
        DO 23  IC1=1,NCPM1
          JIPC=JIPC0+IC1
          IPI=IPC(JIPC)
          DX1=XD(IPI)-X0
          DY1=YD(IPI)-Y0
          DZ1=ZD(IPI)-Z0
          IC2MN=IC1+1
          DO 22  IC2=IC2MN,NCP0
            JIPC=JIPC0+IC2
            IPI=IPC(JIPC)
            DX2=XD(IPI)-X0
            DY2=YD(IPI)-Y0
            DNMZ=DX1*DY2-DY1*DX2
            IF(DNMZ.EQ.0.0d0)    GO TO 22
            DZ2=ZD(IPI)-Z0
            DNMX=DY1*DZ2-DZ1*DY2
            DNMY=DZ1*DX2-DX1*DZ2
            IF(DNMZ.GE.0.0d0)    GO TO 21
            DNMX=-DNMX
            DNMY=-DNMY
            DNMZ=-DNMZ
   21       NMX=NMX+DNMX
            NMY=NMY+DNMY
            NMZ=NMZ+DNMZ
   22     CONTINUE
   23   CONTINUE
        JPD0=5*IP0
        PD(JPD0-4)=-NMX/NMZ
        PD(JPD0-3)=-NMY/NMZ
   24 CONTINUE
C ESTIMATION OF ZXX, ZXY, AND ZYY
   30 DO 34  IP0=1,NDP0
        JPD0=JPD0+5
        X0=XD(IP0)
        JPD0=5*IP0
        Y0=YD(IP0)
        ZX0=PD(JPD0-4)
        ZY0=PD(JPD0-3)
        NMXX=0.0d0
        NMXY=0.0d0
        NMYX=0.0d0
        NMYY=0.0d0
        NMZ =0.0d0
        JIPC0=NCP0*(IP0-1)
        DO 33  IC1=1,NCPM1
          JIPC=JIPC0+IC1
          IPI=IPC(JIPC)
          DX1=XD(IPI)-X0
          DY1=YD(IPI)-Y0
          JPD=5*IPI
          DZX1=PD(JPD-4)-ZX0
          DZY1=PD(JPD-3)-ZY0
          IC2MN=IC1+1
          DO 32  IC2=IC2MN,NCP0
            JIPC=JIPC0+IC2
            IPI=IPC(JIPC)
            DX2=XD(IPI)-X0
            DY2=YD(IPI)-Y0
            DNMZ =DX1*DY2 -DY1*DX2
            IF(DNMZ.EQ.0.0d0)    GO TO 32
            JPD=5*IPI
            DZX2=PD(JPD-4)-ZX0
            DZY2=PD(JPD-3)-ZY0
            DNMXX=DY1*DZX2-DZX1*DY2
            DNMXY=DZX1*DX2-DX1*DZX2
            DNMYX=DY1*DZY2-DZY1*DY2
            DNMYY=DZY1*DX2-DX1*DZY2
            IF(DNMZ.GE.0.0d0)    GO TO 31
            DNMXX=-DNMXX
            DNMXY=-DNMXY
            DNMYX=-DNMYX
            DNMYY=-DNMYY
            DNMZ =-DNMZ
   31       NMXX=NMXX+DNMXX
            NMXY=NMXY+DNMXY
            NMYX=NMYX+DNMYX
            NMYY=NMYY+DNMYY
            NMZ =NMZ +DNMZ
   32     CONTINUE
   33   CONTINUE
        PD(JPD0-2)=-NMXX/NMZ
        PD(JPD0-1)=-(NMXY+NMYX)/(2.0d0*NMZ)
        PD(JPD0)  =-NMYY/NMZ
   34 CONTINUE
      RETURN
      END
