C+
C NAME:
C	MapReadSingle
C PURPOSE:
C	Read files containing Carrington maps into a 2D array
C	All files containing data inside range [XCBeg, XCEnd] are averaged together
C	using a weighted mean.
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
	function MapReadSingle(cWild,BListFnc,BReadFnc,NCoff,XCBeg,XCEnd,nLng,nLat,Z,R0)
C INPUTS:
C	cWild		character	wildcard for locating synoptic map files
C					in format used by href=ForeignFile=.
C
C	BListFnc	integer		external function
C					name of integer function used by href=BListAll=
C	BReadFnc	integer		external function
C					name of integer function used to read individual files
C
C	NCoff		integer		integer Carrington variable offset
C
C	XCBeg		real		start Carrington variable
C	XCEnd		real		stop Carrington variable (usually XCbeg+3)
C
C	nLng		integer		# grid longitudes covering [XCbeg,XCend]
C	nLat		integer		# grid latitudes  covering [-90,90]
C OUTPUTS:
C	Z(nLng,nLat)	real		function values for synoptic map
C	R0		real		source surface distance in AU
C CALLS:
C	ArrR4Zero, ArrR4DivideByArrR4, ArrR4TimesConstant, Say, SplineGridX, SplineGridY
C	BadR4, BListAll, ArrR4Total
C EXPLICIT:
	integer		BListFnc
	integer		BReadFnc
C EXTERNAL:
	external	BListFnc
	external	BReadFnc
C PROCEDURE:
C >	The external integer function BListFnc is called as follows:
C	    I = BListFnc(cFile, NCoff, TFound,XCFoundBeg,XCFoundEnd,R0)
C
C	with input arguments
C		cFile		character	fully-qualified file name
C		NCoff		integer		a Carrington offset to be subtracted from
C						the output time and Carrington variables
C	and output arguments
C		I		integer		return status
C						0 : failure
C						1 : success
C		TFound		real		time assinged to the Carrington map
C		XCFoundBeg	real		start Carrington variable of the synoptic map
C		XCFoundEnd	real		end Carrington variable of the Carrington map
C						(always XCFoundEnd=XCFoundBeg+1.0
C		R0		real		source surface distance in AU
C
C >	The external integer function BReadFnc is called as follows
C >	This user-written function BReadFnc reads each seperate file into the scratch array ZFile
C
C	    I = BReadFnc(cFile,nMax,ZFile,nLng,nLat)
C
C	with input arguments
C		cFile		character	fully-qualified file name
C		nMax		integer		Max. # elements read from file
C	and output arguments
C		I		integer		return status
C						0 : failure
C						1 : success
C		ZFile(nLng,nLat)
C				real		2D array for synoptic map
C		nLng		integer		# grid longitudes in the data file
C		nLat		integer		# grid latitudes  in the data file
C >	nLng = 0 or nLat = 0 signals automatically sets the map 
C	dimensions, matching the resolution of the data files
C >	The output array ZFile covers the full heliographic longitude range [0,360],
C	i.e. ZFile(1,*) corresponds to 0 deg; ZFile(nLng,*) to 360 deg.
C >	The latitude range covers the full latitude range [-90,90],
C	i.e. ZFile(*,1) corresponds to -90 deg; ZFile(*,nLat) to +90 deg.
C >	The longitude grid must be evenly spaced in degrees
C >	The latitude grid must be evenly space in latitude or sin(latitude)
C
C >	MapReadSingle extracts data from ZFile in the range [XCbeg,XCend].
C >	The output array Z will be such that Z(1,*,*) corresponds to
C	Carrington variable XCend and Z(nLng,*,*) to XCbeg.
C >	Currently the following types of synoptic data can be read:
C	- Magnetic source surface map:
C	  Original WSO maps: href=BList_WSO=,  href=BRead_WSO=
C	  NOAAS WSO maps: href=BList_WSO_NOAA=, href=BRead_WSO_NOAA=
C MODIFICATION HISTORY:
C	JUL-1993, Paul Hick (UCSD/CASS)
C	MAY-2002, Paul Hick (UCSD/CASS)
C		Rewritten to support the tomography program. All plotting support
C		for GRPACK programs has been dropped.
C	AUG-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Bug fix. Changed order of two pairs of XCBeg,XCEnd variables in call to
C		SplineGridX (XCEnd needs to come first).
C-
	    character	cWild*(*)
	    integer	NCoff
	    real	XCBeg
	    real	XCEnd
	    integer	nLng
	    integer	nLat
	    real	Z(*)

	character	cSay*13		/'MapReadSingle'/


	parameter	(nLngFileMax =  73)	! Max. expected # longitudes in data files
	parameter	(nLatFileMax =  37)	! Max. expected # latitudes in data files
	parameter	(nLngOutMax  = 109)	! Max. expected # longitudes in output array (3x36+1) 

	parameter	(nBinFileMax = nLngFileMax*nLatFileMax)	! Max. expected # bins in data files
	parameter	(nBinGridMax = (3*(nLngFileMax-1)+1+1)*nLatFileMax)

	real		ZFile(nBinFileMax)
	real		ZGrid(nBinGridMax*2)	! Includes room for accumulating weights

	parameter	(nBinTempMax = nLngOutMax*nLatFileMax)
	real		ZTemp(nBinTempMax)


	parameter	(nFileMax = 100)	! Max. # data files to be read

	character	cFile    (nFileMax)*256
	real		TFile    (nFileMax)
	real		XCFileBeg(nFileMax)
	real		XCFileEnd(nFileMax)
	real		WFile	 (nFileMax)
	integer		Nr	 (nFileMax)

	character	cXCvarFormat*9
	integer		BListAll

	locMap(nLng,iLng,iLat) = (iLat-1)*nLng+iLng

	MapReadSingle = 0

	Bad = BadR4()


	!-------
	! nLngFile, nLatFile refer to the array as stored in the data files,
	! and correspond to a single rotation.

	nFile = BListAll(0,cWild,BListFnc,NCoff,XCBeg,XCEnd,nFileMax,
     &		cFile,TFile,XCFileBeg,XCFileEnd,WFile,Nr,R0,Scale)

	if (nFile .eq. 0) return		! If no file found, return with status=0

	!-------
	! Function values read by BReadFnc into array ZFile are accumulated in
	! the array ZGrid. The first half of ZGrid are used to add the function values.
	! The second half of Z are used to count the number of function values added
	! (in the averaging phase these will be used to calculate average function values).


	do iFile=1,nFile			! Loop over all candidate files

	    if (BReadFnc(cFile(iFile),nBinFileMax,ZFile,nLngFile,nLatFile) .ne. 1) return
						! On read error, return with status=0
	    call Say(cSay,'I','#'//cFile(iFile),'@['//
     &		cXCvarFormat(NCoff, XCFileBeg(iFile))//','//cXCvarFormat(NCoff, XCFileEnd(iFile))//']')

	    !-------
	    ! nLngFile and nLatFile are set in BReadFnc, and are the dimensions
	    ! of the synoptic map as stored in the data files. The first nLngFile*nLatFile
	    ! elements of ZFile contain the array read from file. nLngFile covers
	    ! [0,360] degrees in longitude, nLatFile covers [-90,90] degrees in latitude.

	    if (iFile .eq. 1) then		! Processing first file

		!-------
		! Accumulate array Z that covers range [XCGridBeg, XCGridEnd], a slightly
		! larger range then the input range [XCbeg,XCend]. Z is set up at the same
		! resolution as the arrays read from file.

		dXC = 1.0/(nLngFile-1)		! Resolution in data files

		!-------
		! Set XCGridBeg to longitudinal grid location just before XCbeg
		! Set XCGridEnd to longitudinal grid location just after  XCend

		XCGridBeg = int(XCBeg) + int( (XCBeg-int(XCBeg))/dXC)*dXC
		XCGridEnd = int(XCEnd) + int( (XCEnd-int(XCEnd))/dXC)*dXC
		if (XCGridEnd .lt. XCend) XCGridEnd = XCGridEnd+dXC

		!-------
		! # longitudinal and latitudinal gridpoints in output array Z.

		nLngGrid = nint( (XCGridEnd-XCGridBeg)/dXC ) + 1
		nLatGrid = nLatFile
		nBinGrid = nLngGrid*nLatGrid

		if (nBinGrid .gt. nBinGridMax) call Say(cSay,'E','Array','too small. Increase nBinGridMax')

		call ArrR4Zero(nBinGrid*2,ZGrid)	! Initialize output array

	    end if

	    !-------
	    ! ZGrid has longitude index running from 1 to nLngGrid covering
	    ! XCGridEnd to XCGridBeg.
	    ! Calulate longitude index iMap corresponding to XCFileBeg(iFile).
	    ! Longitude nLngFile (of array ZFile) is dropped into longitude iMap
	    ! of array ZGrid, if it is inside the range 1,nLngGrid

	    WW   = WFile(iFile)
	    iLng = nLngFile
	    iMap = nLngGrid - nint( (XCFileBeg(iFile) - XCGridBeg)/dXC )

	    !-------
	    ! Counting back from nLngFile, find the first longitude of ZFile
	    ! that falls inside ZFile (at iMap).
	
	    do while (iLng .ge. 1 .and. iMap .gt. nLngGrid) ! Find first longitude in ZZ
		iLng = iLng-1				! Index in ZFile
		iMap = iMap-1				! Index in ZGrid
	    end do

	    !-------
	    ! Copy longitude iLng from ZFile to iMap in ZGrid.

	    do while (iLng .ge. 1 .and. iMap .ge. 1)	! Add map iFile to output array Z

		!-------
		! Accumulate into column iMap after multiplying with proper weight WFile.
		! The weighted function values are accumulated in the first nBinGrid elements
		! of ZGrid; the weights are accumulated in the second nBinGrid elements.

		do iLat=1,nLatFile
		    nP = locMap(nLngFile,iLng,iLat)
		    iP = locMap(nLngGrid,iMap,iLat)
		    if (ZFile(nP) .ne. Bad) then
			ZGrid(iP) = ZGrid(iP) + WW*ZFile(nP)
			iP = iP+nBinGrid
			ZGrid(iP) = ZGrid(iP) + WW
		    end if
		end do

		iLng = iLng-1
		iMap = iMap-1

	    end do

	end do

	!--------
	! If the weights array is not empty then at least one map was processed,
	! and presumably we have a sensible result. A completely empty weights array
	! probably means that all files read contained bad values only.

	if (ArrR4Total(nBinGrid,ZGrid(nBinGrid+1),iP) .gt. 0) then

	    !-------
	    ! Calculate weighted average by dividing weighted function values
	    ! by weights.

	    call ArrR4DivideByArrR4(nBinGrid,ZGrid,ZGrid(nBinGrid+1),ZGrid)

	    !-------
	    ! If the input values of nLng and/or nLat were not set, then set them
	    ! to the resolution of the data files.

	    if (nLng .eq. 0) nLng = nLngGrid
	    if (nLat .eq. 0) nLat = nLatGrid

	    if (nLng*nLatGrid .gt. nBinTempMax) call Say(cSay,'E','nBinTempMax','parameter too small')

	    !-------
	    ! Interpolate to exact range [XCBeg,XCEnd], and change to final longitude resolution nLng.
	    ! ZGrid(nLngGrid,nLatGrid) --> ZTemp(nLng,nLatGrid)

	    call SplineGridX(1,XCGridEnd,XCGridBeg,XCEnd,XCBeg,nLngGrid,nLatGrid,ZGrid,nLng,ZTemp)

	    !-------
	    ! Change to final latitude resolution nLat.
	    ! ZTemp(nLng,nLatGrid) --> Z(nLng,nLat)

	    call SplineGridY(1,-90.0,+90.0,-90.0,+90.0, nLng,nLatGrid,ZTemp,nLat,Z)

	    !-------
	    ! Multiply by constant specified in BListFnc

	    call ArrR4TimesConstant(-nLng*nLat,Z,Scale,Z)

	    MapReadSingle = 1

	end if

	return
	end