C+
C NAME:
C	ExtractInsitu
C PURPOSE:
C	Extracts interpolated velocity and density IPS deconvolved maps at
C	a given time and position.
C CALLING SEQUENCE:
	subroutine ExtractInsitu(bExtract,VV3D,DD3D,BR3D,BT3D,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	VV3D(nLng,nLat,nRad) real	3-D velocity
C	DD3D(nLng,nLat,nRad) real	3-D normalized density
C	BR3D(nLnG,nLat,nRad) real	3-D magnetic field radial component
C	BT3D(nLnG,nLat,nRad) real	3-D magnetic field tangential component
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C OUTPUTS:
C INCLUDE:
	include		'openfile.h'
	include		'sun.h'
	include		't3d_array.h'
C CALLS:
C	BadR4, FLINT8, Julian, ExtractPosition, bOpenFile, XMAP_OBS_POS
C 	Get3Dval, iFreeLun, XCvarFormat, cTime2System
C PROCEDURE:
C-
	    real	VV3D(*)
	    real	DD3D(*)
	    real	BR3D(*)	
	    real	BT3D(*)
	    integer	nCar
	    real*8	JDCar(nCar)
	
	parameter	(NTP = 12)
	logical		bExtract(NTP)
	character 	cPrefix(NTP)*2	/'me','ve','ea','ma','ju','sa','ur','ne','pl','ul','h1','h2'/

	real*8		JD
	real*8		JEpoch
	real*8		JDCntr
	real*8		FLINT8

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

	character 	cFile*512
		   !write (iU,'(I5,I4,4F7.3,2F10.3)') iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV
	character	cStr*90	/';Year DOY    Hour     Height    Lat     Long   Density Velocity     Brad   Btang'/
	character	cTime2System*80
	logical		bOpenFile

	parameter	(NT = 80)		! 2*NT+1=# times extracted (covers 30 days)
	real*8		dT		/0.25d0/! Time increment for extracted time series

	integer		XCvarFormat

	Bad = BadR4()

	!-------
	! Retrieve grid information

	call T3D_get_grid(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1)

	call T3D_iget(T3D__NCoff  ,0,NCoff)

	call T3D_get (T3D__XCMAT  ,0,XCbeg)	! Needed for statement fnc XCvar
	call T3D_get (T3D__XCMAT+1,0,XCend)

	call T3D_get (T3D__XCROI  ,0,XCbegROI)
	call T3D_get (T3D__TIME	  ,0,XCcntr  )

	JDCntr = FLINT8(1,nCar,JDCar,dble(XCcntr),0d0)! Time when Earth was at center of rotation
	
	do ITP=1,NTP
	    if (bExtract(ITP)) then
						! Heliographic pos of extraction point at JDcntr
		call Julian(1,iYr,Doy,JDCntr,JEpoch)
		call ExtractPosition(ITP,iYr,Doy,RVect)
		JD = JDcntr+(XCcntr-XMAP_OBS_POS(XCcntr,ObsLng))*SUN__SIDEREALP
		I  = JD
		JD = I+nint((JD-I)/dT)*dT

		cFile = cPrefix(ITP)//'_'
		I = XCvarFormat(NCoff,XCcntr,cFile)

		iRecl = 0
		if (bOpenFile(OPN__TEXT+OPN__NEW+OPN__TRYINPUT+OPN__ONEPASS,iU,cFile,iRecl)) then
		    write (iU,'(3A,A19)') '; ',cFile(:itrim(cFile)),' created on ',cTime2System(' ')
		    write (iU,'(A)' ) cStr(:itrim(cStr))
		    XCinc = 0.0

		    do I=-NT,NT
			call Julian(1,iYr,Doy,JD+I*dT,JEpoch)

			call ExtractPosition(ITP,iYr,Doy,RVect)

			ObsD = Bad
			ObsV = Bad

			ObsXC = XMAP_OBS_POS(XCbegROI,ObsLng)
			if (I .gt. -NT) then
			    dObs = ObsXC-ObsXCL
			else
			    dObs = 0.
			end if
			ObsXCL = ObsXC
			if (abs(dObs) .gt. 0.5) XCinc = XCinc + 1.0
			ObsXC = ObsXC + XCinc

			!print *, 'XCbegROI ObsLng ObsLat ObsDis ObsXCL ObsXC', XCbegROI, ObsLng, ObsLat, ObsDis, ObsXCL, ObsXC

			call Get3Dval(XCbeg,XCend,RR,dRR,nLng,nLat,nRad,DD3D,1,ObsXC,ObsLat,ObsDis,ObsD)
			if (ObsD .ne. Bad) ObsD = ObsD/(ObsDis**2)

			call Get3Dval(XCbeg,XCend,RR,dRR,nLng,nLat,nRad,VV3D,1,ObsXC,ObsLat,ObsDis,ObsV)

			call Get3Dval(XCbeg,XCend,RR,dRR,nLng,nLat,nRad,BR3D,1,ObsXC,ObsLat,ObsDis,BR)
			if (BR .ne. Bad) BR = BR/(ObsDis**2)

			call Get3Dval(XCbeg,XCend,RR,dRR,nLng,nLat,nRad,BT3D,1,ObsXC,ObsLat,ObsDis,BT)
			if (BT .ne. Bad) BT = BT/ObsDis

			ObsD = max(ObsD,-9999.999)		! Should fit output format
			ObsV = max(ObsV,-9999.999)
			BR   = max(BR  , -999.999)
			BT   = max(BT  , -999.999)
			write (iU,'(I5,I4,3F9.3,5F9.3)') iYr,int(Doy),(Doy-int(Doy))*24,ObsDis,ObsLat,ObsLng,ObsD,ObsV,BR,BT

		    end do

		    iU = iFreeLun(iU)

		end if
	    end if
	end do

	return
	end