C+
C NAME:
C	ReadSB8_d.f
C PURPOSE:
C	Read STEREO B density data into arrays
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call ReadSB8_d(cWildSB,nCar,JDCar,XCtest1,XCtest2,Dlimitu,DlimitL,NCoff,Radius,
C     &		NLmax,NLmaxs,NL,N_IN,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
C     &		IYRF,IREC,IYRS,DOYS8,DIST,XLS,XLL,XDL,XCE,XE,XC,YL,DD,NBS)
C
C INPUTS:
C   ReadSA8_d
C	cWildSB		character	location of the STEREO B input file
C	nCar		integer		dimension of JDCar
C	JDCar(nCar)	integer		start times (Julian days) of Carrington rotations
C	XCtest1		real		modified Carrington variable for start of search
C	XCtest2		real		modified Carrington variable for end of search
C	Dlimitu		real		upper density limit
C	DlimitL		real		lower density limit
C	NCoff		integer		JDCar(I) is the start of rotation NCoff+I
C	Radius		real		reference distance (AU)
C
C	XElow		real
C	XEhigh		real
C
C	NLmax		integer		max. # points read into arrays
C	NLmaxs		integer		max. # points read into small arrays
C	NL		integer
C
C OUTPUTS:
C	NL		integer		# valid data points in arrays XC,YL,DD
C	SRC(NLmax)	character*8	Character Source identifier
C	IYRF(NLmax)	integer		year of file from which observation was extracted (should be same as IYRS)
C	IREC(NLmax)	integer		record number on file IYRF
C	DOYS8(NLmax)	real*8		time of observation: day of year (incl. fraction for time of day)
C	XDS (NLmax)	real		Sun-STEREO B distance
C	XLS (NLmax)	real		geocentric ecliptic longitude Sun
C	XLL (NLmax)	real		geocentric ecliptic lng(P)-lng(Sun) (deg)
C	XDL (NLmax)	real		geocentric ecliptic lat(P) (deg)
C	XCE (NLmax)	real		modified Carrington variable of sub-STEREO B point on Sun
C	XE  (NLmax)	real		elongation (deg) (>0: East of Sun; <0: West of Sun)
C	XC  (NLmax)	real		modified Carrington variable point P after traceback to
C					  heliocentric distance Radius at speed VV
C	YL  (NLmax)	real		heliographic latitude point P
C	DD  (NLmax)	real		STEREO B Density (Np)
C	NBS (NLmax)	integer		The type of source - 14 for STEREO-B
C	NSmax		integer		# data points saved in XCsav,YLsav
C	SRCsav(NSmax)	character	scratch arrays (used internally only),
C	IYRFsav(NSmax)	integer		
C	IRECsav(NSmax)	integer
C	IYRSsav(NSmax)	integer
C	DOYSsav(NSmax)	real*8
C	XDSsav (NSmax)	real
C	XLLsav (NSmax)	real
C	XDLsav (NSmax)	real
C	XCEsav (NSmax)	real
C	XEsav  (NSmax)	real
C	XCsav  (NSmax)	real
C	YLsav  (NSmax)	real
C	DDsav  (NSmax)	real
C MODIFICATION HISTORY:
C	April-2009, Bernard Jackson (UCSD)
C-
	subroutine ReadSB8_d(cWildSB,nCar,JDCar,XCtest1,XCtest2,Dlimitu,DlimitL,NCoff,Radius,
     &		NLmax,NLmaxs,NL,N_IN,cvgfiles,SRCVG,SRC,
     &		IYRS,DOYS8,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,DD,NBS)

	real*8		JDCar(nCar),
     &			JEpoch,
     &			Doy8,
     &			DOYS8(NLmax)			! Time: doy of year

	real*8		dLngSun, dLatSun, dDisSun
	
	integer		IYRS(NLmax),
     &			NBS (NLmax)

	real		XDS (NLmax),			! Sun-Earth distance
     &			XLS (NLmax),			! Ecliptic longitude Sun
     &			XLL (NLmax),
     &			XDL (NLmax),
     &			XCE (NLmax),
     &			XE  (NLmax),			! Elongation point P
     &			XC  (NLmax),			! Carrington variable of Point P
     &			YL  (NLmax),			! Heliographic latitude of Point P
     &			DD  (NLmax)			! Solar wind STEREO B density (Np)

	real		RVect(3)
	real		VVect(5)

	character	cvgfiles(*)*11,SRCVG(NLmaxs)*115,SRC(NLmaxs)*10,
     &			cFile*80,
     &			cWildSB*(*),
     &			cWildINSITU*80
     
C	function iReadNagoyan8(iU,iFirst)
	character	cMon*3,
     &			cSrc*10,cSRCVG*66
C	character	cFmt41*41		/'(A9,3I2,F6.2,F5.2,I4,I5,I4,3I5,I4,F8.5)'/
	character	cFmt31*31		/'(A10,3I2,F6.2,27X,I5,4X,F8.5)'/
	character	cFmt35*35		/'(I4,1X,I3,1X,I2,7X,F5.2,6X,E12.6)'/

	external	EARTH8
	include		'sun.h'
	include		'dirspec.h'

	ITPA = 13						! STEREO A
	ITPB = 14						! STEREO B

	call Julian(1,iYr1,Doy,JDCar(int(XCtest1)),JEpoch)
	write (*,'(A,F9.4,A,I5)') 'In readSB8_d. The in-situ analysis begins after Doy',Doy,' of',iYr1

	call Julian(1,iYr2,Doy,JDCar(int(XCtest2+1.0)),JEpoch)
	write (*,'(A,F9.4,A,I5)') 'In readSB8_d. The in-situ analysis ends before Doy',Doy,' of',iYr2
	
	Iyrno = 1
	if(iYr1.ne.iYr2) Iyrno = 2

	I = N_IN
	Ifst = 1
	do N=1,Iyrno
	 if(Ifst.eq.1) iYr = iYr1
	 if(Ifst.eq.2) iYr = iYr2
         write (*,*) cWildSB
	 cWildINSITU = cWildSB

         III = iSetFileSpec(cWildINSITU)
         III = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile)
         iFile = index(cFile, cWildChar(:4))
         III = 1
         call Str2Flt_Int(.TRUE.)
         call Str2Flt(cFile(iFile:iFile+4-1),III,0)
         iFile = index(cWildINSITU,cWildChar(:4))
         cFile = cWildINSITU
         III = Int2Str(iYr,cFile(iFile:iFile+4-1))
         write (*,*) cFile

	 open (13, file=cFile,status='old',recl=120,access='sequential',form='formatted',iostat=iReadSB)
	 if(iReadSB.ne.0.and.Ifst.eq.2) then
	   close(13)
           return	! assume if no file exists on second year, the second year has not yet happened
	 end if

C	 print *, cfile

	 do II=1,10000

	   read (13,cFmt35,iostat=iReadSB) iY,iDoy,iUT,DEN,VEL

	   if(iReadSB.eq.0) then

	     call DATE_DOY(1,iY,cMon,iM,iD,iDoy)

	     write(cSrc(1:10),'(A2,4I2)') 'SB',iY-2000,iM,iD,iUT
	
	     UT   = float(iUT)
	     VEL  = VEL/100.0					! Convert velocity from m/s to km/s 
	     iVEL = nint(VEL)

	     write(cSRCVG(1:66),cFmt31)
     &		cSrc,iY-2000,iM,iD,UT,iVEL,DEN

	     if (DEN.ge.DlimitL.and.DEN.le.Dlimitu) then	! If density inside limits, accept	  
	       nYr      = iY
	       cMon = '   '					! Make sure iM is used
	       call DATE_DOY(0,nYr,cMon,iM,iD,iDoy)
	       Doy8  = iDoy+dble(UT/24.)
 	       XCcntr = XMAP_SC_POS8(EARTH8,nYr,Doy8,-nCar,JDCar)
	       if (XCcntr .ne. BadR4()) then			! Time outside JDCar range
	         XCbeg = XCtest1
	         XCend = XCtest2
	         if (XCcntr.ge.XCbeg.and.XCcntr.le.XCend) then	! Sub-Earth past XCend
	           I = I + 1
C	           if(I.lt.(N_IN+5)) print *, 'XCcntr =',XCcntr,' Sub-Earth after XCbeg =',XCbeg,' or before XCend =', XCend, nYr,DOY8
	           IYRS(I+NL)  = nYr
	           SRCVG(I+NL) = cSRCVG
	           SRC(I+NL)   = cSrc
	           DOYS8(I+NL) = Doy8
	           DD(I+NL)    = DEN				! Density
		   NBS(I+NL)   = ITPB				! Spacecraft number
		   call get_scparams(ITPB,nCar,JDCar,NCoff,NLmax,I+NL,nYR,Doy8,XCcntr,Radius,VEL,			! Stereo B -14
     &		                XDS,XLS,XLL,XDL,XCE,XE,XC,YL,RA,DEC)
		   XLL(I+NL) = XLL(I+NL) - 360.0		! STEREO B
C		   if(I.lt.(N_IN+5)) then
C		     write(*,'(3I10,2A)') I,NL,IYRS(I+NL),SRCVG(I+NL),SRC(I+NL)
C		     write(*,'(I5,12F12.4)') 
C     &			I+NL,DOYS8(I+NL),XDS(I+NL),XLS(I+NL),RVect(1),XLL(I+NL),XDL(I+NL),XCE(I+NL),XE(I+NL),XC(I+NL),YL(I+NL),DD(I+NL),RA,DEC     
C		   end if
		 end if
	       end if
	     end if
	   end if
	 end do
	 Ifst = Ifst + 1
	 close(13)
	end do
	print *, 'Number of STEREO B in-situ density sources =',I-N_IN
	N_IN = I
	print *, ' '
	return
	end