C+
C NAME:
C	GridRan2Reg
C PURPOSE:
C	Generate function values in a regular grid of NX by NY points.
C	The input function values are specified on a set of NPTS random
C	points XP,YP
C CATEGORY:
C	Plotting: contours
C CALLING SEQUENCE:
	subroutine GridRan2Reg(DistIn,NPTS,XP,YP,ZP,XB,YB,XE,YE,NX,NY,Z)
C INPUTS:
C	DistIn		real		data points closer than abs(Dist) grid
C					spacings from a grid point are included
C					in the averaging
C			If Dist < 0 then input fnc-values in arrays Z and ZP of
C			BadR4() values are disregarded in the averaging
C
C	NPTS		integer		number of points in the random set
C			If NPTS < 0 then Z returns standard deviations
C	XP(NPTS)	real		X/Y-coordinates of points in the
C	YP(NPTS)			random set in user-specified units
C	ZP(NPTS)	real		function values in the random set
C	XB,YB		real		X/Y-coordinate of grid point (1,1)
C					in user units
C	XE,YE		real		X/Y-coordinate of grid point (NX,NY)
C					in user units
C
C	NX,NY		integer		dimensions of regular output grid
C OUTPUTS:
C	Z(NX,NY)	real		grid function values.
C					If no function value was calculated for
C					a particular grid point the value
C					BadR4() is returned
C CALLS:
C	Say, BadR4, ArrR4Zero, ArrI4Zero, ArrR4DivideByArrI4
C SEE ALSO:
C	GridReg2Reg
C RESTRICTIONS:
C >	The user units for XB,XE,YB,YE should be the same as for XP and YP
C >	The # elements in the output grid are limited to the
C	value set for parameter NMAX. Currently NMAX=10000.
C PROCEDURE:
C >	BadR4() is used to identify invalid fnc-values in in- and output Z
C >	The output grid defines a regular grid of NYxNX squares.
C	The function values Z are calculated by averaging over points ZP
C	inside a grid square.
C MODIFICATION HISTORY:
C	1989-1990, Paul Hick (SRON,MPAE,UCSD/CASS)
C	1993, Paul Hick (UCSD/CASS; pphick@ucsd.edu); complete revision
C-
	    real	DistIn
	    integer	NPTS
	    real	XP(*)
	    real	YP(*)
	    real	ZP(*)
	    real	XB
	    real	YB
	    real	XE
	    real	YE
	    integer	NX
	    integer	NY
	    real	Z (*)

	parameter	(NMAX = 10000)
	integer		N(NMAX)
	real		S(NMAX)
	logical		bNoFlag
	logical		bSigma

	indx(I,J) = (J-1)*NX+I		! Used to index arrays

	NXY = NX*NY
	if (NXY .gt. NMAX) call Say('GridRan2Reg','E','Parameter','NMAX too small')

	bNoFlag = DistIn .ge. 0
	bSigma  = NPTS .lt. 0

	rBad = BadR4()
	Dist = abs(DistIn)

	DX = (XE-XB)/(NX-1)		! Grid spacings in user units
	DY = (YE-YB)/(NY-1)
	X1 = XB-DX
	Y1 = YB-DY
	XD = abs(DX*Dist)
	YD = abs(DY*Dist)
	IDIS = .5+Dist

	call ArrR4Zero(NXY,Z)		! Clear arrays associated with new grid
	call ArrI4Zero(NXY,N)

	do L=1,abs(NPTS)
	    if (bNoFlag .or. ZP(L) .ne. rBad) then
		I = nint(1+(XP(L)-XB)/DX)
		I = max(1,min(I,NX))
		J = nint(1+(YP(L)-YB)/DY)
		J = max(1,min(J,NY))	! Identify grid point closest to point L
		IX1 = -min(IDIS,I-1)
		IX2 =  min(IDIS,NX-I)
		JY1 = -min(IDIS,J-1)
		JY2 =  min(IDIS,NY-J)
		do J1 = J+JY1,J+JY2
		    Y = Y1+J1*DY
		    do I1 = I+IX1,I+IX2
			X = X1+I1*DX
			if (abs(XP(L)-X) .lt. XD .and. abs(YP(L)-Y) .lt. YD) then
			    nn = indx(I1,J1)
			    Z(nn) = Z(nn)+ZP(L)
			    S(nn) = S(nn)+ZP(L)*ZP(L)
			    N(nn) = N(nn)+1
			end if
		    end do
		end do
	    end if
  	end do

	if (bSigma) then
	    do nn=1,NXY
		S(nn) = rBad
		if (N(nn) .gt. 1) Z(nn) = sqrt(max(0.,(S(nn)-N(nn)*Z(nn)*Z(nn))/(N(nn)-1)))
	    end do
	else
	    call ArrR4DivideByArrI4(NXY,Z,N,Z)
	end if

	return
	end