C+
C NAME:
C	CheckMass
C CALLING SEQUENCE:
	subroutine CheckMass(XCbeg,XCend,V3D,D3D)
C INCLUDE:
	include		'math.h'
	include		't3d_array.h'
C	include		't3d_grid_fnc.h'
C	include		't3d_loc_fnc.h'
C-
	    real	XCbeg
	    real	XCend
	    real	V3D(*)
	    real	D3D(*)

	character	cSay*9			/'CheckMass'/
	integer		Str2Str, Flt2Str
	character 	cStr*200

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

	call T3D_iget(T3D__MODE,0,MODE)

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

	XCrange = XCend-XCbeg

	BaseM   = 0.
	BaseMV  = 0.
	BaseMV2 = 0.
	TopM	= 0.
	TopMV	= 0.
	TopMV2	= 0.

	dPhi	= (XCvar(1)-XCvar(2))*(2.*MATH__Pi)	! Longitudinal resolution (radians)
	dTheta	= (XLdeg(2)-XLdeg(1))/2.		! One half of latitudinal resolution (deg)

	dAPole	= 1.-cosd(dTheta)
	dAElse	= 2.*sind(dTheta)

	do J=1,nLat
	    if (J .eq. 1 .or. J .eq. nLat) then		! Poles
		Omega = dAPole
	    else
		Omega = cosd(XLdeg(J))*dAElse
	    end if
	    Omega = Omega*dPhi				! Area/volume per bin

	    do I=1,nLng
		dOmega = Omega				! Divide equally between 0 and 360 degrees
		if (I .eq. 1 .or. I .eq. nLng) dOmega = dOmega/2.

		L = locVOL(I,J,1,1)

		V = V3D(L)
		D = D3D(L)
		BaseM   = BaseM  +D*    dOmega
		BaseMV  = BaseMV +D*V*  dOmega
		BaseMV2 = BaseMV2+D*V*V*dOmega

		L = locVOL(I,J,nRad,1)

		V = V3D(L)
		D = D3D(L)
		TopM   = TopM  +D*    dOmega
		TopMV  = TopMV +D*V*  dOmega
		TopMV2 = TopMV2+D*V*V*dOmega
	    end do
	end do

	I = 0
	I = I+Str2Str('Source surface vs. outer boundary',cStr(I+1:))
	I = I+Str2Str('#Mass ratio',cStr(I+1:))+1
	I = I+Flt2Str(TopM  /BaseM  ,3,cStr(I+1:))
	I = I+Str2Str('#Mass flux ratio',cStr(I+1:))+1
	I = I+Flt2Str(TopMV /BaseMV ,3,cStr(I+1:))
	I = I+Str2Str('#Momentum flux ratio',cStr(I+1:))+1
	I = I+Flt2Str(TopMV2/BaseMV2,3,cStr(I+1:))

	call Say(cSay,'I','Info',cStr)

	return
	end