SUBROUTINE ITRP(P,NX,NY,BX,BY,BB)
C
C        MAY 1993   GLAHN, CHAMBERS   TDL   HP9000
C
C        PURPOSE
C            TO INTERPOLATE INTO FIELD P(NX,NY) AT POINT BX, BY, AND
C            RETURN VALUE IN BB.  BIQUADRATIC INTERPOLATION IS USED
C            WITHIN GRID EXCEPT BILINEAR IS USED IN OUTSIDE BORDER.
C            LINEAR EXTRAPOLATION IS USED FOR POINTS OUTSIDE THE GRID.
C
C        DATA SET USE
C           NONE.
C
C        VARIABLES
C            SEE PURPOSE.
C
C        NONSYSTEM SUBROUTINES CALLED
C            NONE.
C
      DIMENSION P(NX,NY),B(4)
C
      DATA KFILDO/12/
C
      NBX=BX
      NBY=BY
      NXM1=NX-1
      NYM1=NY-1
      IF(NBX-1)114,120,111
 111  IF(NBX-NXM1)112,120,115
 112  IF(NBY-1)121,130,113
 113  IF(NBY-NYM1)140,130,123
 114  NBX=1
      GO TO 120
 115  NBX=NXM1
 120  IF(NBY-1)121,130,122
 121  NBY=1
      GO TO 130
C
 122  IF(NBY-NY)130,123,123
 123  NBY=NYM1
C
C        START BI-LINEAR INTERPOLATION-EXTRAPOLATION.
C
 130  NBXP1=NBX+1
      NBYP1=NBY+1
      DX=BX-FLOAT(NBX)
      DY=BY-FLOAT(NBY)
C
      BB=P(NBX,NBY)+(P(NBXP1,NBY)-P(NBX,NBY))*DX+(P(NBX,NBYP1)-
     1P(NBX,NBY))*DY+(P(NBX,NBY)+P(NBXP1,NBYP1)-P(NBX,NBYP1)-P(NBXP1,
     2NBY))*DX*DY
      RETURN
C
C        START BI-QUADRATIC INTERPOLATION.
C
 140  DX=BX-FLOAT(NBX)
      DY=BY-FLOAT(NBY)
      NBYP2=NBY+2
      NBYP1=NBY+1
      NBYM1=NBY-1
      FCT=(DY**2-DY)/4.
      FET=(DX**2-DX)/4.
C
      DO 145 J=1,4
      N=NBX-2+J
      B(J)=P(N,NBY)+(P(N,NBYP1)-P(N,NBY))*DY+(P(N,NBYM1)+P(N,NBYP2)-
     1P(N,NBY)-P(N,NBYP1))*FCT
 145  CONTINUE
C
      BB=B(2)+(B(3)-B(2))*DX+(B(1)+B(4)-B(2)-B(3))*FET
      RETURN
      END