C+
C NAME:
C	MkTimes
C PURPOSE:
C	Makes Carrington longitude intervals and year and doy intervals for
C	the time dependent tomography program.
C CALLING SEQUENCE:
	subroutine MkTimes(bForecast,MJDcntr,NmidHR,aNday,NTmax,Alng,nLng,nCar,JDCar,
     &		NT,iYr,Idaybeg,XCintD,XCinT,XCbe,XCtbeg,XCtend,IYRS,DOYS)
C INPUTS:
C	bForecast	logical		.TRUE. Forecast mode
C					.FALSE. Non forecast mode
C	MJDcntr		double precision Modified Julian date at the time of the observation
C	NmidHR		integer		Number of hours before midnight
C	aNday		real		Number of days to combine for digital result (can be fractions of days)
C	NTmax		integer		# of intervals
C	ALng		real		Carrington interval per nLng number
C	nLng		integer		# of longitude resolution intervals in aLng
C	nCar		integer		# Carrington rotations
C	JDCar		double precision Julian Date at beginning of rotations
C OUTPUTS:
C	NT		integer		# of intervals (set by this routine depending on bForecast)
C	iYr		integer		The beginning year
C	Idaybeg		integer		Number of days in the beginning time year (varies with leap year)
C	XCinTD(NTmax)	double precision Time interval values in DOY since beginning value XCtbeg at beginings of each interval
C	XCinT(NTmax)	real		Time interval values in XC since beginning value XCtbeg at beginings of each interval
C	XCbe(2,NTmax)	real		Carrington longitude begining and end values
C       XCtbeg		double precision Start time of sources in doy (Time from beginning)
C       XCtend          double precision End time of sources in doy (Time from beginning)
C	IYRS(NTmax)	integer		Interval year values
C	DOYS(NTmax)	double precision Interval day of year values
C EXTERNAL:
	external	EARTH8
C CALLS:
C	Julian8, XMAP_SC_POS8
C MODIFICATION HISTORY:
C	APR-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Bug fix. Variable Leap was never initialized if XCtbeg >= 1.
C	OCT-2007, ???
C	    Added line forcing nT to be an even number
C-
	real		XCint(NTmax)
     	real            XCbe(2,NTmax)
     	real	        DOYS(NTmax)
     
     	double precision XCintD(NTmax)
     
     	double precision XCtbeg
        double precision XCtend
        double precision dXC
        double precision dXC2
        double precision Doy8

	double precision XCintdd

	integer		IYRS(NTmax)

	logical		bForecast

	double precision JDCar(nCar)
        double precision JD
        double precision JEpoch
        double precision JDXC
        double precision MJDcntr

C	print *, ' '
C	print *, 'xInc, XCtest1, XCtest2, NT, aNday', xInc, XCtest1, XCtest2, NT, aNday
C	Only allow times in multiples of UT days - unlike earlier tomography.

	call Julian8(11,iYrM,Doy8,MJDcntr,JEpoch)		! Determine iYr and Doy8 from MJDcntr

	iDoy  = nint(Doy8)
	Doy8 = dble(iDoy) + dble(NmidHR/24.0 + 0.5)		! UT local noon at time of day at interval midpoint to begin
	ibeg = nint(33./aNday)

	endday = 14.						! Normally go 14 days beyond Carrington rotation middle
	if(bForecast) endday = 5.				! Add five days beyond present in forecast mode
	iend = nint(endday/aNday)

	nT = iend + ibeg + 1
	if(nT.ne.(nT/2)*2) nT = nT + 1				! Make nT an even number
	if(nT.ge.NTmax) then
	    print *, 'Errror! nT =', nT,'  is greater than NTmax-1 =',NTmax
	    NT = NTmax-1
	end if

	XCtbeg = Doy8 - dble(ibeg*aNday)			! Start of interval in DOY at interval midpoint

	iYr = iYrM
	Leap = 0						! Suppose it's not a leap year
	if(XCtbeg.lt.1.0D0) then				! Remember DOY begins with 1.0
	    iYr = iYrM - 1					! Year associated with beginning time
	    if ((iYr)/4*4 .eq. (iYr) .and. ((iYr) .lt. 1582 .or.! It's a leap year
     &         iYr/100*100 .ne. iYr .or. iYr/400*400 .eq. iYr)) Leap = 1
	    XCtbeg = XCtbeg + dble(365+Leap)
	end if

	Idaybeg = 365+Leap					! Number of days in beginning iYr year
	XCtend  = XCTbeg + dble((nT-1)*aNday)
	dXC     = (XCtend - XCtbeg)/dble(nT-1)			! Interval length in days

C	print *, 'XCtbeg, XCtend, dXC, NT', XCtbeg, XCtend, dXC, nT 

	dXC2 = dXC/2.0d0

	ALng2 = ALng/2
	XCM = XMAP_SC_POS8(EARTH8,iYrM,Doy8,nCAR,JDCar)		! Earth location in Carrington coords. at Doy8
	aNdeg = (nLng-1)/alng					! Number of resolution intervals per rotation 
	IXCM = nint(XCM*aNdeg)					! Set midpoint at even rotation resolution
	XCM = IXCM/aNdeg

	do N=1,nT+1
	    XCintD(N) = XCtbeg + dXC*dble(N-1)			! Tomography times in double precision Doys
	    XCintdd = XCtbeg + dXC*dble(N-1) - dXC2		! Interval ends in dble prec.Doys (between tomo. times)
	    iYrr = iYr
	    if(XCintdd.gt.dble(Idaybeg)) then
		iYrr = iYr + 1
		XCintdd = XCintdd - dble(Idaybeg)
	    end if
	    IYRS(N)   = iYrr
	    DOYS(N)   = sngl(XCintdd)				! Interval ends in single precision XCs
	    XCint(N)  = XMAP_SC_POS8(EARTH8,iYrr,XCintdd,nCAR,JDCar)
	    if(N.lt.(NT+1)) then
	        XCbe(1,N) = XCM - ALng2				! Rotation beginnings in single precision XC's
	        XCbe(2,N) = XCM + ALng2				! Rotation ends in single precision XC's
	    end if
	end do
		
	return
	end