C+
C NAME:
C	MapReadt
C PURPOSE:
C	Read files containing a Carrington map as a 2D array. See PROCEDURE.
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call MapRead(MapReadFile,XCbeg,XCend,ICnr,XCadj,XCoff,
C    &			MMAX,mLng,mLat,Z,ZRead,iBlank,Zmin,Zmax)
C INPUTS:
C	XCbeg		real		start Carrington variable
C	XCend		real		stop Carrington variable (usually 
C					XCbeg+1)
C	ICnr		integer		# rotations to be averaged
C	XCadj(3)	real		only XCadj(1) is used as input.
C					If non-zero, a margin is added to 
C					either side of the map.
C					Only used by ReadSRF.
C	XCoff		real		offset to be subtracted from XCbeg and
C					XCend (will be 0 usually).
C					(if the array returned by
C					MapReadFile refers to OFF+[0,360],
C					set XCoff = 1-OFF/360.)
C	MMAX		integer		# elements in arrays Z and ZRead
C	mLng		integer		# grid longitudes, or 0
C	mLat		integer		# grid latitudes, or 0
C	Z(MMAX*2)	real		used as output array
C	ZRead(MMAX)	real		scratch array used by MapReadFile
C	iBlank		integer		= 0 : do not flag any fnc-values
C					= 1 : flag fnc-values BadR4()
C OUTPUTS:	(Z and iBlank are used as input for MAP_PLOT)
C	Z(mLng,mLat,2)	real		function value array for averaged
C					synoptic map
C	iBlank		integer		indicates blanking status
C	Zmin		real		minimun value in Z (BadR4() for no data)
C	Zmax		real		maximum value in Z (BadR4() for no data)
C	XCadj(3)	real		parameter defining extra map margin
C					(see subroutine CvPlot)
C CALLS:
C	nrSqueezeX, ArrR4Zero, ArrR4DivideByArrR4, ArrR4Mask
C	iArrR4ValuePresent, ArrR4GetMinMax, Say
C EXTERNAL:
C	MapReadFile	must be declared EXTERNAL in the calling program
C RESTRICTIONS:
C >	The external integer*4 function MapReadFile is called with the
C	following syntax
C
C	    I = MapReadFile(IC,MMAX,mLng,mLat,ZRead,XCadj)
C
C	with input arguments
C		IC	integer		Carrington rotation number
C		MMAX	integer
C		mLng	integer		# grid longitudes, or 0
C		mLat	integer		# grid latitudes, or 0
C	and output arguments
C		I	integer		return status
C					0 : failure
C					1 : success
C		ZRead(mLng,mLat)
C			real		2D array for synoptic map
C					(BadR4() values will not be used, if
C					iBlank=1 on input)
C >	mLng = 0 or mLat = 0 signals to the MapReadFile to set the map 
C	dimensions itself. Usually this option is used to switch back to a
C	full resolution map.
C >	The output array ZRead must cover the full heliographic longitude
C	range [0,360], i.e. ZRead(1,*) corresponds to 0 deg; ZRead(mLng,*) to
C	360 deg.
C >	If this is not the case for the array on the input file, some
C	interpolation scheme must be used to remap the data array read
C	from file (see e.g. the Fe XIV data)
C >	The latitude range must be the same for all files of a
C	particular data type (otherwise the averaging scheme in MapRead
C	gives wrong results), but can be less than the full [-90,90] range.
C >	The longitude grid must be evenly spaced in degrees
C >	The latitude grid must be evenly space in latitude or sin(latitude)
C PROCEDURE:
C >	MMAX refers to the maximum size of the map.
C >	MapRead creates a file Z ready for input in the MAP_PLOT synoptic map
C	plotting subroutine, in case the synoptic maps are available as
C	seperate files of 2D arrays.
C >	The user-written external function MapReadFile reads each seperate file
C	into the scratch array ZRead (see RESTRICTIONS)
C >	MapRead extracts data in the range [XCbeg,XCend]+i (i=0,ICnr-1)
C	from ZRead and accumulates averages over ICnr rotations in the Z array.
C >	The output array Z will be such that Z(1,*,*) corresponds to
C	Carrington variable XCend and Z(mLng,*,*) to XCbeg.
C <	If XCend=XCbeg+1, Z covers a full 360 deg in heliographic longitude.
C >	Currently the following types of synoptic data can be read:
C	- Stanford source surface maps	(MapReadFile = MapReadSRF)
C	- Greg's Yohkoh synoptic maps	(MapReadFile = MapReadYOH)
C	- Sac Peak Fe XIV maps		(MapReadFile = MapReadSAC)
C	- Mark III maps			(MapReadFile = MapReadMRK)
C >	If iBlank=0 on input then iBlank=0 on output. If iBlank=1 on input,
C	it will be 1 on output, UNLESS all function values are valid (in which
C	case iBlank=0 is returned, i.e. blanking is switched off).
C MODIFICATION HISTORY:
C	JUL-1993, Paul Hick (UCSD)
C-
	subroutine MapReadt(MapReadFile,XCbegIn,XCendIn,ICnrIn,XCadj,XCoff,
     &			MMAX,mLngfull,mLat,Z,Zread,iBlank,Zmin,Zmax)

	real		Z(MMAX*2),
     &			ZRead(MMAX),
     &			XCadj(3)
     	integer			iDay		/100./
	

	Bad   = BadR4()

	mLng  = mLngfull

	ICnr  = max(1,ICnrIn)		! # rotations to be averaged
	XCbeg = XCbegIn-XCoff		! Optional offset
	XCend = XCendIn-XCoff
C-------
C If ICnr > 1 the last rotation for which data are required is ICfinal
C Data are accumulated for IN=1,ICnr: step IN requires data from rotations
C IN-1+[ICfrst,IClast]

	ICfrst = XCbeg			! First Carrington rotation
	IClast = XCend			! Last Carrington rotation
	if (IClast .eq. XCend) IClast = IClast-1
	ICfinal = IClast+ICnr-1

	XCadj1 = XCadj(1)

	call ArrR4Zero(2*MMAX,Z)	! Initialize output array
C-------
C Function values read by MapReadFile into array ZRead are accumulated in
C the array Z. The first half of Z (first mLngfull*mLat elements) are used to
C add the function values. The second half of Z (second mLngfull*mLat elements)
C are used to count the number of function values added (in the averaging phase
C these will be used to calculate average function values).

	MTOT = 0
	do IC=ICfrst,ICfinal		! Loop over Carrington rotations

C-------
C Note that the returned value of array XCadj is set by the last call to the
C MapReadFile.

	    if (XCadj1 .ne. 0.) XCadj(1) = XCadj1
	    	    
c	    if (MapReadFilet(IC,iDay,MMAX,mLng,mLat,ZRead,XCadj) .eq. 1) then
	    if (MapReadSRFt(IC,iDay,MMAX,mLng,mLat,ZRead,XCadj) .eq. 1) then


					! File read succesfully
c*tjd* note that mLng and mLat are reset by MapReadFile					
					
C-------
C Since mLng and mLat may be set in MapReadFile the calculation of Ibeg
C and Iend can NOT be moved in front of the DO-loop.
C The array ZRead returned by MapReadFile has dimension mLngxmLat

C-------
C Start longitude index: Ibeg; End longitude index: Iend.
C Note that the Ibeg and Iend values involve integer roundoff and do not
C correspond exactly to XCbeg and XCend respectively. Hence the array Z
C accumulated here does not correspond exactly to the required range
C [XCbeg,XCend]. This is fixed later by a spline interpolation.


c		print *, 'in mapread_t, checking for bad values in B'
c		print *, 'ZRead(1),ZRead(10), ZRead(30), ZRead(55) ZRead(MMAX)'
c		print *, ZRead(1), ZRead(10), ZRead(30), ZRead(55), ZRead(MMAX)
		
		if (MTOT .eq. 0) then

		    mLng1 = mLng-1
		    Ibeg = nint(1+(1-(XCbeg-ICfrst))*mLng1)
		    Iend = nint(1+(1-(XCend-IClast))*mLng1)


c*tjd* - the following checks stopped perfectly good Iend's

c		    if (Iend .eq. mLng) stop 'MapRead: darn, should not have happened'
c		    if (Iend .eq. mLng) Iend = 1	! Probably redundant


C--------
C The output array Z maintains the same resolution as the mLngxmLat ZRead 
C array. The range [XCbeg,XCend] translates to a Z array of mLngfullxmLat.


		    mLngfull = 1+Ibeg-Iend+mLng1*(IClast-ICfrst)
		    MTOT = mLngfull*mLat
		    if (MTOT .gt. MMAX) call Say('MapRead','E','Output',
     &			'array for map too small. Increase MMAX')

		end if

C-------
C For rotation IC parts of ZRead are copied into at most ICnr sections of Z:
C	maximum range of IN=1,ICnr.
C The actual range of IN=INfrst,INlast is determined as follows:
C For IC=ICfrst only the last IN value used is INlast=1 (by definition).
C	For subsequent rotations this value is incremented.
C For IC=ICfinal only the first IN value used is INfirst=ICnr (by definition).
C	For subsequent rotations this value is decremented.

		INfrst = max(1   ,ICnr+(IC-ICfinal))
		INlast = min(ICnr,1   +(IC-ICfrst ))

		do IN=INfrst,INlast	! Loop to accumulate rotations to be averaged

C-------
C Each column in the ZRead array is identified by IC (rotation number) and
C I=1,mLng (identifies the longitude in [0,360]). The columns in ZRead are
C copied to a column IP in the output array Z. The position IP in Z depends
C on IC, I and IN (identifies a step in the accumulation of maps for averaging).
C The mapping from (IC,I,IN) to IP is given by:
C
C	IP(IC,I,IN) = I+1-Iend+(IClast+IN-1-IC)*(mLng-1)
C		    = IP0(IC,IN)+I
C
C	IP0(IC,IN)  = 1-Iend+(IClast+IN-1-IC)*(mLng-1)
C
C If IP is inside the range [1,mLngfull] it is copied to column IP in Z. If
C it is outside this range IP<1 or IP>mLngfull, it is not copied.
C
C (the sketch below depicts the mapping (IC,I,IN=1) to IP
C
C	rot. IClast		     rot. IC			rot. ICfrst
C
C  IClast+1	  XCend      IC+1		  IC		XCbeg	  ICfrst
C   |		   |	       |		   |		|	     |
C --+-------------------+--..--+-------------------+--..--+------------------+-
C   |		   |	|      |		   |	  |	|	     |
C I= 1		  Iend	mLng    1		  mLng   1     Ibeg	   mLng
C		   |						|
C		   |						|
C		   ----------------------------------------------
C IP=		   1					     mLngfull
C
C If longitude 0 (I=1) [same as longitude 360 (I=mLng)] is included
C in the range [XCbeg,XCend] (but not on the edge, as would happen for integer
C values of XCbeg and XCend), then the corresponding column in Z implicitly
C gets a special treatment:
C * The first column (I=1, longitude 0) of rotation IC is copied into
C	IP=IP0(IC,IN)+1
C * The last column (I=mLng, longitude 360) of rotation IC+1 is copied into
C	IP=IP0(IC+1,IN)+mLng=IP0(IC,IN)-(mLng-1)+mLng=IP0(IC,IN)+1
C i.e. both are copied (added) into the same column.
C
C In most cases the first column for IC is the same as the last column for 
C IC+1, so as a result the column is added twice to the accumulated Z array.
C In the averaging phase this is automatically corrected since the counter
C for the corresponding column in Z will also be twice as big.

		    IP0 = 1-Iend+(IClast+IN-1-IC)*mLng1

C-------
C If IP(mLng)=IP0+mLng >= 1 and IP(1)=IP0+1 <= mLngfull, then at least some of
C the columns in the ZRead array will have to be copied to Z.
C (The test for IP0+1 .le. mLngfull is redundant, I think. Retained it to
C maintain symmetry.)

		    if (.not. (IP0+1 .le. mLngfull .and. IP0+mLng .ge. 1) )
     &			stop 'MapRead: oops'

c		    if (IP0+1 .le. mLngfull .and. IP0+mLng .ge. 1) then

C-------
C Column I=mLngfull-IP0 of ZRead maps to column IP(I)=IP0+I=mLngfull.

			I = min(mLng,mLngfull-IP0)
			do while (I .ge. 1 .and. IP0+I .gt. mLngfull)
			    I = I-1
			end do

			do while (I .ge. 1 .and. IP0+I .ge. 1)
c	print *, I,' into ', IP0+I
			    do J=1,mLat
				JP = (J-1)*mLng+I
				if (iBlank .eq. 0 .or. ZRead(JP) .ne. Bad) then
				    JQ = (J-1)*mLngfull+IP0+I
				    Z(JQ) = Z(JQ)+ZRead(JP)
				    JQ = JQ+MTOT
				    Z(JQ) = Z(JQ)+1
				end if
			    end do
			    I = I-1
			end do
c		    end if

		end do

	    end if

	end do

	if (MTOT .eq. 0) then
	    Zmin = Bad
	    Zmax = Bad
	    call Say('MapRead','W','No files','for this set of rotations')
	    return
	end if
c	print *, 'crash check#3'

	call ArrR4DivideByArrR4(MTOT,Z,Z(MTOT+1),Z)! Average over ICnr rotations
C-------
C Spline interpolation to exact output grid. Since Ibeg and Iend were calculated
C as nearest integers, the differences XCbegTmp-XCbeg and XCendTmp-XCend should
C never exceed half the grid distance 0.5/(mLng-1) in absolute value.
C (not sure about the following piece of code yet)

	call ArrR4GetMinMax(-MTOT,Z,Zmin,Zmax)
	if (Zmax .ne. Bad) then
	    XCbegTmp = ICfrst+1.*(mLng-Ibeg)/mLng1
	    XCendTmp = IClast+1.*(mLng-Iend)/mLng1

	    if (XCend .ne. XCendTmp .or. XCbeg .ne. XCbegTmp) then
		call Say('MapRead','I','Tune','finetuning output grid')
		call nrSqueezeX(XCendTmp,XCbegTmp,XCend,XCbeg,mLngfull,mLat,Z,.TRUE.)
	    end if
	end if

	call ArrR4Mask(MTOT,Z,Bad,0.,0.,-1.,0,Z(MTOT+1))! -1=Valid, 0=Invalid
C-------
C  If iBlank=0 on input then iBlank=0 on output (since all Z(I,J,2) values
C  will be unequal 0). If iBlank=1 on input, it will be 1 on output, UNLESS all
C  function values are valid (all Z(I,J,2) not equal 0).

	iBlank = iArrR4ValuePresent(MTOT,Z,Bad)
	iBlank = min(iBlank,1)

	call ArrR4GetMinMax(-MTOT,Z,Zmin,Zmax)
	if (Zmax .eq. Bad) then
	    call Say('MapRead','W','No data','for this set of rotations')
	    return
	end if

c	dAng = 10.
c	iFillBadZ = 0
c	WThreshold = 2.
c	ClipLng = 0.
c	call AskI4('smooth',ismo)
c	if (ismo .eq. 1) call GridSphere2D(MLngfull,MLat,1,Z,dAng,iFillBadZ,WThreshold,ClipLng)
c		print *, 'exitting mapread'
c		print *, ''

	return
	end