C+
C NAME:
C	Extractd
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.
C CALLING SEQUENCE:
C	call Extractd(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,FALLOFF,
C    &		        VV,DD,XCshift,DVfact,DDfact,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 ending Carrington variable of arrays
C	XCtbeg		real*8		beginning time of observations
C       XCtend          real*8		ending time of observations
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMap		integer		# heights
C	nT		integer		# times
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		Start rotation value
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C OUTPUTS:
C CALLS:
C 	SunNewcomb, ECLIP_EQUATOR,WR2DARR
C	Time2Str, cTime2System, itrim
C-
	subroutine Extractd(bExtract,RR,dRR,XCbe,XCtbeg,XCtend,nLng,nLat,nMap,nT,FALLOFF,
     &		            VV,DD,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)
	parameter	(NTP = 12)
	logical		bExtract(NTP)
	character 	cPrefix(NTP)*2	/'ME','VE','EA','MA','JU','SA','UR','NE','PL','UL','H1','H2'/
	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
	
	real*8		XCtbeg,XCtend,XCtfull,X,Doy8,ObsTim,ObsTimN
	
	real*8		dT		/0.25d0/	! Time increment for extracted time series
	real*8		DoypY		/0.0d0/
	
	real 		ALng		/3.0/
	
	real		RVect(3),XCbe(2,nT)

	equivalence	(ObsLng,RVect(1)),(ObsLat,RVect(2)),(ObsDis,RVect(3))

	character 	cFile*14
	character	cStr*90	/';Year DOY   Hour Height    Lat     Long   Density  Velocity'/
	character       cTime2System*80

	logical		bGetLun

	parameter	(NTT = 60)			! 2*NT+1=# times extracted (covers 30 days)
	
	external	EARTH8

	include		'mapcoordinates.h'
	
C	print *, ' Into extractd.f XCstrt, XCtbeg, XCtend =',XCstrt, XCtbeg, XCtend

	Bad = BadR4()
	dTinc = 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)
	do ITP=1,NTP
	    if (bExtract(ITP)) then
						! Heliographic pos of extraction point at JDcntr
		call Julian8(1,iYr,Doy8,JDCntr,JEpoch)
		call ExtractPositionn8(ITP,iYr,Doy8,RVect)
		JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*SUN__SynodicP
		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')

		write (iU,'(3A,A19)') '; ',cFile(:itrim(cFile)),' created on ',cTime2System(' ')
		write (iU,'(A)' ) cStr(:itrim(cStr))

		XCinc = 0.0
		iYrL  = 0
		do I=-NTT,NTT
		    call Julian8(1,iYr,Doy8,dble(JD+I*dT),JEpoch)
		    call ExtractPositionn8(ITP,iYr,Doy8,RVect)
		    ObsD = Bad
		    ObsV = Bad		    
		    ObsXC = XMAP_OBS_POS(XCstrt,ObsLng)
		    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 = 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))
		    
		    if(iYr.gt.iYrL.and.I.gt.-NTT) then					! Year change
	    		Leap = 0							! Suppose it's not a leap year
	    		if ((iYrL)/4*4 .eq. (iYrL) .and. ((iYrL) .lt. 1582 .or. 	! It's a leap year
     &	        	iYrL/100*100 .ne. iYrL .or. iYrL/400*400 .eq. iYrL)) Leap = 1
     	    		DoypY = dble(365+Leap)
			iYrL = iYr
     		    end if	
		    iYrL = iYr

		    ObsTimN = ObsTim + DoypY
		    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)

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

		    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

C		    print *, 'VF ObsV', VF, ObsV
		    if (VF .ne. Bad .and. ObsV .ne. Bad) ObsV = ObsV*VF
		    ObsD = max(ObsD,-9999.999)			! Should fit output format
		    ObsV = max(ObsV,-9999.999)
		    write (iU,'(I5,I4,3F7.3,3F9.3)') iYr,int(Doy8),(Doy8-int(Doy8))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
		end do
		iU = iFreeLun(iU)
	    end if
	end do
	return
	end