C+
C NAME:
C	nrBCUINT
C PURPOSE:
C	Find function value of point inside grid cell by bicubic interpolation
C CATEGORY:
C	Math: function approximation
C CALLING SEQUENCE:
	subroutine nrBCUINT(Y,Y1,Y2,Y12,X1L,X1U,X2L,X2U,X1,X2,ANSY,ANSY1,ANSY2)
C INPUTS:
C	Given are four gridpoints of a rectangular grid cell, numbered
C	counterclockwise from the lower left)
C	Y(4)	real		function values in grid points
C	Y1(4)	real		partial derivatives in x1 direction
C	Y2(4)	real		partial derivatives in x2 direction
C	Y12(4)	real		cross-derivatives
C	X1L	real		lower corner in x1 direction
C	X1U	real		upper cornner in x1 direction
C	X2L	real		lower corner in x2 direction
C	X2U	real		upper corner in x2 direction
C	X1	real		x1 coordinate and ..
C	X2	real		x2 coordinate of point where fnc-value is required
C OUTPUTS:
C	ANSY	real		interpolated fnc-value in X1,X2
C	ANSY1	real		partial derivative in x1 direction
C	ANSY2	real		partial derivative in x2 direction
C CALLS:
C	nrBCUCOF, BadR4
C PROCEDURE:
C	See Numerical Recipes, par. 3.6, p. 99
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    real	Y(4)
	    real	Y1(4)
	    real	Y2(4)
	    real	Y12(4)
	    real	X1L
	    real	X1U
	    real	X2L
	    real	X2U
	    real	X1
	    real	X2
	    real	ANSY
	    real	ANSY1
	    real	ANSY2

	real		C(4,4)

	ANSY  = BadR4()
	ANSY1 = ANSY
	ANSY2 = ANSY

	D1 = X1U-X1L
	D2 = X2U-X2L
	if (D1*D2 .eq. 0.) return	! Zero-area

	call nrBCUCOF(Y,Y1,Y2,Y12,D1,D2,C)

	T = (X1-X1L)/D1
	U = (X2-X2L)/D2
	ANSY  = 0.
	ANSY1 = 0.
	ANSY2 = 0.
	do I=4,1,-1
	    ANSY  = T*ANSY+((C(I,4)*U+C(I,3))*U+C(I,2))*U+C(I,1)
	    ANSY1 = T*ANSY2+(3.*C(I,4)*U+2.*C(I,3))*U+C(I,2)
	    ANSY2 = U*ANSY1+(3.*C(4,I)*T+2.*C(3,I))*T+C(2,I)
	end do

	ANSY1 = ANSY1/D1
	ANSY2 = ANSY2/D2

	return
	end