C+
C NAME:
C	ReadGIPS
C PURPOSE:
C	Read collection of data from Cambridge daily IPS g-level files
C CALLING SEQUENCE:
	subroutine ReadGIPS(cDat,nCar,JDCar,
     &		iEdt,bAuto,bForeCast,
     &		XCtst1,XCtst2,MJDref,MJDfrst,MJDlast,
     &		Radius,Speed,Power,
     &		XElow,XEhigh,iEorW,XRlim,
     &		NLmax,NL,
     &		iXP,iYP,iMJD,iYRS,DOYS,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG,
     &		NSmax,iXPsav,iYPsav,iMJDsav,iYRSsav,DOYSsav,XDSsav,XLSsav,
     &		XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
C INPUTS:
C	Arguments passes unmodified to ReadG:
C
C	cDat		character*(*)	name of logical for dir containing data files
C	nCar		integer		(read-only)
C					# start times for Carrington rotations
C	JDCar(nCar)	real*8		(read-only)
C					Carrington start times (Julian days)
C	iEdt		integer		(read-only)
C					0 = Read unedited data file
C					1 = Read edited data file
C	Radius		real		Reference distance in AU
C					(not used if Speed=0)
C	Speed		real		trace back velocity (km/s)
C					(if Speed=0 then no traceback is done)
C	Power		real		if non-zero used to scale the measured
C					G-value: G = G(measured)/RP^GPow
C					where RP is the heliocentric distance 
C					of the point P.
C					(i.e. if GPow=0 the unscaled value is 
C					returned)
C
C	Arguments determining selection of data:
C
C	XCtst1,XCtst2
C	MJDref		real*8
C	MJDfrst,MJDlast	integer		modified Julian days for earliest, and latest data
C					files available
C
C	XElow,XEhigh	real		minimum and maximum elongation (degrees)
C	XEhigh		real		Only lines of sight with elongation between XElow and XEhigh
C					are used
C	iEorW		integer		= +1: use only lines of sight east of the Sun
C					= -1: use only lines of sight west of the Sun
C					=  0: use data east and west of the Sun
C	XRlim		real		Only points with Right Ascenscion less than XRlim away 
C					from the Sun are used
C
C
C	Arguments controlling amount of data to be read:
C
C	NLmax		integer		max # data points read into output arrays
C	NSmax		integer		max # data points saved in save arrays
C	NS		integer		# points in save arrays to be re-used
C					(usually this will be set by a previous call;
C					the save arrays can be discarded by setting NS=0;
C					See PROCEDURE)
C OUTPUTS:
C	NL		integer		# data points in output arrays
C	NS		integer		# data points save in save arrays
C	iXP (NL)	integer		location in original daily Cambridge file (hour angle)
C	iYP (NL)	integer		location in original daily Cambridge file (declination)
C	iMJD(NL)	integer		MJD for original daily Cambridge file
C	IYRS(NL)	integer		time: year
C	DOYS(NL)	real		time: doy of year
C	XDS (NL)	real		Earth-Sun distance (AU)
C	XLS (NL)	real		Geocentric ecliptic longitude Sun (degrees)
C	XLL (NL) 	real		Difference between geocentric ecliptic longitudes of
C					line of sight and Sun
C	XDL (NL) 	real		Geocentric ecliptic latitude of line of sight
C	XCE (NL) 	real		Modified Carrington variable of Earth
C	XE  (NL)	real		Elongation of the line of sight used
C					If XE<0	the l.o.s. points west of the Sun
C					If XE>0 the l.o.s. points east of the Sun
C
C	XC  (NL)	real		Modified Carrington variable of 'Point P', except for
C					an offset of NCoff, i.e. the conventional Carrington
C					numbers are NCoff+XC(I)
C	YL  (NL)	real		Heliographic latitudes (degrees) of 'Point P'
C	XV  (NL)	real		velocities (currently set to zero)
C	XG  (NL)	real		G-factors
C
C OUTPUT SAVE ARRAYS:
C	(See PROCEDURE; do not meddle with these arrays !!!!!!!)
C	IYRSsav(NS)	integer		time: year
C	DOYSsav(NS)	real		time: doy of year
C	XDSsav(NL)	real		Earth-Sun distance (AU)
C	XLSsav(NL)	real		Geocentric ecliptic longitude Sun (degrees)
C	XLLsav(NL) 	real		Difference between geocentric ecliptic longitudes of
C					line of sight and Sun
C	XDLsav(NL) 	real		Geocentric ecliptic latitude of line of sight
C	XCEsav(NL) 	real		Modified Carrington variable of Earth
C	XEsav(NL)	real		Elongation of the line of sight used
C					If XE<0	the l.o.s. points west of the Sun
C					If XE>0 the l.o.s. points east of the Sun
C
C	XCsav(NL)	real		Modified Carrington variable of 'Point P', except for
C					an offset of NCoff, i.e. the conventional Carrington
C					numbers are NCoff+XC(I)
C	YLsav(NL)	real		Heliographic latitudes (degrees) of 'Point P'
C	XVsav(NL)	real		velocities (currently set to zero)
C	XGsav(NL)	real		G-factors
C
C CALLS:
C	BadR4, ReadG, Say, itrim, AskWhat
C PROCEDURE:
C >	The save arrays can be used to speed up processing if a new set of data is to be read
C	for nearly the same input specifications as the previous call: most of the required
C	data will be extracted from the save arrays, rather than reading them from file.
C >	If the save arrays are to be used set NSmax ~ NLmax.
C >	If the save arrays are not be used at all, set NS=0. In that case the save arrays
C	are not accessed at all, and some memory can be saved by setting NSmax=1, and declaring
C	the save arrays are scalars in the calling program.
C MODIFICATION HISTORY:
C	JUN-1992, Paul Hick (UCSD)
C	MAR-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		merged synoptic map version with version for the deconvolution program.
C-
	    character	cDat*(*)
	    integer	nCar
	    real*8	JDCar(nCar)
	    integer	iEdt
	    logical	bAuto
	    logical	bForeCast
	    real	XCtst1
	    real	XCtst2
	    real*8	MJDref
	    integer	MJDfrst
	    integer	MJDlast
	    real	Radius
	    real	Speed
	    real	Power
	    real	XElow
	    real	XEhigh
	    integer	iEorW
	    real	XRlim
	    integer	NLmax
	    integer	NL

	    !-------
	    ! Output and save arrays
	    ! If not in auto mode (bAuto=.FALSE.) then NS=0. Then these arrays are
	    ! not accessed. I.e. if the save option is not needed set bAuto=.FALSE.
	    ! and NSmax=1 to save memory.

	    integer	iXP (NLmax)
	    integer	iYP (NLmax)
	    integer	iMJD(NLmax)
	    integer	iYRS(NLmax)
 
	    real	DOYS(NLmax)
	    real	XDS (NLmax)
	    real	XLS (NLmax)
	    real	XLL (NLmax)
	    real	XDL (NLmax)
	    real	XCE (NLmax)
	    real	XE  (NLmax)
	    real	XC  (NLmax)
	    real	YL  (NLmax)
	    real	VV  (NLmax)
	    real	GG  (NLmax)

	    integer	NSmax

	    integer	iXPsav (NSmax)		! Array indices from daily maps
	    integer	iYPsav (NSmax)
	    integer	iMJDsav(NSmax)		! Original MJD of daily maps
	    integer	iYRSsav(NSmax)		! Time: year

	    real	DOYSsav(NSmax)		! Time: doy of year
	    real	XDSsav (NSmax)		! Sun-Earth distance
	    real	XLSsav (NSmax)		! Ecliptic longitude Sun
	    real	XLLsav (NSmax)
	    real	XDLsav (NSmax)
	    real	XCEsav (NSmax)
	    real	XEsav  (NSmax)
	    real	XCsav  (NSmax)		! Carrington variable of Point P
	    real	YLsav  (NSmax)		! Heliographic latitude of Point P
	    real	VVsav  (NSmax)		! Solar wind velocity (km/s)
	    real	GGsav  (NSmax)		! G-factor

	!-------
	! Arrays passed to ReadG to read the Cambridge data files.
	! All, except XR_t, are copied into output arrays. XR_t is used internally only for
	! selection of data points.
 
	parameter	(NTIME	= 72)
	parameter	(NDEC	= 14)

	integer		iYRS_t(NTIME,NDEC)
	real		DOYS_t(NTIME,NDEC)
	real		XDS_t (NTIME,NDEC)
	real		XLS_t (NTIME,NDEC)
	real		XLL_t (NTIME,NDEC)
	real		XDL_t (NTIME,NDEC)
	real		XCE_t (NTIME,NDEC)
	real		XE_t  (NTIME,NDEC)
	real		XC_t  (NTIME,NDEC)
	real		YL_t  (NTIME,NDEC)
	real		XV_t  (NTIME,NDEC)
	real		XG_t  (NTIME,NDEC)
	real		XR_t  (NTIME,NDEC)

	character	cError(-1:0)	/'E','W'/
	character	cStr*120
	character	cSay*8		/'ReadGIPS'/

	integer		NS		/0/
	save		MJDsav, NS

	Bad = BadR4()

	!-------
	! Check the save scratch arrays for valid data points.
	! NS can only be unequal zero in AUTO mode. In AUTO mode only XCtst1
	! and XCtst2 change. All other data selection criteria remain the same.

	NL = 0			! Initialize counter for data points
	if (NS .eq. 0) then		! Not in AUTO mode or first map
	    MJDstop = MJDfrst		! Stop at earliest available MJD

	else

	    !-------
	    ! Scan the save arrays. Only points with XC >= XCtst1 are needed for 
	    ! subsequent maps. Discard points with XC < XCtst1.

	    N = 0
	    do I=1,NS
		if (XCtst1 .le. XCsav(I)) then
		    N = N+1

		    iMJDsav(N) = iMJDsav(I)
		    iXPsav (N) = iXPsav (I)
		    iYPsav (N) = iYPsav (I)

		    iYRSsav(N) = iYRSsav(I)
		    DOYSsav(N) = DOYSsav(I)
 		    XDSsav (N) = XDSsav (I)
 		    XLSsav (N) = XLSsav (I)
 		    XLLsav (N) = XLLsav (I)
 		    XDLsav (N) = XDLsav (I)
 		    XCEsav (N) = XCEsav (I)
 		    XEsav  (N) = XEsav  (I)
 
		    XCsav  (N) = XCsav  (I)
		    YLsav  (N) = YLsav  (I)
		    VVsav  (N) = VVsav  (I)
		    GGsav  (N) = GGsav  (I)

		endif
	    end do
	    NS = N

	    !--------
	    ! Pick up all points in the save arrays which are needed for the current
	    ! map. Since all points with XC < XCtst1 have just been discarded, the 
	    ! points are only tested for XC <= XCtst2. This takes care for all MJD's
	    ! until the current MJDsav value. Only MJD's later than this (i.e. starting
	    ! at MJDstop = MJDsav+1 have to be scanned for additional points; 
	    ! Remember files are scanned backward in time, so MJDstop is the last file
	    ! to be scanned for the current map.

	    do I=1,NS
		if (XCsav(I) .le. XCtst2) then
		    NL = NL+1

		    iMJD(NL) = iMJDsav(I)
		    iXP (NL) = iXPsav (I)
		    iYP (NL) = iYPsav (I)

		    iYRS(NL) = iYRSsav(I)
		    DOYS(NL) = DOYSsav(I)
 		    XDS (NL) = XDSsav (I)
 		    XLS (NL) = XLSsav (I)
 		    XLL (NL) = XLLsav (I)
 		    XDL (NL) = XDLsav (I)
 		    XCE (NL) = XCEsav (I)
 		    XE  (NL) = XEsav  (I)

		    XC  (NL) = XCsav  (I)
		    YL  (NL) = YLsav  (I)
		    VV  (NL) = VVsav  (I)
		    GG  (NL) = GGsav  (I)
		endif
	    end do

	    MJDstop = MJDsav+1			! Earlier files are processed ..
						! .. and stored in save arrays
	endif

	MJD = min(MJDlast,int(MJDref))		! MJDlast is latest day available
	if (MJD .lt. MJDstop) return		! Check whether any files ...
						! .. have to be scanned
	!-------
	! Read first IPS data file

	call ReadG(cDat,nCar,JDCar,MJD,iEdt,Radius,Speed,Power,iOK,XCmin,XCmax,
     &		iYRS_t,DOYS_t, XDS_t,XLS_t,XLL_t,XDL_t,XCE_t, XC_t,YL_t,XE_t,XR_t,XV_t,XG_t,cStr,.TRUE.)

	if (iOK .ne. 0 .and. .not. bForeCast) then

	    !-------
	    ! Make sure not to miss any data.
	    ! Data points are collected going from large to small XC values (i.e. backwards
	    ! in time. Point with XC in [XCtst1,XCtst2] are needed.

	    ! First scan forward in time to find an MJD for which the range
	    ! [XCmin,XCmax] lies entirely above [XCtst1,XCtst2], i.e. XCmin > XCtst2.

	    MJDsav = MJD
	    do while (MJD .lt. MJDlast .and. iOK .ne. 0 .and. (iOK .eq. 3 .or. XCmin .le. XCtst2))
		MJD = MJD+1
		I = iOK
		call ReadG(cDat,nCar,JDCar,MJD,iEdt,Radius,Speed,Power,iOK,XCmin,XCmax,
     &		iYRS_t,DOYS_t, XDS_t,XLS_t,XLL_t,XDL_t,XCE_t, XC_t,YL_t,XE_t,XR_t,XV_t,XG_t,cStr,.TRUE.)
	    end do
	    if (iOK .eq. 0) then
		MJD = MJD-1
		iOK = I
	    end if
	    if (MJD .ne. MJDsav) then
		write (cStr(itrim(cStr)+1:),'(2(A,I5))') '; MJD upward adjusted: ',MJDsav,' --> ',MJD
		call Say(' ',' ',' ',cStr)
	    end if

	    !-------
	    ! Now scan backwards until the first MJD with data inside the range
	    ! [XCtst1,XCtst2] is found, i.e. until XCmin <= XCtst2 (if the previous
	    ! DO WHILE loop was entered this will the day before the current MJD.
	    ! (This loop isn't strictly necessary, but it doesn't hurt either.)

	    MJDsav = MJD
	    do while (MJD .gt. MJDstop .and. iOK .ne. 0 .and. (iOK .eq. 3 .or. XCmin .gt. XCtst2))
		MJD = MJD-1
		call ReadG(cDat,nCar,JDCar,MJD,iEdt,Radius,Speed,Power,iOK,XCmin,XCmax,
     &		iYRS_t,DOYS_t, XDS_t,XLS_t,XLL_t,XDL_t,XCE_t, XC_t,YL_t,XE_t,XR_t,XV_t,XG_t,cStr,.TRUE.)
	    end do
	    if (MJD .ne. MJDsav) then
		write (cStr(itrim(cStr)+1:),'(2(A,I5))') '; MJD downward adjusted: ',MJDsav,' --> ',MJD
		call Say(' ',' ',' ',cStr)
	    end if

 	end if

	!-------
	! If ReadG returned with a false iOK, then stop (if in AUTO mode) or go back
	! for new prompts (if not in AUTO mode)

	if (iOK .eq. 0) then
	    I = 0
	    if (bAuto) I = -1
	    call Say(cSay,cError(I),'RDERR','error reading MJD '//cStr(:Int2Str(MJD,cStr)))
	    NL = 0
	    return
	end if

	MJDsav = MJD			! Used as MJDstop for next map in AUTO mode

	!-------
	! Read into arrays XC,YL,V until 
	! -  low end XCtst1 is reached, or
	! -  no more files available
	! -  arrays are full
	! Only points in range [XCtst1,XCtst2] are accepted
 
	do while (iOK .eq. 3 .or. XCmax .ge. XCtst1)

	    if (iOK .eq. 1) then		! Data points available
		NLtemp = NL
		! do I=1,NTIME
		do I=NTIME,1,-1
		    do J=1,NDEC
			if (iEorW .eq. 0) then	! All XE (East+West)
			    Rdum = abs(XE_t(I,J))
			else
			    Rdum = iEorW*XE_t(I,J) ! Only positive XE (East) or
			end if			! only negative XE (West)

			if (XC_t(I,J) .ne. Bad .and.
     &				XElow .le. Rdum .and. Rdum .le. XEhigh .and.
     &				abs(XR_t(I,J)) .le. XRlim) then

			    if (XCtst1 .le. XC_t(I,J)) then

				!-------
				! In AUTO mode the next map will be for points in a range
				! [XCtst1+XCincr,XCtst2+XCincr], i.e. shifted to higher XC values 
				! relative to the current range [XCtst1,XCtst2].
				! Collect all points that may be needed for subsequent maps (XC_t >= XCtst1)
				! in the save scratch arrays. Note that this includes points with
				! XC_t > XCtst2, which are not used in the current map.

				if (bAuto) then
				    if (NS .ge. NSmax) return
				    NS = NS+1

				    iMJDsav(NS) = MJD
				    iXPsav (NS) = I
				    iYPsav (NS) = J

				    iYRSsav(NS) = iYRS_t(I,J)
				    DOYSsav(NS)	= DOYS_t(I,J)
 				    XDSsav (NS)	= XDS_t (I,J)
 				    XLSsav (NS)	= XLS_t (I,J)
 				    XLLsav (NS) = XLL_t (I,J)
 				    XDLsav (NS) = XDL_t (I,J)
 				    XCEsav (NS) = XCE_t (I,J)
 				    XEsav  (NS) = XE_t  (I,J)

				    XCsav  (NS) = XC_t  (I,J)
				    YLsav  (NS) = YL_t  (I,J)
				    VVsav  (NS) = XV_t  (I,J)
				    GGsav  (NS) = XG_t  (I,J)

				end if

				!-------
				! Pick up all points inside [XCtst1,XCtst2] for the current map.

				if (XC_t(I,J) .le. XCtst2) then
				    if (NL .ge. NLmax) return
				    NL = NL+1

				    iMJD(NL) = MJD
				    iXP (NL) = I
				    iYP (NL) = J

				    iYRS(NL) = iYRS_t(I,J)
				    DOYS(NL) = DOYS_t(I,J)
 				    XDS (NL) = XDS_t (I,J)
 				    XLS (NL) = XLS_t (I,J)
 				    XLL (NL) = XLL_t (I,J)
 				    XDL (NL) = XDL_t (I,J)
 				    XCE (NL) = XCE_t (I,J)
 				    XE  (NL) = XE_t  (I,J)
 
				    XC  (NL) = XC_t  (I,J)
				    YL  (NL) = YL_t  (I,J)
				    VV  (NL) = XV_t  (I,J)
				    GG  (NL) = XG_t  (I,J)

				end if
			    end if
			end if
		    end do
		end do

		write (cStr(itrim(cStr)+1:),'(A,I4,A)') '; ',NL-NLtemp,' used'
		call Say(' ',' ',' ',cStr)

	    end if

	    if (MJD .eq. MJDstop) return		! All files scanned

	    MJD = MJD-1					! Next data file
	    iOK = -1					! Suppress messages
	    call ReadG(cDat,nCar,JDCar,MJD,iEdt,Radius,Speed,Power,iOK,XCmin,XCmax,
     &		iYRS_t,DOYS_t, XDS_t,XLS_t,XLL_t,XDL_t,XCE_t, XC_t,YL_t,XE_t,XR_t,XV_t,XG_t,cStr,.FALSE.)

	    !-------
	    ! If ReadG returned with a false iOK, then stop (if in AUTO mode) or go back
	    ! for new prompts (if not in AUTO mode)

	    if (iOK .eq. 0) then
		I = 0
		if (bAuto) I = -1
	        call Say(cSay,cError(I),'RDERR','error reading MJD '//cStr(:Int2Str(MJD,cStr)))
		call AskWhat('Cancel, Skip to next day, Plot map$0',I)
		if (I .eq. 0) NL = 0
		if (I .eq. 0 .or. I .eq. 2) return
	    end if

	end do

	!-------
	! Number of points read is now NL

	call Say(cSay,'I','nobs',cStr(:Int2Str(NL,cStr))//' source observations')

	return
	end