C+
C NAME:
C	MapReadTimes
C PURPOSE:
C	Read files containing Carrington maps into a 2D array
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
	function MapReadTimes(cWild,BListFnc,BReadFnc,NCoff,nTim,TT,XCBeg,XCEnd,nLng,nLat,nRad,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	nTim		integer		# Carrington times
C	TT		real		Carrington times (with NCoff subtracted)
C
C	XCBeg		real		start Carrington variable (with NCoff subtracted)
C	XCEnd		real		stop  Carrington variable (with NCoff subtracted)
C
C	nLng		integer		# grid longitudes covering [XCbeg,XCend]
C	nLat		integer		# grid latitudes covering [-90,90]
C	nRad		integer		# grid radial distances
C OUTPUTS:
C	Z(nLng,nLat,nRad,nTim)
C			real		function values for synoptic map
C	R0		real		source surface distance in AU
C CALLS:
C	ArrR4Copy, ArrR4Bad, ArrR4GetMinMax, BListAll, Say, iArrR4ValuePresent
C	SplineY, SplineGridX, SplineGridY, ArrR4TimesConstant, BadR4
C INCLUDE:
	include		'filparts.h'
C EXPLICIT:
	integer		BListFnc
	integer		BReadFnc
C EXTERNAL:
	external	BListFnc
	external	BReadFnc
C RESTRICTIONS:
C	At least two synoptic maps need to available for the spline interpolation.
C	This could be a problem if only a single time TT is specified. If no files earlier
C	or later than TT is found than the function BListAll will fail to bracket the
C	time and return only a single file. In this case array Z will be filled with bad
C	values.
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
C
C >	The external integer function BReadFnc is called as follows
C	This user-written function 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		3D 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 >	The output array Z will be such that Z(1,*,*) corresponds to
C	Carrington variable XCend and Z(nLng,*,*) to XCbeg.
C	If XCEnd-XCBeg > 1 (covers more than one Carrington rotation) than heliographic
C	longitudes occur more than once in the output array. These will all be filled
C	with the same value.
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	SEP-2002, Paul Hick (UCSD/CASS)
C		The output array for times that are not bracketed by two files
C		are now flagged as bad (the spline interpolation produces crazy
C		values when it has to extrapolate).
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
	    integer	nTim
	    real	TT(*)
	    real	XCBeg
	    real	XCEnd
	    integer	nLng
	    integer	nLat
	    integer	nRad
	    real	Z(*)

	character	cSay*12		/'MapReadTimes'/

	parameter	(nFileMax = 100)	! Max. # data files to be read
	parameter	(nTimeMax = 200)	! Max. # times for interpolation

	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 

	parameter	(nBinFileMax = nLngFileMax*nLatFileMax)	! Max. expected # bins in data files

	real		ZFile(nBinFileMax*nFileMax)
	real		ZTime(nBinFileMax*nTimeMax)

	real		ZLng (nLngFileMax)

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

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

	real		dTExtra			/0.036664/		! 1.0 day as fraction of CR

	character	cXCvarFormat*9
	character	cStr*(FIL__LENGTH)
	integer		BListAll

	!-------
	! Statement functions

	locMap(iLng,iLat,iTim) = (  (iTim-1)       *nLatFile+iLat-1 )*nLngFile+iLng
	locOut(iLng,iLat,iTim) = ( ((iTim-1)*nRad )*nLat    +iLat-1 )*nLng    +iLng

	!-------

	MapReadTimes = 0

	Bad  = BadR4()
	nRad = max(nRad,1)

	!-------
	! Set the time range included all requested times TT. BListAll will try
	! to bracket this range, i.e. ideally there will be a file with
	! TFile <= TBeg and a file with TFile > TEnd.

	call ArrR4GetMinMax(nTim,TT,TBeg,TEnd)

	!-------
	! Get a list of all the files covering [TBeg,TEnd]. An attempt is made to
	! bracket the range with TFile(1) <= TBeg and TFile(nFile) >= TEnd.
	! There is however no guarantee that this is actually the case.

	nFile = BListAll(1,cWild,BListFnc,NCoff,TBeg,TEnd,nFileMax,
     &		cFile,TFile,XCFileBeg,XCFileEnd,WFile,Nr,R0,Scale)

	if (nFile .eq. 0) then
	    cStr = cWild
	    call Say(cSay,'W','#'//cStr,'no files found')
	    return				! If no file found, return with status=0
	end if

	if (nFile .eq. 1) then			! Would crash SplineY below
	    cStr = cWild
	    call Say(cSay,'W','#'//cStr,'only one file found; need at least two')
	    return				! If one file found, return with status=0
	end if

	!-------
	! Function values read by BReadFnc are accumulated in ZFile.
	! Read the first file. On read error, return with status=0

	if (BReadFnc(cFile(1),nBinFileMax,ZFile,nLngFile,nLatFile) .ne. 1) then
	    call Say(cSay,'W','#'//cFile(1),'read error on first file')
	    return
	end if

	!-------
	! The data files contain a nLngFile x nLatFile array
	! nLngFile covers [0.360] degrees
	! in longitude, nLatFile covers [-90,90] degrees in latitude.

	nBinFile = nLngFile*nLatFile		! # bins in single file
	dXC = 1.0/(nLngFile-1)			! Resolution in data files

	!-------
	! Check scratch space requirements:
	! Size ZFile

	if (nBinFile*nFile .gt. nBinFileMax*nFileMax) call Say(cSay,'E',
     &	    'Not enough','scratch space. Increase nBinFileMax or nFileMax')

	!-------
	! Size ZLng

	if (nLngFile .gt. nLngFileMax) call Say(cSay,'E',
     &	    'Not enough','scratch space. Increase nLngFileMax')

	!-------
	! Size ZTime

	if (nBinFile*nTim .gt. nBinFileMax*nTimeMax) call Say(cSay,'E',
     &	    'Not enough','scratch space. Increase nBinFileMax or nTimeMax')

	call Say(cSay,'I','#'//cFile(1),'@'//cXCvarFormat(NCoff, TFile(1)))

	!-------
	! The maps of nBinFile elements are stored one after the other in ZFile
	! On read error, return with status=0

	do iFile=2,nFile			! Read remaining files

	    if (BReadFnc(cFile(iFile),nBinFileMax,ZFile(locMap(1,1,iFile)),nLngFile,nLatFile) .ne. 1) return
	    call Say(cSay,'I','#'//cFile(iFile),'@'//cXCvarFormat(NCoff, TFile(iFile)))

	end do

	!------
	! The first nBinFile*nFile elements of ZFile now contain the maps for
	! a series of nFile times TFile. Rearrange the file so that they all cover the
	! same range of longitudes as the first file

	call ArrR4GetMinMax(-nBinFile*nFile,ZFile,xmin,xmax)

	do iFile=2,nFile

	    !------
	    ! iLng is the number of elements that each row must be shifted to smaller
	    ! longitudes (larger Carrington variables).
	    ! If XCFileBeg(iFile)-XCFileBeg(1) is not an exact multiple of dXC
	    ! the realignment is accurate only to within 0.5*dXC. This is fixed
	    ! below with a spline interpolation.

	    iLng = nint( mod(XCFileBeg(iFile)-XCFileBeg(1),1.0)/dXC )

	    if (iLng .lt. 0) iLng = nLng-1+iLng

	    if (iLng .lt. 0) print *, 'mapreadtimes, negative iLng=',iLng
	    if (iLng .eq. 0) print *, 'mapreadtimes, zero iLng=',iLng

	    do iLat=1,nLatFile

		iP = locMap(1,iLat,iFile)	! Start of row iLat of map iFile in array ZFile

		!------
		! Copy row iLat of length nLngFile into scratch array ZLng

		call ArrR4Copy(nLngFile       ,ZFile(iP)   , ZLng)

		!------
		! Move the last nLngFile-iLng row elements to the beginning of
		! the row. The last element of the row is discarded (it should
		! have the same value as the first element in the row.

		call ArrR4Copy(nLngFile-iLng-1,ZLng(iLng+1),ZFile(iP))

		!------
		! Move the first iLng+1 row elements to the end of the row
		! The first element is dropped at the spot of the discarded
		! element from the previous copy operation. Element iLng+1
		! ends up at the end of the row; the same element was also
		! put at the beginning of the row in the previous copy
		! operation.

		call ArrR4Copy(iLng+1         ,ZLng        ,ZFile(iP+nLngFile-iLng-1))

	    end do

	    !call SplineGridX(3,XCFileEnd(iFile),XCFileBeg(iFile),XCFileEnd(1),XCFileBeg(1),
     &		!	nLngFile,nLatFile,ZFile(locMap(1,1,iFile)),nLngFile,ZFile(locMap(1,1,iFile)))

	end do

	!-------
	! ZFile is an nLngFile x nLatFile x nFile array, representing maps at
	! times TFile(i), i=1,nFile. The next step is to create an array ZTime
	! which is an nLngFile x nLatFile x nTim array, representing maps at times
	! TT(i), i=1,nTim. If array TFile does not bracket TT then the spline
	! interpolation will put wacky values in the ZTime array for times TT far
	! outside the range covered by TFile. Below we will set ZTime to bad for
	! these times.

	call SplineY(0,nFile,TFile,nTim,TT,nBinFile,ZFile,ZTime)

	!-------
	! 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 = nint( (XCEnd-XCBeg)/dXC ) + 1
	if (nLat .eq. 0) nLat = nLatFile

	!-------
	! Check scratch space requirements:
	! Size ZTemp

	if (nLng*nLatFile .gt. nBinTempMax) call Say(cSay,'E','Not enough','scratch space. Increase nBinTempMax')

	!-------
	! The array ZTime(nLngFile,nLatFile,nTim) contains synoptic maps for the
	! correct times TT. The files cover 360 degrees, corresponding to the
	! Carrington range [XCFileBeg(1),XCFileEnd(1)]. The final step is to
	! change to the grid nLng x nLat x nTim.

	nBin = nLng*nLat

	!-------
	! Only process times TT inside the range covered by TFile.
	! For times outside this range the output array is flagged as bad.
	! TFile is already sorted so TBeg = TFile(1) and TEnd = TFile(nFile)

	call ArrR4GetMinMax(nFile,TFile,TBeg,TEnd)

	TMin = Bad
	Tmax = Bad

	do iTim=1,nTim

	    iP = locOut(1,1,iTim)		! First element in array Z for this time

	    if (TBeg-dTExtra .le. TT(iTim) .and. TT(iTim) .le. TEnd+dTExtra) then

		if (TMin .eq. Bad) TMin = TT(iTim)
		TMax = TT(iTim)

		!-------
		! Interpolate to exact range [XCBeg,XCEnd], and change to final
		! longitude resolution nLng. ZTime(nLngFile,nLatFile) -->
		! ZTemp(nLng,nLatFile). XCFileEnd(1)-XCFileBeg(1)=1. Longitudinal
		! locations in [XCBeg,XCEnd] outside [XCFileEnd(1),XCFileBeg(1)]
		! are mapped back into this range by what amounts to a mod(360)
		! operation.

		call SplineGridX(3,XCFileEnd(1),XCFileBeg(1),XCEnd,XCBeg,nLngFile,nLatFile,
     &			ZTime(locMap(1,1,iTim)),nLng,ZTemp)

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

		call SplineGridY(1,-90.0,90.0,-90.0,90.0, nLng,nLatFile,ZTemp,nLat,Z(iP))

		if (iArrR4ValuePresent(nBin,Z(iP),Bad) .eq. nBin) call Say(cSay,'W',
     &		    'Time',cXCvarFormat(NCoff,TT(iTim))//'contains only bad values')

	    else				! Outside TFile time range: flag bad

		call ArrR4Bad(nBin,Z(iP))
		!call Say(cSay, 'W','Time',cXCvarFormat(NCoff,TT(iTim))//'outside range of magnetic field maps')

	    end if

	end do

	call Say(cSay,'I','Bracketed','['//cXCvarFormat(NCoff,TMin)//','//cXCvarFormat(NCoff,TMax)//']')

	!------
	! Multiply by constant specified in BListFnc
	! (this converts magnetic fields to nano Tesla)

	call ArrR4TimesConstant(-nBin*nTim,Z,Scale,Z)

	MapReadTimes = 1

	return
	end