C+
C NAME:
C	Extractd
3d
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.  This is the same program as Extracttd.for
C 	except that it extracts an interpolated time series from the density 
C	matrix directly (using a 4-D interpolative scheme rather than from 
C	the 3D base density using XCshift and DVfact and DDfact.  The densities 
C	have been extracted at the times XCint a value centered on local midnight 
C	so that the last values (times) of the 4D matrix are no longer the times 
C	iterated upon in the program.  To correct for this, the xcint file is read 
C	and used to determine the location of the time associated with each 4D
C	matrix.
C CALLING SEQUENCE:
C	call Extractd3D(bExtract,RR,dRR,XCbe,XCint,XCtbe,XCten,nLng,nLat,
C    &			nMap,nT,nTn,nTmax,V3DT,D3DT,Rconst,NCoff,XCstrt,nCar,JDCar)
C INPUTS:
C	bExtract(12)	logical		.TRUE.: extract; .FALSE.: do not extract
C		1<=ITP<=9 Planets (1=Mercury, 2=Venus, 3=Earth, etc.)
C		10=Ulysses, 11=Helios 1; 12=Helios 2
C	RR		real		height of deconvolution
C	dRR		real		distance between heights
C	XCbe(2,nT)	real		start and end times
C	XCint(nTmax)	real		start and end times of intervals
C       XCbeg           real		start Carrington variable of arrays
C       XCend           real		end Carrington variable of arrays
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# heights
C	nT		integer		# of times
C	nTn		integer		# of times of extracted matrix
C	nTmax		integer		maximum # of times
C	V3DT(nLng,nLat,nMap,nTn*nT) real  3 D velocity over time	
C	D3DT(nLng,nLat,nMap,nTn*nT) real  3 D density over time
C
C	Rconst		real		Tomography time constant ~25.38/27.275, but specific to the interval.		
C	NCoff		integer
C	XCstrt		real		Start rotation value
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C OUTPUTS:
C FUNCTIONS/SUBROUTINES:
C 	SunNewcomb, ECLIP_EQUATOR,WR2DARR
C-
	subroutine Extractd3D(bExtract,RR,dRR,XCbe,XCint,XCtbe,XCten,nLng,nLat,
     &		nMap,nT,nTn,nTmax,V3DT,D3DT,Rconst,NCoff,XCstrt,nCar,JDCar)
	parameter	(NTP = 12)
	logical		bExtract(NTP)
	character 	cPrefix(NTP)*2	/'ME','VE','E3','MA','JU','SA','UR','NE','PL','UL','H1','H2'/
	real		V3DT(nLng,nLat,nMap,nTn*nT),
     &			D3DT(nLng,nLat,nMap,nTn*nT)
	real*8		JDCar(nCar),JD,JEpoch,JDCntr
	real		RVect(3),XCbe(2,nTn*nT),XCint(nTmax)
	equivalence	(ObsLng,RVect(1)),(ObsLat,RVect(2)),(ObsDis,RVect(3))
	character 	cFile*14,
     &			cStr*90	/';Year DOY   Hour Height    Lat     Long   Density  Velocity'/,
     &			cDate*23
	logical		bGetLun
	parameter	(NTT = 60)		! 2*NT+1=# times extracted (covers 30 days)
	real*8		dT		/0.25d0/! Time increment for extracted time series
	include		'MAPCOORDINATES.H'
	Bad = BadR4()
	ObsTimNL = 0.0
	Ninterp  = nTn
	ANinterp  = Ninterp
	abe2 =  (XCint(2)    - XCint(1)) *0.5/ANinterp
	aen2 = (XCint(nT+1) - XCint(nT))*0.5/ANinterp
	XCtbeg	= XCint(1)    + abe2 - XCstrt + 1.		! "New" beginning interpolative value 
	XCtend	= XCint(nT+1) - aen2 - XCstrt + 1.		! "New" ending interpolative value
	nT1	= nTn*nT - 1					! "New" number of interpolative intervals
	print *, ' INTO EXTRACTD3D XCstrt, XCtbeg, XCtend =',XCstrt, XCtbeg, XCtend
	dTinc = dT
	XCcntr = XCstrt+0.5			! Carrington variable at center of rotation
	I = XCcntr				! Time when Earth was at center of rotation
	JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*(XCcntr-I)
	do ITP=1,NTP
	    if (bExtract(ITP)) then
						! Heliographic pos of extraction point at JDcntr
		call Julian(1,iYr,Doy,JDCntr,JEpoch)
		call ExtractPositionn(ITP,iYr,Doy,RVect)
		JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*SUN__SiderealP
		I  = JD
		JD = I+nint((JD-I)/dT)*dT
		write (cFile,'(A,A,F8.3)') cPrefix(ITP),'_',NCoff+XCstrt
		if (.not. bGetLun(iU,cFile)) call Say('Extract','E','NOUNIT',
     &			'iGetLun unable to assign logical unit number')
		open (iU,file=cFile,status='NEW')
		call OSGetDateAndTime(cDate)
		write (iU,'(4A)') '; ',cFile(:itrim(cFile)),' created on ',cDate
		write (iU,'(A)' ) cStr(:itrim(cStr))
		XCinc = 0.0
		do I=-NTT,NTT
		    call Julian(1,iYr,Doy,JD+I*dT,JEpoch)
		    call ExtractPositionn(ITP,iYr,Doy,RVect)
		    ObsD = Bad
		    ObsV = Bad
		    ObsXC = XMAP_OBS_POS(XCstrt,ObsLng)
		    if(I.eq.-NTT) Xtbeg = ObsXC - XCstrt + 1.0	! 1'st time
		    if(I.gt.-NTT) then				! Other times
			dObs = ObsXC-ObsXCL
		    else
			dObs = 0.
		    end if
		    ObsXCL = ObsXC
		    if(abs(dObs).gt.0.5) XCinc = XCinc + 1.0
		    ObsXC  = ObsXC + XCinc
C	print *, 'ObsXC', ObsXC
		    ObsTim = Xtbeg + dTinc*(I+NTT)/(25.38
/Rconst)
C		    ObsTim = Xtbeg + dTinc*(I+NTT)/25.38
C	print *, 'ObsTim',ObsTim
C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXC', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXC

		    ObsXCN = ObsXC
		    ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN))

		    ObsTimN = ObsTim
		    ObsTimN = min(XCtend,max(XCtbeg,ObstimN))     ! This is the interval time, not the time to interpolate

		    do JJ = 1, nT
		      XCintjp1 = XCint(JJ+1) - XCstrt + 1.
		      XCintj   = XCint(JJ)   - XCstrt + 1.
C	if(JJ.lt.10) print *, JJ, nT, XCintjp1, XCintj, ObsTimN
		      if (ObsTimN .gt. XCintj .and. ObsTimN .le. XCintjp1) then
C	print *, XCintjp1,XCintj
			AI     = (XCintjp1 - XCintj)/ANinterp
			AI2    = AI*0.5
			AnTint = ObsTimN - XCintj
			nTint  = AnTint/AI + 1.				! This is the interval index value at this time
			if(nTint .lt. 1) then
			  nTint = 1
			  print *, ' nTint LESS THAN 1    	- SET TO 1       !!!!!!!!!!!!!!'
			end if
			if(nTint .gt. Ninterp) then
			  nTint = Ninterp
			  print *, ' nTint GREATER THAN Ninterp - SET TO Ninterp !!!!!!!!!!!!!!'
			end if
			LL = Ninterp*JJ + nTint - Ninterp
			ObsTimNN = XCtfull(LL) + AnTint - nTint*AI + AI2 		! Interpolate time at this time
			Obsdiff  = ObsTimN - ObsTimNN					! Correction difference between times
			ObsTimNN = ObsTimN + Obsdiff					! Interpolate time to use
			ObsXCint = XCint(JJ)+AI*(nTint-0.5) - XCstrt + 1.		! Value of TT at ObsTimN as given by JJ
			ObsdiffP = ObsXCint - ObsTimN					! Difference between TT and ObsTimN
			ObsCheck = XCtvar(ObstimNN) 					! What value is given for interpolation
			IObsChec = (ObsCheck - .5)/ANinterp + 1
			nTinttC  = ((ObsCheck- .5)/ANinterp + 1 - IObsChec)*ANinterp + 1. ! What is the associated nTint value
			ObsXCNN = ObsXCN
C	print *, 'I JJ LL ObsTimN ObsTimNN etc.', I,JJ,LL,ObsTimN,ObsTimNN,nTint,ObsXCint,ObsCheck,IObsChec,nTinttC,ObsdiffP,Obsdiff,AI2
		if(nTint .ne. nTinttC) print *, 'nTint NOT EQUAL TO nTintt', I, JJ, nTint, nTinttC,ObsCheck,IObsChec,nTinttC
C	print *, 'ObsTimN',ObsTimN
C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN
			call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			  nMap,nTn*nT,D3DT,1,ObsTimNN,ObsXCNN,ObsLat,ObsDis,DF)
			call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			  nMap,nTn*nT,V3DT,1,ObsTimNN,ObsXCNN,ObsLat,ObsDis,VF)
			ObsD = Bad
			if (DF .ne. Bad) ObsD = DF/(ObsDis*ObsDis)
		        if(ObsTimN .eq. ObsTimNL) ObsD = Bad
			ObsV = Bad
			if (VF .ne. Bad) ObsV = VF		        
			if(ObsTimN .eq. ObsTimNL) ObsV = Bad
			ObsTimNL = ObsTimN
			ObsD = max(ObsD,-9999.999)			! Should fit output format
			ObsV = max(ObsV,-9999.999)
C			print *, iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsXCN,ObsLat,ObsTimN,ObsLng,dXc,DF,ObsD,VF,ObsV
			write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
		      end if
		     end do
		end do
		iU = iFreeLun(iU)
	    end if
	end do
	return
	end