      SUBROUTINE  BVIP2 (MD,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI,
     1                      IWK,WK,IERR)
C THIS SUBROUTINE PERFORMS SMOOTH SURFACE FITTING WHEN THE PRO-
C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
C DISTRIBUTED IN THE PLANE.
C THE INPUT PARAMETERS ARE
C     MD  = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
C         = 1 FOR NEW NCP AND/OR NEW XD-YD,
C         = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
C         = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
C     NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C           MATING PARTIAL DERIVATIVES AT EACH DATA POINT
C           (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
C     NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE X
C           COORDINATES OF THE DATA POINTS,
C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE Y
C           COORDINATES OF THE DATA POINTS,
C     ZD  = ARRAY OF DIMENSION NDP CONTAINING THE Z
C           COORDINATES OF THE DATA POINTS,
C     NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE
C           (MUST BE 1 OR GREATER),
C     NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE
C           (MUST BE 1 OR GREATER),
C     XI  = ARRAY OF DIMENSION NXI CONTAINING THE X
C           COORDINATES OF THE OUTPUT GRID POINTS,
C     YI  = ARRAY OF DIMENSION NYI CONTAINING THE Y
C           COORDINATES OF THE OUTPUT GRID POINTS.
C THE OUTPUT PARAMETERS ARE
C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (NXI,NYI),
C           WHERE THE INTERPOLATED Z VALUES AT THE OUTPUT
C           GRID POINTS ARE TO BE STORED.
C     IERR = ERROR INDICATOR. IERR IS SET TO 0 IF NO ERRORS
C           ARE DETECTED.
C THE OTHER PARAMETERS ARE
C     IWK = INTEGER ARRAY OF DIMENSION
C              MAX0(31,27+NCP)*NDP+NXI*NYI
C           USED INTERNALLY AS A WORK AREA,
C     WK  = ARRAY OF DIMENSION 5*NDP USED INTERNALLY AS A
C           WORK AREA.
C ERROR RETURN
C     IERR = 1      MD IS NOT 1, 2, OR 3.
C     IERR = 2      EITHER  2.LE.NCP.LT.NDP  OR  NCP.LE.NCPMX
C                   IS VIOLATED.
C     IERR = 3      NDP IS LESS THAN 4.
C     IERR = 4      NXI OR NYI IS LESS THAN 1.
C     IERR = 5      NCP OR NDP IS MODIFIED. THIS CANNOT BE
C                   DONE WHEN MD = 2 OR 3.
C     IERR = 6      NXI OR NYI IS MODIFIED. THIS CANNOT BE
C                   DONE WHEN MD = 3.
C     IERR = 7      POINTS (XD(I),YD(I)) AND (XD(J),YD(J))
C                   ARE EQUAL OR ARE TOO CLOSE WHERE
C                   IWK(1) = I AND IWK(2) = J.
C     IERR = 8      THE POINTS IN XD,YD,ZD ARE COLLINEAR OR
C                   ARE ALMOST COLLINEAR.
C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
C YD ARRAYS MUST BE MADE WITH MD=1.  THE CALL WITH MD=2 MUST BE
C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS.  THE CALL WITH
C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
C NXI, AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD, YD,
C XI, AND YI ARRAYS.  BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
C THIS SUBROUTINE CALLS THE IDCLDP, IDGRID, IDPDRV, IDPTIP, AND
C IDTANG SUBROUTINES.
C DECLARATION STATEMENTS
      DIMENSION   XD(NDP),YD(NDP),ZD(NDP),XI(NXI),YI(NYI),
     1            ZI(1),IWK(6),WK(1)
      COMMON/IDPI/ITPV,DMMY(27)
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
C (FOR MD=1,2,3)
   10 MD0=MD
      NCP0=NCP
      NDP0=NDP
      NXI0=NXI
      NYI0=NYI
C ERROR CHECK.  (FOR MD=1,2,3)
   20 IERR=0
      IF(MD0.LT.1.OR.MD0.GT.3)           GO TO 90
      IF(NCP0.LT.2.OR.NCP0.GE.NDP0)      GO TO 91
      IF(NDP0.LT.4)                      GO TO 92
      IF(NXI0.LT.1.OR.NYI0.LT.1)         GO TO 93
      IF(MD0.GE.2)        GO TO 21
      IWK(1)=NCP0
      IWK(2)=NDP0
      GO TO 22
   21 NCPPV=IWK(1)
      NDPPV=IWK(2)
      IF(NCP0.NE.NCPPV)   GO TO 94
      IF(NDP0.NE.NDPPV)   GO TO 94
   22 IF(MD0.GE.3)        GO TO 23
      IWK(3)=NXI0
      IWK(4)=NYI0
      GO TO 30
   23 NXIPV=IWK(3)
      NYIPV=IWK(4)
      IF(NXI0.NE.NXIPV)   GO TO 95
      IF(NYI0.NE.NYIPV)   GO TO 95
C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY.  (FOR MD=1,2,3)
   30 JWIPT=16
      JWIWL=6*NDP0+1
      JWNGP0=JWIWL-1
      JWIPL=24*NDP0+1
      JWIWP=30*NDP0+1
      JWIPC=27*NDP0+1
      JWIGP0=MAX0(31,27+NCP0)*NDP0
C TRIANGULATES THE X-Y PLANE.  (FOR MD=1)
   40 IF(MD0.GT.1)   GO TO 41
      CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
     1            IWK(JWIWL),IWK(JWIWP),WK,IERR)
      IF (IERR.NE.0) GO TO 96
      IWK(5)=NT
      IWK(6)=NL
      GO TO 50
   41 NT=IWK(5)
      NL=IWK(6)
C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT.  (FOR MD=1)
   50 IF(MD0.GT.1)   GO TO 60
      CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC),IERR)
      IF (IERR.NE.0) RETURN
C SORTS OUTPUT GRID POINTS IN ASCENDING ORDER OF THE TRIANGLE
C NUMBER AND THE BORDER LINE SEGMENT NUMBER.  (FOR MD=1,2)
   60 IF(MD0.EQ.3)   GO TO 70
      CALL IDGRID(XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),NXI0,NYI0,
     1            XI,YI,IWK(JWNGP0+1),IWK(JWIGP0+1))
C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
C (FOR MD=1,2,3)
   70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
C INTERPOLATES THE ZI VALUES.  (FOR MD=1,2,3)
   80 ITPV=0
      JIG0MX=0
      JIG1MN=NXI0*NYI0+1
      NNGP=NT+2*NL
      DO 89  JNGP=1,NNGP
        ITI=JNGP
        IF(JNGP.LE.NT)    GO TO 81
        IL1=(JNGP-NT+1)/2
        IL2=(JNGP-NT+2)/2
        IF(IL2.GT.NL)     IL2=1
        ITI=IL1*(NT+NL)+IL2
   81   JWNGP=JWNGP0+JNGP
        NGP0=IWK(JWNGP)
        IF(NGP0.EQ.0)     GO TO 86
        JIG0MN=JIG0MX+1
        JIG0MX=JIG0MX+NGP0
        DO 82  JIGP=JIG0MN,JIG0MX
          JWIGP=JWIGP0+JIGP
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI0+1
          IXI=IZI-NXI0*(IYI-1)
          CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
     1                ITI,XI(IXI),YI(IYI),ZI(IZI))
   82   CONTINUE
   86   JWNGP=JWNGP0+2*NNGP+1-JNGP
        NGP1=IWK(JWNGP)
        IF(NGP1.EQ.0)     GO TO 89
        JIG1MX=JIG1MN-1
        JIG1MN=JIG1MN-NGP1
        DO 87  JIGP=JIG1MN,JIG1MX
          JWIGP=JWIGP0+JIGP
          IZI=IWK(JWIGP)
          IYI=(IZI-1)/NXI0+1
          IXI=IZI-NXI0*(IYI-1)
          CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
     1                ITI,XI(IXI),YI(IYI),ZI(IZI))
   87   CONTINUE
   89 CONTINUE
      RETURN
C ERROR EXIT
   90 IERR=1
      RETURN
   91 IERR=2
      RETURN
   92 IERR=3
      RETURN
   93 IERR=4
      RETURN
   94 IERR=5
      RETURN
   95 IERR=6
      RETURN
   96 IF (IERR.NE.7) RETURN
      IWK(1)=IWK(JWIWP)
      IWK(2)=IWK(JWIWP+1)
      RETURN
      END
