C+
C NAME:
C       write3D_bbtm_HR_3
C PURPOSE:
C  Write magnetic field output to bb3d* file from input magnetic source files at the number and cadence of the tomography program.
C  This writes a 3D file and pertinent information over at times XCint(N) + XCdif/ANinterp of the interval and also  
C  allows input from other than the standard model values (heights RRhght(K) + ddR/ANHt, and spatial dimensions  
C  (XCfull(I+1) - XCfull(I))/ANLngp and (XLdeg(J+1) - XLdeg(J))/ANLatp. This assumes that you have dimensioned the 
C  output arrays VV,DD,V3D,D3D AM3D (input in the calling sequence) appropriately. To help make an array that does
C  not get so big, the radial distances nMap of the array can also be an input parameter. (B Jackson, June, 2008)
C
C  XCint(N) values are the beginnings (and ends) of the time intervals used to select data.  When the data selection
C  used is set to a daily cadence with a one day resolution, these XCint(N) values are the beginning and 
C  ending midnight on Earth at the longitude the data is obtained.  This may be different if one site is 
C  used for density and one is used for velocity.  In any case, the inputing XCint(N) values rule in the 
C  selection of the matrix to be output from this routine. The actual start time of the time matrix XCtbeg 
C  is set to the middle of the first set of these values and proceeds at precise nT interval values from this 
C  value until the end value XCtbeg, which is the middle of the last interval. Thus, the time value for the 
C  extraction of the matrix as near as possible should be the middle of the temporal XCint(N) intervals in
C  order to match the Earth location at noon.  The motive here is to select an extract time when most of the 
C  IPS sources are observed.  If temporal intervals are set at a daily cadence, but at higher temporal 
C  resolutions than one day, then the matrix output, should still be set to the middle of the interval. 
C
C CALLING SEQUENCE:
C
C       call write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintD,XCbe,XCtbeg,XCtend,iYr,RR,dRR,
C			nLng,nLat,nMapb,nMapHR,nMap,nT,nTmax,nCar,JDCar,XCbegin,xInc,NC,
C	      		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
C			ClipLng,Scale,VMAP,XCshift,Vratio,
C			BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRss,BRCs,BTCs,BNCs,BR3D,BT3D,BN3D)
C INPUTS:
C	bForecast	logical		If .true. this is forecast mode and output beginning and ending times different from archival are used.
C	Nit		integer 	interation number
C	NiterT		integer		maximum number of iterations
C	NintLng		integer 	# of intermediate longitude steps
C	NintLat		integer 	# of intermediate latitude steps
C	NintHt		integer 	# of intermediate height steps
C	Ninter		integer		# of intermediate time steps
C	NCoff		integer		Carrington variable offset
C	XCintD(nTmax)	real*8		Beginning and ends of the time intervals used to select data
C	XCbe(2,nT)	real		Beginning (1,nT) and end (2,nT) of Carrington rotation times 
C	XCtbeg		real		Start time
C	XCtend		real		End time
C	iYr		integer		year at beginning time
C	RR		real		radial distance of the source surface
C	dRR		real		radial resolution
C
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nMapb		integer		beginning radial distance number value (usually 1)
C	nMapHR		integer		# radial distances output (original values as nMap)
C	nMap		integer		# radial distances
C	nT		integer		# of time values
C	nTmax		integer		maximum # of time values
C	nCar		integer		# Carrington rotations
C	JDCar		real*8		Julian Date at beginning of rotations
C	XCbegin		real		Beginning Carrington rotation value - xInc
C	xInc		real		increment prior to the beginning of the rotation to begin
C	NC		integer		# of Carrington rotations averaged after the beginning
C
C	DEN1AU				real		Density at 1 AU
C	CONRV				real		Velocity spatial smoothing constant
C	CONRD				real		Density spatial smoothing constant
C	CONSTV				real		Velocity spatial hole-filling constant
C	CONSTG				real		Density spatial hole-filling constant
C	CONVT				real		Velocity temporal smoothing constant
C	CONDT				real		Density temporal smoothing constant
C
C	ClipLng				real		Clip longitude
C	Scale				real		extra scaling factor (set to 1 if no scaling is required)
C	VMAP(nLng,nLat)			real		Velocity map at source surface
C       XCshift(nLng,nLat,nMap,nT,3) 	real 		shift amount from array
C       Vratio(nLng,nLat,nMap,nT)  	real 		velocity difference at height
C
C	BFACTORRR			real		A multiplicative factor used to provide a magnetic field enhancement over the CSSS
C									as suggested by T. Hoeksema (Fall AGU 2013)
C	BFACTORRC			real		A multiplicative factor used to provide a closed r field enhancement over the usual
C	BFACTORTC			real		A multiplicative factor used to provide a closed t field enhancement over the usual
C	BFACTORNC			real		A multiplicative factor used to provide a closed n field enhancement
C	iBFlag				integer		The number of the magnetic maps found as in b3d_param_3.h as B3D__PREFIX(B3D__N) 
C	BRss(nLng,nLat,nT)		real		input array; CSSS radial source surface magnetic field
C	BRCs(nLng,nLat,nT)		real		input array; Closed radial source surface magnetic field
C	BTCs(nLng,nLat,nT)		real		input array; Closed tangential source surface magnetic field
C	BNCs(nLng,nLat,nT)		real		input array; Closed normal source surface magnetic field
C	BR3D(nLng*(NintLng+1)+1,nLat*(NintLat+1)+1,nMap*(NintHt+1)+1)	real	scratch array; radial component of written to file
C	BT3D(nLng*(NintLng+1)+1,nLat*(NintLat+1)+1,nMap*(NintHt+1)+1)	real	scratch array; tangential component written to file
C	BN3D(nLng*(NintLng+1)+1,nLat*(NintLat+1)+1,nMap*(NintHt+1)+1)	real	scratch array; normal component written to file
C
C OUTPUTS:
C
C FUNCTIONS/SUBROUTINES:
C       T3D_fill_global, T3D_fill_time, T3D_write_nv
C-
	subroutine write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintD,XCbe,XCtbeg,XCtend,iYr,RR,dRR,
     &				nLng,nLat,nMapb,nMapHR,nMap,nT,nTmax,nCar,JDCar,XCbegin,xInc,NC,
     &	      			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &				ClipLng,VMAP,XCshift,Vratio,
     &				BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRss,BRCs,BTCs,BNCs,BR3D,BT3D,BN3D)

	character	cStrAdd(1)*120,
     &			cFile*120

	logical		bForecast	! If .true. this is forecast mode and output beginning and ending times different from 
					! archival are used.

	real		XCbe(2,nTmax)

	real	VMAP	(nLng,nLat,nT),			! Velocity source surface map
     &		XCshift	(nLng,nLat,nMap,nT,3),		! shift amount from array
     &		Vratio	(nLng,nLat,nMap,nT),		! velocity difference at height
     &		BRss    (nLng,nLat,nT), 		! input array; CSSS radial source surface magnetic field
     &		BRCs    (nLng,nLat,nT),			! input array; Closed radial source surface magnetic field
     &		BTCs    (nLng,nLat,nT),			! input array; Closed tangential source surface magnetic field
     &		BNCs    (nLng,nLat,nT),			! input array; Closed normal source surface magnetic field
     &		BR3D    (nLng*(NintLng+1)-NintLng,nLat*(NintLat+1)-NintLat,nMap*(NintHt+1)-NintHt), ! scratch array; radial component
     &		BT3D    (nLng*(NintLng+1)-NintLng,nLat*(NintLat+1)-NintLat,nMap*(NintHt+1)-NintHt), ! scratch array; tangential component
     &		BN3D    (nLng*(NintLng+1)-NintLng,nLat*(NintLat+1)-NintLat,nMap*(NintHt+1)-NintHt), ! scratch array; normal component
     &		XCsh	(3)
     
	real		Scale(4)		/0.0,1.0,0.0,1.0/
	real		PwrR (3)		/0.0, 2.0, 1.34/
     
	real*8		XCtbeg,XCtend,XCtfull,X,XCtim,XCt,JDCar(nCar)
	real*8		XCintD(nTmax)	

	external	EARTH8

C	include		't3d_array.h'
C	include		'MAPCOORDINATES.H'

	include		't3d_array_3.h'
	include		'b3d_param_3.h'
C	include		't3d_grid_fnc.h'
	include		'MAPCOORDINATES.H'

	print *, ' '
	write(*,'(A,I4,F10.6,3I5)') 'In write3d_bbtm_HR_3', nT,RR,nMapb,nMapHR,iBflag

	if(iBflag.lt.1.or.iBflag.gt.4) then
	  print *, 'There has been no magnetic field of this type data read in. Return to main.'
	  return
	end if
	nLng1   = nLng-1
	PwrR(3) = FALLOFFBN
	Bad = BadR4()

	VT = 100.0*Sun__Omega*Sun__AU

C	PwrR(3) = Bad

	Scale(2) = Scale(2)*BFACTORR
	Scale(4) = Scale(4)*BFACTORR

C	do L=1,nT
C	    do K=1,nMap
C		do J=1,nLat
C		    do I=1, nLng
C			if(L.eq.22.and.J.eq.4) write (*,'(A,2I3,100F12.5)') L,J,VMAP(i,
C		    end do
C		end do
C	    end do
C	end do

C	print *, ' '
C	write (*,'(A,100(I4,E12.5))') 'BRss', (L,BRss(10,4,L),L=1,nt)
C	print *, ' '
	
	XCdif  = sngl(XCtend-XCtbeg)/(sngl((nT-1)*(Ninter+1)))
	XCbegt = sngl(XCtbeg)

	Ninterp  = Ninter + 1	
	ANinterp = Ninterp
	nHtp     = NintHt + 1			! (B Jackson, June, 2008)
	ANHtp    = nHtp
	AHt      = dRR/ANHtp
	nLngp    = NintLng + 1			! (B Jackson, June, 2008)
	ANLngp   = NLngp
	nLatp    = NintLat + 1			! (B Jackson, June, 2008)
	ANLatp   = NLatp	

C	call T3D_fill_global(t3d,0,Nit,NiterT,NCoff,XCbegt,XCdif,RRhght(nMapb),AHt,(nLng-1)*nLngp+1,(nLat-1)*nLatp+1,
C     &	(nMapHR-1)*nHtp+1,nT,PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,ClipLng,Scale)

	call T3D_fill_global(t3d,0,Nit,NiterT,NCoff,XCbegt,XCdif,RRhght(nMapb),AHt,nLng*nLngp-nLngp+1,nLat*nLatp-nLatp+1,
     &	nMapHR*nHtp-nHtp+1,nT,PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,ClipLng,Scale)

     	call T3D_set(T3D__R_PWR,3,PwrR )
	print *, XCtend,XCtbeg,nT,Ninter
	dTT = (XCtend-XCtbeg)/(sngl((nT-1)*(Ninter+1)))
        call T3D_set (T3D__DTT , 0,dTT )
	print *, dTT
	iset =1

C	do iBB=1,4						! there are four possible component maps for each magnetic field 

	XCdif = sngl(XCintD(2)-XCintD(1))			! Change 1/29/2014 BVJ

	AI    = XCdif/ANinterp					! Change 1/29/2014 BVJ

	XCbeg = XCbe(1,1)					! Change 1/28/2014 BVJ
	XCend = XCbe(2,2)					! Change 1/28/2014 BVJ

	if(NC.gt.0) XCending = XCbegin+xInc+float(NC)		! In case more than one rotation is averaged.
	if(NC.eq.0) XCending = XCend

	BBegin = (XCbegin+xInc)-((0.5*AI)/27.2753)		! Archival begining time (one half rotation before center)
	EEnd   = XCending+((0.5*AI)/27.2753)			! Archival ending time (one half rotation after center)
	if(bforecast) then
	  BBegin = (XCbegin+xInc+0.5-20.0/27.2753)		! Forecast begining time (20 days before the run time)
	  EEnd   = (XCbegin+xInc+0.5+(4.0 + AI)/27.2753)	! Forecast ending time (Four days after the run time)
	end if
	call T3D_set (T3D__TT		, 0, BBegin)
	nTTT = int((EEnd - BBegin)/dTT) 
	LLbeg = 0	
	do N=1,nT

C	  XCbeg = XCbe(1,N)					! Change 1/28/2014 BVJ
C	  XCend = XCbe(2,N)					! Change 1/28/2014 BVJ
	  XCdif = sngl(XCintD(N+1)-XCintD(N))
	  AI    = XCdif/ANinterp

	  do L=1,Ninterp

	   XCtim = XCintD(N) + dble(AI*(L-1))			! Times the matrix is to be extracted this N

C	print *, 'Before XMAP',iyr,XCtim,nCAR

	   XCM = XMAP_SC_POS8(EARTH8,iYr,XCtim,nCAR,JDCar)	! Earth location in Carrington coords. at XCtim
	    
	   XCbegROI = XCM - 0.5 				! Region of interest beginning (-NCoff)
	   XCendROI = XCM + 0.5					! Region of interest ending (-NCoff)
	    
C	   print *, 'In write3D XCtim, ROIB, ROIE', XCTim,XCM,XCbegROI,XCendROI

	   TT = XCM
	   LL = Ninterp*N - Ninterp + L
	   
C	   if(TT .gt. (XCbeg+0.75) .and. TT .lt. (XCbeg+2)) then
C	     XCTT(LL) = TT			! Store this value in XCTT to pass out of the subroutine
C	   end if
	   
C	   if(TT .gt. (XCbeg+1-.2) .and. TT .lt. (XCbeg+2)) then

C	   print *, XCbegin+xInc+0.7,XCending+((0.5*AI)/27.2753), NT, N, TT

C	print *, 'to here 0',TT, BBegin,EEnd

	   if(TT .ge. BBegin .and. TT .le. EEnd) then				! Beginning and end of the time loop

C	   print *, 'Calculating 3D file',N,' at XCtim =',XCtim,' or TT =',TT + float(NCoFF)
	   write (*,'(A,I3,A,F8.3,A,F11.5)') 'Calculating 3D Mag file',N,' at XCtim =',XCtim,' or TT =',TT + float(NCoFF)

	   if(LLbeg.eq.0) then
             call T3D_set (T3D__TT ,  0,TT )
	   end if
C	print*, 'To here 1 Write3d_bbtm_HR_3.f'
	   do iBB=1,4						! there are four possible component maps for each magnetic field 	    	  
	      if(iBB.eq.1) then					! First through at any time, make all bad to start
C		do iRad=1,nRad
		do K=nMapb,nMapHR
		  nHtpp = nHtp
	          if(K.eq.(nMapHR)) nHtpp = 1
	          do KL=1,nHtpp
		    iRad = K*nHtp - nHtp + KL
C		    do iLat=1,nLat
		    do J=1,nLat
		      nLatpp = nLatp
		      if(J.eq.(nLat)) nLatpp = 1
		      do JL=1,nLatpp
			iLat = J*nLatp - nLatp + JL
C			do iLng=1,nLng
			do I=1,nLng
			  nLngpp = nLngp
			  if(I.eq.(nLng)) nLngpp = 1
			  do IL=1,nLngpp
			    iLng = I*nLngp - nLngp + IL
				BR3D(iLng,iLat,iRad) = Bad
				BT3D(iLng,iLat,iRad) = Bad
				BN3D(iLng,iLat,iRad) = Bad
			  end do
		        end do
		      end do
		    end do
		  end do
		end do
	      end if
C	print*, 'To here 2 Write3d_bbtm_HR_3.f'
C	print *, 'To here 1', nMapb,nMapHR,nHtp

	    do K=nMapb,nMapHR
	     nHtpp = nHtp
	     if(K.eq.(nMapHR)) nHtpp = 1
	     do KL=1,nHtpp
	      RRht = RRhght(K) + AHt*(KL-1)
	      iRad = K*nHtp - nHtp + KL
C	print *, 'To here 2',K,RRhght(K),AHt,RRHt,iRad,nLat,nLatp
	      do J=1,nLat
	       nLatpp = nLatp
	       if(J.eq.(nLat)) nLatpp = 1
	       do JL=1,nLatpp
	        ALat = (XLdeg(J+1) - XLdeg(J))/ANLatp
		XLat = XLdeg(J)	+ ALat*(JL-1)			! Latitude the matrix is to be extracted
		XL = cosd(XLat)
		iLat = J*nLatp - nLatp + JL
C	if(J.eq.1) print *, 'To here 3', J,ALat,XLat,XL,iLat,nLng,nLngp
	        do I=1,nLng
		 nLngpp = nLngp
		 if(I.eq.(nLng)) nLngpp = 1
		 do IL=1,nLngpp
		  ALng = (XCfull(I+1) - XCfull(I))/ANLngp
		  XClng = XCfull(I) + ALng*(IL-1)		! Longitude the matrix is to be extracted
		  iLng = I*nLngp - nLngp + IL
C	if(I.le.nLng) print *, 'To here 4', I,ALng,XClng,iLng
		  call Get4Dval(3,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT,
     &			XCshift,1,XCtim,XClng,XLat,RRht,XCsh)
C	if(I.le.nLng) print *, 'To here 5', I,ALng,XClng,iLng
		  call Get4Dval(1,XCbe,XCtbeg,XCtend,RR,dRR,nLng,nLat,nMap,nT,
     &			Vratio,1,XCtim,XClng,XLat,RRht,Vrat)

C	print *, 'In write3d_infotd3dM_HR_3 1 after Get4Dvals'
		  
		  XC	= XClng + XCsh(1)		! To get the source surface projected center XCint spatial values	
		  XC    = max(XCbeg,min(XC,XCend))

		  XLa	= XLat + XCsh(2)		! To get the source surface projected heliographic latitude values

		  XCt	= XCtim + dble(XCsh(3))		! To get the source surface projected center XCint time values
		  XCt   = max(XCtbeg,min(XCt,XCtend))

		  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VMAP,1,XC,XLa,XCt,VEL) ! Get velocity value at this height  
		  if(VEL.ne.Bad.and.Vrat.ne.Bad) then
		    VV = Vrat*VEL
		  else
		    VV = Bad
		  end if

C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10) write(*,'(A,4I4,4F12.6)') 'VV ', iBB,iLng,iLat,iRad,XC,XLa,XCT,VV

				if(iBB.eq.1) then
			call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BRss,1,XC,XLa,XCt,BR)	! Get open radial value at this height
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.1) write(*,'(A,4I4,4F12.6)') 'BR ', iBB,iLng,iLat,iRad,XC,XLa,XCT,BR
				if(iBB.eq.2) then
			call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BRCs,1,XC,XLa,XCt,BRC)	! Get closed radial value at this height
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.2) write(*,'(A,4I4,4F12.6)') 'BRC', iBB,iLng,iLat,iRad,XC,XLa,XCT,BRC
				if(iBB.eq.3) then
			call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BTCs,1,XC,XLa,XCt,BTC)	! Get closed tang. value at this height
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.3) write(*,'(A,4I4,4F12.6)') 'BTC', iBB,iLng,iLat,iRad,XC,XLa,XCT,BTC

				if(iBB.eq.4) then

			call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,BNCs,1,XC,XLa,XCt,BNC)	! Get closed norm. value at this height
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.4) write(*,'(A,4I4,4F12.6)') 'BNC', iBB,iLng,iLat,iRad,XC,XLa,XCT,BNC
				if(iBB .eq. 1) then
				  if(BR .ne. Bad) then
				    BR3D(iLng,iLat,iRad) = BR
				  else
				    BR3D(iLng,iLat,iRad) = Bad
				  end if
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.1) write(*,'(A,4I4,4F12.6)') 'BR3D', iBB,iLng,iLat,iRad,XC,XLa,XCT,
C     &		BR3D(iLng,iLat,iRad)
				if(iBB.eq.2) then
				  if (BR3D(iLng,iLat,iRad).eq.Bad) BR3D(iLng,iLat,iRad) = 0.0
				  if (BRC .ne. Bad .and. BR3D(iLng,iLat,iRad) .ne. Bad) then
				    BR3D(iLng,iLat,iRad) = BR3D(iLng,iLat,iRad) + BRC				! r^2*Br
				  else  
				    BR3D(iLng,iLat,iRad) = Bad
				  end if
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.2) write(*,'(A,4I4,4F12.6)') 'BR3D', iBB,iLng,iLat,iRad,XC,XLa,XCT,
C     &		BR3D(iLng,iLat,iRad)
				if(iBB.eq.3) then
				  if (BTC .ne. Bad) then
				    BT3D(iLng,iLat,iRad) = BTC							! r*Btc
				  else
				    BT3D(iLng,iLat,iRad) = Bad
				  end if
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.3) write(*,'(A,4I4,4F12.6)') 'BT3D', iBB,iLng,iLat,iRad,XC,XLa,XCT,
C     &		BT3D(iLng,iLat,iRad)
				if (iBB .eq. 4) then
				  if(BT3D(iLng,iLat,iRad) .eq. Bad .and. BR3D(iLng,iLat,iRad).ne. Bad) BT3D(iLng,iLat,iRad) = 0.0
 				  if(BT3D(iLng,iLat,iRad) .ne. Bad .and. BR3D(iLng,iLat,iRad).eq. Bad) BR3D(iLng,iLat,iRad) = 0.0
				  if(BR3D(iLng,iLat,iRad) .ne. Bad .and. BT3D(iLng,iLat,iRad).ne. Bad .and. VV .ne. Bad) then
				    BT3D(iLng,iLat,iRad) = -BR3D(iLng,iLat,iRad)*(VT*XL/VV) + BT3D(iLng,iLat,iRad)	! r*Bt + r*Btc
				  else
				    BT3D(iLng,iLat,iRad) = Bad
				  end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.3) write(*,'(A,4I4,4F12.6)') 'BT3D', iBB,iLng,iLat,iRad,XC,XLa,XCT,
C     &		BT3D(iLng,iLat,iRad)
				  if(BNC .ne. Bad) then
				    BN3D(iLng,iLat,iRad) = BNC								
				  else
				    BN3D(iLng,iLat,iRad) = Bad
				  end if
				end if
C	if(iLng.eq.22.and.ilat.eq.18.and.irad.eq.10.and.iBB.eq.4) write(*,'(A,4I4,4F12.6)') 'BN3D', iBB,iLng,iLat,iRad,XC,XLa,XCT,
C     &		BN3D(iLng,iLat,iRad)		  
C	print *, 'In write3d_infotd3dm_HR_3 3 after Get3DTvals'

		 end do
	        end do
	       end do
	      end do
	     end do
	    end do

	    if(iBB.eq.4) call T3D_fill_time(t3d,nTTT,TT,XCbeg,XCend,XCbegROI,XCendROI)

C	print *, 'After iBflag, B3D__PREFIX(iBflag)', iBB, iBflag, B3D__PREFIX(iBflag)
C	write(*,'(A,7I8)') 'After', iLng,iLat,iRad,nLng*nLngp-nLngp+1,nLat*nLatp-nLatp+1,nMapHR*nHtp-nHtp+1,nT

	    if(iBB.eq.4) call T3D_write_bb_3(B3D__PREFIX(iBflag),t3d, BR3D, BT3D, BN3D, 0, ' ', cFile)

	   end do

	  end if

	 end do
	end do

      return
      end