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:
C	call MkTimes(bForecast,MJDcntr,NmidHR,aNday,NTmax,ALng,nLng,nCar,JDCar,
C	   inTn,nTn,iYr,Idaybeg,XCintD,XCinT,XCbe,XCtbeg,XCtend,IYRS,DOYS)
C INPUTS:
C	bForecast	logical		.TRUE. Forecast mode
C					.FALSE. Non forecast mode
C	MJDcntr		real*8		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		real*8		Julian Date at beginning of rotations
C	inTn		integer		0 - set inteval number internally, >0 set interval number to this 
C	nTn		integer		number of intervals to use
C OUTPUTS:
C	nTn		integer		# of intervals (set by this routine depending on bForecast, inTn, NTMmax)
C	iYr		integer		The beginning year
C	Idaybeg		integer		Number of days in the beginning time year (varies with leap year)
C	XCinTD(NTmax)	real*8		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		real*8		Start time of sources in doy (Time from beginning)
C       XCtend          real*8		End time of sources in doy (Time from beginning)
C	IYRS(NTmax)	integer		Interval year values
C	DOYS(NTmax)	real*8		Interval day of year values
C FUNCTIONS/SUBROUTINES:
C 	Julian
C-
	subroutine MkTimes(bForecast,MJDcntr,NmidHR,aNday,NTmax,Alng,nLng,nCar,JDCar,
     &		inTn,nTn,iYr,Idaybeg,XCintD,XCinT,XCbe,XCtbeg,XCtend,IYRS,DOYS)

	real		XCint(NTmax),
     &			XCbe(2,NTmax),
     &			DOYS(NTmax)
     
     	real*8		XCintD(NTmax)
     
     	real*8		XCtbeg,
     &			XCtend,
     &			dXC,dXC2,
     &			Doy8

	real*8		XCintdd

	integer		IYRS(NTmax)

	logical		bForecast

	real*8		JDCar(nCar),JD,JEpoch,JDXC,MJDcntr

	external	EARTH8
	
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.
C
	nTn = nint(float(nTn)/aNday)

	do I=1,NTmax
	  XCint(I)  = 0.0
	  XCbe(1,I) = 0.0
	  XCbe(2,I) = 0.0
	  DOYS(I)   = 0.0
	  XCintD(I) = 0.0
	  IYRS(I)   = 0
	end do

        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
	    
C	ibeg = nint(33./aNday)						! Old version before 11/11/08 BVJ
C	ibeg = nint(25./aNday)						! Number of intervals before midpoint to begin
	ibeg = nint((1.0*27.2753)/aNday)				! Number of intervals before midpoint to begin
	    
	endday = 20.							! 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						! Automatically determines a minimum length
	if(nTn.lt.nT) then
	  nT  = nTn							! Accepts and sets a length lesser than the automatically set limit
	else
	  nTn = nT							! Otherwise returns the automatically set length 
	end if
	if(inTn.gt.0) then						! Placed here 03/05/09 by BVJ work with the intel analysis
	  nT  = inTn-1							! The length inTn is what is used if inTn is non-zero 
	  nTn = inTn-1
	end if

	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
	
C	print *, 'To here 1'
	XCM = XMAP_SC_POS8(EARTH8,iYrM,Doy8,nCAR,JDCar)				! Earth location in Carrington coords. at Doy8
C	print *, 'To here 2'
		
	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)
C	    print *, 'To here 3',iYrr,XCintdd,XCtbeg,dXC,dXC2,N,NT
	    XCint(N)  = XMAP_SC_POS8(EARTH8,iYrr,XCintdd,nCAR,JDCar)		! Interval ends in single precision XC's
C	    print *, 'To here 4',iYrr,XCintdd,N,XCint(N)
	    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