C+
C NAME:
C	GridFill
C PURPOSE:
C	Removes empty bins. An empty bin is given the average over all 
C	neighbours with valid function values.
C CATEGORY:
C	Plotting: contours
C CALLING SEQUENCE:
	subroutine GridFill(nFillIn,nX,nY,Z,NSIDE,Zmin,Zmax)
C INPUTS:
C	nFill		integer	The three least-significant digits of nFill are used:
C				mod(nFill,1000) >= 100: messages to screen are suppressed
C				mod(nFill, 100) >=  10: Z(.,.,2) is used as output array
C					(otherwise Z(*,*,2) is never accessed)
C				mod(nFill,  10) >=   0: threshold on # valid neighbours
C					If zero then all empty bins with extrapolated values
C					are filled in. If larger than zero (and at most 8)
C					only empty bins with valid neighbours equal or larger
C					than the specified threshold are filled in.
C	nX,nY		integer	dimensions of input grid
C	Z(nX,nY) or Z(nX,nY,2)	!!! Z(nX,nY) can be used only if nFill < 10 !!!
C			real	
C		Z(.,.,1)	array of function values (BadR4() indicates bad value)
C		Z(.,.,2)	used as scratch space; content does not matter
C	NSIDE(nX,nY,-1:1,-1:1)
C			integer*1 scratch array
C OUTPUTS:
C	Z(nX,nY,2)	real
C		Z(.,.,1)	array of function values with empty bins removed
C				or (only nFill>0)
C		Z(.,.,2)	array identifying the extrapolated values:
C				= -1  contents of bin is same as valid input value
C				=  1  contents of bin is extrapolated value
C				nFill>0 only:
C				=  0  empty bin with less than nFill neighbours
C				      Z(.,.,1) = BadR4() (same as input).
C	Zmin		real	minimum function value (incl interpolated fncv)
C	Zmax		real	maximum function value (incl interpolated fncv)
C CALLS:
C	ArrR4Mask, iArrR4ValuePresent, ArrR4GetMinMax
C	iGridScan, Int2Str, Str2Str, Say, GridFillWeight
C RESTRICTIONS:
C	If the input array Z is declared in the calling program as Z(nX,nY)
C	then nFill MUST be less than 10 to avoid accessing Z(*,*,2)
C PROCEDURE:
C >	If nFill = 0 then:
C	Step 1: for each empty bin, count the number of non-empty neighbours
C	Step 2: find the subset of empty bins with the maximum number of
C		non-empty neighbours
C	Step 3: for the subset of step 2, calculate the average over the non-
C		empty neighbours and assign this average to the empty bin
C	Step 4: Go to step 1
C	Repeat until there are no empty bins left.
C
C >	nFill > 0:
C	Step 1: for each empty bin, count the number of non-empty neighbours
C	Step 2: find the subset of empty bins with nFill or more non-empty
C		neighbours
C	Step 3: for the subset of step 2, calculate the average over the non-
C		empty neighbours and assign this average to the empty bin
C	Step 4: Return.
C
C >	Neighbouring function values are averaged using the function
C	GridFillWeight to get weighting factors. The default function
C	returns GridFillWeight = 1 (see end of this file).
C	Override the default function with a different version by
C	making sure it precedes this file in the compilatation.
C
C	The function has the form
C		function GridFillWeight(I,J,INX,JNY,ZN)
C	I,J and INX,JNY are the indices of empty bin and neighbour bin.
C	ZN the function value in the neighbour bin.
C
C >	'Bad' bins are usually identified by testing for BadR4().
C	This can be modified using GridFillBad:
C
C	rbad = BadR4Set(0.0)
C	call GridFill(nFillIn,nX,nY,Z,NSIDE,Zmin,Zmax)
C	rbad = BadR4Set(rbad)
C
C	The second BadR4Set call restores the original BadR4 value
C MODIFICATION HISTORY:
C	1990, Paul Hick (UCSD/CASS)
C	APR-1999, Paul Hick (UCSD/CASS)
C	    simplified code by introducing iGridScan function
C	MAY-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    added option to suppress output messages
C	NOV-2004, Paul Hick (UCSD/CASS)
C	    Added entry point GridFillBad to modify the actual
C	    value used to indicate 'bad'.
C	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Removed GridFillBad again. The BadR4 value can
C	    now be set by BadR4Set.
C-
	    integer	nFillIn
	    integer	nX
	    integer	nY
	    real	Z(nX,nY,2)
	    integer*1	NSIDE(nX,nY,-1:1,-1:1)
	    real	Zmin
	    real	Zmax

	character	cSay*8		/'GridFill'/
	logical		bFill
	logical		b2
	logical		bLoud

	character	cStr*80
	integer		Str2Str

	InZ(I) = (1+sign(1,I))/2		! Used to keep array indices for
						! ... Z and NSIDE valid
	rBad = BadR4()

	nXY  = nX*nY
	nBad = iArrR4ValuePresent(nXY,Z,rBad)	! # bad values
	if (nBad .eq. nXY) call Say(cSay,'E','Empty','grid; no valid function values')

	nFill = nFillIn

	nFill = mod(nFill,1000)

	bLoud = nFill .lt. 100			! If .FALSE. then suppress messages
	nFill = mod(nFill,100)

	b2 = nFill .ge. 10			! If .TRUE. then Z(*,*,2) is accessed
	nFill = mod(nFill,10)			! Threshold for # valid neighbours

	!-------
	! Set up Z(*,*,2): -1 for valid elements of Z(*,*,1); 0 for invalid elements.

	if (b2) call ArrR4Mask(nXY,Z,rBad,0.0,0.0,-1.0,0.0,Z(1,1,2)) ! -1=Valid; 0=Invalid

	!-------
	! At this point Z(I,J,2) has value -1 or 0. The bins with value 0 are to
	! be replaced by averaged values.

	nX1  = nX-1
	nY1  = nY-1
	iBad = 0

	!-------
	! If nMax=-1 or nMax=0 then no bins need filling:
	! nMax=-1: no empty bins left;  nMax=0: all bins are empty

	nMax = iGridScan(nX,nY,Z,rBad,NSIDE)
	bFill = nMax .gt. 0

	do while (bFill)

	    do J=1,nY
		JY1 = -InZ(J  -2)
		JY2 =  InZ(nY1-J)

		do I=1,nX
		    IX1 = -InZ(I  -2)
		    IX2 =  InZ(nX1-I)

		    if (Z(I,J,1) .eq. rBad) then

			if (( nFill .eq. 0 .and. NSIDE(I,J,0,0) .eq. nMax ) .or.
     &			    ( nFill .gt. 0 .and. NSIDE(I,J,0,0) .ge. nFill)) then

			    ZZ = 0.0
			    WW = 0.0
			    do JY=JY1,JY2
				do IX=IX1,IX2		! Loop over filled neighbour bins
				    if (IX .ne. 0 .or. JY .ne. 0) then
					INX = I+IX	! Neighbour indices
					JNY = J+JY
							! Replacing this test by ZN .ne. rBad
							! DOES change the result !!!
					if (NSIDE(I,J,IX,JY) .eq. 1) then
					    ZN = Z(INX,JNY,1)
					    W  = GridFillWeight(I,J,INX,JNY,ZN)
					    ZZ = ZZ+W*ZN
					    WW = WW+W
					end if
				    end if
				end do
			    end do

			    Z(I,J,1) = ZZ/WW		! Average over non-empty neighbours
			    if (b2) Z(I,J,2) = 1.0	! Bin with averaged intensity
			    iBad = iBad+1		! Count # bins filled
			end if
		    end if

		end do

	    end do

	    bFill = nFill .eq. 0			! Continue only if iterating and if ..
	    if (bFill) then
		nMax  = iGridScan(nX,nY,Z,rBad,NSIDE)
		bFill = nMax .gt. 0			! .. there are still empty bins left
	    end if

	end do

	call ArrR4GetMinMax(-nXY,Z,Zmin,Zmax)		! Get minimum and maximum

	if (bLoud) then
	    if (nBad .ne. 0 .and. nFill .ne. 0) then
	    	I = 0
		I = I+Int2Str(iBad,cStr(I+1:))
		I = I+Str2Str('/' ,cStr(I+1:))
		I = I+Int2Str(nBad,cStr(I+1:))+1
		I = I+Str2Str('invalid function values',cStr(I+1:))
		call Say(cSay,'I','filled',cStr)
	    end if
	end if

	return
	end