      SUBROUTINE  IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD)
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.0
        NMY=0.0
        NMZ=0.0
        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.0)    GO TO 22
            DZ2=ZD(IPI)-Z0
            DNMX=DY1*DZ2-DZ1*DY2
            DNMY=DZ1*DX2-DX1*DZ2
            IF(DNMZ.GE.0.0)    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.0
        NMXY=0.0
        NMYX=0.0
        NMYY=0.0
        NMZ =0.0
        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.0)    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.0)    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.0*NMZ)
        PD(JPD0)  =-NMYY/NMZ
   34 CONTINUE
      RETURN
      END
