C+
C NAME:
C	ipsg2t
C PURPOSE:
C	Deconvolve IPS observations
C CATEGORY:
	program ipsg2t
C CALLING SEQUENCE:
C	run $exe:ipsgt
C INPUTS:
C	IPS data files
C OUTPUTS:
C	t3d* files in current directory
C CALLS:
C	Say, itrim, LocFirst, ForeignArg, ForeignFile, ReadVIPSValidData, ReadVIPS
C	ReadGIPS, AskYN, AskR4
C	SetGrid, IPSProcessing, IPSConstraints, T3D_set, T3D_iset, T3D_iget, T3D_show,
C	T3D_get_grid, LOSPosition, LOSProjection, LOSReach, LOSClean, LOSWeights
C	LOSIntegralG, LOSIntegralV, LOSTweak, Str2Str, Flt2Str, ArrR4TimesArrR4
C	ArrR4Constant, ArrR4Copy, ArrR4Function, ArrR4Mask, MapGrid, SanityCheck
C	CvG2D, CvD2G, MkD2V, MkV2D, iArrI4Total
C	SW_Model_Kinematic, BuildSourceSurface, GridSphere2D, GridFill, Write3D_nv, Write3D_bb
C	Write3D_nv_UT, Write3D_bb_UT, WriteLOSY
C	WriteLOSD, StopWatch, T3D_marker
C	iReadNagoya, iProcessNagoya, iReadOoty, iProcessOoty
C INCLUDE:
	include		'dirspec.h'
	include		't3d_param.h'
	include		't3d_array.h'
	include		't3d_index.h'
C EXTERNAL:
	external	iReadNagoya
	external	iProcessNagoya		! Passed to ReadVIPS
	external	iReadOoty
	external	iProcessOoty
C PROCEDURE:
C > !!!	Normalized densities are used throughout.
C	These are defined as n*r^2 with n the density in cm^-3 and r the heliocentric
C	distance in AU. Typical values would be in the range 1-10.
C > !!!	Velocities are in km/s
C
C >	The modified Julian day MJD used for the IPS data is related to the
C	regular Julian day (as used in subroutine Julian) by
C		MJD = JD-2400000.5
C	(i.e. the modified Julian day starts at midnight).
C >	Once the logical $CAM is set (during the first run of IPSD), the
C	only way to force IPSD to look in another directory is to delete it.
C
C	II = iand(I__MODE,TOM__MOD):
C
C	LOSReach:
C	uses I__MODE to stop execution when LOS segments fall outside
C	the range [XCbegMAT,XCendMAT] (only if II = 0).
C
C	SW_Model_Kinematic:
C	uses I__MODE three times to deal with mod(XCendMAT-XCbegMAT) issues (if II != 0)
C	and twice to avoid array index errors (if II = 0).
C	The 2nd of these latter is the only place in the NOMOD scheme where a
C	rather ugly kludge is needed.
C
C	LOSIntegralG, LOSIntegralV
C	passes I__MODE to Get3DValue
C
C	LOSProjection
C	passes I__MODE to Get3DValue
C
C	BuildSourceSurface
C	uses I__MODE to map XC values into the range [XCbegMAT,XCendMAT] (II != 0 only).
C
C	POTENTIAL MODIFICATIONS/IMPROVEMENTS:
C
C	The first and last grid point in longitude and latitude now are located
C	on the edges of the interval covered. It makes more sense to put the grid
C	at bin centers instead. E.g. at a 10 degree resolution the latitude grid
C	now is located at -90,-80,...,80,90. It should be located at -85,-75,75,85.
C	This is probably fairly easy to fix in the Fortran tomography program.
C	It may be more difficult to adapt the IDL visualization programs.
C
C	Switch all Carrington variables (time and longitude) to double precision
C MODIFICATION HISTORY:
C	JUN-1992, Paul Hick (UCSD)
C		point-P synoptic map program
C	NOV-1995, Bernie Jackson (STEL)
C		conversion to tomography program
C	APR-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		subtantial rewrite. Program now uses g^2 rather than density.
C-

	!-------
	! Line of sight resolution

	parameter	(NLG0	= 12000)	! Max # G data points
	parameter	(NLV0	=  3000)	! Max # V data points

	parameter	(NLOS0  =   40)		! # resolution elements along each line of sight
	parameter	(dLOS0	= 0.05)		! Resolution along los in AU

	parameter	(NLmax0	= NLG0)		! Use the larger of NLG0 and NLV0 to define NLmax0
	parameter	(NLTot0	= NLOS0*NLmax0)	! Max # los elements

	!-------
	! We set up the grid size separately here. They are used to declare
	! array sizes, and are also the default grid sizes used. The grid size
	! can be modified by specifying the appropriate command line argument.
	! As long as the total # of grid points is not bigger than nSST0 (at
	! the source surface) and n3DT0 (in the 3D volume) everything should
	! work.

	parameter	(nLng0	= 37)		! # Longitudes, resolution is 360/(nLng-1)

	parameter	(nLat0	= 19)		! # Latitudes , resolution is 180/(nLat-1)
	parameter	(nRad0	= 31)		! # Radial heights. Range covered is [0,(nRad-1)*dRR]
	parameter	(nTim0	=  1)
						! Max # long/lat elements in source surface map
	parameter	(nSS0	= (3*(nLng0-1)+1)*nLat0)
	parameter	(n3D0	= nSS0 *nRad0)	! Max # long/lat/rad elements
	parameter	(nSST0	= nSS0 *nTim0)	! Max # long/lat/tim elements in source surface map
	parameter	(n3DT0	= nSST0*nRad0)	! Max # long/lat/rad/tim elements in whole volume

	real		XCbegMAT(nTim0)
	real		XCendMAT(nTim0)
	real		XCbegROI(nTim0)
	real		XCendROI(nTim0)

	!-------
	! Arrays for IPS velocity observations (max. of NLmaxV los observations)
	! The first two array are needed to retrieve the original data from the
	! input file and write them to a separate file with the model IPS data
	! attached.

	integer		IYRFV	(NLV0)		! Year of file where los observation is located
	integer		IRECV	(NLV0)		! Record number in file

	integer		IYRSV	(NLV0)		! Year of los observation
	integer		NBSV	(NLV0)		! Initialized by LOSReach

	! IPS V data arrays

	real		DOYSV	(NLV0)		! Day of year (incl. fraction for time of day) of observation
	real		DISTV	(NLV0)		! Sun-Earth distance at time of obs.
	real		XLSV	(NLV0)		! Geocentric ecliptic longitude of Sun
	real		XLLV	(NLV0)		! Geocentric ecliptic longitude diff.lng(los)-lng(Sun)
	real		XDLV	(NLV0)		! Geocentric ecliptic latitude los
	real		XCEV	(NLV0)		! Carrington variable of sub-Earth point
	real		XEV	(NLV0)		! Los elongation (deg) (>0: East of Sun; <0: West of Sun)
	real		XCV	(NLV0)		! Carrington variable of point-P traced back to source surface at speed VOBS
	real		YLV	(NLV0)		! Heliographic latitude of point-P (deg)
	real		VOBS	(NLV0)		! Observed IPS velocities

	real		VMDL	(NLV0)		! Model IPS velocities at end
	real		VMDL0	(NLV0)		! Model IPS velocities after iteration 0
	real		PPlosV	(LOS__N,NLOS0,NLV0)
	real		VPlosV	(NLOS0,NLV0)	! sin(Chi) = Vperp/Vrad
	real		WWlosV	(NLOS0,NLV0)

	!-------
	! Arrays for IPS G-level observations (max. of NLmaxG los observations)
	! The first two array are needed to retrieve the original data from the
	! input file and write them to a separate file with the model IPS data
	! attached.

	integer		iXPG	(NLG0)		! Array indices from daily G maps
	integer		iYPG	(NLG0)
	integer		iMJDG	(NLG0)		! Original MJD of daily G maps

	integer		IYRFG	(NLG0)		! Year of yearly IPS file
	integer		IRECG	(NLG0)		! Record # in yearly file

	equivalence	(iMJDG, IYRFG)
	equivalence	(iXPG , IRECG)

	! IPS G data arrays

	integer		IYRSG	(NLG0)		! Year of source observation
	integer		NBSG	(NLG0)		! Initialized by LOSReach

	real		DOYSG	(NLG0)		! Day of year (incl. fraction for day of year) of obs.
	real		DISTG	(NLG0)		! Sun-Earth distance at time of obs.
	real		XLSG	(NLG0)		! Geocentric ecliptic longitude of Sun
	real		XLLG	(NLG0)		! Geocentric ecliptic longitude diff.lng(los)-lng(Sun)
	real		XDLG	(NLG0)		! Geocentric ecliptic latitude los
	real		XCEG	(NLG0)		! Carrington variable of sub-Earth point
	real		XEG	(NLG0)		! Los elongation (deg) (>0: East of Sun; <0: West of Sun)
	real		XCG	(NLG0)		! Carrington variable of point-P traced back to source surface at speed V0
	real		YLG	(NLG0)		! Heliographic latitude of point-P (deg)
	real		GOBS	(NLG0)		! Observed IPS g-levels
	real		G2OBS	(NLG0)		! Observed IPS g^2-levels

	real		G2MDL	(NLG0)		! Model IPS G-levels at end
	real		G2MDL0	(NLG0)		! Model IPS G-levels after iteration 0
	real		PPlosG	(LOS__N,NLOS0,NLG0)
	real		WWlosG	(NLOS0,NLG0)

	real		M2Mean	(NLG0)

	!-------
	! Scratch space for line of sight manipulations.
	! Note that some of the arrays are equivalenced, so be careful!!

	integer		iLOSTmp	(NLTot0,3)
	real		rLOSTmp	(NLTot0,3)
	equivalence	(rLOSTmp, iLOSTmp)

	real		GMDL	(NLmax0)
	real		VSUM	(NLmax0)
     	equivalence	(GMDL,VSUM)

	! Arrays shared by IPS V and G

	real		PPprj	(PRJ__N,NLTot0)	! Lng, lat, time of projected point on los
	real		Weight	(NLTot0)
	real		FIX	(NLmax0)	! Fix to model to make OK

	!-------
	! Source surface arrays (includes scratch arrays)

	real		VMap	(nSST0)		! Source surface velocity map
	real		VFin	(nSST0)
	real		GMap	(nSST0)		! Source surface density map
	real		GFin	(nSST0)

	real		BINXV	(nSST0)		! # los crossings in V (scratch array)
	real		BINXG	(nSST0)		! # los crossings in G (scratch array)

	real		TmpSST	(nSST0*3)	! Scratch array for SW_Model_Kinematic 
	real		TmpSST1 (nSST0)
	real		TmpSST2 (nSST0)

	equivalence	(TmpSST1, TmpSST(      1))
	equivalence	(TmpSST2, TmpSST(nSST0+1))

	real		Gtmp	(nSS0,2)	! Scratch arrays for SW_Model_Kinematic, BuildSourceSurface
	real		Vtmp	(nSS0,2)

 	byte		Btmp	(nSS0,9)	! Scratch array
 	equivalence	(Btmp,Vtmp)


	!-------
	! Volume array (includes scratch arrays)

	real		XC3D	(n3DT0*3)	! 3D 'traceback' array
	real		V3D	(n3DT0)		! Velocities
	real		G3D	(n3DT0)		! g^2
	real		Tmp3D1	(n3D0 )		! Scratch array
	real		Tmp3D2	(n3D0 )		! Scratch array



	logical		bGcon
	logical		bVcon
	logical		bVNago
	logical		bVOoty
	logical		bGNago
	logical		bGOoty
	logical		bGCamb
	logical		bMV2				! Use mv^2 = constant
	logical		bBoundV		/.TRUE./	! Limit V-map boundaries?
	logical		bBoundG		/.TRUE./	! Limit D-map boundaries?
	logical		bForecast
	logical		bRemove
	logical		bPointP
	logical		bMod
	logical		bDiff
	logical		bTrackG
 
	integer		iEorWG		/ 0/
	integer		iEorWV		/ 0/
	integer		NAV		/ 0/

	integer		NitTotal	/18/
	integer		NitBadV		/-1/
	integer		NitBadG		/-1/
	integer		nFill		/ 4/		! Min. # neighbours for fill
	integer		NT		/ 0/

	real		BINXVmin	/  4.0/		! Minimum required # los crossings at source surface
	real		BINXGmin	/ 10.0/		! (passed to BuildSourceSurface)
	real		YLbeg		/-90.0/		! Heliographic latitude limits (MapGrid)
	real		YLend		/ 90.0/
	real		XElowG		/ 30.0/		! Limits on los elongations (ReadVIPS,ReadGIPS)
	real		XElowV		/  5.0/
	real		XEhighG		/ 80.0/
	real		XEhighV		/180.0/
	real		XRlimG		/ 90.0/		! Limits on hour angle difference with Sun (ReadVIPS,ReadGIPS)
	real		XRlimV		/ 90.0/
	real		GPOWER		/  0.0/		! Determines radial dependence of G (ReadGIPS)
	real		ClipLng		/  0.0/		! Needs to be a non-zero for NOMOD version (90?)

	real		V0		/400.0/		! Passed to ReadGIPS; constant V map if no IPS V data
	real		PWN		/  0.3/		! Power for normalized density dependence
	real		PWR		/  1.4/		! Power for radial dependence (deviation from r^-2)

	real		SMOOTH_D	/  6.0/		! Gaussian width (deg) for angular smoothing (GridSphere2D)
	real		SMOOTH_V	/  7.5/
	real		SMOOTH_TIME	/  0.0/
	real		D1AU		/  1.0/		! Density at 1 AU
	real		GSIGB		/  3.0/
	real		VSIGB		/  3.0/
	real		DAV		/  1.0/		! Passed to MapGrid

	parameter	(nCar	= 1000)
	real*8		JDCar(nCar)
	real*8		MJDref

	character	cCam*4
	character	cStr*80
	character	cSay*6		/'ipsg2t'/
	character	cVar(2)*100
	character	cArg*100
	character	cWildNagoya*50
	character	cWildOoty  *50

	intrinsic	sqrt

	integer		Flt2Str
	integer		Str2Str

	real		Scale(4)	/0.0,1.0,0.0,1.0/

	iVar = 2
	call ForeignArg(' ',iVar,cVar,cArg)

	Bad = BadR4()

	I__MODE	= TOM__NORATIO

	!-------
	! Check command line for switches (the grid spacings will be checked later)

	if (itrim(cArg) .gt. 0) then
	    call Say(cSay,'I','CmdLine',cArg)
	    if (LocFirst(cSwitch//'pointp',cArg) .ne. 0) I__MODE = I__MODE+TOM__POINTP	! Point-P map as initial condition
	    if (LocFirst(cSwitch//'mod'   ,cArg) .ne. 0) I__MODE = I__MODE+TOM__MOD	! Mod vs nomod version
	    if (LocFirst(cSwitch//'diff'  ,cArg) .ne. 0) I__MODE = I__MODE+TOM__DIFF	! Differences instead of ratios
	    if (LocFirst(cSwitch//'mass'  ,cArg) .ne. 0) I__MODE = I__MODE+TOM__MASS	! Goes to SW_Model_Kinematic 
	    if (LocFirst(cSwitch//'inclg' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__INCLG	! Goes to LOSIntegralG
	    if (LocFirst(cSwitch//'inclv' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__INCLV	! Goes to LOSIntegralV
	    if (LocFirst(cSwitch//'trackg',cArg) .ne. 0) I__MODE = I__MODE+TOM__TRACKG	! Tracks G-convergence in addition to G^2
	end if

	bPointP = iand(I__MODE,TOM__POINTP) .ne. 0
	bMod    = iand(I__MODE,TOM__MOD   ) .ne. 0
	bDiff   = iand(I__MODE,TOM__DIFF  ) .ne. 0
	bTrackG = iand(I__MODE,TOM__TRACKG) .ne. 0

	I = 0
	if (.not. bMod ) I = I+Str2Str('in quasi-time-dependent mode',cStr(I+1:))
	if (	  bMod ) I = I+Str2Str('in time-independent mode'    ,cStr(I+1:))
	if (.not. bDiff) I = I+Str2Str('#using ratios of model and observed los values#.',cStr(I+1:))
	if (	  bDiff) I = I+Str2Str('#using differences of model and observed los values#.',cStr(I+1:))
	call Say(cSay,'I','Running',cStr)


	call IPSProcessing(bForecast,NitTotal,bVcon,bVNago,bVOoty,bGcon,bGNago,bGCamb,bGOoty,bMV2)

	if (bVNago .or. bGNago) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya)
	if (bVOoty .or. bGOoty) call ForeignFile(iVar,cVar,'ooty'  ,cWildOoty  )
	
	! If bMV2=.TRUE. then either bGcon or bVcon is .FALSE. and the other is .TRUE.

	if (bMV2) then			! Either bGcon or bVcon is .FALSE.
	    if (bVcon) SMOOTH_D = SMOOTH_V
	    if (bGcon) SMOOTH_V = SMOOTH_D
	else				! Both bGcon and bVcon are .TRUE.
	    if (bVNago .and. bGNago) SMOOTH_D = SMOOTH_V ! G and V from same source
	end if

	call AskR4('Power of g-dependence on normalized density'   ,PWN)
	call AskR4('Power of g-dependence on heliocentric distance',PWR)

	if (bVcon) call IPSConstraints(bForecast,TOM__V,NitTotal,bBoundV,BINXVmin,NitBadV,iEorWV,XElowV,XEhighV,XRlimV)

	if (bGcon) then

	    if (bVNago .and. bGNago) then
		bBoundG  = bBoundV
		BINXGmin = BINXVmin
		NitBadG  = NitBadV
		iEorWG   = iEorWV
		XElowG   = XElowV
		XEhighG  = XEhighV
		XRlimG   = XRlimV

	    else
		call IPSConstraints(bForecast,TOM__G,NitTotal,bBoundG,BINXGmin,NitBadG,iEorWG,XElowG,XEhighG,XRlimG)

	    end if

	    call AskR4('Solar wind speed for G-level traceback',V0)

	end if


	!-------
	! Establish defaults for numerical grid (updated in SetGrid)

	nLng = nLng0
	nLat = nLat0
	nRad = nRad0
	nTim = nTim0

	!-------
	! Set up the t3d array. The array is saved internally. From here on access to
	! its elements should be done only using T3D_iget, T3D_iset, T3D_get, T3D_set.

	call T3D_iset(T3D__MODE		,0,I__MODE)

	!-------
	! Store grid dimensions.
	! Then check whether arrays are big enough to handle them.


	call T3D_iset(T3D__NLNG		,0,nLng	)
	call T3D_iset(T3D__NLAT		,0,nLat	)
	call T3D_iset(T3D__NRAD		,0,nRad	)
	call T3D_iset(T3D__NTIM		,0,nTim	)


	call T3D_set (T3D__PWN_V	,0,PWN	)
	call T3D_set (T3D__PWN_G	,0,PWN	)
	call T3D_set (T3D__PWR_V	,0,PWR	)
	call T3D_set (T3D__PWR_G	,0,PWR	)


	call T3D_set (T3D__D1AU		,0,D1AU	)

	call T3D_set (T3D__SMOOTH_V	,0,SMOOTH_V)
	call T3D_set (T3D__SMOOTH_D	,0,SMOOTH_D)
	call T3D_set (T3D__FILL_V	,0,SMOOTH_V)
	call T3D_set (T3D__FILL_D	,0,SMOOTH_D)
	call T3D_set (T3D__SMOOTH_TIME_V,0,SMOOTH_TIME)
	call T3D_set (T3D__SMOOTH_TIME_D,0,SMOOTH_TIME)

	call T3D_iset(T3D__ITERATION	,0,0	)
	call T3D_iset(T3D__ITERATIONS	,0,NitTotal)

	call T3D_set (T3D__BINX_V	,0,BINXVmin)
	call T3D_set (T3D__BINX_D	,0,BINXGmin)

	call T3D_set (T3D__CLIPLNG	,0,ClipLng)
	call T3D_set (T3D__SCALE	,4,Scale)


 					! Set time limits on data selection
	call SetGrid(bForecast, nCar, JDCar, XCbegMAT, XCendMAT, XCbegROI, XCendROI,
     &		XCtst1, XCtst2, MJDref,  bGCamb, MJDfrst, MJDlast, cCam, iEdt, NT)



	XCrange = XCendMAT(1)-XCbegMAT(1)

	!-------
	! Pick up constants set by SetGrid

	call T3D_iget(T3D__MODE ,0,I__MODE)
	call T3D_iget(T3D__NCOFF,0,NCoff )
	call T3D_get_grid(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1)

	call T3D_show()


	nSS1	= nLng*nLat		! # points at source surface at given time
	nSST	= nSS1*nTim		! # points at source surface at all times combined
	n3DT	= nSST*nRad		! # point in reconstruction space (spatial + temporal)

	if (nSS1 .gt. nSS0 ) call Say(cSay,'F','Too','many bins at source surface; nLng*nLat too large')
	if (nSST .gt. nSST0) call Say(cSay,'F','Too','many bins at source surface; nLng*nLat*nTim too large')
	if (n3DT .gt. n3DT0) call Say(cSay,'F','Too','many voxels in volume; nLng*nLat*nRad*nTim too large')

	!-------
	! The range of the Carrington variable needed is [XCbegMAT,XCendMAT]. The  target range
	! [XCtst1,XCtst2] includes an extra safety margin to make sure that all relevant
	! data are read from file.

	NLmaxG	= NLG0
	NLmaxV	= NLV0
	NLOS	= NLOS0
	dLOS	= dLOS0

	call T3D_set (T3D__DLOS,0,dLOS	)
	call T3D_iset(T3D__NLOS,0,NLOS	)

	! Only one of (bVNago, bVOoty) is .TRUE.

	call ReadVIPSValidData(LOS__CHECKV,Bad,Bad,Bad,Bad)	! Check for valid V data

	if (bVNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref,
     &		RR,XElowV,XEhighV,iEorWV,XRlimV,NLmaxV,NLV,
     &		IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bVOoty) call ReadVIPS(iReadOoty  ,iProcessOoty  ,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref,
     &		RR,XElowV,XEhighV,iEorWV,XRlimV,NLmaxV,NLV,
     &		IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav)

	! Only one of (bGCamb, bGOoty, bGNago) is .TRUE.
	! NOTE: the dummy array rLOSTmp is needed to avoid overwriting the array VOBS

	call ReadVIPSValidData(LOS__CHECKG,Bad,Bad,Bad,Bad)	! Check for valid G data

	if (bGCamb) call ReadGIPS(cCam,nCar,JDCar,
     &		iEdt,.FALSE.,bForecast,
     &		XCtst1,XCtst2,MJDref,MJDfrst,MJDlast,
     &		RR,V0,GPOWER,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG,
     &		iXPG,iYPG,iMJDG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS,
     &		1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav,XDSsav,XLSsav,
     &		XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bGNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref,
     &		RR,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG,
     &		IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bGOoty) call ReadVIPS(iReadOoty  ,iProcessOoty  ,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref,
     &		RR,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG,
     &		IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	call T3D_iset(T3D__NL_V,0,NLV	)
	call T3D_iset(T3D__NL_G,0,NLG	)



	call StopWatch(' ',' ')			! Start stopwatch


	!-------
	! The nRad entry in the t3d array may be modified by LOSReach.
	! Save the value here so we can restore it later.

	nRadFull = nRad


	!-------
	! Process the line-of-sight data:
	! - LOSPosition	: Calculate positions of all segments
	! - LOSReach	: Calculate range of positions of all segments; lines of sight
	!		  that fall outside the numerical grid are flagged as bad in NBS
	! - LOSClean	: iscard lines-of-sight flagged by LOSReach
	! - LOSWeights	: Calculate line-of-sight weights for all segments

	if (bVcon) then
	    call LOSPosition(TOM__V,IYRSV,DOYSV,DISTV,XLSV,XCEV,XLLV,XDLV,PPlosV,VPlosV)

	    call LOSReach(TOM__V,XCbegMAT,XCendMAT,PPlosV,NBSV)

	    call LOSClean(TOM__V,.FALSE.,NBSV,iXPsav,iYPsav,iMJDsav,
     &		  IYRFV,IRECV, IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,PPlosV)

	    IPSV = 1
	    if (bVNago) IPSV = 1		! Calculate weights, WWlosV, along all lines of sight
	    call LOSWeights(TOM__V,IPSV,XEV,DISTV,PPlosV,WWlosV,rLOSTmp)

	    call Say(cSay,'I',' ','.#.')
	end if

        if(bGcon) then
	    call LOSPosition(TOM__G,IYRSG,DOYSG,DISTG,XLSG,XCEG,XLLG,XDLG,PPlosG,rLOSTmp)

	    call LOSReach(TOM__G,XCbegMAT,XCendMAT,PPlosG,NBSG)

	    call LOSClean(TOM__G,bGCamb,NBSG,iXPG,iYPG,iMJDG, 
     &		  IYRFG,IRECG, IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS,PPlosG)

	    IPSG = 2				! Default: Cambridge
	    if (bGNago) IPSG = 1		! Calculate weights, WWlosG, along all lines of sight
	    call LOSWeights(TOM__G,IPSG,XEG,DISTG,PPlosG,WWlosG,M2Mean)

	    call Say(cSay,'I',' ','.#.')
	end if

	call T3D_iget(T3D__NL_V,0,NLV )		! Updated by LOSClean
	call T3D_iget(T3D__NL_G,0,NLG )		! Updated by LOSClean
	call T3D_iget(T3D__NRAD,0,nRad)		! Updated by LOSReach

	if (bGcon) call ArrR4TimesArrR4(NLG,GOBS,GOBS,G2OBS)	! Switch to g^2

	!-------
	! At this point all lines of sight to be used in the reconstruction have
	! been collected. Before starting the reconstruction set up a standard
	! for the 'goodness of fit', based on a 'flat' source surface velocity
	! and density map (these are saved internally by LOSTweak)

	call ArrR4Constant(nSST,V0,VFin)
	call ArrR4Constant(nSST,1.,GFin)

	call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VFin,GFin,XC3D,V3D,G3D,TmpSST)

	if (NLV .gt. 0) then
	    call LOSIntegralV(TOM__V,XCbegMAT,XCendMAT,V3D,G3D,PPlosV,WWlosV,VPlosV,VMDL,Weight,VSUM,rLOSTmp)
	    call LOSTweak(TOM__V,VOBS,VMDL,NBSV,FIX,FIXMEAN,FIXSTDV)
	end if

	if (NLG .gt. 0) then
	    call LOSIntegralG(TOM__G,XCbegMAT,XCendMAT,G3D,M2Mean,PPlosG,WWlosG,G2MDL,Weight,rLOSTmp)

	    if (bTrackG) then
		call ArrR4Function(NLG,sqrt,G2MDL,GMDL)
		call LOSTweak(TOM__G,GOBS,GMDL,NBSG,FIX,FIXMEAN,FIXSTDV)
	    end if
	    call LOSTweak(TOM__G,G2OBS,G2MDL,NBSG,FIX,FIXMEAN,FIXSTDV)
	end if

	!-------
	! Set up the initial maps at the source surface.
	! The default is constant V, constant normalized density

	call ArrR4Constant(nSST,1.,GMap)
	call ArrR4Constant(nSST,V0,VMap)

	!-------
	! Replace the default by a 'point-P' map. The point-P map may have holes in it
	! so call GridSphere2D to fill the holes (it also smooths the map)

	if (bPointP) then
	    NC = 1
	    if (bVcon) then
		call MapGrid(XCbegMAT(1),XCendMAT(1),YLbeg,YLend,NC, NLV,XCV,YLV,VOBS,
     &			NAV,DAV,nLng,nLat,VMap,Gtmp,Vtmp,Zmin,Zmax,ZDmin,ZDmax)
		call GridSphere2D(XCrange,nLng,nLat,nTim,VMap,SMOOTH_V,3,0.0,ClipLng)

		N = 1
		do I=2,nTim		! Use same V source surface map for all times
		    N = N+nSS1
		    call ArrR4Copy(nSS1,VMap,VMap(N))
		end do
	    end if

	    if (bGcon) then
		call MapGrid(XCbegMAT(1),XCendMAT(1),YLbeg,YLend,NC,NLG,XCG,YLG,G2OBS,
     &			NAV,DAV,nLng,nLat,GMap,Gtmp,Vtmp,Zmin,Zmax,ZDmin,ZDmax)
		call GridSphere2D(XCrange,nLng,nLat,nTim,GMap,SMOOTH_D,3,0.0,ClipLng)

		N = 1
		do I=2,nTim		! Use same G source surface map for all times
		    N = N+nSS1
		    call ArrR4Copy(nSS1,GMap,VMap(N))
		end do

	    end if
        end if

	call LOSSanityCheck(bVcon,bGcon)

	!-------
	! Iteration loop

	do Nit=1,NitTotal
	    call Say(cSay,'I',' ','.#Iteration '//cStr(:Int2Str(Nit,cStr)))

	    call T3D_iset(T3D__ITERATION,0,Nit)

	    if (bVcon) then				! V deconvolution
							! bMV2=.TRUE., implies bGcon=.FALSE.
		if (bMV2) then
		    call MkV2D(nSST,VMap,GMap,V0)	! Calculate nr^2
		    call CvD2G(PWN,-nSST,GMap,GMap)	! nr^2 -> g^2
		end if

		call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST)

		call LOSIntegralV(TOM__V,XCbegMAT,XCendMAT,V3D,G3D,PPlosV,WWlosV,VPlosV,VMDL,Weight,VSUM,rLOSTmp)

		if (Nit .eq. 1) call ArrR4Copy(NLV,VMDL,VMDL0)

		call LOSProjection(TOM__V,XCbegMAT,XCendMAT,XC3D,PPlosV,PPprj)

		call LOSTweak(TOM__V,VOBS,VMDL,NBSV,FIX,FIXMEAN,FIXSTDV)

		bRemove = Nit .eq. NitBadV
c		bRemove = Nit/50*50 .eq. Nit

		call BuildSourceSurface(TOM__V,bRemove,bBoundV,BINXVmin,NBSV,FIX,FIXMEAN,FIXSTDV,
     &			PPprj,Weight,XCbegMAT,XCendMAT,VMap,VSIGB,Vtmp,Gtmp,BINXV,iLOSTmp)

		if (Nit .eq. NitTotal) call ArrR4Copy(nSST,VMap,VFin)

		call GridSphere2D(XCrange,nLng,nLat,nTim,VMap,SMOOTH_V,3,0.0,ClipLng)

		call Say(cSay,'I',' ','.')
	    end if

	    if (bGcon) then				! G deconvolution
							! bMV2=.TRUE., implies bVcon=.FALSE.
		if (bMV2) then
		    call CvG2D(PWN,-nSST,GMap,VMap)	! In: GMap = g^2 -> Out: VMap = nr^2
		    call MkD2V(nSST,VMap,VMap,V0)	! nr^2 --> Velocity
		end if

		call SW_Model_Kinematic(TOM__G,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST)

		call LOSIntegralG(TOM__G,XCbegMAT,XCendMAT,G3D,M2Mean,PPlosG,WWlosG,G2MDL,Weight,rLOSTmp)

		if (Nit .eq. 1) call ArrR4Copy(NLG,G2MDL,G2MDL0)

		call LOSProjection(TOM__G,XCbegMAT,XCendMAT,XC3D,PPlosG,PPprj)

		if (bTrackG) then
		    call ArrR4Function(NLG,sqrt,G2MDL,GMDL)
		    call LOSTweak(TOM__G,GOBS,GMDL,NBSG,FIX,FIXMEAN,FIXSTDV)
		end if
		call LOSTweak(TOM__G,G2OBS,G2MDL,NBSG,FIX,FIXMEAN,FIXSTDV)
		
		bRemove = Nit .eq. NitBadG
c		bRemove = Nit/50*50 .eq. Nit

		call BuildSourceSurface(TOM__G,bRemove,bBoundG,BINXGmin,NBSG,FIX,FIXMEAN,FIXSTDV,
     &			PPprj,Weight,XCbegMAT,XCendMAT,GMap,GSIGB,Vtmp,Gtmp,BINXG,iLOSTmp)

		if (Nit .eq. NitTotal) call ArrR4Copy(nSST,GMap,GFin)

		call GridSphere2D(XCrange,nLng,nLat,nTim,GMap,SMOOTH_D,3,0.0,ClipLng)! Fill gaps and smooth

		call Say(cSay,'I',' ','.')
	    end if

	end do

	call Say(cSay,'I',' ','.#.#Finished iterating#.#.')


	nRad = nRadFull
	call T3D_iset(T3D__NRAD,0,nRad)
	call T3D_iset(T3D__ITERATION,0,Nit)

	!-------
	! The velocity map VFin and density map GFin (as produced by the last
	! BuildSourceSurface calls) contain bad values indicating where the source surface
	! maps are invalid. GMap and VMap have gone to one more GridSphere2D,
	! filling in all bad values and smoothing the array.

	! Use VMap and GMap to create the final shift array XC3D. Output from
	! SW_Model_Kinematic (XC3D, V3D, G3D) does not contain any bad values (the 'bad
	! value' information is stored in the source surface maps VMap, GMap).

	call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST)

	! In the original program there was an inconsistency here: XC3D, which
	! is based on VMap and GMap, was combined with VFin and GFin to find
	! the final velocity array. Here I'm using VMap and GMap after copying
	! the 'holes' from VFin and GFin. To go back to the old treatment,
	! comment the next 5 lines

	!=======

	call ArrR4Mask(nSST,VFin,Bad,Bad,Bad,0.,1.,VMap)! Copy bad values of VFin into VMap
	call ArrR4Mask(nSST,GFin,Bad,Bad,Bad,0.,1.,GMap)! Copy bad values of GFin into GMap

	call ArrR4Copy(nSST,VMap,VFin)			! Copy VMap into VFin
	call ArrR4Copy(nSST,GMap,GFin)			! Copy GMap into GFin

	!=======

	N = -nSS1+1
	do I=1,nTim					! Fill in small holes
	    N = N+nSS1
	    call GridFill(nFill,nLng,nLat,VFin(N),Btmp,Zmin,Zmax)
	    call GridFill(nFill,nLng,nLat,GFin(N),Btmp,Zmin,Zmax)
	end do


	!-------
	! Write3D_nv_UT does not modify any of its input arguments

	call Write3D_nv_UT(NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,VFin,GFin,XC3D,V3D,G3D,Tmp3D1,Tmp3D2)


	!-------
	! Write the t3d* file.
	! Write3D_nv returns velocity in V3D and normalized density in GFin and G3D.
	! Output arrays contain bad values based on the bad values in the
	! input arrays VFin and GFin.

	call Write3D_nv(XCbegMAT,XCendMAT,XCbegROI,XCendROI,VFin,GFin,XC3D,V3D,G3D,BINXV,BINXG)

	!-------
	! Write the b3d* file. Note that V3D is the full 3D velocity matrix

	call Write3D_bb   (		 XCbegMAT,XCendMAT,XCbegROI,XCendROI,XC3D,V3D,Tmp3D1,Tmp3D2,Gtmp)
	call Write3D_bb_UT(NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,XC3D,V3D,Tmp3D1,Tmp3D2,TmpSST1,TmpSST2,Btmp)

	if (bVcon) then
	    I = Str2Str('los'//char(TOM__V),cStr)
	    I = Flt2Str(NCoff+TT,-4,cStr(I+1:))
	    if (dTT .ne. 0) I = T3D_marker(cStr,icount)
	    if (bVNago) call WriteLOSY(cStr,cWildNagoya,NLV,IYRFV,IRECV,VOBS,VMDL0,VMDL,XEV,NBSV)
	    if (bVOoty) call WriteLOSY(cStr,cWildOoty  ,NLV,IYRFV,IRECV,VOBS,VMDL0,VMDL,XEV,NBSV)
	end if

	if (bGcon) then
	    I = Str2Str('los'//char(TOM__G),cStr)
	    I = Flt2Str(NCoff+TT,-4,cStr(I+1:))
	    if (dTT .ne. 0) I = T3D_marker(cStr, icount)
	    if (bGNago) call WriteLOSY(cStr,cWildNagoya,NLG,IYRFG,IRECG,G2OBS,G2MDL0,G2MDL,XEG,NBSG)
	    if (bGOoty) call WriteLOSY(cStr,cWildOoty  ,NLG,IYRFG,IRECG,G2OBS,G2MDL0,G2MDL,XEG,NBSG)
	    if (bGCamb) call WriteLOSD(cStr,cCam,iEdt,NLG,iMJDG,iXPG,iYPG,G2OBS,G2MDL0,G2MDL,XEG,NBSG)
	end if

	call Say(cSay,'I','Removed',cStr(:Int2Str(NLV-iArrI4Total(NLV,NBSV,I),cStr))//' bad V sources')
	call Say(cSay,'I','Removed',cStr(:Int2Str(NLG-iArrI4Total(NLG,NBSG,I),cStr))//' bad G sources')

	call StopWatch(' ',' ')
	call Say(cSay,'S','Stop','Program successfully completed')

	end
