C      ALGORITHM 677, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 15, NO. 4, PP. 365-374.
      PROGRAM MAIN
C  IT IS A TEST PROGRAM FOR THE MASUB SUBPROGRAM PACKAGE.
C  THE INPUT CONSISTS OF THE SCATTERED DATA POINTS
C            (XD(I),YD(I),ZD(I)),I=1,..,30.
C  X-Y COORDINATES ARE CONTAINED IN DATA STATEMENTS;
C  Z COORDINATES ARE CALCULATED AS VALUES OF THE GAUSSIAN FUNCTION
C                    F(X,Y) = 8./3.14*EXP(-8.*(X**2+Y**2))
C  AT THESE POINTS.
C  RECONSTRUCTION IS OBTAINED AT 10X10 RECTANGULAR GRID POINTS
C  USING EXTRAPOLATION (IEX=1) AND THE  TENSION PARAMETER TP=10.
C  THE OUTPUT OF THE PROGRAM CONSISTS OF THE DOUBLY DIMENSIONED
C  ARRAYS OF THE ANALYTIC AND RECONSTRUCTED VALUES AT THE GRID POINTS.
C  THE NOUT CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS
C  THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C  THEREFORE, SYSTEM DEPENDENT.
C
C  DECLARATION STATEMENTS.
      DIMENSION XD(40),YD(40),ZD(40),
     *          XI(10),YI(10),ZI(10,10),ZI1(10,10),
     *          IWK(1111),WK(664)
      DATA ANS1/0.06695094/,ANS2/0.6059994/,TOL/1E-6/
      DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6),
     *     XD(7), XD(8), XD(9), XD(10),XD(11),XD(12),
     *     XD(13),XD(14),XD(15),XD(16),XD(17),XD(18),
     *     XD(19),XD(20),XD(21),XD(22),XD(23),XD(24),
     *     XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/
     * -1.0,     1.0,    -0.93572,-0.24358,   -0.59047,-1.0,
     * -0.33608, 0.69231,-0.31579, 0.23818,   -0.13344, 0.30510,
     *  0.30086, 0.70214, 0.26879,-3.5654E-02,-0.41440,-0.57115,
     * -0.71437,-0.59515,-0.84836, 1.0,       -0.41212, 0.66413,
     * -0.66319, 0.26784, 0.53973, 1.9985E-02, 0.64104,-7.993E-02/
      DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6),
     *     YD(7), YD(8), YD(9), YD(10),YD(11),YD(12),
     *     YD(13),YD(14),YD(15),YD(16),YD(17),YD(18),
     *     YD(19),YD(20),YD(21),YD(22),YD(23),YD(24),
     *     YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/
     * -1.0,    -1.0,       -0.77752,-0.46379, 0.63127, 1.0,
     * -0.38873,-0.34749,   -0.76732, 0.33503, 0.18401, 0.17446,
     *  0.23498, 9.8021E-02, 0.73058,-0.78921, 0.61654,-0.97579,
     *  0.49590,-3.4130E-02,-0.78300, 1.0,    -0.50750, 0.55233,
     * -0.95014, 0.15838,   -0.18706,-0.19651,-0.38508, 0.98614/
      DATA ND/30/,NXI/10/,NYI/10/,A/-1./,B/1./,C/-1./,D/1./
      DATA NOUT/6/,NDIM/10/
C  STATEMENT FUNCTION.
      FUN(X,Y)=8./3.14*EXP(-8.*(X**2+Y**2))
C  CALCULATION.
      CALL MACEPS(SRELPR)
      DO 10 I=1,ND
        ZD(I)=FUN(XD(I),YD(I))
10    CONTINUE
      HX=(B-A)/FLOAT(NXI-1)
      DO 20 I=1,NXI
        XI(I)=A+FLOAT(I-1)*HX
20    CONTINUE
      HY=(D-C)/FLOAT(NYI-1)
      DO 30 I=1,NYI
        YI(I)=C+FLOAT(I-1)*HY
30    CONTINUE
      DO 40 J=1,NYI
      DO 40 I=1,NXI
        ZI(I,J)=FUN(XI(I),YI(J))
40    CONTINUE
      IEX=1
      IC=1
      TP=10.
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *           NOUT,IWK,WK)
C  PRINTING OF INPUT DATA.
      WRITE(NOUT,1000)ND
      DO 50 I=1,ND
        WRITE(NOUT,2000)I,XD(I),YD(I),ZD(I)
50    CONTINUE
      WRITE(NOUT,3000)
      WRITE(NOUT,4000)XI
      DO 60 J=1,NYI
        WRITE(NOUT,5000)YI(J),(ZI(I,J),I=1,NXI)
60    CONTINUE
C  PRINTING OF OUTPUT RESULTS.
      WRITE(NOUT,6000)
      WRITE(NOUT,4000)XI
        DO 70 J=1,NYI
      WRITE(NOUT,5000)YI(J),(ZI1(I,J),I=1,NXI)
70    CONTINUE
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      IF ((ABS(ERMEAN-ANS1)/ANS1.GT.TOL).OR.
     *    (ABS(ERMAX-ANS2)/ANS2.GT.TOL)) THEN
        WRITE(NOUT,7000)
                                         ELSE
        WRITE(NOUT,8000)
      ENDIF
      WRITE(NOUT,9000)ERMAX,ANS2,ERMEAN,ANS1
      STOP
1000  FORMAT(1H1,3X,8HTEST RUN,///3X,10HINPUT DATA,8X,4HND =,I3//,
     *       30H      I      XD     YD     ZD /)
2000  FORMAT(5X,I2,2X,2F7.3,E14.7)
3000  FORMAT(//3X,8HANALYTIC)
4000  FORMAT(//7X,2HYI,3X,3HXI=/12X,10F7.3/)
5000  FORMAT(/1X,F9.3,2X,10F7.3)
6000  FORMAT(//3X,14HRECONSTRUCTION)
7000  FORMAT(//3X,29HTHE CODE HAS NOT SUCCESSFULLY,
     *        /3X,25HEXECUTED ON YOUR COMPUTER)
8000  FORMAT(//3X,25HTHE CODE HAS SUCCESSFULLY,
     *        /3X,25HEXECUTED ON YOUR COMPUTER)
9000  FORMAT(//3X,29HMAX ERROR ON YOUR COMPUTER = ,E14.7,
     *        /3X,21HMAX ERROR BUILT IN = ,E14.7,
     *       //3X,30HMEAN ERROR ON YOUR COMPUTER = ,E14.7,
     *        /3X,22HMEAN ERROR BUILT IN = ,E14.7)
      END
      SUBROUTINE MACEPS(SRELPR)
C IT CARRIES OUT AN APPROXIMATION OF THE SINGLE RELATIVE PRECISION
C CONSTANT
      SRELPR=1.0
10    SRELPR=0.5*SRELPR
C THE DO LOOP STATEMENT IS NECESSARY WHEN THE ARITHMETIC UNIT HAS MORE
C BITS THAN IN STORAGE (INTEL 8087 FAMILY OF ARITHMETIC UNITS IS OF
C THIS TYPE), INFACT A DO LOOP INVOLVES A STORE FROM REGISTER TO MEMORY
C OF THE VALUE FOR SRELPR.
      DO 20 I=1,2
20    CONTINUE
      IF (SRELPR+1.0.GT.1.0) GOTO 10
      SRELPR=2.0*SRELPR
      RETURN
      END
      SUBROUTINE MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
C IT CARRIES OUT MAX ERROR AND MEAN ERROR AT PRESCRIBED RECTANGULAR
C GRID POINTS OF THE RECONSTRUCTED SURFACE.
C
C THE INPUT PARAMETERS ARE
C     NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE (MUST BE
C           1 OR GREATER),
C     NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE (MUST BE
C           1 OR GREATER),
C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           ANALYTIC  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     ZI1 = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           INTERPOLATED  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     NDIM= DECLARED ROW DIMENSION OF THE ARRAYS ZI AND ZI1.
C
C THE OUTPUT PARAMETERS ARE
C     ERMAX = MAX ERROR
C     ERMEAN= MEAN ERROR
C
      DIMENSION ZI(NDIM,*),ZI1(NDIM,*)
      ERMAX=0.0
      SUM=0.0
      N=0
      DO 10 I=1,NXI
        DO 10 J=1,NYI
          IF (ZI1(I,J).EQ.1.E6) GOTO 10
          N=N+1
          ERR=ABS(ZI(I,J)-ZI1(I,J))
          SUM=SUM+ERR
          IF (ERR.GT.ERMAX) ERMAX=ERR
10    CONTINUE
      ERMEAN=SUM/FLOAT(N)
      RETURN
      END
      PROGRAM MAIN
C  IT IS A SOPHISTICATED TEST TO VERIFY SOME OF THE MASUB
C  SUBPROGRAM PACKAGE CAPABILITIES.
C  THE INPUT CONSISTS OF THE SCATTERED DATA POINTS
C            (XD(I),YD(I),ZD(I)),I=1,..,100.
C  X-Y COORDINATES ARE GENERATED BY A RANDOM NUMBER GENERATOR;
C  Z COORDINATES ARE CALCULATED AS VALUES OF THE ANALYITIC
C  EXPONENTIAL FUNCTION (F(X,Y)) AND OF THE ANALYTIC SADDLE
C  FUNCTION (G(X,Y)) AT THE GIVEN POINTS.
C  RECONSTRUCTIONS ARE OBTAINED BOTH WITH AND WITHOUT EXTRAPOLATION,
C  WITH DIFFERENT TENSION PARAMETERS AND AT DIFFERENT GRID POINTS.
C  THE OUTPUT OF THE PROGRAM CONSISTS OF A GREY LEVEL MAP OF THE
C  ANALYTIC FUNCTION AND OF SOME MAPS OF THE VARIOUS RECONSTRUCTIONS
C  FOR EACH TEST FUNCTION.
C  THE  NOUT  CONSTANT  IN  THE  DATA  INITIALIZATION  STATEMENT
C  IS THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C  THEREFORE , SYSTEM DEPENDENT.
C  PROGRAM MAIN2 CALLS THE RAND, GGRAF AND MASUB SUBROUTINES.
C
C  DECLARATION STATEMENTS
      DIMENSION XD(112),YD(112),ZD(112)
      DIMENSION XI(40),YI(40),ZI(40,40),ZI1(40,40)
      DIMENSION IWK(4145),WK(2108)
      DATA A,B,C,D/1.,10.,1.,10./
      DATA ND/100/,NOUT/6/,NDIM/40/
C  ANALYTIC EXPRESSION OF THE EXPONENTIAL FUNCTION
      F(X,Y)=9.*(3./4.*EXP((-(X-3.)**2-(Y-3.)**2)/4.)+EXP(-(X/7.)**2-(Y/
     .10.))-1./5.*EXP(-(X-5.)**2-(Y-8.)**2)+1./2.*EXP((-(X-8.)**2-(Y-4.)
     .**2)/4.))
C  ANALYTIC EXPRESSION OF THE SADDLE FUNCTION
      G(X,Y)=3./2.*(COS(3./5.*(Y-1.))+5./4.)/(1.+((X-4.)/3.)**2)
      CALL MACEPS(SRELPR)
C  RANDOM CHOICE OF THE ND DATA POINTS
      BMA=B-A
      DMC=D-C
      IX=1
      DO 10 I=1,ND
      XCASO=BMA*RAND(IX)
      XD(I)=XCASO+A
      YCASO=DMC*RAND(IX)
      YD(I)=YCASO+C
      ZD(I)=F(XD(I),YD(I))
  10  CONTINUE
C  CONSTRUCTION OF THE REGULAR 30*30 GRID POINTS
      NXI=30
      NYI=30
      HX=BMA/FLOAT(NXI-1)
      HY=DMC/FLOAT(NYI-1)
      DO 20 I=1,NXI
      XI(I)=A+FLOAT(I-1)*HX
  20  CONTINUE
      DO 30 J=1,NYI
      YI(J)=C+FLOAT(J-1)*HY
  30  CONTINUE
C  PRINT OF THE EXPONENTIAL ANALYTIC MAP
      DO 40 I=1,NXI
      DO 40 J=1,NYI
      ZI(I,J)=F(XI(I),YI(J))
  40  CONTINUE
      WRITE(NOUT,1000)
1000  FORMAT(/,15X,'EXPONENTIAL: ANALYTIC MAP',//)
      ZMIN=-0.5
      ZMAX=12.5
      CALL GGRAF(ZI,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
C  FIRST CALL WITH EXTRAPOLATION
      IEX=1
      IC=1
      TP=3.
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *NOUT,IWK,WK)
      WRITE(NOUT,2000)TP
2000  FORMAT (///,15X,'EXPONENTIAL: RECONSTRUCTED MAP WITH EXTRAPOLATION
     *',/,29X,'AND TENSION PARAMETER TP=',1F5.2,//)
      CALL GGRAF(ZI1,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      WRITE(NOUT,7000)ERMAX,ERMEAN
7000  FORMAT(//3X,12HMAX ERROR = ,E14.7,/3X,13HMEAN ERROR = ,E14.7)
C  FIRST CALL WITHOUT EXTRAPOLATION
      IEX=0
      IC=1
      TP=0.
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *NOUT,IWK,WK)
      WRITE(NOUT,3000)TP
3000  FORMAT (///,15X,'EXPONENTIAL: RECONSTRUCTED MAP WITHOUT EXTRAPOLAT
     *ION',/,29X,'AND TENSION PARAMETER TP=',1F5.2,//)
      CALL GGRAF(ZI1,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      WRITE(NOUT,7000)ERMAX,ERMEAN
C  CALL WITH NEW TENSION PARAMETER
      IC=3
      TP=3
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *NOUT,IWK,WK)
      WRITE(NOUT,3000)TP
      CALL GGRAF(ZI1,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      WRITE(NOUT,7000)ERMAX,ERMEAN
C  NEW REGULAR 40*40 GRID POINTS
      NXI=40
      NYI=40
      HX=BMA/FLOAT(NXI-1)
      HY=DMC/FLOAT(NYI-1)
      DO 60 I=1,NXI
      XI(I)=A+FLOAT(I-1)*HX
60    CONTINUE
      DO 70 J=1,NYI
      YI(J)=C+FLOAT(J-1)*HY
70    CONTINUE
      DO 80 I=1,NXI
      DO 80 J=1,NYI
      ZI(I,J)=F(XI(I),YI(J))
80    CONTINUE
      WRITE(NOUT,1000)
      CALL GGRAF(ZI,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
C  CALL WITH NEW GRID POINTS
      IC=4
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *NOUT,IWK,WK)
      WRITE(NOUT,3000)TP
      CALL GGRAF(ZI1,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      WRITE(NOUT,7000)ERMAX,ERMEAN
C  NEW ZD VALUES
      DO 90 I=1,ND
      ZD(I)=G(XD(I),YD(I))
90    CONTINUE
      DO 100 I=1,NXI
      DO 100 J=1,NYI
      ZI(I,J)=G(XI(I),YI(J))
100   CONTINUE
C  PRINT OF THE SADDLE ANALYTIC MAP
      WRITE(NOUT,5000)
5000  FORMAT(///,15X,'SADDLE: ANALYTIC MAP',//)
      ZMIN=0.
      ZMAX=3.5
      CALL GGRAF(ZI,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
C  CALL WITH NEW ZD VALUES
      IC=2
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *NOUT,IWK,WK)
      WRITE(NOUT,6000)TP
6000  FORMAT(///,15X,'SADDLE: RECONSTRUCTED MAP WITHOUT EXTRAPOLATION'
     *,/,24X,'WITH TENSION PARAMETER TP=',1F5.2)
      CALL GGRAF(ZI1,NDIM,NXI,NYI,ZMIN,ZMAX,NOUT)
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      WRITE(NOUT,7000)ERMAX,ERMEAN
      STOP
      END
      FUNCTION RAND(IX)
C  PORTABLE RANDOM NUMBER GENERATOR
C   USING THE RECURSION
C    IX = IX*A MOD P
C
C  SOME COMPILERS, E.G., THE HP3000, REQUIRE THE FOLLOWING
C  DECLARATION TO BE INTEGER *4
      INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K
C
C  7**5, 2**15, 2**16, 2**31-1
      DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/
C
C  GET 15 HI ORDER BITS OF IX
      XHI=IX/B16
C  GET 16 LO BITS OF IX AND FORM LO PRODUCT
      XALO=(IX-XHI*B16)*A
C  GET 15 HI ORDER BITS OF LO PRODUCT
      LEFTLO=XALO/B16
C  FORM THE 31 HIGHEST BITS OF FULL PRODUCT
      FHI=XHI*A+LEFTLO
C  GET OVERFLOW PAST 31ST BIT OF FULL PRODUCT
      K=FHI/B15
C  ASSEMBLE ALL THE PARTS AND PRESUBTRACT P
C  THE PARENTHESES ARE ESSENTIAL
      IX=(((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16)+K
C  ADD P BACK IN IF NECESSARY
      IF(IX.LT.0)IX=IX+P
C  MULTIPLY BY 1/(2**31-1)
      RAND=FLOAT(IX)*4.656612875E-10
      RETURN
      END
      SUBROUTINE GGRAF(A,NDIM,NX,NY,AMIN,AMAX,NOUT)
C
C  IT PRINTS A SIX GREY-LEVEL MAP OF THE  BIVARIATE
C  FUNCTION, WHOSE VALUES ARE STORED IN A DOUBLY DIMENSIONED ARRAY.
C
C  THE INPUT PARAMETERS ARE
C      A   = DOUBLY DIMENSIONED ARRAY OF DIMENSION (NX X NY)
C            CONTAINING THE VALUES OF THE FUNCTION.
C      NDIM= DECLARED ROW DIMENSION OF THE ARRAY A.
C      NX  = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE.
C      NY  = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE.
C      AMIN= MINIMUM FUNCTION'S VALUE.
C      AMAX= MAXIMUM FUNCTION'S VALUE.
C      NOUT= LOGICAL UNIT NUMBER FOR THE STANDARD OUTPUT UNIT OF THE
C            SYSTEM.
C
C  ATTENTION: IN THE CASE OF AMIN=AMAX, THE SUBROUTINE ESTIMATES
C  THE MINIMUM AND MAXIMUM VALUE OF THE FUNCTION AS THE MINIMUM
C  AND MAXIMUM VALUE BETWEEN THE ELEMENTS OF THE ARRAY A,
C  THE VALUES OF THE FUNCTION MUST BE STORED IN THE ARRAY A IN THE
C  FOLLOWING FORM:
C        A(I,J) = F(XI(I),YI(J))  I=1,NX, J=1,NY
C  WHERE (XI(I),YI(J)) ARE THE COORDINATES X AND Y OF THE OUTPUT
C  GRID POINTS, IN ASCENDING ORDER.
C  THE MAX DIMENSION OF THE MAP IS 50X50.
      DIMENSION A(NDIM,*)
      CHARACTER *2 GREY(11),LINE(52)
      DATA GREY/'..','//','++','OO','$$','**',
     *          '**','  ','..','  ','. '/
      IF(AMIN.NE.AMAX)GO TO 20
      AMIN=1.E7
      AMAX=-1.E7
      DO 10 J=1,NY
      DO 10 I=1,NX
      IF(A(I,J).EQ.1.E6)GO TO 10
      IF(A(I,J).GT.AMAX)AMAX=A(I,J)
      IF(A(I,J).LT.AMIN)AMIN=A(I,J)
   10 CONTINUE
   20 H=(AMAX-AMIN)/6.
      BMIN=AMIN+H/2.
      BMAX=AMAX-H/2.
      FRM=5./(BMAX-BMIN)
      WRITE(NOUT,1)
      DO 60 J=1,NY
      L=NY-J+1
      DO 50 I=1,NX
      NL=8
      IF(A(I,L).GE.AMIN.AND.A(I,L).LE.AMAX)NL=(A(I,L)-BMIN)*FRM+1.51
      LINE(I)=GREY(NL)
   50 CONTINUE
      WRITE(NOUT,2)(LINE(I),I=1,NX)
   60 CONTINUE
      K=NX
      DO 80 I=1,3
      L=I+7
      DO 70 J=1,K
      LINE(J)=GREY(L)
   70 CONTINUE
      K=K+1-2*(I/3)
      LINE(K)=GREY(11)
      WRITE(NOUT,3)(LINE(J),J=1,K)
   80 CONTINUE
      WRITE(NOUT,4)
      Q2=AMIN
      DO 90 I=1,6
      Q1=Q2
      Q2=Q1+H
      WRITE(NOUT,5)GREY(I),Q1,Q2
   90 CONTINUE
300   RETURN
    1 FORMAT(11X,1H./10X,3H...)
    2 FORMAT(11X,1H.,3X,50A2)
    3 FORMAT(15X,52A2)
    4 FORMAT(10X,6HMARKER,13X,6HHEIGHT,/)
    5 FORMAT(12X,A2,8X,1H(,F7.2,3H , ,F7.2,1H))
      END
      SUBROUTINE MACEPS(SRELPR)
C IT CARRIES OUT AN APPROXIMATION OF MACHINE PRECISION
      SRELPR=1.0
10    SRELPR=0.5*SRELPR
C THE DO LOOP STATEMENT IS NECESSARY WHEN THE ARITHMETIC UNIT HAS MORE
C BITS THAN IN STORAGE (INTEL 8087 FAMILY OF ARITHMETIC UNITS IS OF
C THIS TYPE), INFACT A DO LOOP INVOLVES A STORE FROM REGISTER TO MEMORY
C OF THE VALUE FOR SRELPR.
      DO 20 I=1,2
20    CONTINUE
      IF (SRELPR+1.0.GT.1.0) GOTO 10
      SRELPR=2.0*SRELPR
      RETURN
      END
      SUBROUTINE MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
C IT CARRIES OUT MAX ERROR AND MEAN ERROR AT PRESCRIBED RECTANGULAR
C GRID POINTS OF THE RECONSTRUCTED SURFACE.
C
C THE INPUT PARAMETERS ARE
C     NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE (MUST BE
C           1 OR GREATER),
C     NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE (MUST BE
C           1 OR GREATER),
C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           ANALYTIC  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     ZI1 = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           INTERPOLATED  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     NDIM= DECLARED ROW DIMENSION OF THE ARRAYS ZI AND ZI1.
C
C THE OUTPUT PARAMETERS ARE
C     ERMAX = MAX ERROR
C     ERMEAN= MEAN ERROR
C
      DIMENSION ZI(NDIM,*),ZI1(NDIM,*)
      ERMAX=0.0
      SUM=0.0
      N=0
      DO 10 I=1,NXI
        DO 10 J=1,NYI
          IF (ZI1(I,J).EQ.1.E6) GOTO 10
          N=N+1
          ERR=ABS(ZI(I,J)-ZI1(I,J))
          SUM=SUM+ERR
          IF (ERR.GT.ERMAX) ERMAX=ERR
10    CONTINUE
      ERMEAN=SUM/FLOAT(N)
      RETURN
      END
      PROGRAM MAIN
C  IT IS A TEST PROGRAM FOR THE MASUB SUBPROGRAM PACKAGE.
C  THE INPUT CONSISTS OF THE DATA POINTS OF LAWSON'S EXAMPLE (REFERENCE
C  5 OF THE MANUSCRIPT) (XD(I),YD(I),ZD(I)),I=1,..,26.
C  X-Y COORDINATES ARE CONTAINED IN DATA STATEMENTS;
C  Z COORDINATES ARE CALCULATED AS VALUES OF THE GAUSSIAN FUNCTION
C                    F(X,Y) = 8./3.14*EXP(-8.*(X**2+Y**2))
C  AT THESE POINTS.
C  RECONSTRUCTION IS OBTAINED AT 40X40 RECTANGULAR GRID POINTS
C  USING EXTRAPOLATION (IEX=1) AND THE  TENSION PARAMETER TP=10.
C  THE OUTPUT OF THE PROGRAM CONSISTS OF THE MAXIMUM AND MEAN ABSOLUTE
C  ERRORS ON THE DOUBLY DIMENSIONED ARRAYS OF THE ANALYTIC AND
C  RECONSTRUCTED VALUES AT THE GRID POINTS.
C  THE NOUT CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS
C  THE LOGICAL UNIT NUMBER FOR THE STANDARD OUTPUT UNIT AND IS,
C  THEREFORE, SYSTEM DEPENDENT.
C
C  DECLARATION STATEMENTS.
      DIMENSION XD(32),YD(32),ZD(32),
     *          XI(40),YI(40),ZI(40,40),ZI1(40,40),
     *          IWK(2400),WK(588)
      DATA ND/26/,NXI/40/,NYI/40/,NDIM/40/
      DATA NOUT/6/,ANS1/0.04544656/,ANS2/0.5812947/,TOL/1.0E-6/
      DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6),
     *     XD(7), XD(8), XD(9), XD(10),XD(11),XD(12),
     *     XD(13),XD(14),XD(15),XD(16),XD(17),XD(18),
     *     XD(19),XD(20),XD(21),XD(22),XD(23),XD(24),
     *     XD(25),XD(26)/
     * -0.95, -0.89, -0.87, -0.84, -0.82, -0.79,
     * -0.77, -0.70, -0.58, -0.52, -0.47, -0.41,
     * -0.22, -0.18, -0.08,  0.00,  0.03,  0.03,
     *  0.05,  0.08,  0.10,  0.15,  0.50,  0.70,
     *  0.90,  0.95/
      DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6),
     *     YD(7), YD(8), YD(9), YD(10),YD(11),YD(12),
     *     YD(13),YD(14),YD(15),YD(16),YD(17),YD(18),
     *     YD(19),YD(20),YD(21),YD(22),YD(23),YD(24),
     *     YD(25),YD(26)/
     * -0.42, -0.10, -0.65, -0.55,  0.12, -0.88,
     * -0.32,  0.32, -0.67, -0.13,  0.13,  0.42,
     *  0.12,  0.95, -0.68,  0.75, -0.79,  0.48,
     *  0.31,  0.10,  0.28,  0.41,  0.11,  0.28,
     * -0.60,  0.68/
C  STATEMENT FUNCTION.
      FUN(X,Y)=8./3.14*EXP(-8.*(X*X+Y*Y))
C  CALCULATION.
      CALL MACEPS(SRELPR)
      DO 10 I=1,ND
        ZD(I)=FUN(XD(I),YD(I))
10    CONTINUE
      A=XD(1)
      B=XD(1)
      C=YD(1)
      D=YD(1)
      DO 15 I=2,ND
        IF (XD(I).GT.B) THEN
          B=XD(I)
                        ELSE
          IF (XD(I).LT.A) A=XD(I)
        ENDIF
        IF (YD(I).GT.D) THEN
          D=YD(I)
                        ELSE
          IF (YD(I).LT.C) C=YD(I)
        ENDIF
15    CONTINUE
      HX=(B-A)/FLOAT(NXI-1)
      XI(1)=A
      DO 20 I=2,NXI-1
        XI(I)=XI(I-1)+HX
20    CONTINUE
      XI(NXI)=B
      HY=(D-C)/FLOAT(NYI-1)
      YI(1)=C
      DO 30 I=2,NYI-1
        YI(I)=YI(I-1)+HY
30    CONTINUE
      YI(NYI)=D
      DO 40 J=1,NYI
      DO 40 I=1,NXI
        ZI(I,J)=FUN(XI(I),YI(J))
40    CONTINUE
      IEX=1
      IC=1
      TP=10.
      CALL MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI1,NDIM,
     *           NOUT,IWK,WK)
C  PRINTING OF INPUT DATA.
      WRITE(NOUT,1000)ND
      DO 50 I=1,ND
        WRITE(NOUT,2000)I,XD(I),YD(I),ZD(I)
50    CONTINUE
C  PRINTING OF OUTPUT RESULTS.
      CALL MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
      IF ((ABS(ERMEAN-ANS1)/ANS1.GT.TOL).OR.
     *    (ABS(ERMAX-ANS2)/ANS2.GT.TOL)) THEN
        WRITE(NOUT,3000)
                                         ELSE
        WRITE(NOUT,4000)
      ENDIF
      WRITE(NOUT,5000)ERMAX,ANS2,ERMEAN,ANS1
      STOP
1000  FORMAT(1H1,3X,8HTEST RUN,///3X,10HINPUT DATA,8X,4HND =,I3//,
     *       30H      I      XD     YD     ZD /)
2000  FORMAT(5X,I2,2X,2F7.3,e14.7)
3000  FORMAT(//3X,29HTHE CODE HAS NOT SUCCESSFULLY,
     *        /3X,25HEXECUTED ON YOUR COMPUTER)
4000  FORMAT(//3X,25HTHE CODE HAS SUCCESSFULLY,
     *        /3X,25HEXECUTED ON YOUR COMPUTER)
5000  FORMAT(//3X,29HMAX ERROR ON YOUR COMPUTER = ,E14.7,
     *        /3X,21HMAX ERROR BUILT IN = ,E14.7,
     *       //3X,30HMEAN ERROR ON YOUR COMPUTER = ,E14.7,
     *        /3X,22HMEAN ERROR BUILT IN = ,E14.7)
      END
      SUBROUTINE MACEPS(SRELPR)
C IT CARRIES OUT AN APPROXIMATION OF THE SINGLE RELATIVE PRECISION
C CONSTANT
      SRELPR=1.0
10    SRELPR=0.5*SRELPR
C THE DO LOOP STATEMENT IS NECESSARY WHEN THE ARITHMETIC UNIT HAS MORE
C BITS THAN IN STORAGE (INTEL 8087 FAMILY OF ARITHMETIC UNITS IS OF
C THIS TYPE), INFACT A DO LOOP INVOLVES A STORE FROM REGISTER TO MEMORY
C OF THE VALUE FOR SRELPR.
      DO 20 I=1,2
20    CONTINUE
      IF (SRELPR+1.0.GT.1.0) GOTO 10
      SRELPR=2.0*SRELPR
      RETURN
      END
      SUBROUTINE MMERR(NXI,NYI,ZI,ZI1,NDIM,ERMAX,ERMEAN)
C IT CARRIES OUT MAX ERROR AND MEAN ERROR AT PRESCRIBED RECTANGULAR
C GRID POINTS OF THE RECONSTRUCTED SURFACE.
C
C THE INPUT PARAMETERS ARE
C     NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE (MUST BE
C           1 OR GREATER),
C     NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE (MUST BE
C           1 OR GREATER),
C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           ANALYTIC  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     ZI1 = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           INTERPOLATED  Z  VALUES AT THE OUTPUT GRID POINTS ARE TO
C           BE STORED.
C     NDIM= DECLARED ROW DIMENSION OF THE ARRAYS ZI AND ZI1.
C
C THE OUTPUT PARAMETERS ARE
C     ERMAX = MAX ERROR
C     ERMEAN= MEAN ERROR
C
      DIMENSION ZI(NDIM,*),ZI1(NDIM,*)
      ERMAX=0.0
      SUM=0.0
      N=0
      DO 10 I=1,NXI
        DO 10 J=1,NYI
          IF (ZI1(I,J).EQ.1.E6) GOTO 10
          N=N+1
          ERR=ABS(ZI(I,J)-ZI1(I,J))
          SUM=SUM+ERR
          IF (ERR.GT.ERMAX) ERMAX=ERR
10    CONTINUE
      ERMEAN=SUM/FLOAT(N)
      RETURN
      END
      SUBROUTINE MASUB(SRELPR,IC,IEX,ND,XD,YD,ZD,TP,NXI,NYI,XI,YI,ZI,
     *                 NDIM,NOUT,IWK,WK)
C IT CARRIES OUT A SMOOTH BIVARIATE INTERPOLATION WHEN THE DATA
C POINTS PROJECTION IN THE  X-Y  PLANE IS IRREGULARLY DISTRIBUTED
C IN THE PLANE: IT ESTIMATES THE INTERPOLANT FUNCTION AT PRESCRIBED
C RECTANGULAR GRID POINTS BY MEANS OF THE 9-PARAMETERS DISCRETIZED
C VERSION OF THE NIELSON TRIANGULAR INTERPOLANT. THE REQUIRED DERIVATIVE
C VALUES ARE CALCULATED BY MINIMIZING SUITABLE TENSION'S FUNCTIONAL
C AND THE TRIANGULATION OF THE CONVEX HULL OF THE X-Y DATA IS REALIZED
C BY USING LAWSON'S LOCAL OPTIMIZATION PROCEDURE.
C IF EXTRAPOLATION IS REQUESTED (IEX PARAMETER) IT ADDS TO THE DATA
C POINTS SOME ADDITIONAL DATA, EVALUATED BY MEANS OF THE SHEPARD'S
C METHOD, SO THAT THE NEW CONVEX HULL OF THE DATA POINTS CONTAINS
C THE RECTANGULAR DOMAIN WHERE THE SURFACE HAS TO BE RECONSTRUCTED.
C
C
C THE INPUT PARAMETERS ARE
C     SRELPR = SINGLE RELATIVE PRECISION,
C     IC  = FLAG OF COMPUTATION (MUST BE 1,2,3 OR 4),
C         = 1 FOR THE FIRST CALL AND FOR NEW IEX,ND,XD-YD,
C         = 2 FOR NEW ZD,TP,
C         = 3 FOR NEW TP,
C         = 4 FOR NEW NXI,NYI,XI-YI,
C     IEX = FLAG OF COMPUTATION (MUST BE 0 OR 1),
C         = 1 (ONE) ALLOWS THE EXTRAPOLATION, 0 (ZERO) DOES NOT;
C           THIS PARAMETER MUST BE THE SAME IN THE FOLLOWING CALL,
C     ND  = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER);
C           GIVEN IEX AND ND WE DEFINE N AS
C                       N = ND + IEX*(2*INT(ND/25)+4),
C     XD  = ARRAY OF DIMENSION N CONTAINING THE  X  COORDINATES OF
C           THE ND DATA POINTS,
C     YD  = ARRAY OF DIMENSION N CONTAINING THE  Y  COORDINATES OF
C           THE ND DATA POINTS,
C     ZD  = ARRAY OF DIMENSION N CONTAINING THE  Z  COORDINATES OF
C           THE ND DATA POINTS,
C     TP  = TENSION PARAMETER (MUST BE GREATER THAN OR EQUAL TO ZERO).
C     NXI = NUMBER OF RECTANGULAR GRID POINTS IN THE X COORDINATE (MUST
C           BE 1 OR GREATER),
C     NYI = NUMBER OF RECTANGULAR GRID POINTS IN THE Y COORDINATE (MUST
C           BE 1 OR GREATER),
C     XI  = ARRAY OF DIMENSION NXI CONTAINING,IN ASCENDING ORDER, THE X
C           COORDINATES OF THE RECTANGULAR GRID POINTS,  WHERE THE
C           SURFACE HAS TO BE RECONSTRUCTED,
C     YI  = ARRAY OF DIMENSION NYI CONTAINING,IN ASCENDING ORDER, THE Y
C           COORDINATES OF THE RECTANGULAR GRID POINTS,  WHERE THE
C           SURFACE HAS TO BE RECONSTRUCTED,
C     NDIM= DECLARED ROW DIMENSION OF THE ARRAY CONTAINING ZI,
C     NOUT= LOGICAL UNIT NUMBER FOR THE STANDARD OUTPUT UNIT OF THE
C           SYSTEM.
C
C THE OUTPUT PARAMETERS ARE
C     XD,YD,ZD = ARRAYS OF DIMENSION N, CONTAINING THE X, Y AND Z
C                COORDINATES OF THE ND DATA POINTS AND,IF EXTRAPOLATION
C                IS REQUESTED (IEX=1), THE X, Y AND Z COORDINATES OF
C                N-ND ADDITIONAL POINTS EVALUATED BY MEANS OF SHEPARD'S
C                METHOD. IF EXTRAPOLATION IS NOT REQUESTED (IEX=0) THEN
C                N=ND AND XD, YD, ZD ARE NOT CHANGED,
C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION NXI*NYI, WHERE THE
C           INTERPOLATED  Z  VALUES AT THE RECTANGULAR GRID POINTS ARE
C           TO BE STORED.
C           IF EXTRAPOLATION IS NOT REQUESTED (IEX=0), THE Z VALUES AT
C           THE RECTANGULAR GRID POINTS OUTSIDE THE CONVEX HULL OF THE
C           DATA POINTS, ARE SET TO 1.0E6.
C
C THE OTHER PARAMETERS ARE
C     IWK = INTEGER ARRAY OF DIMENSION  23*N-31+MAX0(8*N+25,NXI*NYI)
C           USED INTERNALLY AS WORK AREA,
C     WK  = ARRAY OF DIMENSION  19*N-20  USED INTERNALLY AS WORK AREA.
C
C THE FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW ND VALUE
C AND/OR NEW CONTENTS OF THE XD AND YD ARRAYS MUST BE MADE WITH IC=1.
C THE CALL WITH IC=2 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME
C IEX,ND,NXI AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD,YD,
C XI AND YI ARRAYS.
C THE CALL WITH IC=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME
C IEX,ND,NXI AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD,YD,
C ZD,XI AND YI ARRAYS.
C THE CALL WITH IC=4 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME
C IEX AND ND VALUES AND WITH THE SAME CONTENTS OF THE XD,YD AND ZD
C ARRAYS.
C IWK AND WK ARRAYS MUST NOT BE DISTURBED BETWEEN THE CALL WITH IC NOT
C EQUAL TO  1  AND THE PRECEDING ONE.
C SUBROUTINE MASUB CALLS THE EXTRP,CTANG,ADJAC,PDSTE,PDMIN,ORDGR AND
C INTRP SUBROUTINES.
C
C DECLARATION STATEMENT.
      DIMENSION XD(*),YD(*),ZD(*),XI(*),YI(*),ZI(NDIM,*),IWK(*),WK(*)
      TOLL=1.0-SRELPR
C  ERROR CHECK.
8     IF(IC.LT.1.OR.IC.GT.4)GO TO 1000
      IF(IEX.LT.0.OR.IEX.GT.1)GO TO 1000
      IF(ND.LT.4)GO TO 1000
      IF(NXI.LT.1.OR.NYI.LT.1)GO TO 1000
      IF(IC.GT.1)GO TO 10
      IWK(1)=ND
      IWK(2)=IEX
      IWK(3)=NXI
      IWK(4)=NYI
      IWK(6)=ND
        IF(IEX.EQ.0)GO TO 50
        WK(1)=XI(1)
        WK(2)=XI(NXI)
        WK(3)=YI(1)
        WK(4)=YI(NYI)
        GO TO 40
10    IF(IWK(1).NE.ND) GO TO 1000
      IF(IWK(2).NE.IEX)GO TO 1000
      IF(IC.GT.3) GO TO 20
      IF(IWK(3).NE.NXI)GO TO 1000
      IF(IWK(4).NE.NYI)GO TO 1000
        IF(IC.EQ.2.AND.IEX.EQ.1)GO TO 40
        GO TO 50
20    IWK(3)=NXI
      IWK(4)=NYI
        IF(IEX.EQ.0)GO TO 50
        IF(XI(1).LT.WK(1).OR.XI(NXI).GT.WK(2))GO TO 1000
        IF(YI(1).LT.WK(3).OR.YI(NYI).GT.WK(4))GO TO 1000
        GO TO 50
C ADDS SOME EXTERNAL POINTS FOR THE EXTRAPOLATION.
40    CALL EXTRP(ND,XD,YD,ZD,IC,IWK,WK(1),WK(2),WK(3),WK(4))
50    NDP=IWK(6)
C ALLOCATION STORAGE AREAS IN THE IWK AND WK ARRAYS.
      JIPT=8
      JIPL=6*NDP-7
      JIWP=12*NDP-7
      JIND=13*NDP-7
      JNGP=19*NDP-22
      JIGP=23*NDP-32
      IPD=5
      IAL=2*NDP+5
      IBE=5*NDP-1
      IGA=8*NDP-7
      IEI=11*NDP-13
      IALS=14*NDP-19
      IBES=IALS+NDP
      IGAS=IBES+NDP
      IZX=IGAS+NDP
      IZY=IZX+NDP
      IF(IC.GT.1)GO TO 60
C TRIANGULATES THE X-Y PLANE.
      CALL CTANG(NDP,XD,YD,NOUT,NT,IWK(JIPT),IWK(JIPL),IWK(JIND),
     *           IWK(JIWP),WK(IPD))
      IWK(5)=NT
C CONSTRUCTS THE ADJACENCIES MONODIMENSIONAL ARRAY.
      CALL ADJAC(NT,IWK(JIPT),NDP,IWK(JIPL),IWK(JIWP-1))
60    NT=IWK(5)
      IF(IC.GT.3)GO TO 70
      IF(IC.EQ.3)GO TO 65
C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
      CALL PDSTE(NDP,XD,YD,ZD,NT,IWK(JIPT),WK(IPD),IWK(JIND))
65    CALL PDMIN(SRELPR,NDP,XD,YD,ZD,NOUT,IWK(JIPL),IWK(JIWP-1),
     *     IWK(JIND),TP,WK(IPD),WK(IAL),WK(IBE),WK(IGA),WK(IEI),
     *     WK(IALS),WK(IBES),WK(IGAS),WK(IZX),WK(IZY))
      IF(IC.GT.1)GO TO 80
C SORTS RECTANGULAR GRID POINTS ACCORDING TO THEIR BELONGING TO THE
C TRIANGLES.
70    CALL ORDGR(XD,YD,NT,IWK(JIPT),NXI,NYI,XI,YI,IWK(JNGP),
     *           IWK(JIGP))
80    DO 85 I=1,NXI
      DO 85 J=1,NYI
85    ZI(I,J)=1.E6
      ITPV=0
      IMAX=0
      I1MIN=NXI*NYI+1
      DO 120 KNGP=1,NT
        ITI=KNGP
        JWNGP=JNGP-1+KNGP
        NGP0=IWK(JWNGP)
        IF(NGP0.EQ.0)GO TO 100
        IMIN=IMAX+1
        IMAX=IMAX+NGP0
        DO 90 KIGP=IMIN,IMAX
          JWIGP=JIGP+KIGP-1
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI+1
          IXI=IZI-NXI*(IYI-1)
          CALL INTRP(SRELPR,TOLL,XD,YD,ZD,IWK(JIPT),WK(IPD),ITI,ITPV,
     *               XI(IXI),YI(IYI),ZI(IXI,IYI))
90      CONTINUE
100     JWNGP=JNGP+2*NT-KNGP
        NGP1=IWK(JWNGP)
        IF(NGP1.EQ.0)GO TO 120
        I1MAX=I1MIN-1
        I1MIN=I1MIN-NGP1
        DO 110 KIGP=I1MIN,I1MAX
          JWIGP=JIGP+KIGP-1
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI+1
          IXI=IZI-NXI*(IYI-1)
          CALL INTRP(SRELPR,TOLL,XD,YD,ZD,IWK(JIPT),WK(IPD),ITI,ITPV,
     *               XI(IXI),YI(IYI),ZI(IXI,IYI))
110     CONTINUE
120   CONTINUE
      RETURN
C ERROR EXIT.
1000  WRITE(NOUT,130)IEX,IC,ND
130   FORMAT(1X,34HIMPROPER INPUT PARAMETER VALUE(S)./
     *       4HIEX=,I4,4H IC=,I4,4H ND=,I4)
      RETURN
      END
      SUBROUTINE EXTRP(ND,X,Y,Z,KC,IWK,A,B,C,D)
C IT ADDS SOME EXTERNAL POINTS TO THE DATA POINT SET AND ESTIMATES
C THE Z COORDINATE AT THESE POINTS FOLLOWING THE SHEPARD METHOD.
C
C THE INPUT PARAMETERS ARE
C     ND  = NUMBER OF DATA POINTS,
C     X,Y,Z = ARRAYS OF DIMENSION  ND  CONTAINING THE  X,Y  AND  Z
C             COORDINATES OF THE DATA POINTS,
C     KC  = FLAG OF COMPUTATION,
C     A,B,C,D = EXTREME OF THE RECTANGULAR GRID POINTS.
C
C THE OUTPUT PARAMETERS ARE
C     X,Y,Z = ARRAYS OF DIMENSION  NDP  CONTAINING THE  X,Y  AND  Z
C             COORDINATES OF THE NEW SET OF DATA POINTS WHERE
C                 NDP = ND+2*INT(ND/25)+4,
C     A,B,C,D = EXTREME OF THE RECTANGULAR REGION CONTAINING THE DATA
C               POINT SET AND THE RECTANGULAR GRID POINTS.
C
C THE OTHER PARAMETER IS
C     IWK = INTEGER ARRAY OF DIMENSION 4*ND+7 USED AS WORK AREA.
C
C DECLARATION STATEMENT.
      DIMENSION X(*),Y(*),Z(*),IWK(*),DIST(5),IPC0(5)
C STATEMENT FUNCTIONS.
      DINF(U1,V1,U2,V2)=AMAX1(ABS(U1-U2),ABS(V1-V2))
      DSQF(U1,V1,U2,V2)=((U2-U1)**2+(V2-V1)**2)**2
      NDP=IWK(6)
      NCP=IWK(7)
      IF(KC.EQ.2)GO TO 200
C  ESTIMATES THE SMALLEST RECTANGLE CONTAINING THE DATA POINT SET
C  AND THE RECTANGULAR GRID POINTS.
      JA=0
      IA=7
      JB=0
      IB=IA+ND
      JC=0
      IC=IB+ND
      JD=0
      ID=IC+ND
      DO 50 I=1,ND
        IF(X(I)-A)16,18,20
16        JA=0
18        JA=JA+1
          IWK(IA+JA)=I
          A=X(I)
          GO TO 30
20      IF(B-X(I))26,28,30
26        JB=0
28        JB=JB+1
          IWK(IB+JB)=I
          B=X(I)
30      IF(Y(I)-C)36,38,40
36        JC=0
38        JC=JC+1
          IWK(IC+JC)=I
          C=Y(I)
          GO TO 50
40      IF(D-Y(I))46,48,50
46        JD=0
48        JD=JD+1
          IWK(ID+JD)=I
          D=Y(I)
50    CONTINUE
C  ESTIMATES THE NUMBER OF POINTS AND WHERE THEY HAVE TO BE ADJOINTED.
      N=ND/25
      BMA=B-A
      DMC=D-C
      HX=BMA
      HY=DMC
      NX=1
      NY=1
      IF(N.EQ.0)GO TO 75
        DO 70 I=1,N
          IF(HX.GT.HY)GO TO 60
            NY=NY+1
            HY=DMC/FLOAT(NY)
          GO TO 70
60          NX=NX+1
            HX=BMA/FLOAT(NX)
70      CONTINUE
      HX=BMA/FLOAT(NX)
      HY=DMC/FLOAT(NY)
C  ADDS THE NEW EXTERNAL POINTS AND CHECKS THAT THEY ARE NOT
C  COINCIDENT WITH THE OLD ONES.
75    NDP=ND+1
      YP=C
      DO 100 I=1,NY
        IF(JA.EQ.0)GO TO 80
          DO 77 J=1,JA
            IF(YP.EQ.Y(IWK(IA+J)))GO TO 90
77        CONTINUE
80      X(NDP)=A
        Y(NDP)=YP
        NDP=NDP+1
90      YP=YP+HY
100   CONTINUE
      XP=A
      DO 130 I=1,NX
        IF(JD.EQ.0)GO TO 110
          DO 105 J=1,JD
            IF(XP.EQ.X(IWK(ID+J)))GO TO 120
105       CONTINUE
110     X(NDP)=XP
        Y(NDP)=D
        NDP=NDP+1
120     XP=XP+HX
130   CONTINUE
      YP=D
      DO 160 I=1,NY
        IF(JB.EQ.0)GO TO 140
          DO 135 J=1,JB
            IF(YP.EQ.Y(IWK(IB+J)))GO TO 150
135       CONTINUE
140     X(NDP)=B
        Y(NDP)=YP
        NDP=NDP+1
150     YP=YP-HY
160   CONTINUE
      XP=B
      DO 190 I=1,NX
        IF(JC.EQ.0)GO TO 170
          DO 165 J=1,JC
            IF(XP.EQ.X(IWK(IC+J)))GO TO 180
165       CONTINUE
170     X(NDP)=XP
        Y(NDP)=C
        NDP=NDP+1
180     XP=XP-HX
190   CONTINUE
      NDP=NDP-1
      IWK(6)=NDP
      NCP=5
      IF(ND.LE.5)NCP=3
      IWK(7)=NCP
C  ESTIMATES THE FUNCTION VALUE AT THE NEW EXTERNAL POINTS.
C
200   N1=ND+1
        DO 270 IP1=N1,NDP
          X0=X(IP1)
          Y0=Y(IP1)
          DMX=0.
            DO 220 IP2=1,ND
              DM=DINF(X0,Y0,X(IP2),Y(IP2))
              DIST(IP2)=DM
              IPC0(IP2)=IP2
                IF(DM.LE.DMX)GO TO 210
                  DMX=DM
                  JMX=IP2
210             IF(IP2.GE.NCP)GO TO 230
220         CONTINUE
230       IP2P1=IP2+1
            DO 250 IP2=IP2P1,ND
              DM=DINF(X0,Y0,X(IP2),Y(IP2))
                IF(DM.GE.DMX)GO TO 250
                  DIST(JMX)=DM
                  IPC0(JMX)=IP2
                  DMX=0.
                    DO 240 J1=1,NCP
                      IF(DIST(J1).LE.DMX)GO TO 240
                        DMX=DIST(J1)
                        JMX=J1
240                 CONTINUE
250         CONTINUE
        ANUM=0.
        ADEN=0.
          DO 260 J2=1,NCP
            IP2=IPC0(J2)
            R4=DSQF(X0,Y0,X(IP2),Y(IP2))
              IF(R4.EQ.0)GO TO 260
                ANUM=ANUM+Z(IP2)/R4
                ADEN=ADEN+1/R4
260       CONTINUE
        Z(IP1)=ANUM/ADEN
270   CONTINUE
      RETURN
      END
      SUBROUTINE CTANG(NDP,XD,YD,NOUT,NT,IPT,IPL,IWL,IWP,WK)
C IT CARRIES OUT TRIANGULATION BY DIVIDING THE X-Y PLANE INTO A
C NUMBER OF TRIANGLES ACCORDING TO THE GIVEN DATA POINTS IN THE PLANE.
C AT THE END OF THE OPERATIONS,THE INDICES OF THE VERTEXES OF THE
C TRIANGLES ARE LISTED COUNTER-CLOCKWISE.
C SUBROUTINE CTANG CALLS THE MAXMN FUNCTION.
C
C
C THE INPUT PARAMETERS ARE
C     NDP = NUMBER OF DATA POINTS,
C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE X COORDINATES OF THE
C           DATA POINTS,
C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE Y COORDINATES OF THE
C           DATA POINTS,
C     NOUT= LOGICAL UNIT NUMBER FOR THE STANDARD OUTPUT UNIT OF THE
C           SYSTEM.
C
C THE OUTPUT PARAMETERS ARE
C     NT  = NUMBER OF TRIANGLES,
C     IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE INDICES OF
C           THE VERTEXES OF THE (IT)TH TRIANGLE ARE TO BE STORED AS THE
C           (3*IT-2)ND, (3*IT-1)ST AND (3*IT)TH ELEMENTS, IT=1,2,..,NT.
C
C THE OTHER PARAMETERS ARE
C     IPL = INTEGER ARRAY OF DIMENSION 6*NDP USED INTERNALLY AS WORK
C           AREA,
C     IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED INTERNALLY AS WORK
C           AREA,
C     IWP = INTEGER ARRAY OF DIMENSION NDP USED INTERNALLY AS WORK
C           AREA,
C     WK  = ARRAY OF DIMENSION NDP USED INTERNALLY AS WORK AREA.
C
C DECLARATION STATEMENTS.
      DIMENSION XD(*),YD(*),IPT(*),IPL(*),IWL(*),IWP(*),WK(*)
      DIMENSION ITF(2)
      DATA  RATIO/1.0E-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 MIDPOINTS.
   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.0)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.0
      YDMP=(YD(IPMN1)+YD(IPMN2))/2.0
C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF DISTANCE
C FROM THE MIDPOINTS AND STORES THE STORED DATA POINTS NUMBERS IN THE
C 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 SO 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 VERTEXES OF
C THE TRIANGLES IN THE IPT ARRAY, AND STORES POINT NUMBERS OF THE
C BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN 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.0)       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 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 SO THAT THE INVISIBLE BORDER LINE
C SEGMENTS ARE 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 SEGMENTS IN
C THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER LINE SEGMENTS TO BE
C 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 THE 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 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 WHETHER THE EXCHANGE IS NECESSARY.
   63     IF(MAXMN(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 THE 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 THE
C 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 ON THE
C 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 WHETHER THE EXCHANGE IS NECESSARY.
   74 IF(MAXMN(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 ARE
C 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.0)     GO TO 81
        IPT(ITT3-2)=IP2
        IPT(ITT3-1)=IP1
   81 CONTINUE
      NT=NT0
      RETURN
C
C ERROR EXIT.
   90 WRITE (NOUT,2090) NDP0
      GO TO 93
   91 WRITE (NOUT,2091) NDP0,IP1,IP2,X1,Y1
      GO TO 93
   92 WRITE (NOUT,2092) NDP0
   93 WRITE (NOUT,2093)
      NT=0
      STOP
C FORMAT STATEMENTS.
 2090 FORMAT(1X/23H ***   NDP LESS THAN 4./7H  NDP =,I5)
 2091 FORMAT(1X/29H ***   IDENTICAL DATA POINTS./
     1       7H  NDP =,I5,5X,5HIP1 =,I5,5X,5HIP2 =,
     2       I5,5X,4HXD =,E12.4,5X,4HYD =,E12.4)
 2092 FORMAT(1X/33H ***   ALL COLLINEAR DATA POINTS./
     1       7H  NDP =,I5)
 2093 FORMAT(35H ERROR DETECTED IN ROUTINE    CTANG)
      END
      FUNCTION MAXMN(X,Y,I1,I2,I3,I4)
C IT DETERMINES WHETHER THE EXCHANGE OF TWO TRIANGLES IS NECESSARY
C OR NOT ON THE BASIS OF THE MAX-MIN-ANGLE CRITERION BY C.LAWSON.
C
C THE INPUT PARAMETERS ARE
C     X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA POINTS,
C     I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1,P2,P3 AND P4
C                   FORMING A QUADRILATERAL WITH P3 AND P4
C                   DIAGONALLY CONNECTED.
C
C FUNCTION MAXMN RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EXCHANGE IS
C NECESSARY, OTHERWISE 0 (ZERO).
C
C DECLARATION STATEMENT.
      DIMENSION X(*),Y(*)
      EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ),
     1            (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ)
C PRELIMINARY PROCESSING.
      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.
      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.0)GO TO 10
      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*AMAX1(A1SQ,B1SQ))
      S2SQ=U2*U2/(C2SQ*AMAX1(A2SQ,B2SQ))
      S3SQ=U3*U3/(C3SQ*AMAX1(A3SQ,B3SQ))
      S4SQ=U4*U4/(C4SQ*AMAX1(A4SQ,B4SQ))
      IF(AMIN1(S1SQ,S2SQ).LT.AMIN1(S3SQ,S4SQ))IDX=1
   10 MAXMN=IDX
      RETURN
      END
      SUBROUTINE ADJAC(NT,IPT,N,IADVE,NADVE)
C IT ESTIMATES THE ADJACENCIES MONODIMENSIONAL ARRAY CONTAINING
C FOR EACH VERTEX THE INDICES OF THE VERTEXES ADJACENT IN THE
C TRIANGULATION.
C
C THE INPUT PARAMETERS ARE
C     NT  = NUMBER OF TRIANGLES,
C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE INDICES
C           OF THE VERTEXES OF THE TRIANGLES,
C     N   = NUMBER OF DATA POINTS.
C
C THE OUTPUT PARAMETERS ARE
C     IADVE = INTEGER ARRAY OF DIMENSION  6*N-12  CONTAINING FOR EACH
C             VERTEX THE INDICES OF THE VERTEXES ADJACENT IN THE
C             TRIANGULATION,
C     NADVE = INTEGER ARRAY OF DIMENSION  N+1  CONTAINING FOR EACH
C             VERTEX THE NUMBER OF THE VERTEXES ADJACENT IN THE
C             TRIANGULATION.
C
C  DECLARATION STATEMENT.
      DIMENSION IPT(*),IADVE(*),NADVE(*),ITEM(30)
      NADVE(1)=0
      KIN=0
      NT3=3*NT
        DO 70 I=1,N
          I2=0
C  STORES THE INDICES OF THE ADJACENT VERTEXES.
          DO 30 J1=1,NT3,3
            J2=J1+2
            DO 10 J=J1,J2
              IF(I.EQ.IPT(J))GO TO 20
   10       CONTINUE
          GO TO 30
   20       I1=I2+1
            I2=I2+2
            ITEM(I1)=IPT(J1+MOD(J,3))
            ITEM(I2)=IPT(J1+MOD(J+1,3))
   30     CONTINUE
C  DISCARDS THE INDICES THAT HAVE BEEN STORED TWICE.
          JIN=KIN+1
          KIN=KIN+2
          JFIN=KIN
          IADVE(JIN)=ITEM(1)
          IADVE(JFIN)=ITEM(2)
            IF(I2.EQ.2)GO TO 60
              DO 50 J=3,I2
                DO 40 L=JIN,JFIN
                  IF(ITEM(J).EQ.IADVE(L))GO TO 50
   40           CONTINUE
                KIN=KIN+1
                IADVE(KIN)=ITEM(J)
                JFIN=KIN
   50         CONTINUE
   60       NADVE(I+1)=KIN
   70   CONTINUE
      RETURN
      END
      SUBROUTINE PDSTE(N,X,Y,Z,NT,IPT,PD,IPD)
C IT ESTIMATES THE FIRST ORDER PARTIAL DERIVATIVE VALUES
C AT THE DATA POINTS FOLLOWING THE KLUCEWICZ METHOD.
C
C THE INPUT PARAMETERS ARE
C     N   = NUMBER OF DATA POINTS,
C     X,Y,Z = ARRAYS OF DIMENSION N CONTAINING THE X,Y AND Z
C             COORDINATES OF THE DATA POINTS,
C     NT  = NUMBER OF TRIANGLES,
C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE INDICES
C           OF THE VERTEXES OF THE TRIANGLES.
C
C THE OUTPUT PARAMETER IS
C     PD  = ARRAY OF DIMENSION 2*N CONTAINING THE PARTIAL DERIVATIVE
C           VALUES AT THE DATA POINTS.
C
C THE OTHER PARAMETER IS
C     IPD = INTEGER ARRAY OF DIMENSION N USED INTERNALLY AS WORK AREA.
C
C DECLARATION STATEMENT.
      DIMENSION X(*),Y(*),Z(*),IPT(*),PD(*),IPD(*)
C PRELIMINARY PROCESSING.
      N2=N+N
        DO 10 I=1,N2
          PD(I)=0.
          IPD(I)=0
10      CONTINUE
C ESTIMATES FOR EACH TRIANGLE THE SLOPES OF THE PLANE THROUGH THE
C FUNCTION'S VALUE AT THE VERTEXES.
      DO 20 I=1,NT
        L=3*I
        I1=IPT(L-2)
        I2=IPT(L-1)
        I3=IPT(L)
        X21=X(I2)-X(I1)
        X31=X(I3)-X(I1)
        Y21=Y(I2)-Y(I1)
        Y31=Y(I3)-Y(I1)
        Z21=Z(I2)-Z(I1)
        Z31=Z(I3)-Z(I1)
        C=Y21*X31-X21*Y31
        DX=(Y21*Z31-Z21*Y31)/C
        DY=(Z21*X31-X21*Z31)/C
C  UPDATES THE IPD AND PD ARRAYS.
        IPD(I1)=IPD(I1)+1
        I1=I1+I1
        PD(I1-1)=PD(I1-1)+DX
        PD(I1)=PD(I1)+DY
        IPD(I2)=IPD(I2)+1
        I2=I2+I2
        PD(I2-1)=PD(I2-1)+DX
        PD(I2)=PD(I2)+DY
        IPD(I3)=IPD(I3)+1
        I3=I3+I3
        PD(I3-1)=PD(I3-1)+DX
        PD(I3)=PD(I3)+DY
20    CONTINUE
C  AVERAGES THE DERIVATIVE VALUES STORED IN THE PD ARRAY.
      J=0
        DO 30 I=2,N2,2
          J=J+1
          DEN=IPD(J)
          PD(I-1)=PD(I-1)/DEN
          PD(I)=PD(I)/DEN
30      CONTINUE
      RETURN
      END
      SUBROUTINE PDMIN(SRELPR,N,X,Y,Z,NOUT,IADVE,NADVE,INDEX,TP,PD,
     *                 ALFA,BETA,GAMMA,EIJQ,ALFAS,BETAS,GAMMAS,ZX,ZY)
C IT ESTIMATES THE FIRST ORDER PARTIAL DERIVATIVE VALUES AT
C THE DATA POINTS BY MEANS OF A GLOBAL METHOD BASED ON A MINIMUM
C NORM UNDER TENSION NETWORK .
C
C THE INPUT PARAMETERS ARE
C     SRELPR = SINGLE RELATIVE PRECISION
C     N   = NUMBER OF DATA POINTS,
C     X,Y,Z = ARRAY OF DIMENSION N CONTAINING THE X,Y AND Z
C             COORDINATES OF THE DATA POINTS,
C     NOUT = LOGICAL UNIT NUMBER FOR THE STANDARD OUTPUT UNIT OF THE
C            SYSTEM,
C     IADVE = INTEGER ARRAY OF DIMENSION 6*N-12 CONTAINING THE INDICES
C             OF THE VERTEXES ADJACENT TO EACH VERTEX IN THE
C             TRIANGULATION,
C     NADVE = INTEGER ARRAY OF DIMENSION N+1 CONTAINING THE NUMBER OF
C             THE VERTEXES ADJACENT TO EACH VERTEX IN THE TRIANGULATION
C     TP  = TENSION PARAMETER,
C     PD  = ARRAY OF DIMENSION 2*N CONTAINING AN INITIAL EVALUATION
C           OF THE PARTIAL DERIVATIVE VALUES AT THE DATA POINTS.
C
C THE OUTPUT PARAMETER IS
C     PD  = ARRAY OF DIMENSION 2*N CONTAINING THE PARTIAL DERIVATIVE
C           VALUES AT THE DATA POINTS.
C
C THE OTHER PARAMETERS ARE
C     INDEX = INTEGER ARRAY OF DIMENSION 6*N-15 USED INTERNALLY AS
C             WORK AREA,
C     ALFA,BETA,GAMMA,EIJQ = ARRAYS OF DIMENSION 3*N-6 USED INTERNALLY
C                            AS WORK AREAS,
C     ALFAS,BETAS,GAMMAS,ZX,ZY = ARRAYS OF DIMENSION N USED INTERNALLY
C                                AS WORK AREAS.
C
C  THE  RELER  CONSTANT IN THE DATA  INITIALIZATION STATEMENT  IS A
C  RELATIVE ERROR TOLERANCE TO STOP THE ITERATIVE METHOD.
C  THEREFORE IT IS MACHINE DEPENDENT;
C  THE ABSOLUTE ERROR TOLERANCE  TAU  IS THEN OBTAINED BY
C         TAU=RELER*AMAX1(ABS(PD(I)),I=1,2*N)+2*N*SRELPR.
C
C DECLARATION STATEMENTS.
      DIMENSION X(*),Y(*),Z(*),IADVE(*),NADVE(*),PD(*)
      DIMENSION INDEX(*),ALFA(*),BETA(*),GAMMA(*),EIJQ(*),
     *          ALFAS(*),BETAS(*),GAMMAS(*),ZX(*),ZY(*)
      DATA RELER/1.E-5/
C  CALCULATES THE PART OF MATRIX COEFFICIENTS INDEPENDENT
C  FROM THE TENSION PARAMETER TP.
      K=0
      PDM=0.0
        DO 60 I=1,N
          J=I+I
          PDM=AMAX1(PDM,ABS(PD(J)),ABS(PD(J-1)))
          ZX(I)=0.
          ZY(I)=0.
          ALFAS(I)=0.
          BETAS(I)=0.
          GAMMAS(I)=0.
          JIN=NADVE(I)+1
          JFIN=NADVE(I+1)
            DO 50 J=JIN,JFIN
              IND=IADVE(J)
              DX=X(I)-X(IND)
              DY=Y(I)-Y(IND)
              DZ=Z(I)-Z(IND)
              DXQ=DX*DX
              DYQ=DY*DY
              ALUNQ=DXQ+DYQ
              AL3=ALUNQ*SQRT(ALUNQ)
              ZX(I)=ZX(I)+DZ*DX/AL3
              ZY(I)=ZY(I)+DZ*DY/AL3
                IF(IND.GT.I)GO TO 30
                  LIN=NADVE(IND)+1
                  LFIN=NADVE(IND+1)
                    DO 10 L=LIN,LFIN
                      IF(I.EQ.IADVE(L))GO TO 20
   10               CONTINUE
   20             INDEX(J)=INDEX(L)
                  GO TO 40
   30         K=K+1
              INDEX(J)=K
              AL3P2=AL3+AL3
              EIJQ(K)=ALUNQ
              ALFA(K)=DXQ/AL3P2
              BETA(K)=DX*DY/AL3P2
              GAMMA(K)=DYQ/AL3P2
   40       ALFAS(I)=ALFAS(I)+ALFA(INDEX(J))
            BETAS(I)=BETAS(I)+BETA(INDEX(J))
            GAMMAS(I)=GAMMAS(I)+GAMMA(INDEX(J))
   50       CONTINUE
          ZX(I)=3.*ZX(I)/2.
          ZY(I)=3.*ZY(I)/2.
          ALFAS(I)=2.*ALFAS(I)
          BETAS(I)=2.*BETAS(I)
          GAMMAS(I)=2.*GAMMAS(I)
   60   CONTINUE
      IF(TP.EQ.0)GO TO 100
C  CALCULATES THE PART OF MATRIX COEFFICIENTS DEPENDING FROM
C  THE TENSION PARAMETER TP.
      TPQ=TP*TP
      TPQ60=TPQ/60.
      TPQ40=TPQ/40.
      TPQ15=TPQ/15.
        DO 80 I=1,N
          AL=0.
          BE=0.
          GA=0.
          Z1=0.
          Z2=0.
          JIN=NADVE(I)+1
          JFIN=NADVE(I+1)
            DO 70 J=JIN,JFIN
              K=INDEX(J)
              IND=IADVE(J)
              DZ=Z(I)-Z(IND)
              AL=AL+ALFA(K)*EIJQ(K)
              BE=BE+BETA(K)*EIJQ(K)
              GA=GA+GAMMA(K)*EIJQ(K)
              SQEIJQ=SQRT(EIJQ(K))
              Z1=Z1+(X(I)-X(IND))*DZ/SQEIJQ
              Z2=Z2+(Y(I)-Y(IND))*DZ/SQEIJQ
   70       CONTINUE
          ALFAS(I)=ALFAS(I)+TPQ15*AL
          BETAS(I)=BETAS(I)+TPQ15*BE
          GAMMAS(I)=GAMMAS(I)+TPQ15*GA
          ZX(I)=ZX(I)+TPQ40*Z1
          ZY(I)=ZY(I)+TPQ40*Z2
   80   CONTINUE
        M=NADVE(N+1)/2
          DO 90 I=1,M
            ALFA(I)=ALFA(I)-TPQ60*ALFA(I)*EIJQ(I)
            BETA(I)=BETA(I)-TPQ60*BETA(I)*EIJQ(I)
            GAMMA(I)=GAMMA(I)-TPQ60*GAMMA(I)*EIJQ(I)
   90     CONTINUE
C  CALCULATES THE SOLUTIONS OF THE SYSTEM FOLLOWING THE GAUSS-SIEDEL
C  METHOD.
  100 ITER=1
      TAU=RELER*PDM+2*N*SRELPR
  110 ERQ=0.
        DO 130 I=1,N
        ADX=0.
        BDX=0.
        BDY=0.
        GDY=0.
        JIN=NADVE(I)+1
        JFIN=NADVE(I+1)
          DO 120 J=JIN,JFIN
            IND=2*IADVE(J)
            K=INDEX(J)
            ADX=ADX+ALFA(K)*PD(IND-1)
            BDY=BDY+BETA(K)*PD(IND)
            BDX=BDX+BETA(K)*PD(IND-1)
            GDY=GDY+GAMMA(K)*PD(IND)
  120     CONTINUE
        DET=ALFAS(I)*GAMMAS(I)-BETAS(I)*BETAS(I)
        S1=ADX+BDY-ZX(I)
        S2=BDX+GDY-ZY(I)
        PDXN=(-GAMMAS(I)*S1+BETAS(I)*S2)/DET
        PDYN=(BETAS(I)*S1-ALFAS(I)*S2)/DET
        IPI=I+I
        ERX=PDXN-PD(IPI-1)
        ERY=PDYN-PD(IPI)
        ERQ=ERQ+ERX*ERX+ERY*ERY
        PD(IPI-1)=PDXN
        PD(IPI)=PDYN
  130 CONTINUE
C  CHECKS WHETHER CONVERGENCE IS REACHED WITH THE PRESCRIBED TOLERANCE.
      IF(ERQ.LT.TAU)RETURN
        IF(ITER.EQ.20)GO TO 150
          ITER=ITER+1
          GO TO 110
C
C  ERROR EXIT.
  150 WRITE(NOUT,1)
      RETURN
    1 FORMAT(1X,26HMINIMIZATION NOT COMPLETED)
      END
      SUBROUTINE ORDGR(XD,YD,NT,IPT,NXI,NYI,XI,YI,NGP,IGP)
C IT ORGANIZES GRID POINTS FOR SURFACE RECONSTRUCTION BY
C SORTING THEM ACCORDING TO THEIR BELONGING TO THE TRIANGLES.
C
C THE INPUT PARAMETERS ARE
C     XD,YD = ARRAY OF DIMENSION N CONTAINING THE X AND Y COORDINATES
C             OF THE DATA POINTS, WHERE N IS THE NUMBER OF THE DATA
C             POINTS,
C     NT  = NUMBER OF TRIANGLES,
C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE INDICES OF
C           THE VERTEXES OF THE TRIANGLES,
C     NXI = NUMBER OF GRID POINTS IN THE X COORDINATES,
C     NYI = NUMBER OF GRID POINTS IN THE Y COORDINATES,
C     XI,YI = ARRAY OF DIMENSION NXI AND NYI CONTAINING THE X AND Y
C             COORDINATES OF THE GRID POINTS,RESPECTIVELY.
C
C THE OUTPUT PARAMETERS ARE
C     NGP = INTEGER ARRAY OF DIMENSION 2*NT WHERE THE NUMBER OF GRID
C           POINTS BELONGING TO EACH TRIANGLE IS TO BE STORED,
C     IGP = INTEGER ARRAY OF DIMENSION NXI*NYI WHERE THE INDICES OF THE
C           GRID POINTS ARE TO BE STORED ACCORDING TO THEIR BELONGING
C           TO THE TRIANGLES CONSIDERED IN ASCENDING ORDER NUMBERS.
C
C DECLARATION STATEMENT.
C
      DIMENSION XD(*),YD(*),IPT(*),XI(*),YI(*),NGP(*),IGP(*)
C
C STATEMENT FUNCTION.
      SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3)-(V1-V3)*(U2-U3)
      NT0 = NT
C
C PRELIMINARY PROCESSING.
      NXI0 = NXI
      NYI0 = NYI
      NXINYI = NXI0*NYI0
C
C DETERMINES GRID POINTS INSIDE THE DATA AREA.
      JNGP0 = 0
      JNGP1 = 2*NT0 + 1
      JIGP0 = 0
      JIGP1 = NXINYI + 1
      DO 160 IT0=1,NT0
        NGP0 = 0
        NGP1 = 0
        IT0T3 = IT0*3
        IP1 = IPT(IT0T3-2)
        IP2 = IPT(IT0T3-1)
        IP3 = IPT(IT0T3)
        X1 = XD(IP1)
        Y1 = YD(IP1)
        X2 = XD(IP2)
        Y2 = YD(IP2)
        X3 = XD(IP3)
        Y3 = YD(IP3)
        XMN = AMIN1(X1,X2,X3)
        XMX = AMAX1(X1,X2,X3)
        YMN = AMIN1(Y1,Y2,Y3)
        YMX = AMAX1(Y1,Y2,Y3)
        INSD = 0
        DO 20 IXI=1,NXI0
          IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 10
          IF (INSD.EQ.0) GO TO 20
          IXIMX = IXI - 1
          GO TO 30
   10     IF (INSD.EQ.1) GO TO 20
          INSD = 1
          IXIMN = IXI
   20   CONTINUE
        IF (INSD.EQ.0) GO TO 150
        IXIMX = NXI0
   30   DO 140 IYI=1,NYI0
          YII = YI(IYI)
          IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 140
          DO 130 IXI=IXIMN,IXIMX
            XII = XI(IXI)
            L = 0
            IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 130, 40, 50
   40       L = 1
   50       IF (SIDE(X2,Y2,X3,Y3,XII,YII)) 130, 60, 70
   60       L = 1
   70       IF (SIDE(X3,Y3,X1,Y1,XII,YII)) 130, 80, 90
   80       L = 1
   90       IZI = NXI0*(IYI-1) + IXI
            IF (L.EQ.1) GO TO 100
            NGP0 = NGP0 + 1
            JIGP0 = JIGP0 + 1
            IGP(JIGP0) = IZI
            GO TO 130
  100       IF (JIGP1.GT.NXINYI) GO TO 120
            DO 110 JIGP1I=JIGP1,NXINYI
              IF (IZI.EQ.IGP(JIGP1I)) GO TO 130
  110       CONTINUE
  120       NGP1 = NGP1 + 1
            JIGP1 = JIGP1 - 1
            IGP(JIGP1) = IZI
  130     CONTINUE
  140   CONTINUE
  150   JNGP0 = JNGP0 + 1
        NGP(JNGP0) = NGP0
        JNGP1 = JNGP1 - 1
        NGP(JNGP1) = NGP1
  160 CONTINUE
      RETURN
      END
      SUBROUTINE INTRP(SRELPR,TOL,X,Y,Z,IPT,PD,ITI,ITPV,XII,YII,ZII)
C IT CARRIES OUT PUNCTUAL INTERPOLATION, I.E., IT DETERMINES
C THE Z VALUE AT A GIVEN POINT IN A TRIANGLE BY MEANS OF THE
C 9-PARAMETER DISCRETIZED VERSION OF NIELSON'S SCHEME.
C
C THE INPUT PARAMETERS ARE
C    SRELPR,TOL = SINGLE RELATIVE PRECISION AND TOLERANCE,
C    X,Y,Z= ARRAYS OF DIMENSION N CONTAINING THE X,Y AND Z COORDINATES
C         OF THE DATA POINTS, WHERE N IS THE NUMBER OF THE DATA POINTS,
C    IPT = INTEGER ARRAY OF DIMENSION 3*NT, WHERE NT IS THE NUMBER OF
C          TRIANGLES, CONTAINING THE INDICES OF THE VERTEXES OF THE
C          TRIANGLES-THEMSELVES,
C    PD  = ARRAY OF DIMENSION 2*N CONTAINING THE PARTIAL DERIVATIVE
C          VALUES AT THE DATA POINTS,
C    ITI = INDEX OF THE TRIANGLE WHERE THE POINT FOR WHICH
C          INTERPOLATION HAS TO BE PERFORMED, LIES,
C    ITPV= INDEX OF THE TRIANGLE CONSIDERED IN THE PREVIOUS CALL,
C    XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH INTER-
C              POLATION HAS TO BE PERFORMED.
C
C THE OUTPUT PARAMETER IS
C    ZII = INTERPOLATED Z VALUE.
C
C DECLARATION STATEMENTS.
      DIMENSION X(*),Y(*),Z(*),IPT(*),PD(*)
      SAVE
      ITO=ITI
        IF(ITO.EQ.ITPV)GO TO 10
C  SELECTS THE TRIANGLE CONTAINING THE POINT (XII,YII).
          ITPV=ITO
          ITO3=3*ITO
          IND1=IPT(ITO3-2)
          IND2=IPT(ITO3-1)
          IND3=IPT(ITO3)
C  CALCULATES THE BASIC QUANTITIES RELATIVES TO THE SELECTED TRIANGLE.
          X12=X(IND1)-X(IND2)
          X13=X(IND1)-X(IND3)
          X23=X(IND2)-X(IND3)
          Y12=Y(IND1)-Y(IND2)
          Y13=Y(IND1)-Y(IND3)
          Y23=Y(IND2)-Y(IND3)
          D=X13*Y23-X23*Y13
          INDEX=2*IND1
          D3V1=-X13*PD(INDEX-1)-Y13*PD(INDEX)
          D2V1=-X12*PD(INDEX-1)-Y12*PD(INDEX)
          INDEX=2*IND2
          D1V2=X12*PD(INDEX-1)+Y12*PD(INDEX)
          D3V2=-X23*PD(INDEX-1)-Y23*PD(INDEX)
          INDEX=2*IND3
          D2V3=X23*PD(INDEX-1)+Y23*PD(INDEX)
          D1V3=X13*PD(INDEX-1)+Y13*PD(INDEX)
          E1=X23*X23+Y23*Y23
          E2=X13*X13+Y13*Y13
          E3=X12*X12+Y12*Y12
          E=2.*E1
          ALFA21=(E2+E1-E3)/E
          ALFA31=(E3+E1-E2)/E
          E=2.*E2
          ALFA12=(E1+E2-E3)/E
          ALFA32=(E3+E2-E1)/E
          E=2.*E3
          ALFA13=(E1+E3-E2)/E
          ALFA23=(E2+E3-E1)/E
C  CALCULATES THE REMAINING QUANTITIES NECESSARY FOR THE ZI EVALUATION
C  DEPENDING FROM THE INTERPOLATION POINT.
   10   B1=((XII-X(IND3))*Y23-(YII-Y(IND3))*X23)/D
        B2=((YII-Y(IND1))*X13-(XII-X(IND1))*Y13)/D
        B3=1.-B1-B2
          IF(B1.GE.TOL)GO TO 30
            IF(B2.GE.TOL)GO TO 40
              IF(B1+B2.LE.SRELPR)GO TO 50
        W=B1*B2*B3/(B1*B2+B1*B3+B2*B3)
        QB=B1*B1
        WB=W*B1
        ZII=Z(IND1)*(QB*(3.-2.*B1)+6.*WB*(B3*ALFA12+B2*ALFA13))+
     *      D3V1*(QB*B3+WB*(3.*B3*ALFA12+B2-B3))+
     *      D2V1*(QB*B2+WB*(3.*B2*ALFA13+B3-B2))
        QB=B2*B2
        WB=W*B2
        ZII=ZII+Z(IND2)*(QB*(3.-2.*B2)+6.*WB*(B1*ALFA23+B3*ALFA21))+
     *      D1V2*(QB*B1+WB*(3.*B1*ALFA23+B3-B1))+
     *      D3V2*(QB*B3+WB*(3.*B3*ALFA21+B1-B3))
        QB=B3*B3
        WB=W*B3
        ZII=ZII+Z(IND3)*(QB*(3.-2.*B3)+6.*WB*(B2*ALFA31+B1*ALFA32))+
     *      D2V3*(QB*B2+WB*(3.*B2*ALFA31+B1-B2))+
     *      D1V3*(QB*B1+WB*(3.*B1*ALFA32+B2-B1))
      RETURN
   30 ZII=Z(IND1)
      RETURN
   40 ZII=Z(IND2)
      RETURN
   50 ZII=Z(IND3)
      RETURN
      END