C+
C NAME:
C	MapGrid
C PURPOSE:
C	Calculates the function values on an equidistant grid covering
C	the full heliographic longitude/latitude range of a synoptic map
C	from a set of 'random' observations scattered across the map.
C CATEGORY:
C	Plotting
C CALLING SEQUENCE:
	subroutine MapGrid(XCbeg,XCend,YLbeg,YLend,NC,
     &		 NL,XC,YL,ZL,NAV,DAV,NX,NY,Z,ZD,ZW,Zmin,Zmax,ZDmin,ZDmax)
C INPUTS:
C	XCbeg	real		Carrington variable at Z(NX,*)
C	XCend	real		Carrington variable at Z(1,*)
C				(XCrange=XCend-XCbeg must be less/equal 1)
C				(XCrange=1 corresponds to 360 deg in longitude)
C	YLbeg	real		heliographic latitude at Z(*,1) (degrees)
C	YLend	real		heliographic latitude at Z(*,NY) (degrees)
C				(YLbeg and YLend must be in [-90,90])
C	NC	integer		# Carrington rotations to be averaged
C	NL	integer		# observations
C				If NL is set to a negative value then values
C				BadR4() in array ZL are ignored.
C	XC(NL)	real		values of 'modified Carrington' variable
C	YL(NL)  real		heliographic latitude (degrees)
C	ZL(NL)  real		function values
C	NAV	integer		determines how the averaging is done
C				= 0 : averaging over square bins (fast!)
C				= 1 : use the weight function SphereWeight
C					(slow, but more sound method)
C	DAV	real		half-width of weighting function in units of 
C				the horizontal grid spacing (used in
C				calculating the grid function values Z)
C	NX,NY	integer		X,Y dimensions of equidistant map grid
C OUTPUTS:
C	Z (NX,NY)	real	function values in grid points
C				BadR4() if no function value available
C	ZD(NX,NY)	real	standard deviations
C	ZW(NX,NY)	real	weighting factors (scratch array)
C	Zmin		real	minimum value in Z
C	Zmax		real	maximum value in Z
C	ZDmin		real	minimum value in ZD
C	ZDmax		real	maximum value in ZD
C CALLS:
C	BadR4, SphereWeight, ArrR4GetMinMax, ArrR4Zero
C RESTRICTIONS:
C >	If XCend .lt. XCbeg their values are interchanged
C >	XCrange = XCend-XCbeg must be less/equal 1
C >	YLbeg and YLend must be in the range [-90,90]
C >	NC must be larger/equal 1
C >	If any of these restrictions is violated the input values of XCend,
C	YLbeg,YLend and NC will be changed
C PROCEDURE:
C >	The 'modified Carrington variable' contains information about the 
C	Carrington rotation to which the observations belongs (integer part)
C	and the heliographic longitude [(1.-decimal fraction)*360.]
C >	I.e. XCbeg lies in Carrington rotation N0 = int(XCbeg) and
C	corresponds to heliographic longitude XLbeg = (1-(XCbeg-N0))*360
C	degrees
C >	The first and last grid point in the X-direction [i.e. Z(1,*) and 
C	Z(NX,*)] have heliographic longitudes XCend*360. and XCbeg*360.
C	(grid spacing XCrange*360./(NX-1) degrees).
C >	The first and last grid point in the Y-direction [i.e. Z(*,1) and 
C	Z(*,NY) correspond to latitude YLbeg and YLend
C	(grid spacing (YLend-YLbeg)/(NY-1))
C >	I.e. XCrange=1,YLbeg=-90,YLend=90,NX=37,NY=19 would result
C	in an equidistant grid with grid spacing of 10 degrees in X and Y
C	direction covering both hemispheres.
C >	XC values in [XCbeg,XCend] are used to construct the synoptic
C	map. Points outside this range are used only for obtaining the
C	function values of grid points near the edge of the range.
C >	If the grid should run from 0 [Z(1,*)] to 360 degrees [Z(NX,*)]
C	XCbeg must be an integer.
C >	The weight function SphereWeight is a decreasing function of angular
C	distance from a given grid point (i.e. a function of the two 
C	independent variables in any spherical coordinate system defined on the
C	unit sphere). The grid function values are calculated as weighted 
C	averages over all data points (with weight SphereWeight). The method
C	is more sound than the averaging over square bins, but is much more
C	computationally intensive.
C > 	The code for SphereWeight is appended to this subroutine
C MODIFICATION HISTORY:
C	JAN-1992, Paul Hick (UCSD)
C	SEP-2001, Paul Hick (UCSD), added option to skip bad values in the
C		input array ZL.
C-
	    real	XCbeg
	    real	XCend
	    real	YLbeg
	    real	YLend
	    integer	NC
	    integer	NL

	    real	XC(*)
	    real	YL(*)
	    real	ZL(*)

	    integer	NAV
	    real	DAV

	    integer	NX
	    integer	NY

	    real	Z(NX,NY)
	    real	ZD(NX,NY)
	    real	ZW(NX,NY)

	    real	Zmin
	    real	Zmax
	    real	ZDmin
	    real	ZDmax

	if (NL .eq. 0) return		! No data points available

	Bad = BadR4()

	!-------
	!  Make sure the latitude range is valid. If not, change input values

	XCrange = XCend-XCbeg

	YLbeg = max(min(YLbeg,90.),-90.)
	YLend = max(min(YLend,90.),-90.)
	YLrange = YLend-YLbeg

	NC = max(1,NC)			! At least one Carrington rotation

	!-------
	! Range in Carrington variable; grid spacings

	dXC = -XCrange/(NX-1)		! Longitudinal grid spacing (1 unit=360 deg)
	dXL = -360.*dXC			! Longitudinal grid spacing (degrees)

	dYL = YLrange/(NY-1)

	IL = NX*NY
	call ArrR4Zero(IL,Z )		! Initialize arrays
	call ArrR4Zero(IL,ZD)
	call ArrR4Zero(IL,ZW)

	if (NAV .ne. 0) then

	!-------
	! AVERAGING IN SPHERICAL COORDINATES (USING ANGULAR DISTANCE IN FUNCTION 
	!	SphereWeight
	! Loop over all grid points IX=1,..,NX; IY=1,..,NY.
	! For each grid point loop over all data points IL=1,NL; include only points
	! with longitude less then 180 degrees away from the grid point (< 0.5 
	! difference in Carrington variable).

!$omp parallel do private(IX,IY,IL,IC,XCIX,YLIY,XCDIF,W)
!$omp&  lastprivate(WHALF)
	    do IY=1,NY
		YLIY = YLbeg+(IY-1)*dYL		! Heliographic latitude in grid

		do IX=1,NX
		    XCIX = XCbeg-(NX-IX)*dXC	! Carrington variable in grid

		    !-------
		    ! Calculate weighted sums Z and ZW. The weight is calculated by the 
		    ! function SphereWeight

		    do IL=1,abs(NL)		! Loop over all data points
			if (NL .gt. 0 .or. ZL(IL) .ne. Bad) then
			    do IC=0,NC-1	! Loop over Carrington rotations

				XCDIF = XC(IL)-IC-XCIX
				if (abs(XCDIF) .lt. .5) then! < 180 deg difference
				    W = SphereWeight(XCDIF*360.,YLIY,YL(IL),DAV*dXL,NAV,3.,WHALF)
				    Z (IX,IY) = Z (IX,IY)+W*ZL(IL)
				    ZD(IX,IY) = ZD(IX,IY)+W*ZL(IL)*ZL(IL)
				    ZW(IX,IY) = ZW(IX,IY)+W
				end if

			    end do
			end if
		    end do

		end do
	    end do
!$omp end parallel do

	else

	    WHALF = 1.

	    !-------
	    ! The averaging is done over square bins in longitude and latitude centered
	    ! on the grid point with width 2*DAV

	    XDIS = abs(DAV*dXC)
	    YDIS = abs(DAV*dYL)
	    IDIS = .5+DAV

	    !-------
	    ! IX0,IY0 are the array indices of the grid point closest to data point IL
	    ! For the range of neighbours of IX0,IY0 given by IX1,...,IX2 and IY1,...,IY2
	    ! we check whether IL is to be included in the weighted average for these grid points.

!$omp parallel private(IX,IY,IL,IC,IY0,IY1,IY2,IX0,IX1,IX2,XCIX,YLIY)
	    do IL=1,abs(NL)
		if (NL .gt. 0 .or. ZL(IL) .ne. Bad) then
		    IY0 = nint(1+(YL(IL)-YLbeg)/dYL)
		    IY1 = max(IY0-IDIS, 1)
		    IY2 = min(IY0+IDIS,NY)

		    do IC=0,NC-1
			IX0 = nint(NX+(XC(IL)-IC-XCbeg)/dXC)
			IX1 = max(IX0-IDIS, 1)
			IX2 = min(IX0+IDIS,NX)

!$omp do
			do IY=IY1,IY2
			    YLIY = YLbeg+(IY-1)*dYL	! Heliographic latitude in grid

			    do IX=IX1,IX2
				XCIX = XCbeg-(NX-IX)*dXC! Carrington variable in grid

				if (abs(XC(IL)-IC-XCIX) .lt. XDIS .and. abs(YL(IL)-YLIY) .lt. YDIS) then
c		print *, il, ix,xcix,xdis,xc(il),iy,yliy,ydis,yl(il)
				    Z (IX,IY) = Z (IX,IY)+ZL(IL)
				    ZD(IX,IY) = ZD(IX,IY)+ZL(IL)*ZL(IL)
				    ZW(IX,IY) = ZW(IX,IY)+1
				end if

			    end do

			end do
!$omp end do

		    end do
		end if
	    end do
!$omp end parallel

	end if

	!-------
	! The grid function is a weighted average or is flagged (if the total weight
	! ZW(IX,IY) is too small)

!$omp parallel do private(IX,IY)
	do IY=1,NY
	    do IX=1,NX
		if (ZW(IX,IY) .lt. WHALF) then
		    Z (IX,IY) = Bad
		    ZD(IX,IY) = Bad
		else
 		    Z (IX,IY) = Z(IX,IY)/ZW(IX,IY)
		    ZD(IX,IY) = ZD(IX,IY)/ZW(IX,IY)-Z(IX,IY)*Z(IX,IY)
		    ZD(IX,IY) = sqrt( max(ZD(IX,IY),0.) )
		end if
	    end do
	end do
!$omp end parallel do

	call ArrR4GetMinMax(-NX*NY,Z ,Zmin ,Zmax)
	call ArrR4GetMinMax(-NX*NY,ZD,ZDmin,ZDmax)

	return	
	end