C+
C NAME:
C	Write3D_bbt
C PURPOSE:
C	Write magnetic field output to b3d* file
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
	subroutine Write3D_bbt(I__MODE,kIter,nIter,nInterD,NCoff,XCint,XCbegG,XCtbeg,XCtend,RR,dRR,nLng,nLat,nRad,nT,nTmaxV,
     &		XCbegMAT,XCendMAT,XCbegROI,XCendROI, XC3D,V2D,V3D,BR3D,BT3D,
     &		PWNV,PWNG,PWRV,PWRG,N1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONGT,ClipLng,Scale,Btmp,BR2DT,BRcheck)

C INPUTS:
C	I__MODE		integer		dummy
C
C	kIter		integer		iteration number
C	nIter		integer		max. # iterations
C	nInterD		integer         # interpolations. not used yet.
C
C	NCoff		integer		Carrington variable offset
C					This value needs to be added to XCbeg, XCend, XCbegROI,
C					XCendROI and XCTime to get the actual Carrington variable.
C
C	TT		real		start time as Carrington variable (- NCoff)
C	dTT		real		time resolution as Carrington variable
C					(if nT=1 then the values of TT and dTT don't matter)
C
C	RR		real		radial distance of source surface (AU)
C	dRR		real		radial resolution (AU)
C
C	nLng		integer		# longitudes covering [0,360] degrees
C	nLat		integer		# latitudes covering [-90,90] degrees
C	nRad		integer		# radial distance covering [RR,RR+(nRad-1)*dRR]		
C	nT		integer		# times
C
C	XCbegMAT(nT)	real		start Carrington variable of matrix (- NCoff)
C	XCendMAT(nT)	real		end Carrington variable of matrix (- NCoff)
C					(range of matrix typically is three Carrington rotations)
C	XCbegROI(nT)	real		start Carrington variable of region of interest
C	XCendROI(nT)	real		end Carrington variable of region of interest
C					(the region of interest typically covers one whole
C					Carrington rotation
C	V3D(nLng,nLat,nRad,nT)
C			real		3D velocity array
C	XC3D(nLng,nLat,nRad,nT)
C			real		shift from source surface (in Carrington variable)
C	PWNG		real		power determining density dependence of g-level
C	PWR		real		power determining dependenc on radial dependence of g-level
C	N1AU		real		density at 1 AU (cm^-3)
C
C	CONV		real		smoothing width for velocity
C	COND		real		smoothing width for density
C	CONT		real		time smoothing???
C	ClipLng		real		clip longitude (argument for GridSphere)
C
C	BR3D(nLng,nLat,nRad,nT)	real	scratch array
C	BT3D(nLng,nLat,nRad,nT)	real	scratch array
C OUTPUTS:
C	(to file)
C RESTRICTIONS:
C	nRad must be larger/equal 3.
C CALLS:
C	Str2Str, Int2Str, Flt2Str, WR2DARR, CvG2D, GridSphere2D, T3D_fill_global
C	T3D_fill_time
C INCLUDE:
	include		'sun.h'
	include		't3d_param.h'
	include		't3d_array.h'

C	include		't3d_time.h'
C	include		't3d_radial.h'
C	include		't3d_angle.h'
C	include		't3d_init_time.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	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	The Stanford source surface maps at RR = 2.5 Rsun give the radial magnetic field BB in muT.
C	We are mostly interested in magnetic field at Earth, i.e. about r0=1 AU. We expect magnetic
C	fields of order (RR/r0)^2*BB ~ 10^-4*BB ~ 0.1 nT. This extra factor of 10^4 is included in the
C	output i.e. output is in units 0f 0.1 nT.
C	(In situ magnetic fields at Earth are typically in the nT range, at least one order of
C	magnitude larger than predicted from the Stanford maps, probably indicating the the Stanford
C	potential field maps underestimate the magnetic field near the current sheet by that much.)
C
C	*TJD*'s unit troubles
C	When using Xuepu's source surface maps, which are in Gauss,we wished to have an output of 
C	nanotesla. Before, we had ~1mT * (2.5/225)^2 = 10^-4 mT => multiply by scale => .1nT
C       Now, we have units ~1 G * (15/200)^2 = 1 * (3/4)^2 * (1/10)^2 = 1 * .5 * 10^-2 ~= 10^-2 G
C       = 10^-2 G* (10^-4T / G ) (10^9 nT / T) = 10^-6 * 10 ^9 = 1000nT
C	or ~1000nT.  Since output is on order of unity with input (ie, .1) , we multiply it by 
C	10^4 to get it from .1 (1000 nT) to 1 (1 nT). Thus, we leave the scale factor as it is,
C	and it appears to work. 
C
C MODIFICATION HISTORY:
C	SEP-1999, Paul Hick (UCSD; pphick@ucsd.edu)
C-
	real		XCbegMAT(nT),
     &			XCendMAT(nT),
     &			XCbegROI(nT),
     &			XCendROI(nT),
     &			XCbegG(2,nT),
     &			XCint(nT+1)

	real		V3D	(nLng,nLat,nRad,nTmaxV),
     &			XC3D	(nLng,nLat,nRad,nTmaxV),
     &			V2D	(nLng,nLat,nT)

	real		BR3D	(nLng,nLat,nRad,nT),
     &			BT3D	(nLng,nLat,nRad,nT),
     &			BR2DT	(nLng,nLat,nT),
     &			Btmp	(nLng,nLat,3),
     &			TT,		!start time in car. var, without XCoff
     &			dTT,		!time increments, in car variable.
     &			Small		/0.001/,
c    &			nCar		/3.0/,
     &			eps		/0.00002/,
     &			XCtend,
     &			XCtbeg,
     &			BRcheck	(nLng,nLat),
     &			LL

	integer		ND(3),
     &			ND4(4),
     &			nLng,
     &			nLat,
     &			nRad,
     &			nT,
     &			nLng1,
     &			nLat1,
c    &			XCbegR,
     &			NinterD,
     &			nIter,
     &			Ninterp,
     &			ANinterp
     
     
	real		PP4(4)
	real 		PP(3)

	real		Scale(2), N1AU,PWNV,PWNG,PWRV,PWRG,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONGT,ClipLng
	character	cFile*120

	logical		bMod
	logical		T3D_ReadBt


	include		't3d_grid_fnc.h'!takes the place of t3d_time.h,t3d_radial.h,t3d_angle.h

c	include		't3d_init_time.h'

C-------
C Get source surface map. The source surface must be at 2.5 solar radii. If this is
C not the case T3D_ReadB returns .FALSE.

	print *, 'in write bb'

	bMod = iand(I__MODE,TOM__MOD ) .ne. 0
	
	Bad = BadR4()

	Scale(1) = 0.0
	Scale(2) = 1E-4				! Conversion to desired units

	print *, 'nInterD in the beginning',nInterD

C-------
C Interpolations between time points.
c	nInterD  = 0				!We will make nInterD = 0 for no interpolations.
	Ninterp  = nInterD + 1
	ANinterp = Ninterp
	
	print *, 'nInterD later',nInterD
		
	nLng1 = nLng-1
	nLat1 = nLat-1 						
	TT    =  XCbegG(1,1) + XCtbeg 					!TT is start time 
	dTT   = (XCtend-XCtbeg)/(nT-1)					!Time increments
	dTTi  =  1/dTT

	do N=1, nT
	  do L = 1, Ninterp
		XCdiff = XCint(N+1) - XCint(N)
		AI     = XCdiff/ANinterp
		LL     = Ninterp*N + L - Ninterp

c		print *, 'TTint', TTint 
		
		
		TTint = XCint(N) + AI*(L - 0.5)

c		print *, 'TTint', TTint 
		
		XCbegROI(N) = TTint -.5 					!region of interest
		XCendROI(N) = TTint +.5					
		XCbegMAT(N) = XCbegG(1,1)
		XCendMAT(N) = XCbegG(2,1)

		if (XCbegROI(N) .lt. XCbegMAT(N)) then
			XCbegROI(N) = XCbegMAT(N)
 			XCendROI(N) = XCbegMAT(N) + 1.0
 		end if
 
		if (XCendROI(N) .gt. XCendMAT(N)) then
 			XCendROI(N) = XCendMAT(N)
 			XCbegROI(N) = XCendMAT(N) - 1.0
		end if

	
	
		end do
	  end do

	print *, 'XCbegMAT', XCbegMAT
	print *, 'XCendMAT', XCendMAT
	print *, 'XCbegROI', XCbegROI
	print *, 'XCendROI', XCendROI

	call T3D_fill_global(t3d, I__MODE,kIter,nIter,NCoff,TT,dTT,RR,dRR,
     &		nLng,nLat,nRad,nT,PWNV,PWNG,PWRV,PWRG, N1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONGT,ClipLng, Scale)


c*t3d is now filled. readBt fills source surfaces
	if (.not. T3D_ReadBt(t3d,XCbegMAT,XCendMAT,XCbegROI,XCendROI,BR3D,BT3D,Btmp)) return

	ND(1) = nLng
	ND(2) = nLat
	ND(3) = nT

	ND4(1) = nLng
	ND4(2) = nLat
	ND4(3) = nRad
	ND4(4) = nT

	nLngLat    = nLng*nLat
C-------
C For the magnetic field extraction we assume that the tomography source surface is the same
C as the magnetic field source surface (i.e. 2.5 solar radii for the Stanford maps).

	RR2 = RR*RR/Scale(2)
	VT  = 100*SUN__Omega*SUN__AU

c multiply each value in SS matrix BR2DT by RR2 since we will always need to in later equations

	do I=1,nLng
	 do J=1,nLat
	  do N=1,nT
	
	  BR2DT(I,J,N)=RR2*BR3D(I,J,1,N)
c	 if ( BR2DT(I,J,N) .eq. NAN ) print *, 'BR2DT nan', I,J,1,N
c	 if ( BR2DT(I,J,N) .eq. -INF ) print *, 'BR2DT -INF', I,J,1,N
	if ( BR3D(I,J,1,N) .eq. BAD ) BR2DT(I,J,N) = BAD
	  end do
	 end do
	end do
	
	print *, 'now in extrapolation loop, interpolations:', Ninterp
	
		
	do N=1,nT
	
		do L = 1, Ninterp
		XCdiff = XCint(N+1) - XCint(N)
		AI     = XCdiff/ANinterp
		LL     = Ninterp*N + L -Ninterp
		TTint = XCint(N) + AI*(L - 0.5)
	
c	    call T3D_fill_time(t3d, N, TTtime(N), XCbegMAT(N),XCendMAT(N),XCbegROI(N),XCendROI(N))
		call T3D_fill_time(t3d, N, TTint, XCbegMAT(N),XCendMAT(N),TTint-0.5,TTint+0.5)
C------
C Multiply Br at source surface with RR^2.
c because where ever Br ss is used in extrapolation, it is always squared.
c so we just do it now.

c	    call ArrR4TimesConstant(-nLngLat,BR3D(1,1,1,L),RR2,BR3D(1,1,1,L))
	    TTind = TTindx(TTint)	


	    do K=1,nRad
		PP4(3) = float(K)
		do J=1,nLat
		
		    PP(2) = J				 ! Index for latitude
		PP4(2) = float(J)
		    XL = cosd(XLdeg(J))

		    do I=1,nLng	
		    
		    PP4(1) = float(I)	  
 		    PP4(4) = TTind	  
C-------
C Find the Carrington variable of origin at the source surface of grid point I,J,K,L
C Note that XCbeg and XCrange are set to the values at time L.

			XCbeg   = XCbegMAT(N)		! Needed for statement fnc XCvar
			XCend   = XCendMAT(N)

			XCrange = XCend-XCbeg
			VVint = FLINT(-4,ND4,V3D,PP4,eps)
			PP(1)   = XCvar(I)+ FLINT(4,ND4,XC3D,PP4,eps)
							! Carrington variable at source surface
							! XCvar comes from t3d angle
							! depends on XCbeg and XCrange
C-------
C Find the time of origin at the source surface of grid point I,J,K,L,
			PP(3)	= TTint+FLINT(4,ND4,XC3D,PP4,eps)! Time at source surface
							
			PP(3)	= TTindx(PP(3)) 	! Index for time in source surface array

				
C-------
C To convert the Carrington longitude to an array index, the values of XCbeg and XCrange need
C to be set to the proper values at time PP(3) (the time of origin at the source surface).
C This is done by interpolation on the XCbegMAT and XCendMAT arrays.

			XCbeg = XCbegMAT(1)
			XCend = XCendMAT(1)! Needed for statement fncs XCindx, XCmod


			if (XCbeg .eq. Bad .or. XCend .eq. Bad) then
			    BR = Bad
			    VV = Bad
			    
c			    print *, 'XCbeg or XCend Bad!'
c			    print *, 'I J K N', I, J, K, N
			else
			    XCrange = XCend-XCbeg
			    if (bMod) then
			    
			    PP(1) = XCmod(PP(1))	! Map to [XCbeg, XCend] if bMod is set
			    
	   	       end if

			    
			    PP(1) = XCindx(PP(1))		! Index for longitude in source surface array
								! ...	PP(1)	= nLng - (PP(1) - XCbeg)/(XCrange - nLng1)
C------
C The velocity at the source surface is calculated by 4-dim linear interpolation at
C longitude P(1), latitude P(2)=J, and time P(3)
c			

c*TJD knocked down a dimension 4 -> 3
c			    BR = FLINT(-4,ND,BR3D,PP) 		! r^2*Br at source surface (ie, we divided out by r^2 for physics then normalized by * r^2 computational reasons 
C			    VV = FLINT(-4,ND,V3D ,PP)		! Velocity at source surface

	
				BR = FLINT(-3,ND,BR2DT,PP,eps) 		! r^2*Br at source surface (ie, we divided out by r^2 for physics then normalized by * r^2 computational reasons 
				if (BR .eq. NaN) BR = Bad
			    	VV = FLINT(-3,ND,V2D,PP,eps)		! Velocity at source surface
				if (VV .eq. NaN) VV = Bad

			    if (BR .ne. Bad .and. VV .ne. Bad) then
			    	
			    
c				note, its OK if	V3D(I,J,K,N) is occaisionally bad,
c				since one expected VMAPD to have a few holes

				VV = 0.5*(VV+VVint)											
				if (VV .eq. 0.0) VV = Small	!tjd got rid of zero-divisor
				VV = -BR*(VT*XL/VV)		! r*Bt ...did we loose an 1/RR
				if (K .eq. 1) VV = Small
				
			    end if

			end if	
	
c			if (BR .le. 10E-6) BR = BAD	
c			if (BR .ge. 10E+6) BR = BAD		
			
c			if (VV .le. 10E-6) BR = BAD				
c			if (VV .ge. 10E+6) BR = BAD			
					
			BR3D(I,J,K,N) = BR			
			BT3D(I,J,K,N) = VV
		    end do
		end do
		end do

c*tjd* check some values, print		

c	print *, 'BRs and BTs going into T3D write bb for N=',N

c	print *, 'after extrapolation:'
c	print *, 'Source'

c	print *, 'BRtD(20,5,1,',N,')',	BT3D(20,5,1,N)
c	print *, 'BRtD(20,8,1,',N,')',	BT3D(20,8,1,N)
c	print *, 'BRtD(35,3,1,',N,')',	BT3D(35,3,1,N)
c	print *, 'BRtD(35,5,1,',N,')',	BT3D(35,5,1,N)
c	print *, 'BRtD(35,8,1,',N,')',	BT3D(35,8,1,N)

c	print *, 'BRtD(20,3,10,1)',	BT3D(20,3,10,N)
c	print *, 'BRtD(20,5,10,1)',	BT3D(20,5,10,N)
c	print *, 'BRtD(20,8,10,1)',	BT3D(20,8,10,N)
c	print *, 'BRtD(35,3,10,1)',	BT3D(35,3,10,N)
c	print *, 'BRtD(35,5,10,1)',	BT3D(35,5,10,N)
c	print *, 'BRtD(35,8,10,1)',	BT3D(35,8,10,N)
	
c	print *, 'BRrD(1,1,1,',N,')',	BR3D(1,1,1,N)
c	print *, 'BRrD(2,1,1,',N,')',	BR3D(2,1,1,N)
c	print *, 'BRrD(3,1,1,',N,')',	BR3D(3,1,1,N)
c	print *, 'BRrD(1,2,1,',N,')',	BR3D(1,2,1,N)
c	print *, 'BRrD(1,3,1,',N,')',	BR3D(1,3,1,N)
c	print *, 'BRrD(20,3,1,',N,')',	BR3D(20,3,1,N)
c	print *, 'BRrD(20,5,1,',N,')',	BR3D(20,5,1,N)
c	print *, 'BRrD(20,7,1,',N,')',	BR3D(20,7,1,N)
c	print *, 'BRrD(20,9,1,',N,')',	BR3D(20,9,1,N)
c	print *, 'BRrD(20,10,1,',N,')',	BR3D(20,10,1,N)	
c	print *, 'BRrD(19,8,1,',N,')',	BR3D(19,8,1,N)
c	print *, 'BRrD(20,8,1,',N,')',	BR3D(20,8,1,N)
c	print *, 'BRrD(21,5,1,',N,')',	BR3D(19,8,1,N)
c	print *, 'BRrD(35,3,1,',N,')',	BR3D(35,3,1,N)
c	print *, 'BRrD(35,4,1,',N,')',	BR3D(35,4,1,N)
c	print *, 'BRrD(35,5,1,',N,')',	BR3D(35,5,1,N)
c	print *, 'BRrD(35,6,1,1)',	BR3D(35,6,1,N)
c	print *, 'BRrD(35,7,1,1)',	BR3D(35,7,1,N)
c	print *, 'BRrD(35,9,1,',N,')',	BR3D(35,9,1,N)
c	print *, 'BRrD(35,10,1,',N,')',	BR3D(35,10,1,N)	
c	print *, 'BRrD(34,3,1,1)',	BR3D(34,3,1,N)
c	print *, 'BRrD(34,4,1,1)',	BR3D(34,4,1,N)
c	print *, 'BRrD(34,5,1,1)',	BR3D(34,5,1,N)
c	print *, 'BRrD(34,6,1,1)',	BR3D(34,6,1,N)
c	print *, 'BRrD(34,7,1,1)',	BR3D(34,7,1,N)
c	print *, 'BRrD(34,9,1,1)',	BR3D(34,9,1,N)
c	print *, 'BRrD(34,10,1,1)',	BR3D(34,10,1,N)
c	print *, 'BRrD(34,8,1,1)',	BR3D(34,8,1,N)	
c	print *, 'BRrD(35,8,1,1)',	BR3D(35,8,1,N)
c	print *, 'BRrD(36,8,1,1)',	BR3D(36,8,1,N)
		
c	print *, 'Xtrap:'

c	print *, 'BRrD(19,8,1,',N,')',	BR3D(19,8,1,N)
c	print *, 'BRrD(19,9,1,',N,')',	BR3D(19,9,1,N)
c	print *, 'BRrD(19,10,1,',N,')',	BR3D(19,10,1,N)
c	print *, 'BRrD(19,8,2,',N,')',	BR3D(19,8,2,N)
c	print *, 'BRrD(19,9,2,',N,')',	BR3D(19,9,2,N)
c	print *, 'BRrD(19,10,2,',N,')',	BR3D(19,10,2,N)
c	print *, 'BRrD(19,8,3,',N,')',	BR3D(19,8,3,N)
c	print *, 'BRrD(19,9,3,',N,')',	BR3D(19,9,3,N)
c	print *, 'BRrD(19,10,3,',N,')',	BR3D(19,10,3,N)
c


c	print *, 'BRrD(1,1,10,',N,')',	BR3D(1,1,10,N)
c	print *, 'BRrD(2,1,10,',N,')',	BR3D(2,1,10,N)
c	print *, 'BRrD(1,2,10,',N,')',	BR3D(1,2,10,N)
c	print *, 'BRrD(1,3,10,',N,')',	BR3D(1,3,10,N)
c	print *, 'BRrD(20,3,10,',N,')',	BR3D(20,3,10,N)
c	print *, 'BRrD(20,5,10,',N,')',	BR3D(20,5,10,N)
cc	print *, 'BRrD(20,7,10,',N,')',	BR3D(20,7,10,N)
c	print *, 'BRrD(20,9,10,',N,')',	BR3D(20,9,10,N)
c	print *, 'BRrD(20,10,10,1)',	BR3D(20,10,10,N)	
c	print *, 'BRrD(19,8,10,1)',	BR3D(19,8,10,N)
c	print *, 'BRrD(20,8,10,1)',	BR3D(20,8,10,N)
c	print *, 'BRrD(21,5,10,1)',	BR3D(19,8,10,N)
c	print *, 'BRrD(35,3,10,1)',	BR3D(35,3,10,N)
c	print *, 'BRrD(35,4,10,1)',	BR3D(35,4,10,N)
c	print *, 'BRrD(35,5,10,1)',	BR3D(35,5,10,N)
c	print *, 'BRrD(35,6,10,1)',	BR3D(35,6,10,N)
c	print *, 'BRrD(35,7,10,1)',	BR3D(35,7,10,N)
c	print *, 'BRrD(35,9,10,1)',	BR3D(35,9,10,N)
c	print *, 'BRrD(35,10,10,1)',	BR3D(35,10,10,N)
c	print *, 'BRrD(34,3,10,1)',	BR3D(34,3,10,N)
c	print *, 'BRrD(34,4,10,1)',	BR3D(34,4,10,N)
c	print *, 'BRrD(34,5,10,1)',	BR3D(34,5,10,N)
c	print *, 'BRrD(34,6,10,1)',	BR3D(34,6,10,N)
c	print *, 'BRrD(34,7,10,1)',	BR3D(34,7,10,N)
c	print *, 'BRrD(34,9,10,1)',	BR3D(34,9,10,N)
c	print *, 'BRrD(34,10,10,1)',	BR3D(34,9,10,N)
c	print *, 'BRrD(34,8,10,1)',	BR3D(34,8,10,N)	
c	print *, 'BRrD(35,8,10,1)',	BR3D(35,8,10,N)
c	print *, 'BRrD(36,8,10,1)',	BR3D(36,8,10,N)
c*should do nothing more than write the file and clean out bad values.

c
c	print *, 'N =', N
c        do I=1, nLng
c         do J=1, nLat
c  		BRcheck(nLng,nLat) = BR3D(I,J,25,N)
c	 end do
c	end do
	
c	call WhatisR4(nLat*nLng, BRcheck,'BR3D check ar rad 25 info as follows')
c	print *, 'writing cFile=',cFile	
      if (N .ge. 10) then
       if (N .le. 46)  call T3D_write_bb('bb3d',t3d, BR3D(1,1,1,N), BT3D(1,1,1,N), 0, ' ', cFile)
      end if	
	print *, 'written out file at N=',N

	end do
    	end do

    	return
    	end