CACK_OOPS_!@#$%^&*()*&^%$#@!_SPOO_KCAC WRITE3D_BBT.FOR
C+
C NAME:
C	Write3D_bbt
C PURPOSE:
C	Write magnetic field output to bb3d* file
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
	subroutine Write3D_bbt(XCbeGG,iT3D,nT3D, TT3D,XC3D,VV3D,BRss,BR3D,BT3D)
C INPUTS:
C	XCbeGG(2,nTim)			real		Somehow related to the tomography times
C							Only XCbeGG(1,1) and XCbeGG(2,1) are used
C	iT3D				integer		1st time index in TT3D processed
C							TT3D(iT3D) must be a valid time.
C							(the last time index processed is determined
C							internally as the last valid time in TT3D)
C	nT3D				integer		size of time dimension in arrays
C							TT3D,XC3D,VV3D and BRss.
C							(essential to control access to array
C							XC3D because of the extra trailing dimension)
C
C	The following three arrays are set up by Write3D_info.
C
C	TT3D(nT3D)			real		Carrington times
C							These most likely will be the times spaced at
C							regular intervals in UT converted to Carrington times
C	XC3D(nLng,nLat,nRad,nT3D,3)	real		Shift needed to go back to source surface
C							(in Carrington variable). Will be negative.
C							Only XC3D(*,*,*,iT3D:jT3D,*) is used
C	VV3D(nLng,nLat,nRad,nT3D)	real		Velocity array
C							Only VV3D(*,*,*,iT3D:jT3D) is used
C	BRss(nLng,nLat,jT3D-kT3D+1)	real		scratch array; source surface magnetic field
C
C	BR3D(nLng,nLat,nRad)		real		scratch array; radial component written to file
C	BT3D(nLng,nLat,nRad)		real		scratch array; tangential component written to file
C OUTPUTS:
C	(to file)
C INCLUDE:
	include		'sun.h'
	include		't3d_array.h'
	include		'b3d_param.h'

C	include		't3d_grid_fnc.h'
C CALLS:
C	FLINT, T3D_write_bb, T3D_get, T3D_set, T3D_iget, T3D_iset, T3D_get_grid
C	ArrR4Step, ArrR4Constant, ArrR4Interpolate, BadR4, BField_Get, Say
C	ArrR4TimesConstant
C PROCEDURE:
C >	TT3D may be only partially filled with valid times, e.g TT3D(kT3D:jT3D).
C	kT3D and jT3D are determined internally as the 1st and last valid time
C	(.ne. BadR4()) in TT3D. The input value iT3D must be greater/equal kT3D
C	and less/equal jT3D). Output files are produced only for times iT3D to jT3D.
C	kT3D usually is smaller than jT3D. Magnetic source surface data are collected
C	between time indices kT3D and jT3D. This ensures that heliospheric grid elements
C	at early times (just above iT3D) after traceback will still have a valid
C	time associated with them (above kT3D) at the source surface.
C
C >	The radial dependence of Br follows from flux conservation:
C		Br = BB*(RR/r)^2
C	The tangential component follows from geometrical considerations
C	by assuming that the magnetic field lies along the Parker spiral:
C	(Omega is angular velocity of Sun):
C		Bt/Br = r*cos(Lat)*Omega/V
C		Bt = BB*(Omega*RR*cos(Lat)/V)*(RR/r)
C	Output are 'normalized' magnetic field components (making the
C	values independent of radial distance):
C		r^2*Br = RR^2*BB
C		r  *Bt = RR^2*BB *Omega*cos(Lat)/V
C	Expressing radial distances in units of r0 = 1 AU):
C		(r/r0)^2*Br = (RR/r0)^2*BB
C		(r/r0)  *Bt = (RR/r0)^2*BB * r0*Omega*cos(Lat)/V
C	In terms of quantities defined in include file sun.h:
C		r0*Omega/V = Sun__AU*10^13*Sun__Omega*10^-6/(V[km/s]*10^5) = 100*Sun__AU*Sun__Omega/V[km/s]
C
C	We are mostly interested in magnetic field at Earth, i.e. about r0=1 AU.
C	Insitu magnetic field are usually expressed in nT, so we use the same units here.
C	We write 'normalized magnetic fields' in units of nT to file: (RR/r0)^2*Br for the radial
C	and (RR/r0)*Bt for the tangential component. These quantities should be more or less
C	constant across the heliosphere, since the radial dropoff has been taken out. At 1 AU
C	(RR/r0 = 1) it should reproduce the insitu values.
C
C	*TJD*'s unit troubles
C	X's output is in mT or 10^-2G. to bring to nT:
C		(15/215)^2 * (10^-6 T / mT) (10^9 nT/ T)
C
C MODIFICATION HISTORY:
C	SEP-1999, Paul Hick (UCSD/CASS)
C	???-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) Made time dep., fleshed out.
C	JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Valid magnetic fields are returned only in grid points where the
C		source surface Br and V are valid. Previously Br would be valid
C		if only the source surface Br was valid (even if V was bad).
C	AUG-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) 
C		Altered Bt component to depend entirely on V3D. If V3D is bad, use VSS.
C		This revealed bug below: 
C	SEP-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Fixed big error in use of XC3D (did not take into account that the time dimension
C		was different from the V3D array). Now three pieces of information are expected
C		that are calculated by Write3D_info: TT3D, times at which magnetic fields are
C		required; XC3D and V3D: shifts and velocity at times TT3D. XC3D and V3D now have
C		the same dimensions.
C	JUN-2003, Paul Hick (UCSD/CASS)
C		Added processing of multiple magnetic field types
C	OCT-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Added extra arguments iT3D and jT3D
C-
	    real	XCbeGG(*)	! (2,nTim)
	    integer	iT3D
	    integer	nT3D
	    real	TT3D(*)
	    real	XC3D(*)
	    real	VV3D(*)
	    real	BRss(*)
	    real	BR3D(*)
	    real	BT3D(*)

	character	cSay*11			/'Write3D_bbt'/
	real		eps			/0.00002/

	integer		Nss(3)				! Array dimensions at source surface
	real 		Pss(3)				! Location at source surface

	equivalence	(Nss(1),nLng)
	equivalence	(Nss(2),nLat)
	equivalence	(Nss(3),nTim)

	real		Scale(4)		/0.0,1.0,0.0,1.0/
	real		PwrR(2)			/2.0,1.0/

	character	cFile*256
	logical		BField_Get
	integer		iB

	parameter	(nMaxT = 250)

	real		XCbegMAT(nMaxT)
	real		XCendMAT(nMaxT)

	include		't3d_grid_fnc.h'

	!-------
	! loc3D converts from 3D location (iLng,iLat,iRad,iTim,iPos) to linear
	! offset relative to start of array.

	loc3D(iLng,iLat,iRad,iTim,iPos) = ((((iPos-1)*nT3D+iTim-1)*nRad+iRad-1)*nLat+iLat-1)*nLng+iLng

	Bad = BadR4()
	VT = 100.0*Sun__Omega*Sun__AU

	if (TT3D(iT3D) .eq. Bad) then			! Make sure TT3D(iT3D) is valid
	    call Say(cSay,'W','time','at index iT3D is bad. No magnetic field calculated')
	    return
	end if

	kT3D = iT3D-1					! Search for first good time preceding iT3D

	do while (kT3D .ge. 1 .and. TT3D(kT3D) .ne. Bad)
	    kT3D = kT3D-1
	end do
	kT3D = kT3D+1					! Index of earliest good time, or 1

	jT3D = iT3D+1					! Search for last good time after iT3D
	do while (jT3D .le. nT3D .and. TT3D(jT3D) .ne. Bad)
	    jT3D = jT3D+1
	end do
	jT3D = jT3D-1					! Index of last good time, or nT3D

	print *, cSay, kT3D,iT3D,jT3D,nT3D

	!-------
	! Setting the iteration counter to one more than the total number of
	! iterations will suppress adding the iteration number to the output file name.

	call T3D_iget(T3D__ITERATIONS, 0, L  )
	call T3D_iset(T3D__ITERATION , 0, L+1)

	call T3D_set(T3D__SCALE,4,Scale)		! No scale change, but set anyways
	call T3D_set(T3D__R_PWR,2,PwrR )		! Radial normalization powers

	!-------
	! Update T3D structure. Storing nT3D is essential because BField_Get will use this
	! to determine the dimension of the TT3D array.

	call T3D_set (T3D__TT	, 0, 0.0 )
	call T3D_set (T3D__DTT	, 0, 0.0 )

	call T3D_iset(T3D__NTIM	, 0, jT3D-kT3D+1)	! Must stay before T3D_get_grid call !!

	!-------
	! Extract grid information. Note that nTim = jT3D. TT, dTT and dTTi (not used) are zero.

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

	if (nTim .gt. nMaxT) call Say(cSay,'E','nMaxT','parameter too small')

	call ArrR4Constant(nTim,XCbeGG(1),XCbegMAT)
	call ArrR4Constant(nTim,XCbeGG(2),XCendMAT)

	do iB=1,B3D__N

	    !-------
	    ! The magnetic field is returned in nano Tesla (nT).

	    if (BField_Get(iB,TT3D(kT3D),XCbegMAT,XCendMAT,BRss)) then

		!-------
		! Distance normalization: multiply radial component BR3D by RR*RR.
		! The resulting r^2*Br should be in the range 1-10 nT.

		call ArrR4TimesConstant(-nLng*nLat*nTim,BRss,RR*RR,BRss)

		do iTim=iT3D,jT3D

		    TTi = TT3D(iTim)

		    XCbeg   = XCbegMAT(iTim-kT3D+1)		! Needed for statement fnc XCvar
		    XCend   = XCendMAT(iTim-kT3D+1)
		    XCrange = XCend-XCbeg

		    XCb = TTi-0.5 				! Define region of interest (ROI)
		    XCe = TTi+0.5	

		    if (XCb .lt. XCbeg) then			! Keep ROI inside [XCbeg,XCend]
			XCb = XCbeg
			XCe = XCb+1.0
		    end if

		    if (XCe .gt. XCend) then
			XCe = XCend
			XCb = XCe-1.0
		    end if

		    call T3D_set (T3D__TIME   , 0, TTi	)	! Time
		    call T3D_iset(T3D__ITIME  , 0, iTim-kT3D+1)	! Time index
		    call T3D_set (T3D__XCMAT  , 0, XCbeg)	! Begin of matrix
		    call T3D_set (T3D__XCMAT+1, 0, XCend)	! End of matrix
		    call T3D_set (T3D__XCROI  , 0, XCb	)	! Begin ROI
		    call T3D_set (T3D__XCROI+1, 0, XCe	)	! End ROI

		    do iRad=1,nRad

			do iLat=1,nLat

			    XL = cosd(XLdeg(iLat))

			    do iLng=1,nLng

				!-------
				! Trace heliospheric location back to the source surface.
				! There are no bad values in XC3D. Statement functions can be used
				! for longitude and latitude but NOT for the time.

				Pss(1) = XCindx( XCvar(iLng) + XC3D( loc3D(iLng,iLat,iRad,iTim,1) ) )

				!-------
				! We assume that the latitude shift is stored in units of grid spacings.
				! This may be wrong. Currently the latitude shifts are zero so it doesn't
				! matter.

				Pss(2) =         float(iLat) + XC3D( loc3D(iLng,iLat,iRad,iTim,2) )

				!-------
				! To get the time index at the source surface we use the ArrR4Interpolate
				! function. DO NOT USE the statement functions TTindx and TTtime here!
				! Pss(3) will be bad if the traced-back time is outside the range of TT3D.
				! Setting it to -1.0 forces FLINT to return a bad BR value.

				Pss(3) = ArrR4Interpolate(nTim, TT3D(iT3D), TTi+XC3D( loc3D(iLng,iLat,iRad,iTim,3) ))
				if (Pss(3) .eq. Bad) Pss(3) = -1.0

				!------
				! BR is the normalized radial magnetic field at source surface location Pss
				! (interpolate on source surface array BRss)

				BR = FLINT(-3,Nss,BRss,Pss,eps)		! r^2*Br at source surface

C	if (iRad .eq. 10 .and. iLat .eq. nLat/2 .and. 133 .le. iTim .and. iTim .le. 143 )
C     &    print *, iLng, iLat, iTim, Pss(1), Pss(2), Pss(3),
C     &	  XC3D(loc3D(iLng,iLat,iRad,iTim,1)), TTi,BR

				VV = VV3D(loc3D(iLng,iLat,iRad,iTim,1))	! Velocity in heliosphere

				if (BR .ne. Bad .and. VV .ne. Bad) then
				    BT = -BR*(VT*XL/VV)			! r*Bt 
				else
				    BR = Bad				! Erase good BR if VV is bad
				    BT = Bad
				end if

				iPos = loc3D(iLng,iLat,iRad,1,1)
				BR3D(iPos) = BR
				BT3D(iPos) = BT 

			    end do

			end do

		    end do

		    call T3D_write_bb(B3D__PREFIX(iB),t3d, BR3D, BT3D, 0, ' ', cFile)

		end do

	    end if

	end do

	return
	end