C+
C NAME:
C	Write3D_bb
C PURPOSE:
C	Write magnetic field output to bb3d* file
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
	subroutine Write3D_bb(XC3D,V3D,BR3D,BT3D)
C INPUTS:
C	cWild		character	wild card specification for locating magnetic field
C					files in the format processed by href=FileSelection=.
C	V3D (nLng,nLat,nRad)	real	3D velocity array
C					(contains bad values)
C	XC3D(nLng,nLat,nRad)	real	shift from source surface (in Carrington variable)
C					(does not contain any bad values)
C
C	(additional inputs, such as the dimensions nLng,nLat,nRad are pulled
C	out of the internal t3d array)
C OUTPUTS:
C	BR3D(nLng,nLat,nRad)	real	scratch array (content does not matter)
C	BT3D(nLng,nLat,nRad)	real	scratch array (content does not matter)
C	(these magnetic field components are also written to file)
C RESTRICTIONS:
C	nRad must be larger/equal 3.		(??? Why ???)
C CALLS:
C	T3D_get_grid, T3D_get, T3D_set, ArrR4Bad, T3D_Read_B, Say
C	ArrR4TimesConstant, FLINT, T3D_write_bb, T3D_get_mode
C INCLUDE:
	include		'sun.h'
	include		't3d_array.h'
	include		'b3d_param.h'

C	include		't3d_grid_fnc.h'
C PROCEDURE:
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 >	We use 'normalized' magnetic field components, which are as much as
C	possible independent of radial distance:
C		r^2*Br = RR^2*BB
C		r  *Bt = RR^2*BB *Omega*cos(Lat)/V
C	Since we are mostly interested in magnetic field at Earth, we express
C	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	These are the quantities written to file and returned in the output
C	arrays. The magnetic units are nT (the same units are usually used
C	for insitu values at 1 AU).
C >	In terms of the quantities defined in include file sun.h:
C		r0*Omega/V = SUN__AU*10^13*SUN__OMEGA*10^-6/(V[km/s]*10^5)
C			   = 100*SUN__AU*SUN__OMEGA/V[km/s]
C
C >	For each grid point XC3D is used to trace the point back to the source surface.
C	If both the velocity V3D and the radial magnetic field BR3D at the source
C	surface are valid then valid BR3D and BT3D in the grid point are returned.
C	The velocity needed for BT3D is taken as the average between the source
C	surface velocity and the velocity in the grid point (if it is valid).
C
C >	*TJD*'s unit reminders
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 MODIFICATION HISTORY:
C	SEP-1999, Paul Hick (UCSD/CASS)
C	NOV-2001, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu)
C	JUL-2002, Paul Hick (UCSD/CASS)
C		Valid magnetic fields are returned only in grid points where the
C		source surface Br and V are both valid. Previously Br would be valid 
C		if only the source surface Br was valid (even if V was bad).
C	AUG-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Added check for type of source surface data (WSO or WSO_NOAA maps)
C	AUG-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) 
C		Altered Bt component to depend entirely on V3D. If V3D doesnt exist, uses VV.
C-
	    real	V3D (*)
	    real	XC3D(*)
	    real	BR3D(*)
	    real	BT3D(*)

	character	cSay*10		/'Write3D_bb'/
	real		Tiny		/2.0e-5/
	real		Scale(4)	/0.0,1.0,0.0,1.0/
	real		PwrR(2)		/2.0,1.0/

	character	cFile*512
	logical		BField_Get

	include		't3d_grid_fnc.h'

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

	locP(iLng,iLat,iRad) = ((iRad-1)*nLat+iLat-1)*nLng+iLng

	Bad = BadR4()

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

	call T3D_get_grid(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1)
	nMap = nLng*nLat			! # bins in map at fixed distance
	nVol = nMap*nRad			! # bins in whole volume

	call T3D_get(T3D__XCMAT  ,0,XCbeg)	! Next 3 lines: needed for statement fnc XCvar
	call T3D_get(T3D__XCMAT+1,0,XCend)
	XCrange = XCend-XCbeg

	VT = 100.0*SUN__OMEGA*SUN__AU		! Rotation velocity at solar surface
	call T3D_set(T3D__SCALE,4,Scale)	! No scale change, but set anyways
	call T3D_set(T3D__R_PWR,2,PwrR )	! Radial normalization powers


	do iB=1,B3D__N				! Loop over all types of source surface maps

	    call ArrR4Bad(nVol, BR3D)		! Initialize magnetic field to bad
	    call ArrR4Bad(nVol, BT3D)

	    !-------
	    ! The magnetic field is returned in nano Tesla (nT).
	    ! If source surface magnetic field found, normalize by multiplying
	    ! the radial component BR3D by RR*RR.

	    if (BField_Get(iB,TT,XCbeg,XCend,BR3D)) then

		call ArrR4TimesConstant(-nMap,BR3D,RR*RR,BR3D) ! Change to r^2*Br

		do iRad=1,nRad

		    do iLat=1,nLat

			XL = cosd(XLdeg(iLat))

			do iLng=1,nLng	
 
			    iP  = locP(iLng,iLat,iRad)	! Location grid point in 3D array
			    iP0 = locP(1   ,iLat,   1)	! Location of start of row at source surface

			    !-------
			    ! Find the Carrington variable of origin at the source surface of grid point

			    PP = XCvar(iLng)+XC3D(iP)	! Carrington variable at source surface
			    PP = XCindx(PP)		! Index for longitude at source surface

			    !------
			    ! The magnetic field at the source surface is calculated by 1-dim
			    ! linear interpolation across rows

			    BR = FLINT(-1,nLng,BR3D(iP0),PP,Tiny)! r^2*BR at source surface
			    VV = V3D(iP)		! Velocity in grid point

			    if (BR .ne. Bad .and. VV .ne. Bad) then
				BT = -BR*(VT*XL/VV)	! r^BT
			    else
				BR = Bad
				BT = Bad
			    end if

			    BR3D(iP) = BR
			    BT3D(iP) = BT

			end do

		    end do

		end do

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

	    end if

	end do

	return
	end