C+
C NAME:
C	nrPolInt2
C PURPOSE:
C	Approximate function value in point (X1,X2) by 2D polynomial
C	interpolation on a 2D array
C CATEGORY:
C	Math: 2D polynomial interpolation 
C CALLING SEQUENCE:
	subroutine nrPolInt2(M,N,X1A,X2A,YA,X1,X2,Y,DY)
C INPUTS:
C	M	integer		dimension in X1 direction
C	N	integer		dimension in X2 direction
C	X1A(M)	real		coordinate values in X1 direction
C	X2A(N)	real		coordinate values in X2 direction
C	YA(M,N)	real		2D array of function values
C	X1	real		X1 coordinate and ..
C	X2	real		.. X2 coordinate of point where interpolated
C				.. fnc-value is required
C OUTPUTS:
C	Y	real		interpolated fnc-value in (X1,X2)
C	DY	real		error estimate (based on X1 direction)
C CALLS:
C	nrPolInt
C RESTRICTIONS:
C	The dimensions M,N must be less then MMAX and NMAX respectively.
C	MMAX and NMAX are define as parameters (both equal 10)
C PROCEDURE:
C	See Numerical Recipes, par. 3.6, p. 97
C	First do M 1-dimensional interpolations in the x2 direction at
C	x1 value X1A(J) to find the fnc-values at X1A(J),X2. Then do one
C	final 1-dimensional interpolation in the x1 direction at X2 to find
C	the fnc-value at X1,X2
C MODIFICATION HISTORY:
C	JUN-1993, Paul Hick (UCSD)
C-
	    integer	M
	    integer	N
	    real	X1A(M)
	    real	X2A(N)
	    real	YA(M,N)
	    real	X1
	    real	X2
	    real	Y
	    real	DY

	parameter	(NMAX=10)
	parameter	(MMAX=10)
	real		YNTMP(NMAX)
	real		YMTMP(MMAX)

	if (M .gt. MMAX .or. N .gt. NMAX) call Say('nrPolInt2','E','Stop','M and N must be <= MMAX')

	do J=1,M
	    do K=1,N
		YNTMP(K) = YA(J,K)
	    end do
	    call nrPolInt(N,X2A,YNTMP,X2,YMTMP(J),DY)
	end do
	call nrPolInt(M,X1A,YMTMP,X1,Y,DY)

	return
	end