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, a small correction 
C	factor is determined based on the dirrerent between the XCint values and 
C	Julian date values, and this is applied to the times when the values are
C	extracted from the 4D matricies in this program.  
C	
C CALLING SEQUENCE:
C	call Extractd3D(bExtract,RR,dRR,XCbe,XCint,XCtbeg,XCtend,nLng,nLat,
C    &			nMap,nT,nTn,nTmax,V3DT,D3DT,
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,nT,nTn) real  3 D velocity over time	
C	D3DT(nLng,nLat,nMap,nT,nTn) real  3 D density over time
C		
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,XCtbeg,XCtend,nLng,nLat,
     &		nMap,nT,nTn,nTmax,V3DT,D3DT,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,nT,nTn),
     &			D3DT(nLng,nLat,nMap,nT,nTn)

	real*8		JDCar(nCar),JD,JEpoch,JDCntr

	real		RVect(3),XCbe(2,nT)
	equivalence	(ObsLng,RVect(1)),(ObsLat,RVect(2)),(ObsDis,RVect(3))

	character 	cFile*14,
c		    write (iU,'(I5,I4,4F7.3,2F10.3)') iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
     &			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

	real		XCcorr(2*NTT+1),
     &			XCint(nTmax)

	include		'MAPCOORDINATES.H'

	Bad = BadR4()

	dTinc = dT

	Ninterp  = nTn
	ANinterp  = Ninterp

	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

		if (.not. bGetLun(iU)) call Say('Extract','E','NOUNIT',
     &			'iGetLun unable to assign logical unit number')

		write (cFile,'(A,A,F8.3)') cPrefix(ITP),'_',NCoff+XCstrt
		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
		    if(I.eq.-NTT) then				! Determine differences from times of written files

		        call Julian(1,iYrr,Doyy,JD,JEpoch)	! Determines yr and doy at midpoint JD
		        call ExtractPositionn(ITP,iYrr,Doyy,RVect)! Determines longitude, latitude and distance from iYr and Doy
		        ObsXCM = XMAP_OBS_POS(XCstrt,ObsLng)	! Determines XC value from middle longitude
			AXCcorrL = 10000.
			do J = 1, nT+1
			  AXCcorr = abs(XCint(J) - ObsXCM - 1.)
			  if(AXCcorr .lt. AXCcorrL) then		! Find minimum difference
			    AXCcorrL = AXCcorr
			    Jsave    = J 
			  end if
			end do
			AXCcorr = (XCint(Jsave) - ObsXCM - 1.)		! Change to the real value
C			print *, 'Jsave, AXCcorr', Jsave, AXCcorr
		    end if
		    call Julian(1,iYr,Doy,JD+I*dT,JEpoch)
	! Advances yr and doy from +- midpoint JD

		    call ExtractPositionn(ITP,iYr,Doy,RVect)
	! Determines longitude, latitude and distance from iYr and Doy

		    ObsD = Bad
		    ObsV = Bad

		    ObsXC = XMAP_OBS_POS(XCstrt,ObsLng)
	 ! Determines XC value from longitude relative to start 

		    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

		    JJ = Jsave + I*dTinc			! Calculates X correction values
		    if(JJ.lt.1) JJ = 1
		    if(JJ.gt.nT+1) JJ = nT+1
		    Anterp = (XCint(JJ+1) - XCint(JJ))*(I - JJ/dTinc + Jsave/dTinc)*dTinc
		    XCcorr(I+NTT+1) = ObsXC - XCint(JJ) + AXCcorr - Anterp

		    ObsTim = Xtbeg + dTinc*(I+NTT)/25.38

C	print *, '/ I XCstrt ObsLng ObsLat ObsXCL ObsXC XCint(JJ) JJ NTT', I, XCstrt, ObsLng, ObsLat, OBsXCL, ObsXC
, XCint(JJ), JJ, NTT
		    ObsXCN = ObsXC
		    ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN))

		    ObsTimN = ObsTim
 - XCcorr(I+NTT+1)			 ! This corrects for the incorrect time extraction
	 	    AI    = (XCint(JJJ+1) - XCint(JJJ))/ANinterp

		    JJJ   = XCtvar(ObsTimN + AI)			! This calculates the time index value of time
		    nTint = (XCtvar(ObsTimN + AI)-JJJ)*Ninterp + 1.0	! This calculates the interval value at this time
		    if(nTint.lt.1) then
			nTint = 1
			print *, ' nTint LESS THAN 1 !!!!!!!!!!!!!!'
		    end if
		    if(nTint.gt.4) then
			nTint = 4
			print *, ' nTint GREATER THAN 4 !!!!!!!!!!!!!!'
		    end if

		    ObsTimN = min(XCtend,max(XCtbeg,ObstimN))

C	print *, 'ObsTim, ObsTimN, Anterp, XCcorr(I+NTT+1), JJJ, XCtvar(ObsTim), nTint', ObsTim, ObsTimN, Anterp
, XCcorr(I+NTT+1), JJJ, XCtvar(ObsTim), nTint

C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN

	print *, 'XCtbeg, ETC', XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT
	print *, 'XCstrt PARAMS', XCstrt,ObsTimN,ObsXCN,ObsLat,ObsDis

C		    call Get4Dval(XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
C     &			nMap,nT,XCshift,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXc)

		    if (dXC .eq. Bad) stop 'Extract: oops'

		    call Get4Dval(XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,D3DT(1,1,1,1,nTint),1,ObsTimN,ObsXCN,ObsLat,ObsDis,DF)

		    call Get4Dval(XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,V3DT(1,1,1,1,nTint),1,ObsTimN,ObsXCN,ObsLat,ObsDis,VF)

C		    ObsXCN = ObsXC+dXC
C		    ObsXCN = min(XCbe(2,nT),max(XCbe(1,1),ObsXCN))
C		    ObsTimN = ObsTim+dXC
C		    ObsTimN = min(XCtend,max(XCtbeg,ObstimN))

C		    call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DD,1,
C     &			ObsXCN,ObsLat,ObsTimN,ObsD)

C		    print *, 'DF ObsD', DF, ObsD

		    ObsD = Bad
		    if (DF .ne. Bad) ObsD = DF
/(ObsDis*ObsDis)

C		    call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1,
C     &			ObsXCN,ObsLat,ObsTimN,ObsV)

C		    print *, 'VF ObsV', VF, ObsV

		    ObsV = Bad
		    if (VF .ne. Bad) ObsV = VF

		    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 do

		iU = iFreeLun(iU)

	    end if
	end do

	return
	end