NAME:
	GAMMALN
 PURPOSE:
	Calculate natural logarithm of gamma function for X>1
 CATEGORY:
	Math: special functions
 CALLING SEQUENCE:
	function GAMMALN(XX)
 INPUTS:
	XX	real*4	(read-only)	x-value where gamma fnc is evaluated
 OUTPUTS:
	GAMMALN real*4	logarithm of gamma function in X
 SIDE EFFECTS:
 RESTRICTIONS:
 PROCEDURE:
	Lanczos formulae. Calculations are performed in double precision.
 MODIFICATION HISTORY:
	Numerical Recipes, p. 157.

 NAME:
	GaussLookup
 PURPOSE:
	Calculates exp(-x*x-y*y) using lookup table
 CATEGORY:
	Math
 CALLING SEQUENCE:
	function GaussLookup(X,Y)
 INPUTS:
	X		real		X-coordinate
	Y		real		Y-coordinate
 OUTPUTS:
	F		real		fnc value
 CALLS: ***
	FLINT
 CALLED BY:
	GridSphere3D, GridSphereRange, GridSphereWeight, RebinSphere
 SEE ALSO:
	CosLookup, SinLookup
 MODIFICATION HISTORY:
	SEP-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	get_compiler
 CALLING SEQUENCE:
	subroutine get_compiler(f77id)
 CALLED BY:
	Time2System

 NAME:
	get_compiler
 CALLING SEQUENCE:
	subroutine get_compiler(f77id)
 CALLED BY:
	Time2System

 NAME:
	GinBar	
 PURPOSE:
	Duummy routine for Windows and Linux
	for functions containing calls to GRPACK functions
 PROCEDURE:
	Needed by: MapReadSrf
 MODIFICATION HISTORY:
	JAN-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	GridFill
 PURPOSE:
	Removes empty bins. An empty bin is given the average over all 
	neighbours with valid function values.
 CATEGORY:
	Plotting: contours
 CALLING SEQUENCE:
	subroutine GridFill(nFillIn,nX,nY,Z,NSIDE,Zmin,Zmax)
 INPUTS:
	nFill		integer	The three least-significant digits of nFill are used:
				mod(nFill,1000) >= 100: messages to screen are suppressed
				mod(nFill, 100) >=  10: Z(.,.,2) is used as output array
					(otherwise Z(*,*,2) is never accessed)
				mod(nFill,  10) >=   0: threshold on # valid neighbours
					If zero then all empty bins with extrapolated values
					are filled in. If larger than zero (and at most 8)
					only empty bins with valid neighbours equal or larger
					than the specified threshold are filled in.
	nX,nY		integer	dimensions of input grid
	Z(nX,nY) or Z(nX,nY,2)	!!! Z(nX,nY) can be used only if nFill < 10 !!!
			real	
		Z(.,.,1)	array of function values (BadR4() indicates bad value)
		Z(.,.,2)	used as scratch space; content does not matter
	NSIDE(nX,nY,-1:1,-1:1)
			integer*1 scratch array
 OUTPUTS:
	Z(nX,nY,2)	real
		Z(.,.,1)	array of function values with empty bins removed
				or (only nFill>0)
		Z(.,.,2)	array identifying the extrapolated values:
				= -1  contents of bin is same as valid input value
				=  1  contents of bin is extrapolated value
				nFill>0 only:
				=  0  empty bin with less than nFill neighbours
				      Z(.,.,1) = BadR4() (same as input).
	Zmin		real	minimum function value (incl interpolated fncv)
	Zmax		real	maximum function value (incl interpolated fncv)
 CALLS: ***
	ArrR4GetMinMax, ArrR4Mask, BadR4, GridFillWeight, Int2Str, Say, Str2Str
	iArrR4ValuePresent, iGridScan
 CALLED BY:
	GridSphere3D, GridSphereRange, GridSphereWeight, RebinSphere, smei_get_lsff
 RESTRICTIONS:
	If the input array Z is declared in the calling program as Z(nX,nY)
	then nFill MUST be less than 10 to avoid accessing Z(*,*,2)
 PROCEDURE:
 >	If nFill = 0 then:
	Step 1: for each empty bin, count the number of non-empty neighbours
	Step 2: find the subset of empty bins with the maximum number of
		non-empty neighbours
	Step 3: for the subset of step 2, calculate the average over the non-
		empty neighbours and assign this average to the empty bin
	Step 4: Go to step 1
	Repeat until there are no empty bins left.

 >	nFill > 0:
	Step 1: for each empty bin, count the number of non-empty neighbours
	Step 2: find the subset of empty bins with nFill or more non-empty
		neighbours
	Step 3: for the subset of step 2, calculate the average over the non-
		empty neighbours and assign this average to the empty bin
	Step 4: Return.

 >	Neighbouring function values are averaged using the function
	GridFillWeight to get weighting factors. The default function
	returns GridFillWeight = 1 (see end of this file).
	Override the default function with a different version by
	making sure it precedes this file in the compilatation.

	The function has the form
		function GridFillWeight(I,J,INX,JNY,ZN)
	I,J and INX,JNY are the indices of empty bin and neighbour bin.
	ZN the function value in the neighbour bin.

 >	'Bad' bins are usually identified by testing for BadR4().
	This can be modified using GridFillBad:

	rbad = BadR4Set(0.0)
	call GridFill(nFillIn,nX,nY,Z,NSIDE,Zmin,Zmax)
	rbad = BadR4Set(rbad)

	The second BadR4Set call restores the original BadR4 value
 MODIFICATION HISTORY:
	1990, Paul Hick (UCSD/CASS)
	APR-1999, Paul Hick (UCSD/CASS)
	    simplified code by introducing iGridScan function
	MAY-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    added option to suppress output messages
	NOV-2004, Paul Hick (UCSD/CASS)
	    Added entry point GridFillBad to modify the actual
	    value used to indicate 'bad'.
	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
	    Removed GridFillBad again. The BadR4 value can
	    now be set by BadR4Set.

 NAME:
	GridFillWeight
 PURPOSE:
	Sample for use in GridFill
 CALLING SEQUENCE:
	function GridFillWeight(I,J,INX,JNY,ZN)
 CALLED BY:
	GridFill
 MODIFICATION HISTORY:
	SEP-2007, Paul Hick (UCSD/CASS)
	    Extracted from gridfill.f

 NAME:
	GridRan2Reg
 PURPOSE:
	Generate function values in a regular grid of NX by NY points.
	The input function values are specified on a set of NPTS random
	points XP,YP
 CATEGORY:
	Plotting: contours
 CALLING SEQUENCE:
	subroutine GridRan2Reg(DistIn,NPTS,XP,YP,ZP,XB,YB,XE,YE,NX,NY,Z)
 INPUTS:
	DistIn		real		data points closer than abs(Dist) grid
					spacings from a grid point are included
					in the averaging
			If Dist < 0 then input fnc-values in arrays Z and ZP of
			BadR4() values are disregarded in the averaging

	NPTS		integer		number of points in the random set
			If NPTS < 0 then Z returns standard deviations
	XP(NPTS)	real		X/Y-coordinates of points in the
	YP(NPTS)			random set in user-specified units
	ZP(NPTS)	real		function values in the random set
	XB,YB		real		X/Y-coordinate of grid point (1,1)
					in user units
	XE,YE		real		X/Y-coordinate of grid point (NX,NY)
					in user units

	NX,NY		integer		dimensions of regular output grid
 OUTPUTS:
	Z(NX,NY)	real		grid function values.
					If no function value was calculated for
					a particular grid point the value
					BadR4() is returned
 CALLS: ***
	ArrI4Zero, ArrR4DivideByArrI4, ArrR4Zero, BadR4, Say
 SEE ALSO:
	GridReg2Reg
 RESTRICTIONS:
 >	The user units for XB,XE,YB,YE should be the same as for XP and YP
 >	The # elements in the output grid are limited to the
	value set for parameter NMAX. Currently NMAX=10000.
 PROCEDURE:
 >	BadR4() is used to identify invalid fnc-values in in- and output Z
 >	The output grid defines a regular grid of NYxNX squares.
	The function values Z are calculated by averaging over points ZP
	inside a grid square.
 MODIFICATION HISTORY:
	1989-1990, Paul Hick (SRON,MPAE,UCSD/CASS)
	1993, Paul Hick (UCSD/CASS; pphick@ucsd.edu); complete revision

 NAME:
	GridReg2Reg
 PURPOSE:
	Generate function values in a regular grid of NX by NY points. The
	input function values are specified on a regular grid of MX by MY points
 CATEGORY:
	Plotting: contours
 CALLING SEQUENCE:
	subroutine GridReg2Reg(Dist,MX,MY,ZP,NX,NY,Z)
 INPUTS:
	Dist		real		data points closer than abs(Dist) grid
					spacings from a grid point are included
					in the averaging
			If Dist < 0 then input fnc-values in arrays Z and ZP of
			BadR4() values are disregarded in the averaging

	MX,MY		integer		dimensions of regular input grid
	ZP(MX,MY)	real		function values in input grid

	NX,NY		integer		dimensions of regular output grid
 OUTPUTS:
	Z(NX,NY)	real		grid function values.
					If no function value was calculated for
					a particular grid point the value
					BadR4() is returned
 CALLS: ***
	BadR4
 SEE ALSO:
	GridRan2Reg
 RESTRICTIONS:
	The user units for XB,XE,YB,YE should be the same as for XP and YP
 PROCEDURE:
 >	BadR4() is used to identify invalid fnc-values in in- and output Z
 >	The output grid defines a regular grid of NYxNX squares.
	The function values Z are calculated by averaging over points ZP
	inside a grid square.
 MODIFICATION HISTORY:
	1989-1990, Paul Hick (SRON,MPAE,UCSD/CASS)
	1993, Paul Hick (UCSD/CASS; pphick@ucsd.edu); complete revision

 NAME:
	GridSphere2D
 PURPOSE:
	Fills 'holes' and/or smooths a 2D array using a Gaussian weighting function.
	The array dimensions represent a regular grid in spherical coordinates
	(longitude and latitude). The third is the grid in radial distance.
 CALLING SEQUENCE:
	subroutine GridSphere2D(dXC,nLng,nLat,nR,Z,WidDeg,iBadZ,WThreshold,ClipLng)
 INPUTS:
	    real	dXC		! range of Carrington variable covered by nLng
					! i.e. number of 360 deg rotations.
	    integer	nLng		! # longitudes
	    integer	nLat		! # latitudes
	    integer	nR		! # radial distances
	    real	Z(nLng,nLat,nR)	! 3D array of function values
	    real	WidDeg		! Width of Gaussian used for angular
					! smoothing (in degrees) (WidDeg<0 has
					! special meaning; see PROCEDURE)
	    integer	iBadZ		! 0: Invalid elements remain invalid;
					!    valid elements are replaced by a
					!    smoothed value
					! 1: All elements (valid and invalid) are
					!    replaced by smoothed values.
					! 2: Invalid elements are replaced by
					!    smoothed values; valid elements remain
					!    untouched.
					!
					! The following two options are useful if all
					! bad bins have to be filled in with something.
					! Use them at your own risk (see PROCEDURE).
					!
					! 3: Same as iBadZ=1, but if any invalid
					!    elements remain at the end, these are
					!    all filled in by a call to GridFill
					! 4: Same as iBadZ=2, but if any invalid
					!    elements remain at the end, these are
					!    all filled in by a call to GridFill
					!
					! >=10: use open grid for longitude and latitude
	    real	WThreshold	! The replacement by a smoothed value is
					! made only if total weight is 
					! larger than/equal to the threshold value.
					! (see PROCEDURE).
	    real	ClipLng		! Longitude difference (see PROCEDURE)
					! (usually set to zero).
 OUTPUTS:
	    real	Z(nLng,nLat,nR)	  smoothed array of function values
 CALLS: ***
	GridSphere3D
 CALLED BY:
	smei_skyd_fts, smei_skyd_init, smei_skyd_make, smei_skyd_sky
 PROCEDURE:
	GridSphere2D is a one-line wrapper around GridSphere3D
 MODIFICATION HISTORY:
	FEB-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	GridSphere3D
 PURPOSE:
	Fills 'holes' and/or smooths a 3D array using a Gaussian weighting function.
	Two of the array dimensions represent a regular grid in spherical
	coordinates (longitude and latitude). The third represents a regular
	grid in radial distance.
 CALLING SEQUENCE:
	subroutine GridSphere3D(dXC,nLng,nLat,nR,Z,WidDeg,R1,iBadZ,WThreshold,ClipLng)
 INPUTS:
	dXC		real		range of Carrington variable covered by nLng
					i.e. number of 360 deg rotations.
	nLng		integer		# longitudes
	nLat		integer		# latitudes
	nR		integer		# radial distances
	Z(nLng,nLat,nR)	real		3D array of function values
	WidDeg		real		Width of Gaussian used for angular
					smoothing (in degrees) (WidDeg<0 has
					special meaning; see PROCEDURE)
	R1		real		(3D version only)
				R1 > 0:Radial distance for inner boundary (in
					grid spacings), i.e. location of Z(*,*,1)
					R1=0 IS INVALID; MUST BE R1>0
				R1 < 0:(value does not matter)
					suppresses the radial smoothing, i.e.
					each subarray Z(*,*,i),i=1,nR is
					2D smoothed over longitude and latitude only.
	iBadZ		integer	     0: Invalid elements remain invalid;
					valid elements are replaced by a
					smoothed value
				     1: All elements (valid and invalid) are
					replaced by smoothed values.
				     2: Invalid elements are replaced by
					smoothed values; valid elements remain
					untouched.

					The following two options are useful if all
					bad bins have to be filled in with something.
					Use them at your own risk (see PROCEDURE).

				     3: Same as iBadZ=1, but if any invalid
					elements remain at the end, these are
					all filled in by a call to GridFill
				     4: Same as iBadZ=2, but if any invalid
					elements remain at the end, these are
					all filled in by a call to GridFill

				  >=10: Use open grid for longitude and latitude (?)
				 >=100: Use lookup table for Gaussian

	WThreshold	real		The replacement by a smoothed value is
					made only if total weight is 
					larger than/equal to the threshold value.
					(see PROCEDURE).
	ClipLng		real		Longitude difference (see PROCEDURE)
					(usually set to zero).
 OUTPUTS:
	Z(nLng,nLat,nR)	real		smoothed array of function values
 CALLED BY:
	GridSphere2D
 INCLUDE:
	include		'math.h'
 CALLS: ***
	ArrR4Copy, BadR4, GaussLookup, GridFill, Say, iArrR4ValuePresent
 RESTRICTIONS:
 >	GridSphere3D uses internal scratch arrays, which are
	dimensioned using 3 parameters, LNG,LAT and LNR. If the scratch
	arrays are too small, program execution stops.
 >	iBadZ=3 and iBadZ=4:
	GridFill uses a 2D algorithm to fill in holes. If nR>1 then GridFill
	is applied to each level k=1,nR separately.
 PROCEDURE:
 >	Invalid or 'Bad' array elements are indicated by the value BadR4()
 >	The 1st dimension represents a regular grid in longitude covering
	[0,360] degrees, i.e. Lng(i) = 360*(i-1)/(nLng-1), i=1,nLng.
	The 2nd dimension represents a regular grid in latitude covering
	[-90,+90] degrees, i.e. Lat(j) = -90*[-1+2*(j-1)/(nLat-1)], j=1,nLat.
	The 3rd dimension represents a regular grid in radial distance,
	i.e. R(k) = R1+k-1
 >	Obviously, 2D arrays can be smoothed by setting nR=1.
 >	GridSphere2D is just a special call to GridSphere3D with R1=-1 (the
	value does not matter as long as it's negative).
 >	The weighting is Gaussian with weight=1 at the origin. To avoid loosing
	valid elements for iBadZ=0 or 1, set the threshold weight to a
	value less than 1.
 >	The angular weighting factor involves the angular distance between
	two locations on a sphere. The points are connected by two arcs,
	one is the minimum angular distance A (0<A<180), the other is the
	complementary angle AC=360-A (180<AC<360).
	If ClipLng=0 the angle A is always used. As a result, the smoothed
	function value for small/large longitudes will depend on the input
	function values at large/small longitudes with fairly large weighting
	factors (since the angular distance is small).
	This effect can be suppressed by setting ClipLng to a non-zero value.
	If the longitude difference abs(Lng(i1)-Lng(i2)) is larger than
	ClipLng, then the corresponding element is excluded.
	BE CAREFUL WHEN USING THIS OPTION. In particular, setting
	ClipLng to values around 180 degrees affects the smoothing for
	points near the poles (where points with a wide range of longitude
	differences are close together). As a guideline 
	ClipLng > 180+3*WidDeg is probably not a bad place to start.
 >	If WidDeg > 0, the smoothing procedure is speeded up by looping
	over only part of the array when calculating the smoothed function 
	value. Excluded are all elements which are more than 3 times the
	Gaussian width away from the element for which the smoothed value is
	calculated (the weight would be less than exp(-9)). The exclusion
	is based on an approximate calculation of the distance between points.
	Near the poles there may be a problem (see the explanation in the code),
	but I think it's under control now.
	In case you don't trust the result set WidDeg to a negative value.
	If WidDeg <0, then the absolute value is used as Gaussian width, and
	the whole array is used in the calculation of the smoothed average.
	The procedure will be a LOT slower though.
 >	If ClipLng=0 then the averaging algorithm is mod(XCend-XCbeg), i.e. if
	in the input array column 1 (XC = XCend) and column nLng (XC = XCbeg)
	are the same this is also true for the output array.
	HOWEVER this symmetry is broken if iBadZ=3 or 4 because of the extra
	call to GridFill.
 MODIFICATION HISTORY:
	FEB-1997, Paul Hick (UCSD/CASS; pphick@ucsd.edu)

 NAME:
	GridSphereRange
 PURPOSE:
 CALLING SEQUENCE:
	entry GridSphereRange(range_lo_,range_hi_)
 CALLS: ***
	ArrR4Copy, BadR4, GaussLookup, GridFill, Say, iArrR4ValuePresent
 CALLED BY:
	smei_skyd_fts, smei_skyd_init, smei_skyd_make, smei_skyd_sky
 MODIFICATION HISTORY:
	MAY-2007, Paul Hick (pphick@ucsd.edu)

 NAME:
	GridSphereWeight
 PURPOSE:
 CALLING SEQUENCE:
	entry GridSphereWeight(nLng,nLat,nR,Z,ZThreshold_)
 CALLS: ***
	ArrR4Copy, BadR4, GaussLookup, GridFill, Say, iArrR4ValuePresent
 CALLED BY:
	smei_skyd_fts, smei_skyd_init, smei_skyd_make, smei_skyd_sky
 MODIFICATION HISTORY:
	MAY-2007, Paul Hick (pphick@ucsd.edu)