C+
C NAME:
C	FLINT8
C PURPOSE:
C	Linear interpolation in N dimensions.
C CALLING SEQUENCE:
	double precision function FLINT8(NN,ND,F,P,EPS)
C CALLS:
C	BadR8, ArrI4Copy, ArrR8Copy
C SEE ALSO:
C	FLINT
C PROCEDURE:
C	Double precision version of href=FLINT=.
C	Same syntax as FLINT, but uses double precision
C	instead of real*4 arguments.
C MODIFICATION HISTORY:
C	SEP-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		NN
	    integer		ND(*)
	    double precision	F(*)
	    double precision	P(*)
	    double precision	EPS

	logical			bCheckBad

	parameter		(N1MAX=10)		! Max # of degenerate dimensions
	integer			ND1(N1MAX)
	double precision	P1 (N1MAX)

	double precision	Bad
	double precision	BadR8
	double precision	FI
	double precision	FP
	double precision	FF

	NP(I) = max(1,min(ND(I)-1,int(P(I))))	! Index I of principal corner

	N = abs(NN)
	bCheckBad = NN .lt. 0

	Bad = BadR8()

	!-------
	! Don't extrapolate more than one unit. Make sure that
	! 0 <= P(I) <= ND(I)+1 or that the dimension I is degenerate
	! (if ND(I) = 1 then the coordinate P(I) is not used).

	I = 1			! Don't extrapolate more than one unit
	do while (I .le. N .and. ((0 .le. P(I) .and. P(I) .le. ND(I)+1) .or. ND(I) .eq. 1))
	    I = I+1
	end do

	if (I .le. N) then	! Invalid coordinate P(I)
	    FF = Bad

	else
	    N1 = 0
	    do I=1,N		! Check for degenerate dimensions (of 1 element)
		if (ND(I) .eq. 1) then
		    N1 = N1+1	! Count # degenerate dimensions
		    if (N1 .le. N1MAX) then
			ND1(N1) = I	! Position of degenerate dimension
			P1 (N1) = P(I)	! Degenerate coordinate
		    end if
		end if
	    end do


	    if (N1 .gt. N1MAX) then	! Too many degenerate dimensions
		FF = Bad

	    else if (N .eq. N1) then	! All dimensions degenerate
		FF = F(1)		! .. only 1 element in F

	    else

		!-----
		! Remove the degenerate dimensions from the ND and P arrays, by
		! shifting the non-degenerate dimensions to a lower index.

		do I=N1,1,-1
		    IP = ND1(I)
		    if (IP .ne. N) then
			call ArrI4Copy(N-IP,ND(IP+1),ND(IP))
			call ArrR8Copy(N-IP,P (IP+1),P (IP))
		    end if
		    ND(N) = 0	! Not really necessary
		    P (N) = 0.0d0
		end do

		N = N-N1	! # non-degenerate dimensions

		!-----
		! ND(1:N) now contains the non-degenerate dimensions and NP(1:N) the
		! the matching coordinates. ND(N+1:*) and P(N+1:*) will not be used.

		!-----
		! In the calling program F is declared as an n-dimensional array.
		! Internally F is addressed as a one dimensional array.
		!
		! If F is an n-dimensional array F(ND1,ND2,..,NDn) then element F(I1,I2,..,In)
		! is located at one-dimensional index IP:
		!
		! (1-dim array)  IP = 1+ I1-1
		! (2-dim array)  IP = 1+ I1-1 +ND1*( I2-1 )
		! (3-dim array)  IP = 1+ I1-1 +ND1*( I2-1 +ND2*(I3-1) )
		! (n-dim array)  IP = 1+ I1-1 +ND1*( I2-1 +ND2*( + ... ND[n-1]*(In-1))...)

		IP = NP(N)-1	! Calculate offset for the principal corner:
		do I=N-1,1,-1
		    IW = NP(I)-1
		    IP = IW+ND(I)*IP
		end do

		IP0 = IP+1	! One-dimensional index of principal corner

		!-----
		! The n-dim cube has 2^n corners, identified by ICN=0,2^n-1.
		! The array index of the principal corner in the i-th dimension: NP(I)
		! The index for the other corners in the i-th dimension  : NP(I)+0 or NP(I)+1
		! The offset 0 or 1 is calculated according to the following scheme
		! (example is for a 3-dim array).
		!
		!           0               1          x = 1st dim
		!         /   \           /   \
		!       0       1       0       1      y = 2nd dim
		!     /   \   /   \   /   \   /   \
		!     0   1   0   1   0   1   0   1    z = 3rd dim
		!
		!     0   1   2   3   4   5   6   7      = ICN
		!
		! E.g. corner ICN=5 has array indices NP(1)+1,NP(2)+0,NP(3)+1
		! The principal corner ICN=0 has indices NP(1),NP(2),NP(3)

		FF  = 0.0d0
		ICN = -1
		NCN = 2**N-1

		do while (ICN .lt. NCN)	! Loop over all ICN corners of the N-dim cube around P
		    ICN = ICN+1

		    IC = 2*ICN
		    FI = 1.0d0
		    IP = 0

		    do I=N,1,-1		! Loop over dimensions N-1 to 1
			IC = IC/2
			IW = mod(IC,2)	! Offset 0/1 for I-th dim
					! Distance to principal corner =
			DI = P(I)-NP(I)	! .. weight 1-DI for offset 0; DI for offset 1

			if (abs(DI  ) .lt. EPS) DI = 0.0d0
			if (abs(DI-1) .lt. EPS) DI = 1.0d0

			FI = FI*((1-IW)*(1-DI)+IW*DI)	! Total weight for corner ICN

			IP = IW+ND(I)*IP! Total 1-dim of corner ICN relative
					! .. to the principal corner ICN=0
		    end do

		    !------
		    ! If weight for corner ICN is zero, then the function value
		    ! does not matter and we can skip this corner. We check for
		    ! FI = 0.0 before we check for bad values.

		    if (FI .ne. 0.0d0) then
			FP = F(IP0+IP)	! If bad value, return with Bad
			if (bCheckBad .and. FP .eq. Bad) then
			    ICN = NCN	! Break do while loop
			    FF  = Bad	! Return bad value
			else
			    FF  = FF+FI*FP	! Accumulate interpolation
			end if
		    end if
		end do

		!-------
		! If there were degenerate dimensions we still need to restore
		! the arrays ND and P to their input values, using the arrays
		! ND1 and P1.

		N = N+N1		! Original # dimensions

		do I=1,N1		! Loop over degenerate dimensions
		    IP = ND1(I)
		    if (IP .ne. N) then	! Make room for insertion by shifting right
			call ArrI4Copy(N-IP,ND(IP),ND(IP+1))
			call ArrR8Copy(N-IP,P (IP),P (IP+1))
		    end if
		    ND(IP) = 1		! Restore input values
		    P (IP) = P1(I)
		end do

	    end if

	end if

	FLINT8 = FF			! Interpolated value (maybe bad)

	return
	end