C+
C NAME:
C	Write3D_nv_XC
C PURPOSE:
C	Writes output from tomographic reconstruction to ascii t3d* files.
C	One file is written for each time in the tomographic reconstruction.
C	These time are equally separated in Carrington time.
C	Ouput are the 3D velocity and normalized density arrays, and, optionally,
C	the source surface maps, and los crossings maps.
C
C	As a secondary function the input arrays G2S and G3D are converted
C	to normalized density.
C CATEGORY:
C	Tomography: I/O
C CALLING SEQUENCE:
	subroutine Write3D_nv_XC(XCbegMAT,XCendMAT,XCbegROI,XCendROI,VSS,G2S,XC3D,V3D,G3D,BINXV,BINXG)
C INPUTS:
C	(except for the scratch arrays all input is read-only)
C
C	XCbegMAT(nTim)	real		start Carrington variable of matrix (- NCoff)
C	XCendMAT(nTim)	real		end Carrington variable of matrix (- NCoff)
C					(range of matrix typically is three Carrington rotations)
C	XCbegROI(nTim)	real		start Carrington variable of region of interest
C	XCendROI(nTim)	real		end Carrington variable of region of interest
C					(the region of interest typically covers one whole
C					Carrington rotation
C	VSS(nLng,nLat,nTim)
C			real		velocity at source surface RR (with bad values)
C	G2S(nLng,nLat,nTim)
C			real		nr^2 (normalized density) or g^2 at source surface at RR
C					(with bad values)
C	XC3D(3,nLng,nLat,nRad,nTim)
C			real		shift from source surface (in Carrington variable (- NCoff))
C	V3D(nLng,nLat,nRad,nTim)
C			real		velocity ratio or velocity (no bad values)
C	G3D(nLng,nLat,nRad,nTim)
C			real		nr^2 / g^2 ratio or nr^2 / g^2 (no bad values)
C
C	BINXV(nLng,nLat,nTim)	real	# line-of-sight crossings on source surface for velocity
C	BINXG(nLng,nLat,nTim)	real	# line-of-sight crossings on source surface for density
C OUTPUTS:
C	(to t3d file)
C	G2S(nLng,nLat,nTim)
C			real		nr^2 (normalized density) at source surface (with bad values)
C	V3D(nLng,nLat,nRad,nTim)
C			real		velocity matrix (with bad values)
C	G3D(nLng,nLat,nRad,nTim)
C			real		nr^2 (normalized density) matrix (with bad values)
C CALLS:
C	T3D_set, T3D_get, T3D_iset, T3D_get_Grid, Str2Str, FLINT, CvG2D
C	BadR4, WR2DARR, T3D_write_nv
C INCLUDE:
	include		't3d_param.h'
	include		't3d_array.h'
	include		't3d_index.h'

C	include		't3d_grid_fnc.h'
C	include		't3d_loc_fnc.h'
C PROCEDURE:
C	G2S and V2S contain the source surface maps, which may contain bad values.
C	V3D and G3D contain velocity, and g^2 or nr^2, information without bad values.
C	The V3D and G3D matrices are updated using the source surface maps and the
C	shift array XC3D (if TOM__RATIO is not set then the source surface data are used
C	only to flag array elements as bad).
C
C	Several bits in MODE to control several aspects of the calculations.
C
C	TOM__MOD	Carrington variables at the source surface are shifted into
C			the range [XCbeg,XCend] by applying a mod(XCend-XCbeg) operation.
C
C	TOM__G2		If set then input is in g^2 rather than normalized density nr^2
C			If set:
C				G2S is g^2 at source surface (with bad values)
C				G3D is 3D g^2 matrix (if TOM__RATIO not set) or 
C					g^2 ratio (if TOM__RATIO set)
C			If not set:
C				G2S is nr^2 at source surface (with bad values)
C				G3D is 3D nr^2 matrix (if TOM__RATIO not set) or 
C					nr^2 ratio (if TOM__RATIO set)
C
C	TOM__NORATIO	If not set then V3D and G3D are ratios
C			If set:
C				V3D is the 3D velocity matrix
C				G3D is the g^2 matrix (if TOM__G2 set) or the
C					nr^2 matrix (if TOM__G2 not set)
C
C	TOM__BINX	Write los crossings array to t3d file
C	TOM__SHIFT	Write shift arrays XC3D to t3d file
C MODIFICATION HISTORY:
C	SEP-1999, Paul Hick (UCSD/CASS)
C	DEC-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu), major rewrite
C-
	    real	XCbegMAT(*)
	    real	XCendMAT(*)
	    real	XCbegROI(*)
	    real	XCendROI(*)

	    real	VSS	(*)		! Velocity at source surface
	    real	G2S	(*)		! G^2 or nr^2 at source surface
	    real	XC3D	(*)
	    real	V3D	(*)
	    real	G3D	(*)

	    real	BINXV	(*)
	    real	BINXG	(*)

	integer		NDSS(3)
	real		PPSS(3)

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

	character	cFile *120
	character	cFmt2D*7	/'(F10.4)'/
	character	cUser *80
	character	cBinX *80	/'Velocity/normalized density line of sight crossings at source surface'/
	character	cShift*80	/'Shift in Carrington variable (add to go back to source surface)'/

	integer		Str2Str

	logical		bG2
	logical		bRatio
	logical		bBinX
	logical		bShift
	logical		bMod
	logical		bPointP
	logical		bInclV
	logical		bInclG

	real		Tiny		/0.0/

	include		't3d_grid_fnc.h'
	include		't3d_loc_fnc.h'

	call T3D_iget(T3D__MODE,0,MODE)

	!-------
	! Affect calculation

	bMod	= iand(MODE,TOM__MOD	) .ne. 0	! Enforce mod(XCend-XCbeg)
	bG2	= iand(MODE,TOM__G2	) .ne. 0	! Input is in G2 (not normalized density)
	bRatio	= iand(MODE,TOM__NORATIO) .eq. 0	! V3D, G3D are ratios

	!-------
	! Only affect output quantities

	bBinX	= iand(MODE,TOM__BINX	) .ne. 0	! Add los crossings array to t3d file
	bShift	= iand(MODE,TOM__SHIFT	) .ne. 0

	!-------
	! Only used to build user string

	bPointP = iand(MODE,TOM__POINTP	) .ne. 0	! Identifies initial source surface maps (point-P map or not)
	bInclV	= iand(MODE,TOM__INCLV	) .ne. 0
	bInclG	= iand(MODE,TOM__INCLG	) .ne. 0

	!-------
	! Summary of various logical settings is translated into a line
	! for the header of the output t3d* file

	I = 0
	I = I+Str2Str('IPSG2',cUser(I+1:))+2
	if (      bPointP) I = I+Str2Str('POINTP',cUser(I+1:))+4
	if (.not. bPointP) I = I+Str2Str('FLAT'  ,cUser(I+1:))+4

	if (.not. bInclV ) I = I+Str2Str('NO',cUser(I+1:))+1
	I = I+Str2Str('INCLV',cUser(I+1:))+4
	if (.not. bInclG ) I = I+Str2Str('NO',cUser(I+1:))+1
	I = I+Str2Str('INCLG',cUser(I+1:))


	call T3D_set(T3D__SCALE,4,Scale)
	call T3D_set(T3D__R_PWR,2,R_pwr)

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

	call T3D_get(T3D__PWN_V,0,PWN)


	NDSS(1) = nLng
	NDSS(2) = nLat
	NDSS(3) = nTim					! May be degenerate (=1)

	Bad   = BadR4()

	do L=1,nTim

	    !-------
	    ! Fill in time-dependent parameters in t3d array

	    call T3D_iset(T3D__ITIME,0,L )
	    call T3D_set (T3D__TIME ,0,TTtime(L))

	    call T3D_set (T3D__XCMAT  ,0,XCbegMAT(L))
	    call T3D_set (T3D__XCMAT+1,0,XCendMAT(L))

	    call T3D_set (T3D__XCROI  ,0,XCbegROI(L))
	    call T3D_set (T3D__XCROI+1,0,XCendROI(L))

	    do K=1,nRad

		do J=1,nLat

		    do I=1,nLng

			!-------
			! Find the Carrington variable of origin at the source surface
			! of grid point I,J,K,L Note that XCbeg and XCrange are set to
			! the values at time L.

			XCbeg	= XCbegMAT(L)			! Needed for statement fnc XCvar
			XCend	= XCendMAT(L)
			XCrange	= XCend-XCbeg

			iloc = locSHFT(0,I,J,K,L)

			PPSS(1)	= XCvar(I)+XC3D(iloc+PRJ__XC)	! Carrington variable at source surface

			!-------
			! Find the time of origin at the source surface of grid point I,J,K,L,

			PPSS(3)	= TTtime(L)+XC3D(iloc+PRJ__TT)	! Time at source surface
			PPSS(3)	= TTindx(PPSS(3)) 		! Index for time in source surface array

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

			XCbeg   = FLINT(1,nTim,XCbegMAT,PPSS(3),Tiny)! Needed for statement fncs XCindx, XCmod
			XCend	= FLINT(1,nTim,XCendMAT,PPSS(3),Tiny)

			iloc = locVOL(I,J,K,L)

			if (XCbeg .eq. Bad .or. XCend .eq. Bad) then
			    V = Bad
			    D = Bad

			else
			    XCrange = XCend-XCbeg

			    if (bMod) PPSS(1) = XCmod(PPSS(1))	! Map to [XCbeg, XCend] if bMod is set
			    PPSS(1) = XCindx(PPSS(1))		! Index for longitude in source surface array

			    PPSS(2) = float(J)			! Latitude index

			    !------
			    ! Interpolate at longitude PPSS(1), latitude PPSS(2) and time PPSS(3)
			    ! in the source surface

			    V = FLINT(-3,NDSS,VSS,PPSS,Tiny)	! Velocity at source surface
			    D = FLINT(-3,NDSS,G2S,PPSS,Tiny)	! nr^2 or g^2 at source surface

			    if (bRatio) then			! V3D, G3D contain ratios
				if (V .ne. Bad) V = V3D(iloc)*V
				if (D .ne. Bad) D = G3D(iloc)*D
			    else				! V3D, G3D do not contain ratios
				if (V .ne. Bad) V = V3D(iloc)
				if (D .ne. Bad) D = G3D(iloc)
			    end if

			end if

			V3D(iloc) = V				! Velocity
			G3D(iloc) = D				! nr^2 or g^2
		    end do
		end do

	    end do
	end do

	!-------
	! At this point V3D is the 3D velocity matrix. G3D is the nr^2 or g^2 matrix.

	if (bG2) then
	    I = nLng*nLat*nTim
	    call CvG2D(PWN,-I     ,G2S,G2S)			! Convert to nr^2
	    call CvG2D(PWN,-I*nRad,G3D,G3D)
	end if

	!-------
	! Write velocity and nr^2

	do L=1,nTim
	    Lss  = locSS  (  1,1,  L)
	    Lvol = locVOL (  1,1,1,L)
	    Lprj = locSHFT(1,1,1,1,L)

	    call T3D_write_nv('nv3d',t3d,VSS(Lss),G2S(Lss),V3D(Lvol),G3D(Lvol),1,cUser,cFile)

	    if (bBinX) then
		call WR2DARR(2,nLng,nLat,BINXV(Lss),cFile,cFmt2D,.FALSE.,1,cBinX)
		call WR2DARR(2,nLng,nLat,BINXG(Lss),cFile,cFmt2D,.FALSE.,0,cBinX)
	    end if

	    if (bShift) call WR2DARR(2,PRJ__N*nLng,nLat*nRad,XC3D(Lprj),cFile,cFmt2D,.FALSE.,1,cShift)
	end do

	return
	end