C+
C NAME:
C	Extractdvd
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.
C CALLING SEQUENCE:
C	call Extractdvd(bForecast,bExtract,cPrefix,NTP,RR,dRR,XCbe,XCtbeg,XCtend,iYr,nLng,nLat,
C    	            nMap,nT,NinterD,FALLOFF,VV,DD,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)
C INPUTS:
C	bForecast	logical		If .TRUE. this is forecast mode and then add a few days to the output
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	cPrefix		character	Prefix for the outparameter file
C	NTP		integer		Total number of prefixes
C	RR		real		height of deconvolution
C	dRR		real		distance between heights
C       XCbe(2,nT)      real		start and ending Carrington variable of arrays
C	XCtbeg		real*8		beginning time of tomography (in days)
C       XCtend          real*8		ending time of tomography (in days)
C	iYr		integer		year at the beginning time of the tomography
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# heights
C	nT		integer		# times
C	NinterD		integer		# of intermediate data times interpolated
C	FALLOFF		real		Density falloff with solar distance
C	VV(nLng,nLat,nT)real		velocity map at height RR
C	DD(nLng,nLat,nT)real		density map at height RR
C	XCshift(nLng,nLat,nMap,nT,3) real shift amount from array
C	DVfact(nLng,nLat,nMap,nT) real  velocity difference at height		
C	DDfact(nLng,nLat,nMap,nT) real  density difference at height		
C	NCoff		integer
C	XCstrt		real		Output data start value (in rot). For archived data - rotation beginning. Forecast -0.5 rotation fm cntr
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 Extractdvd(bForecast,bExtract,cPrefix,NTP,RR,dRR,XCbe,XCtbeg,XCtend,iYr,nLng,nLat,
     &			nMap,nT,NinterD,FALLOFF,VV,DD,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)

	logical		bExtract(NTP)
	logical		bForecast
	character 	cPrefix(NTP)*2
	real		VV(nLng,nLat,nT),
     &			DD(nLng,nLat,nT),
     &			XCshift(nLng,nLat,nMap,nT,3),
     &			DDfact (nLng,nLat,nMap,nT),
     &			DVfact (nLng,nLat,nMap,nT),
     &			dXC(3)

	real*8		JDCar(nCar),JD,JEpoch,JDCntr,JDStrt,JDBeg
	
	real*8		XCtbeg,XCtend,XCtF,XCtE,XCtfull,X,Doy8,ObsTim,ObsTimN
	
	real*8		dT				! Time increment for extracted time series
	
	real 		ALng		/3.0/
	
	real		RVect(3),XCbe(2,nT)

	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
	
	external	EARTH8

	include		'MAPCOORDINATES.H'
	
	print *, ' Into extractdvd.f nT, XCstrt, XCtbeg, XCtend =',nT,XCstrt, XCtbeg, XCtend

	XCtF = 0.0d0						! number of days past the tomography end time
	XCtE = 0.5d0*dble(Sun__SynodicP)			! Number of days past the end Carrington rotation center time	
	if(bForecast) then
	  XCtF = 1.0d0						! number of days past the tomography end time time
	  XCtE = 4.0d0						! Number of days past the forecast time					
	end if

	Bad = BadR4()
	
	dT = (XCtend - XCtbeg)/(dble((NinterD+1)*(nT-1)))	! extract interval length in days - typically 0.25 days or 6 hours

	NTTT = (XCtend+XCtF - XCtbeg)/dT + 1			! total number of time steps possible from beginning of the tomography

C	print *, 'Extractdvd Interval length in days', dT
	
	XCcntr = XCstrt + 0.5					! Carrington variable at center of rotation
	I = XCcntr						! Carrington interval before Earth was at center of rotation
	JDCntr = JDCar(I)+(JDCar(I+1)-JDCar(I))*dble(XCcntr-I)  		 ! Julian date in days at center of rotation
	call Julian8(0,iYr,XCtbeg,JD,JEpoch) 			! Julian date at XCtbeg
	I = int(JD)
	JDStrt = I + nint((JD-I)/dT)*dT				! Julian date at XCtbeg to nearest value of dT
	I = int(JDCntr) 
	JDBeg  = I + nint((JDCntr-I)/dT)*dT - 15.0d0		! Beginning Julian date 15 days before center of rotation	

	print *, iYr, XCcntr, JDCntr, JDStrt, JDBeg

	do ITP=1,NTP
	    if (bExtract(ITP)) then
		print *, 'The location extracted is: ', cPrefix(ITP)
						! Heliographic pos of extraction point at JDcntr
		call Julian8(1,iYrc,Doy8,JDCntr,JEpoch)		
		call ExtractPositionn8(ITP,iYrc,Doy8,RVect)
		JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*Sun__SynodicP
		I  = int(JD)
		JD = I+nint((JD-I)/dT)*dT					! Julian date in days at observer location to nearest value of 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

		XCinc = 0.0
		ObsXCL = -100.							! last observer 
		xxinc = 0.0	

		do I=1,NTTT
	
		  dTT = (dT * dble(I-1))
		  JD   = JDStrt + dTT
		  Doy8 = XCtbeg + dTT

C		print *, 'to here 1', JD, JDBeg, Doy8, XCtbeg, XCtend+XCtF, JDCntr+XCtE

		  if(JD.ge.JDBeg.and.Doy8.gt.XCtbeg.and.Doy8.le.XCtend+XCtF.and.JD.le.JDCntr+XCtE) then		! Limits on extracted data

C		print *, 'to here 2', JD, JDBeg, Doy8, XCtbeg, XCtend+XCtF, JDCntr+XCtE 

		    call ExtractPositionn8(ITP,iYr,Doy8,RVect)
		    ObsD = Bad
		    ObsV = Bad		    
		    ObsXC = XMAP_OBS_POS(XCstrt,ObsLng)

		    if(ObsXCL.gt.-99.) then						! Changed from previous on 20 April 2011 BVJ
			dObs = ObsXC-ObsXCL
		    else
			dObs = 0.
		    end if
		    ObsXCL = ObsXC
		    if(abs(dObs).gt.0.5) XCinc = XCinc + 1.0 + xxcinc
		    ObsXC  = ObsXC + XCinc
C	print *, 'ObsXC', ObsXC
		    ObsTim = Doy8
		    
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))

C	print *, 'ObsTimN',ObsTimN

		    call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,XCshift,1,ObsTimN,ObsXCN,ObsLat,ObsDis,dXC)
     
C	print *, 'XCstrt ObsLng ObsLat ObsXCL ObsXCN dXC ', XCstrt, ObsLng, ObsLat, OBsXCL, ObsXCN, dXC(1), dXC(2), dXC(3)

C		    if (dXC(1) .eq. Bad) stop 'Extract: oops 1'
C		    if (dXC(2) .eq. Bad) stop 'Extract: oops 2'
C		    if (dXC(3) .eq. Bad) stop 'Extract: oops 3'

C	write(*,'(A,3F12.4)') 'XCstrt ObsXCN ObsLng', float(NCOFF)+XCstrt, float(NCOFF)+ObsXCN, ObsLng
C	write(*,'(A,3F12.4)') 'Doy8, ObsTimN, dXC(1) ', Doy8, ObsTimN, dXC(1)

		    if (dXC(1).ne.Bad.and.dXC(2).ne.Bad.and.dXC(3).ne.Bad) then

		      call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nT,DDfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,DF)
		      call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,
     &			nMap,nt,DVfact,1,ObsTimN,ObsXCN,ObsLat,ObsDis,VF)

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

		      ObsLa  = ObsLat  + dXC(2)

		      ObsTimN = ObsTim + dXC(3)
		      ObsTimN = min(XCtend,max(XCtbeg,ObstimN))

		      call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DD,1,
     &			ObsXCN,ObsLa,ObsTimN,ObsD)
C		      print *, 'DF ObsD', DF, ObsD
		      if (DF .ne. Bad .and. ObsD .ne. Bad) ObsD = ObsD*DF*((RR/ObsDis)**FALLOFF)
		      call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VV,1,
     &			ObsXCN,ObsLa,ObsTimN,ObsV)

C		      print *, 'After = ', JD+I*DT, iYr, DF, VF, Doy8, ObsXCN, ObsLa, ObsTimN, ObsDis, ObsD, ObsV

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

		      print *, 'DF VF, ObsD ObsV', DF, VF, ObsD, ObsV

		      ObsD = max(ObsD,-9999.999)			! Should fit output format
		      ObsV = max(ObsV,-9999.999)
		    else
C		      print *, cPrefix(ITP),' is out of range of the 3D matrix' 
		      ObsD = -9999.999
		      ObsV = -9999.999
		    end if

		    write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
		  end if
		end do
		iU = iFreeLun(iU)
	    end if
	end do
	return
	end