C+
C NAME:
C	ipshtd77 
C
C	The latest 9/15/00 work - tried a copydtodv here and this filled
C	the density matrix a lot better when it was placed just before write3dinfo - 
C	There was no apparent change in extractd from the ea_file but maybe something 
C	was different with ulysses.  At this stage, this program does not give exactly
C	the same thing as the current version on the vax, but it is close.
C	Found an error 10/15/00 (A gridsphere2d was placed outside a do loop after 
C	the first MKshiftd. This version should now have no seam in each temporal value.
C
C	Modifications to the current program until BVJ 03/01/02
C
C	To compile the program with the g77 compler type:
C 	
C	./mk_ipshtd77
C
C	To run the program type:
C
C	$myfor/IPSHTD/ipshtd77 ucsd=SANIPS,$dat/ucsd/ helios1=A77B4_318_"(4)".AVE,$HOME/dat/Helios/
C	helios2=B77B4_318_"(4)".AVE,$HOME/dat/Helios/
C
C	Program defaults should allow the November 24 1977 data files to be run to give density and 
C	magnetic field.
C	 
C PURPOSE:
C	Construct deconvolved synoptic maps for the IPS observations using
C	a "time-dependent" technique that projects all observations to a
C	single lower reference surface.  Note - uneven time cadences are
C 	allowed for velocity and density data sets.
C       
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	run ipshtd
C INPUTS:
C	Data file $CAM:G*.
C OUTPUTS:
C	Rectangular grid of data points
C FUNCTIONS/SUBROUTINES:
C	Local2UT,MAP_TZERO,READ_GO,
C	SetGipsy
C INCLUDE:
C	include		'sun.h'
C COMMON BLOCKS:
C SIDE EFFECTS:
C RESTRICTIONS:
C PROCEDURE:
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 ipshtd), the
C	only way to force ipshtd to look in another directory is to delete the
C	logical $CAM.
C MODIFICATION HISTORY:
C	JUN-1992, Paul Hick (UCSD)
C       NOV-1995, Bernie Jackson (STEL)
C	March-1999 Bernie Jackson (UCSD)
C	June -1999 Bernie Jackson (UCSD)
C	August -1999 Bernie Jackson (UCSD)
C	Jan/Feb 2001  Bernie Jackson (UCSD)
C	Mar 2001 Bernie Jackson (UCSD)
C-
	program ipshtd
	parameter      (NLmaxG	=  22000,	! Max # G data points   (13000) (400,000)
     &			NLmaxV	=  500,		! Max # V data points   (2000)
     &			NGREST	=  1,		! G temporal resolution factor better than 1 day
     &			NVREST	=  1,		! V temporal resolution factor better than 1 day
     &			NGRESS	=  2,		! Spatial G resolution factor better than 20 degrees
     &			NVRESS	=  2,		! Spatial V resolution factor better than 20 degrees
     &			NF	=  1,		! LOS resolution factor (Should be 1 if Spatial resolution is >= 2)
     &			nTmaxG	=  50*NGREST,	! Max # G times 	(50) (100) (150)
     &			nTmaxV	=  50*NVREST)	! Max # V times 	(50) (100) (150)
	parameter      (NLOSG   =  20*NF*NGRESS, ! # resolution elements along each G line of sight  (40) (80) (125)
     &			NLOSV	=  20*NF*NVRESS, ! # resolution elements along each V line of sight (40) (80) (125)
     &			NLOSGP1 =  NLOSG+1,
     &			NLOSVP1 =  NLOSV+1,
     &			dLOSG	=  0.1/(NGRESS*NF), ! LOS G Resolution in units of observer-Sun distance  (0.05) (0.025) (0.016)
     &			dLOSV	=  0.1/(NVRESS*NF), ! LOS V Resolution in units of observer-Sun distance (0.05) (0.025) (0.016)
     &			xInc	=  0.5)		! Carrington map "extra" increment
	parameter	(nLng	= 54*NGRESS+1,	! # Longitudes, resolution is 360/(nLng-1) 	(55) (109) (217) (325)
     &			 iLng	= 18*NGRESS+1,  ! # Longitudes to use in extracted map     	(19)  (37)  (73) (109)
     &			 nLat	=  9*NGRESS+1,	! # Latitudes , resolution is 180/(nLat-1)   (9)(10)  (19)  (37)  (55)
     &			 nMap	=  31,		! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31)
     &			 nLngLat = nLng*nLat)
	parameter	(dRR	= 0.1)		! Resolution (AU) in radial direction 		(0.1)
	integer		iXPG	(NLmaxG),	! Array indices from daily G maps
     &			iYPG	(NLmaxG),
     &			iMJDG	(NLmaxG),	! Original MJD of daily G maps
     &			IYRFG	(NLmaxG),
     &			IRECG	(NLmaxG),
     &			IYRSG	(NLmaxG),	! Year of source observation
     &			IYRSGG	(NLmaxG),	! Year of source observation
     &			NBSG	(NLmaxG),
     &			LON	(NLOSGP1,NLmaxG),
     &			LAT	(NLOSGP1,NLmaxG),
     &			ITIM	(NLOSGP1,NLmaxG)

	real		DOYSG	(NLmaxG),	! DOY (and fraction) of source obs.
     &			DOYSGG	(NLmaxG),	! DOY (and fraction) of source obs.
     &			DISTG	(NLmaxG),	! G distance from Sun to Earth
     &			DISTGG	(NLmaxG),	! G distance from Sun to Earth
     &			XLSG	(NLmaxG),	! Ecliptic longitude of Sun
     &			XLSGG	(NLmaxG),	! Ecliptic longitude of Sun
     &			XLLG	(NLmaxG),	! Longitude diff. (Earth-Sun)
     &			XLLGG	(NLmaxG),	! Longitude diff. (Earth-Sun)
     &			XDLG	(NLmaxG),	! Declination of Source
     &			XDLGG	(NLmaxG),	! Declination of Source
     &			XCEG	(NLmaxG),	! Carrington variable of Earth
     &			XCEGG	(NLmaxG),	! Carrington variable of Earth
     &			 XCEARG	 (NLmaxG),	 ! Helios Carrington variable of Earth
     &			XCEARGG	(NLmaxG),	! Helios Carrington variable of Earth
     &			XEG	(NLmaxG),	! Source Elongation
     &			XEGG	(NLmaxG),	! Source Elongation
     &			XCG	(NLmaxG),	! Carrington variable
     &			XCGG	(NLmaxG),	! Carrington variable
     &			YLG	(NLmaxG),	! Heliographic latitude (degree)
     &			YLGG	(NLmaxG),	! Heliographic latitude (degree)
     &			GOBS	(NLmaxG),	! Solar wind observed G or brightness values
 - final
     &			GGOBS	(NLmaxG),	! Solar wind observed G or brigntness values
 - initial and modified
     &			BBASE	(NLmaxG),	! Solar wind observed base brightness values
     &			ZZ	(NLmaxG),	! Thomson scattering density values for MapGRID
     &			ZZZ	(NLmaxG),	! Thomson scattering density values for MapGRID
c    &			VOBS	(NLmaxV),	! Solar wind observed velocities
     &			VOBS	(NLmaxG),	! Solar wind observed velocities
     &			GMOD2	(NLmaxG),	! Model G values
     &			FIX	(NLmaxG),	! Fix to model to make OK
     &			GSSIG	(NLmaxG),	! G-level source sigma
     &			WTSG	(NLOSG ,NLmaxG),
     &			GWTij	(NLOSG ,NLmaxG),
     &			XLON	(NLOSGP1,NLmaxG),
     &			XLONP	(NLOSGP1,NLmaxG),
     &			XLAT	(NLOSGP1,NLmaxG),
     &			RPX	(NLOSGP1,NLmaxG),
     &			XCtpr	(NLOSGP1,NLMAXG),
     &			XLproj	(NLOSGP1,NLMAXG),
     &			XCtim	(NLOSGP1,NLMAXG),
     &			VFAC	(NLOSGP1,NLmaxG),	! New array (3/99)
     &			DFAC	(NLOSGP1,NLmaxG)

	character	SRCV	(NLmaxV)*9,		! Velocity source name
c    &			SRCVS*9,			! TJD, never referenced
     &			SRCVsav(NLmaxV)*9,
     &			SRCGG	(NLmaxG)*9,		! G-level source name
     &			SRCG	(NLmaxG)*9,
     &			SRCGS*9,SRCGsav(NLmaxG)*9
	integer		IYRFV	(NLmaxV),
     &			IRECV	(NLmaxV),
     &			NSRV	(NLmaxV),
     &			NSRG	(NLmaxG),
     &			IYRSV	(NLmaxV),		! Year of source observation
     &			NBSV	(NLmaxV)
	real		DOYSV	(NLmaxV),	! DOY (and fraction) of source obs.
     &			DISTV	(NLmaxV),	! V distance from Sun to Earth
     &			XLSV	(NLmaxV),	! Ecliptic longitude of Sun
     &			XLLV	(NLmaxV),	! Longitude diff. (Earth-Sun)
     &			XDLV	(NLmaxV),	! Declination of Source
     &			XCEV	(NLmaxV),	! Carrington variable of Earth
     &			XEV	(NLmaxV),	! Source Elongation
     &			XCV	(NLmaxV),	! PLH
     &			YLV	(NLmaxV),	! PLH
     &			VVOBS	(NLmaxG),
     &			VMOD	(NLmaxV),	! Model velocities
     &			VSSIG	(NLmaxV),	! Velocity source sigma	
     &			VWTij	(NLOSV,NLmaxV),
     &			WTSV	(NLOSV,NLmaxV)
	real		XCinTV	(nTmaxV),
     &			XCinTG	(nTmaxG),
     &			DOYV	(nTmaxV),
     &			DOYG	(nTmaxG),
     &			XCbeV	(2,nTmaxV),
     &			XCbeGG	(2,nTmaxG)
	integer		IYRV	(nTmaxV),
     &			IYRG	(nTmaxG) 
c    &			NMAPD	(nLng,nLat,nTmaxG),	TJD, never ref.		! Scratch integer arrays
c    &			NMAPV	(nLng,nLat,nTmaxV)	TJD, never ref.	 
	real		VDEV	(nLng,nLat,nTmaxV),
     &			VPNT	(nLng,nLat,nTmaxV),
     &			PPGMAP	(nLng,nLat,nTmaxG),	! Point-P G map (with holes)
     &			VMAPHOLE(nLng,nLat,nTmaxV),
     &			DM	(nLng,nLat),
     &			VM	(nLng,nLat),
     &			DMAPHOLE(nLng,nLat,nTmaxG),
     &			DMAP	(nLng,nLat,nTmaxG),	! Density map (no holes)
     &			DBASE	(nLng,nLat,nTmaxG),	! Base map to use in Thomson scattering
     &			VMAP	(nLng,nLat,nTmaxV),
     &			VMAPD	(nLng,nLat,nTmaxG),
     & 			DMAPV	(nLng,nLat,nTmaxV),
     &			GDEV	(nLng,nLat,nTmaxG),
     &			GPNT	(nLng,nLat,nTmaxG),
     &			DDfact	(nLng,nLat,nMap,nTmaxG),
     &			DVfact	(nLng,nLat,nMap,nTmaxG),	! New array (3/99)
     &			XCshift	(nLng,nLat,nMap,nTmaxG,3),
     &			V3D	(nLng,nLat,nMap),
     &			D3D	(nLng,nLat,nMap),
     &			Vtmp 	(nLng,nLat,nTmaxG,3),		! Scratch real arrays
     &			VtmpV	(nLng,nLat,nTmaxV,3),
     &			VZtmp	(nLng,nLat,nTmaxV),
     &			Dtmp 	(nLng,nLat,nTmaxG),
     &			XLTtmp	(nLng,nTmaxG,2),
     &			XLtmp	(nLng,nTmaxG),
     &			XLtmpt	(nLng,nTmaxG),
     &			Dtmpt	(nLng,nTmaxG),
c    &			Vtp	(nLng,nLat),			TJD never ref
c    &			Dtp	(nLng,nLat),			TJD never ref
     &			ANMAPD	(nLng,nLat,nTmaxG),
     &			DMAD	(nLng,nLat,nTmaxG),
     &			FIXMD	(nLng,nLat,nTmaxG),
     &			WTSMD	(nLng,nLat,nTmaxG),
     &			DMEANFD	(nLng,nLat,nTmaxG),
     &			DFIXM2D	(nLng,nLat,nTmaxG),
     &			WTPD	( NLOSG,NLmaxG),
     &			GMWTi	(NLmaxG),
     &			ANMAPV	(nLng,nLat,nTmaxV),
     &			VMAV	(nLng,nLat,nTmaxV),
     &			FIXMV	(nLng,nLat,nTmaxV),
     &			WTSMV	(nLng,nLat,nTmaxV),
     &			VMEANFV	(nLng,nLat,nTmaxV),
     &			VFIXM2V	(nLng,nLat,nTmaxV),
     &			WTPV	(NLOSV,NLmaxV),
     &			VVV	(nLng),
     &			DDD	(nLng),
     &			VVT	(nTmaxV),
     &			DDT	(nTmaxG),		
     &			VRATION	(200),
     &			GRATION	(200),
     &			Scale	(2)
	byte		NSIDE(9*nLng*nLat)
	parameter	(nCar	  = 2000)
	parameter	(nCarHOS1 = 2000)
	parameter	(nCarHOS2 = 2000)
	real*8		JDCar(nCar),
     &			JD, JEpoch,
     &			dLngSun,dLatSun,dDisSun ,
     &			JDCarHOS1(nCarHOS1),
     &			JDCarHOS2(nCarHOS2),
     &			JDCX, JDCX2, JDCP,
     &			MJD_to_JD	/2400000.5d0/,	! JD = MJD+MJD_to_JD
     &			MJDref,
     &			MJDcntr,
c     &			MJDend,				! AUTO mode: last MJD (TJD, never ref)
     &			MJDinc		/.757639d0/,	! AUTO mode: increment
     &			DLONG,
     &			DLAT,
     &			DDIST
C-------
C String written to output files. First three chars are also used in file name construction
        character	cStrPPV(1)*80	/'PPV: Point P velocity map with holes'/,
     &			cStrPPG(1)*80	/'PPG: Point P G-level map with holes'/,
     &			cStrPPD(1)*80	/'PPD: 1 AU Point P density map with holes'/,
     &			cStrDV (1)*80	/'DVV: Deconvolved velocity map'/,
     &			cStrDD (1)*80	/'DDD: Deconvolved density map'/,
     &			cStrDVs(1)*80	/'DVS: Deconvolved, smoothed velocity map'/,
     &			cStrDDs(1)*80	/'DDS: Deconvolved, smoothed density map'/,
     &			cStrVname*14	/'vts1941_00.txt'/,
     &			cStrGname*14	/'gts1941_00.txt'/,
     &                  cSym*17         /'LIB__SIGNAL_QUIET'/
	logical		bGcon,	bVcon,  bTcon,
     &			bVNago, bVOoty,  bVUCSD,
     &			bGcamb, bHelio, bGOoty, bGNago, bh1, bh2,
     &			bVmv2		/.FALSE./,	! Use mv^2 = constant 	
     &			bGmv2		/.FALSE./,	! Use mv^2 = constant
     &			bDaydV, bDaydG,
     &			bOnly1,				! If .TRUE., use only one shift matrix (on density matrix) 
     &                  bVproxy         /.FALSE./,      ! Velocity proxy map used for deconvolution
     &                  bDproxy         /.FALSE./,      ! Density proxy map used for deconvolution
     &			bMagnetic	/.FALSE./,	! Magnetic map available
     &			bGWRITE		/.FALSE./,	! Write out a G data file
     &			bVWRITE		/.FALSE./,	! Write out a V data file
     &			bMirror		/.FALSE./,
     &			bjMirror	/.FALSE./,
     &			bWrite		/.FALSE./,
     &			bWrite3D	/.FALSE./,
     &			bWriteMS	/.TRUE./,
     &			bStopw		/.FALSE./,
     &			bForecast,
     &			bExtract(12)	/2*.FALSE.,.TRUE.,9*.FALSE./,
     &			bEWRITE,			! Write out an Earth file
     &			bMWRITE,			! Write out a Mercury file
     &			bUWRITE,			! Write out a Ulysses file
     &                  bH1WRITE,        		! Write out an Helios 1 file
     &                  bH2WRITE        		! Write out an Helios 2 file
	equivalence	(bExtract( 3), bEWRITE),
     &			(bExtract(1), bMWRITE),
     &			(bExtract(10), bUWRITE) ,
     &			(bExtract(11), bH1WRITE),
     &			(bExtract(12), bH2WRITE)
 
	integer		iEorWG		/0/,
     &			iEorWV		/0/,
c     &			MJD,
c     &			MJDfrst		/48180/, ! MJD for earliest data file
c    &			MJDlast		/49622/, ! MJD for latest data file
     &			MJDfrst		/ 44100/,
     &			MJDlast		/ 51800/,
     &			LimitoV		/-1/,
     &			LimitoG		/-1/,
     &			NLOSWV		/1/,
     &			NLOSWG		/1/,
     &			NTV		/44/,
     &			NTG		/44/
	integer		NiterT		/18/,
     &			iOffUT		/ 8/,
     &			NAV		/ 0/,
     &			NC		/ 1/,		! # rotations averaged
     &			nDUMP		/ 0/,		! Threshold for #pnts/bin
     &			nFILL		/ 3/,		! Min. # neighbours for fill
     &			NmidHRV		/0/,
     &			NmidHRG		/0/
	real		YLbeg		/-90./,
     &			YLend		/ 90./,
     &			XCadj(3)	/3*0./,
     &			XCadjB(3)	/3*0./,
     &			XCmargin 	/.5/,
     &			XElowG		/30./,
     &			XElowV		/5./,
     &			XEhighG		/80./,
     &			XEhighV		/180./,
     &			XRlimG		/90./,
     &			XRlimV		/90./,
     &                  GPOWER          /0./,
     &			R1AU		/1.0/,
C     &			DEN1AU		/6.8/ ,
     &			DEN1AU		/5.0/,
     &			VEL1AU		/425.0/,
     &			ANLIMITV	/1.5/,
     &			ANLIMITG	/1.5/,
     &			aNdayV		/1.0/,
     &			aNdayG		/1.0/
	real		VMAPSIG		/1.0/,
     &			DMAPSIG		/1.0/,
     &			GSIG		/1.0/,
     &			GSIGB		/3.0/,
     &			VSIG		/1.0/,
     &			VSIGB		/3.0/,
     &			FCONSV		/1.0/,
     &			FCONSG		/1.0/,
     &			CONRV		/15./,
     &			CONRD		/6./,
     &			CONVT		/0.55/,
     &			CONDT		/1.3/,
     &			ConsTV		/1.0/,
     &			ConsTG		/1.0/,
     &			VLIM		/1200./,	! Velocity upper limit for sources
     &			VLIL		/100./,		! Velocity lower limit for sources
     &			GLIMP		/4.0/,		! G-level upper limit for sources
     &			GLIMF		/1.0/,		! G-level input multiplier
     &			GLIMM		/0.2/,		! G-level lower limit for sources
     &			HELIOmax	/10000./,	! Maximum Helios brightness allowed
     &			HELIOmin	/-10000./	! Minimum Helios brightness allowed
 	real		RRS		/    -15. /,	! Reference distance for deconvolution 	
     &			VUL		/ 1500.0 /,	! Velocity upper limit for maps	
     &			VLL		/  100.0 /,	! Velocity lower limit for maps
     &			DUL		/  500.0 /,	! Density upper limit for maps
     &			DLL		/  0.0 /	! Density lower limit for maps
	real		DAV		/1./,
     &			Speed   	/400./,
     &			PWG		/0.35/,	!/0.5//1.0//2.0/	! Power of density
     &			PWV		/0.35/,	!/0.5//1.0//2.0/	! Power of density
     &			PWRG		/0.3/ ,				 ! Radial power of density
     &			PWRV		/0.3/,				! Radial power of velocity
     &			FALLOFF		/2.20/				! Density falloff with solar distance (>2.00 Implies
									! an acceleration of the solar wind with height.)
	character	cCam*4,
     &			cMon*3,
c    &			cStrID*11,		TJD, never ref
     &			cVar(3)*100,	
     &			cArg*100,
     &			cStr*80,
     &			cPre*4
	character	cWildNagoya*80,
     &			cWildOoty*80,
     &			cWildUCSD*80,
     &			cWildHelios1*80,
     &			cWildHelios2*580
C Below specific to Helios
	character	cHOS(10)*80
	integer		iSectBeg	/-15/,
     &			iSectEnd	/ 16/
	integer		iDayHOS		/ 0/,	! Time filter on Helios data
     &			nDevHOS		/10/
	 ! Glitch removal on Helios data (doesn't matter if iDayHOS is 0)
		
	real		dSigHOS		/0./
	 ! Offset of Helios data from mean (doesn't matter if iDayHOS is 0)
        character       STRING*80
	external	iReadNagoyan, iProcessNagoya n
	external	iReadOotyn, iProcessOotyn
	external	ireadUCSDn, iprocessUCSDn
 	external 	EARTH , HELIOS_1, HELIOS_2
	include		'dirspec.h'
	include		'sun.h'
	integer		Flt2Str
 
 
c**TJD**added for magnetic time dependant version
	real 		Btmp(nLng,nLat,3),		!scratch
     &			XCbegMAT(nTmaxG),		!3-Car time frame
     &			XCendMAT(nTmaxG),
     &			XCbegROI(nTmaxG),		!ROI time frame
     &			XCendROI(nTmaxG),
     &			TrueVel(nLng,nLat,nMap,nTmaxG),	!real velocity array
     &			BR3D(nLng,nLat,nMap,nTmaxG),	!scratch
     & 			BR2DT(nLng,nLat,nTmaxG),	!scratch
     &			Bt3D(nLng,nLat,nMap,nTmaxG)
 
 
 
c  C:\smei\exe\win\tipshtd.exe ucsd=SANIPS,C:\dat\ucsd\ helios1=A77B4_318_(4).R08,C:\dat\Helios\ helios2=B77B4_318_(4).R08,C:\dat\Helios\
	iVar = 3
	call ForeignArg(' ',iVar,cVar,cArg)
C-------
C Calculate start times for Carrington rotations between 1985 and now.
C Rotation NCoff+I starts at time JDCar(I).
	aNdayV	= aNdayV/NVREST
	aNdayG	= aNdayG/NGREST
	ConsTMV = ConsTV*nTV
	ConsTMD = ConsTG*nTG
	ConstL = nLng/3
	I = iDeleteSymbol(cSym,2)
	iYr = 1974
	Doy = 330.
	call MAP_TZERO(EARTH,iYr,Doy,.01,nCar,JDCar)
	NCoff = N_CARRINGTON(iYr,Doy)
	NCoffHOS1 = 0
	NCoffHOS2 = 0
	call ArrR4Constant(NLmaxV,5000.,XCEV)	! Initialize XCEV array 
	call ArrR4Constant(NLmaxG,5000.,XCEG)	! Initialize XCEGG array 
	call ArrR4Constant(NTmaxV,5000.,XCintV)	! Initialize XCintV array 
	call ArrR4Constant(NTmaxG,5000.,XCintG)	! Initialize XCintG array
C-------
C Set range of Helios Carrington variables to work with: [XClowHOS,XChighHOS]
	    XClowHOS1  =   1
C	    XClowHOS1  =  16
C	    XClowHOS1  =  28
C	    XClowHOS1  =  30
	    XChighHOS1 = 129
	    XClowHOS2  =  1
C	    XClowHOS2  =  16
C	    XClowHOS2  =  29
C	    XClowHOS2  =  30
	    XChighHOS2 =  65
	XCbegHOS1 = 38.5
	XCbegHOS2 = 38.5 
C------
C iOffUT is an offset time (in hours) between local time and UT. This is used only if the program
C is operated in forecast mode. Note that the local time is assumed to be the same as the
C computer time. The default offset is 8 hours (the difference between PST and GST).
	call AskYN('Forecast mode?$no',bForecast)		! Regular or forecast mode
	if (bForecast) call AskI4('Offset between local time and UT (hours)?',iOffUT)
	call AskI4('Number of iterations?',NiterT)		! # iterations
								
        call AskR4('Deconvolution distance (AU; 0=solar surface)?',RRS)
        if (RRS .lt. 0) RRS = -RRS*Sun__RAu			! Neg. RRS are in units of solar radii
        RRS = max(RRS,sngl(Sun__RAu))				! Keep RRS > 1 solar radius
								! Deconvolve velocity, G-level or both?
	call AskWhat('Deconvolve: velocity, density, both$1$3',I)
C-------
C A proxy map for density and/or velocity can be used to replace the point-P
C map determined from the observations.
C The map should be in files dcrot.xxx and vcrot.xxx. They must be ascii files of dimension
C nLng, nLat.
C
C If bVmv2=TRUE then bVcon=FALSE
C If bGmv2=TRUE then bGcon=FALSE
	bVcon = I .eq. 1 .or. I .eq. 3
	bGcon = I .eq. 2 .or. I .eq. 3
	if (bVcon) then			! V deconvolution: select IPS station
	    call AskWhat('Use velocity data from: Nagoya, Ooty, UCSD$1$3',I)
	    bVNago = I .eq. 1
	    bVOoty = I .eq. 2
	    bVUCSD = I .eq. 3
	    bVmv2  = .FALSE.
	    if(bVNago) NLOSWV = 1
	    if(bVOoty) NLOSWV = 2
	    if(bVUCSD) NLOSWV = 3
	else				! No V deconvolution
	    bVNago = .FALSE.
	    bVOoty = .FALSE.
	    bVUCSD = .FALSE.
	    call AskYN('Do you want to make a velocity map using mv^2 = con$no',bVmv2)
	end if
	if (bGcon) then			! G or density deconvolution: select IPS station  or Thomson scattering source
	    call AskWhat('Use density or G-level data from: Helios, Camb, Ooty, Nagoya$1$1',I)
	    bHelio = I .eq. 1
	    bh1 =   .TRUE.
	    bh2 =   .TRUE.
	    bTcon = .FALSE.
	    if(bHelio) call ASKYN('HELIOS 1',bh1)
	    if(bHelio) call ASKYN('HELIOS 2',bh2)
	    if(bHelio) GLIMP = HELIOmax
	    if(bHelio) GLIMM = HELIOmin
	    if (bh1) call MAP_TZERO(HELIOS_1,iYr,Doy,.01,nCarHOS1,JDCarHOS1)
	    if (bh2) call MAP_TZERO(HELIOS_2,iYr,Doy,.01,nCarHOS2,JDCarHOS2)
	    if(bHelio) bTcon = .TRUE.
	    bGCamb = I .eq. 2
	    bGOoty = I .eq. 3
	    bGNago = I .eq. 4
	    bGmv2  = .FALSE.
	    cCam = cEnvi//'CAM'
	    if(bh1)    NLOSWG = 5
	    if(bh2)    NLOSWG = 5
	    if(bGCamb) NLOSWG = 2
	    if(bGOoty) NLOSWG = 1
	    if(bGNago) NLOSWG = 1
	else				! No g deconvolution
	    bHelio = .FALSE.
	    bGCamb = .FALSE.
	    bGOoty = .FALSE.
	    bGNago = .FALSE.
	    call AskYN('Do you want to make a density map using mv^2 = con',bGmv2)
	end if
	    if (bGCamb) call SetGipsy(cCam,iOffUT,iEdt,MJDfrst,MJDlast)
	    if (bVNago .or. bGNago ) call ForeignFile(iVar,cVar,'nagoya' ,cWildNagoya)
	    if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty'   ,cWildOoty)
	    if (bVUCSD)		     call ForeignFile(iVar,cVar,'ucsd'   ,cWildUCSD)
	    if (bh1)		     call ForeignFile(iVar,cVar,'helios1',cWildHelios1)
	    if (bh2)		     call ForeignFile(iVar,cVar,'helios2',cWildHelios2)
	if(bVcon) then
	  bDaydV = .TRUE.
	  call AskYN('Deconvolve V from one day to the next$yes',bDaydV)
	  if(bDaydV) then
	      if(bVUCSD) aNdayV = 12.
	      call AskR4('Number of days to combine for velocity?',aNdayV)
	      if(bVNago) NmidHRV = -10
	      call AskI4('Number of hours from V UT midnight?',NmidHRV)
	  end if
	  call AskI4('Number of velocity time values?',NTV) ! # NTV time itervals
	  if(NTV.gt.NTmaxV-2) then
	      NTV = NTmaxV-2
	      call AskI4('Number of velocity time values?',NTV) ! Ask again
	  end if
	end if
	if(bGcon) then
	  bDaydG = .TRUE.
	  call AskYN('Deconvolve density or G from one day to the next$yes',bDaydG)
	  if(bDaydG) then
	      call AskR4('Number of days to combine for density or G-level?',aNdayG)
	      if(bGNago) NmidHRG = -10
	      call AskI4('Number of hours from density or G UT midnight?',NmidHRG)
	  end if
	  call AskI4('Number of density or G-level time values?',NTG) ! # NTG time itervals
	  if(NTG.gt.NTmaxG-2) then
	      NTG = NTmaxG-2
	      call AskI4('Number of density or G-level time values?',NTG) ! Ask again
	  end if
	end if
C	I = iFilePath('/home/bjackson/dat/nagoya/',1,'yearly','nagoya.'//cWildChar(:4),cWildNagoya)
C	I = iFilePath(cEnvi//'DAT',1,'Ooty'  ,'Ooty.'  //cWildChar(:4),cWildOoty  )
C-------
C Select the rotation for Helios 1.
C
	if(bHelio) then
          if(bh1) then
	    XCbegHOS1 = max(XClowHOS1,NCoffHOS1+XCbegHOS1)
	    write (STRING,'(2(A,F7.2),A)') '$',XClowHOS1,'$',XChighHOS1,'$0'
	    call ASKR4('Helios 1 rotation (0 = STOP)'//STRING,XCbegHOS1)
	    if (XCbegHOS1 .le. 0) STOP 'Stop in IPSH12Dtd'
	    XCbegHOS1 = min(max(XCbegHOS1,XClowHOS1),XChighHOS1)
	    XCbegHOS1 = XCbegHOS1-NCoffHOS1
C------
C For the moment use only a single Helios rotation
c	    call ASKI4('# Rotations to be averaged',NCHOS)
c	    NCHOS = min(NCHOS,1+int(XChighHIGH-(NCoffHOS1+XCbegHOS1))
	    NCHOS = 1
C
C-------	Info only: start time of Carrington rotation
C
	    I = int(XCbegHOS1)
	    JDCX = JDCarHOS1(I)
	    JDCP = JDCarHOS1(I+NCHOS)
	    NCoffP = NCoffHOS1 + I
C
	    call JULIAN(1,IYR,DOY,JDCX,JEPOCH)
	    IDOY = DOY
	    call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY)
	    JDCX2 = JDCarHOS1(I+1)
	    call JULIAN(1,IYR2,DOY2,JDCX2,JEPOCH)
	    IDOY2 = DOY2
	    call DATE_DOY(1,IYR2,CMON2,IMON2,IDAY2,IDOY2)
	    write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/,30X,A,I2,3A,I4,A,I3,A)') 'Helios 1 rotation ',NCoffP,
     &		' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )',
     &		' and ends on date ',IDAY2,'-',CMON2,'-',IYR2,' ( DOY ',IDOY2,' )'
	    iYrHOS1 = IYR
	end if
C-------
C Select the rotation for Helios 2.
C
          if(bh2) then
	    XCbegHOS2 = max(XClowHOS2,NCoffHOS2+XCbegHOS2)
	    write (STRING,'(2(A,F7.2),A)') '$',XClowHOS2,'$',XChighHOS2,'$0'
	    call ASKR4('Helios 2 rotation (0 = STOP)'//STRING,XCbegHOS2)
	    if (XCbegHOS2 .le. 0) stop 'Stop in IPSH12Dtd'
	    XCbegHOS2 = min(max(XCbegHOS2,XClowHOS2),XChighHOS2)
	    XCbegHOS2 = XCbegHOS2-NCoffHOS2
C------
C For the moment use only a single Helios rotation
c	    call ASKI4('# Rotations to be averaged',NCHOS)
c	    NCHOS = min(NCHOS,1+int(XChighHIGH-(NCoffHOS2+XCbegHOS2))
	    NCHOS = 1
C
C-------	Info only: start time of Carrington rotation
C
	    I = int(XCbegHOS2)
            JDCX = JDCarHOS2(I)
	    JDCP = JDCarHOS2(I+NCHOS)
	    NCoffP = NCoffHOS2 + I
	    call JULIAN(1,IYR,DOY,JDCX,JEPOCH)
	    IDOY = DOY
	    call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY)
	    JDCX2 = JDCarHOS2(I+1)
	    call JULIAN(1,IYR2,DOY2,JDCX2,JEPOCH)
	    IDOY2 = DOY2
	    call DATE_DOY(1,IYR2,CMON2,IMON2,IDAY2,IDOY2)
	  write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/,30X,A,I2,3A,I4,A,I3,A)') 'Helios 2 rotation ',NCoffP,
     &		' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )',
     &		' and ends on date ',IDAY2,'-',CMON2,'-',IYR2,' ( DOY ',IDOY2,' )'
	    iYrHOS2 = IYR
	  end if
C-------
C NEarth is the Earth-based Carrington rotation closest to the Helios rotation
C XCbeg.
	    JD = 0.5d0*(JDCP+JDCX)				! Halfway Helios data range
	    call JULIAN(1,IYR,DOY,JD,JEPOCH)
	    NEarth = N_CARRINGTON(IYR,DOY)			! Nearest Carrington rotation
	    IDOY = DOY
	    call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY)
	    write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A)') 'Nearest IPS rotation ',NEarth,
     &		' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )'
	    XCbeg = NEarth-NCoff-NC/2		! Center IPS range on Helios range
	    write (*,'(/,10X,2(A,F6.1))') 'Carrington range for IPS data ',NCoff+XCbeg,'-',NCoff+XCbeg+NC
C
C Set up file parameters for different data sets output
C
	    IROT = NEarth + 0.00001
	    RDUM = IROT
	    IHGT = RRS*10.0 + 0.00001
	    NCOFH1 = int(XCbegHOS1)
	    NCOFH2 = int(XCbegHOS2)
	    NCOFI  = int(XCbeg)
	    PLUS = 0.0
C
C-------	Get modified Carrington variable for center of plot
C		Get heliographic coordinates for center of plot (Earth position)
	    XCEarth = XCbeg+.5		! Center of plot
	    I = XCEarth			! Rotation containing plot center
	    MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I))
						! Julian day for plot center
	    call JULIAN(1,IYR,DOY,MJDcntr,JEPOCH) ! IYR,DOY for plot center
	    MJDcntr = MJDcntr-MJD_to_JD	! Modified Julian day for plot center
	    call SunNewcomb(0,IYR,DOY,DLONG,DLAT,DDIST)
	    EARTH_LNG = mod(sngl(DLONG)+180.,360.)
	    EARTH_LAT = DLAT
	    call ECLIPTIC_HELIOGRAPHIC(0,IYR,DOY,EARTH_LNG,EARTH_LAT)
						! Heliographic coord. for plot center
	    XCtest1 = XCbeg-XCmargin		! Search limits
	    XCtest2 = XCbeg+NC+XCmargin
	    I = XCtest2
 	    MJDref = JDCar(I)-MJD_to_JD+(XCtest2-I)*(JDCar(I+1)-JDCar(I))
						! Start MJD for data search
	else						! Non-Helios mode
C Calculate range of IPS Carrington variables to work with: [XClo,XChi]
C XClo is based on the first data file available (MJDfrst).
C XChi is based on the current time.
	    if(bVUCSD.and..not.bGcon) then
		Doy = Doy + 30.
	    else
 	         call Julian(11,iYr,Doy,dble(MJDfrst),JEpoch)
	    end if
	print *, 'before JD', iYr,DOY,XClo,xchi 
	    XClo = NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)
	    XClo = XClo-.5
	    call Local2UT(iOffUT,iYr,Doy)
	    XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo)
 	print *, 'before JD', XClo,xchi 
	end if
	RRV = RRS
	RRG = RRS
	if (bVcon) then
	    ANLIMITV = ANLIMITV/(CONRV/CONRD)
	    call AskR4('# of V-points per bin that constitute a deconvolution',ANLIMITV)
	    if(NiterT.gt.8) LimitoV = NiterT/2
	    call AskI4('Iteration at which to limit V-source deviations, -1=none',LimitoV)
C-------
C  Fix the reference distance used for the V map formed.
C  Get E or W or EW; Set the elongation angles
	    call AskWhat('All velocity data, East data, West data ?$0',iEorWV)
	    if (iEorWV .eq. 2) iEorWV = -1
	    call AskR4('Minimum elongation for V-data$0$180$',XElowV)
	    call AskR4('Maximum elongation for V-data$'//cStr(:Flt2Str(XElowV,1,cStr))//'$180$',XEhighV)
	    call AskR4('Maximum RA relative to Sun for V-data$0$180$',XRlimV)
	    RRV = RRS
	    call AskR4('V-map reference distance in AU (0=solar surface)?',RRV)
	    if (RRV .lt. 0) RRV = -RRV*Sun__RAu
	    RRV = max(RRV,sngl(Sun__RAu))
	    call AskYN('Do you want to use a proxy map for velocity?',bVproxy)
	    call AskR4('Maximum limit for V deconvolution?',VLIM)  
	end if
	if (bGcon) then
	    ANLIMITG = ANLIMITG/aNdayG
	    call AskR4('# of D or G-points per bin that constitute a deconvolution',ANLIMITG)
	    if(NiterT.gt.8) LimitoG = NiterT/2
	    call AskI4('Iteration at which to limit D or G-source deviations, -1=none',LimitoG)
C  Get E or W or EW; Set the elongation angles
	    call AskWhat('All D or G data, East data, West data ?$0',iEorWG)
	    if (iEorWG .eq. 2) iEorWG = -1
	    if(bHelio) XElowG = 15.
	    if(bGNago) XElowG = 5.
	    call AskR4('Minimum elongation for D or G-data$0$180$',XElowG)
	    if(bHelio) XEhighG = 180.
	    if(bGNago) XEhighG = 180.
	    call AskR4('Maximum elongation for D or G-data$'//cStr(:Flt2Str(XElowG,1,cStr))//'$180$',XEhighG)
	    if(bHelio) XRlimG = 180.
	    call AskR4('Maximum RA relative to Sun for D or G-data$0$180$',XRlimG)
	    RRG = RRS
	    call AskR4('D or G-map reference distance in AU (0=solar surface)?',RRG)
	    if (RRG .lt. 0) RRG = -RRG*Sun__RAu
	    RRG = max(RRG,sngl(Sun__RAu))
	    call AskYN('Do you want to use a proxy map for density?',bDproxy)
	    call AskR4('Maximum limit for D or G deconvolution?',GLIMP)
	    call AskR4('Minimum limit for D or G deconvolution?',GLIMM)
  	    call AskR4('Factor by which to modify D or G-level?',GLIMF)
	end if
  	    call AskR4('Density time hole filter?',ConsTG)
  	    call AskR4('Velocity time hole filter?',ConsTV)
 
  	    call AskR4('Density spatial filter?',CONRD)
		if(bVNago.and.bGNago) CONRV = CONRD 
  	    call AskR4('Velocity spatial filter?',CONRV)
  	    call AskR4('Density temporal filter?',CONDT)
		if(bVNago.and.bGNago) CONVT = CONDT
  	    call AskR4('Velocity temporal filter?',CONVT)
 		 if(bGCamb) PWG = 0.15
  	    if(.not.bTcon) call AskR4(' G-level power to fit density?',PWG)
 	     if(bTcon) call AskR4('1 AU average density to fit base?',DEN1AU)
		if(bVNago.and.bGNago) PWV = PWG
  	    call AskR4('G-level power to fit velocity?',PWV)
  	    if(.not.bTcon) call AskR4('Radial G-level falloff to fit G density?',PWRG)
	    if(bTcon) call AskR4('Radial Density falloff?',FALLOFF)
 		if(bVNago.and.bGNago) PWRV = PWRG
	    call AskR4('Radial G-level falloff to fit V density?',PWRV)
	    call AskR4('Solar wind speed for G-level traceback?',Speed)
	RRD = RRG
C-------	Prompts for time
	if (bForecast) then			! Forecast mode
	    call Local2UT(iOffUT,iYr,Doy)		! Doy = universal time
	    call Julian(10,iYr,Doy,JD,JEpoch)		! Current time in MJD
	    write (*,'(/,10X,A,I6,A,I6,A,F11.3)') 'First MJD =',MJDfrst,
     &			'   Last MJD =',MJDlast,'   Current MJD =',JD
	    MJDref = JD					! MJD at center of plot (Earth position)
	    write (cStr,'(A,I5,A,F7.1,A)') '$',MJDfrst,'$',JD,'$0$'
	    call AskR8('Modified Julian Day (0 = STOP)'//cStr,MJDref)
 	    if (MJDref .le. 0) call Say('ipshtd','I','StopPlay','StopPlay')
C-------	Make sure iYr,Doy match MJDref
C		Get heliographic coordinates for center of plot (Earth position)
	    call Julian(11,iYr,Doy,MJDref,JEpoch)
	    call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun)
	    EarthLng = mod(sngl(dLngSun)+180.,360.)
	    EarthLat = -dLatSun
	    call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat)
C-------	Get modified Carrington variable for center of plot
C		[Note that EarthLng=(1-frac(XCEarth))*360]
	    XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)
	    XCbeg   = XCEarth-.5
	    XCtest1 = XCbeg-XCmargin		! Search limits
	    XCtest2 = XCbeg+NC+XCmargin
	    XCend= XCbeg + NC
	    XCstrt = XCbeg
	    XCbeg= XCbeg - xInc
	    XCend= XCend + xInc
	    MJDcntr = MJDref
	else						! Normal mode
	    if(bVUCSD.and..not.bHelio) XCbeg = 1661.6
	    if(bHelio) then
		XClo  = NCoff+XCbeg-1
		XChi  = NCoff+XCbeg+1
	    end if
	    write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$'
	    if(bHelio) XCbeg = 1661.6
	    if(bVNago) XCbeg = 1941
	    if(bGCamb) XCbeg = 1884
	    call AskR4('Carrington rotation (0 = STOP)'//cStr,XCbeg)
	    if (XCbeg .le. 0) call Say('ipshtd','I','StopPlay','StopPlay')
	    NintF = 1
	    if(bForecast) then
	        nTF = nTG - 1
	        NintF = 0
	    end if
	    call AskI4('# Intermediate rotations interpolated',NintF)	
	    call AskI4('# Rotations averaged',NC)
	    NC = min(NC,1+int(XChi-XCbeg) )	! Stay within available data range
	    XCbeg = XCbeg-NCoff
C-------	Info only: start time of Carrington rotation
	    call Julian(1,iYr,Doy,JDCar(int(XCbeg)),JEpoch)
	    iDoy = Doy
	    call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy)
	    write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+int(XCbeg),
     &		' starts on date ',iDay,'-',cMon,'-',iYr,' ( Doy ',iDoy,' )'
C-------	Get modified Carrington variable for center of plot
C		Get heliographic coordinates for center of plot (Earth position)
	    XCEarth = XCbeg+.5			! Center of plot
	    I = XCEarth			! Rotation containing plot center
	    MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I))
						! Julian day for plot center
	    call Julian(1,iYr,Doy,MJDcntr,JEpoch) ! iYr,Doy for plot center
	    MJDcntr = MJDcntr-Sun__MJDtoJD	! Modified Julian day for plot center
	    call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun)
	    EarthLng = mod(sngl(dLngSun)+180.,360.)
	    EarthLat = -dLatSun
	    call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat)
						! Heliographic coord. for plot center
	    XCtest1 = XCbeg-XCmargin		! Search limits
	    XCtest2 = XCbeg+NC+XCmargin
	    XCend= XCbeg + NC
	    XCstrt = XCbeg
	    XCbeg= XCbeg - xInc
	    XCend= XCend + xInc
	    I = XCtest2
	    MJDref = JDCar(I)-Sun__MJDtoJD+(XCtest2-I)*(JDCar(I+1)-JDCar(I))
						! Start MJD for data search
 	end if
C	write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCbeg
	call AskYN('Do you want to output model-source files$no',bWriteMS)
	call AskYN('Do you want one shift matrix per iteration$no',bOnly1)
        call AskYN('Output a Three D file$yes',bWrite3D) 	  	! 3D file write
        call AskYN('Output an Earth parameter file$yes',bEWRITE)	! Earth extraction
        call AskYN('Output a Mercury parameter file$no',bMWRITE)	! Mercury extraction
        call AskYN('Output a Ulysses parameter file$no',bUWRITE)	! Ulysses extraction
        call AskYN('Output a Helios 1 parameter file$yes',bH1WRITE)	! Helios 1 extraction
        call AskYN('Output a Helios 2 parameter file$yes',bH2WRITE)	! Helios 2 extraction
C Do you want to include a magnetic map at the reference surface height?
C If yes, then read the source surface map (at 2.5Rs) at
C the resolution of the reference arrays.
        call AskYN('Do you want to make a magnetic map?$yes',bMagnetic)
C-------
C The range of the Carrington variable plotted is [XCbeg,XCend]. The  target range
C [XCtest1,XCtest2] includes an extra safety margin to make sure that all relevant
C data are read from file.
 	
	if (bVNago) call ReadVIPSn(iReadNagoyan,iProcessNagoyan,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,SRCV,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
	if (bVOoty) call ReadVIPSn(iReadOotyn, iProcessOotyn,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,SRCV,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav)
	if (bVUCSD) call ReadVIPSn(iReadUCSDn,iProcessUCSDn,cWildUCSD,
     &          nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &          RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &          NLmaxV,NLV,SRCV,SRCVsav,
     &          IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &          1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &          XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav)
	if (bGCamb) call ReadGIPS('$CAM',nCar,JDCar,
     &		iEdt,.FALSE.,bForecast,
     &		XCtest1,XCtest2,MJDref,MJDfrst,MJDlast,
     &		RRG,Speed,GPOWER,XElowG,XEhighG,iEorWG,XRlimG,
     &		 NLmaxG,NLG,
C     &		NLmaxG,NLG,SRCGG,SRCGsav,
     &		iXPG,iYPG,iMJDG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,GGOBS,
     &		1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav,XDSsav,XLSsav,
     &		XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
	if (bGOoty) call ReadVIPSn(iReadOotyn, iProcessOotyn,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
c     &		RRV,XElowG,XEhighG,iEorWG,XRlimG,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxV,NLG,SRCGG,SRCGsav,
c     &		NLmaxV,NLV,SRCGG,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS,
c     &		IYRFG,IRECG,IYRSGG,DOYSGG,DISTV,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
	if (bGNago) call ReadVIPSn(iReadNagoyan,iProcessNagoyan,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxV,NLG,SRCGG,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
C
C  Set up the next two limits for the deconvolution just in case they are used
C
	if(bh1) then
		XCbegHOS = XCbegHOS1
	else
		XCbegHOS = XCbegHOS2
	end if
C-------
C If HELIOS 1 data is to be deconvolved, read the HELIOS 1 data points.
C
	if(bHelio) then
	NS = 0
	iHOS = 10
	  if(bh1) then
	    iSc = 1
	    LT1 = 2
	    NLT1 = 0
	    call READ_HOSF(cWildHelios1,nCarHOS1,JDCarHOS1,XCbegHOS1,NCHOS,iYrHOS1,
     &		RRS,Speed,0.,                       ! Radius,Speed,Eref
     &		XElowG,XEhighG,iEorWG,XRlimG,
     &		iSc,iHOS,cHOS,iSectBeg,iSectEnd,LT1,iDayHOS,nDevHOS,dSigHOS,
     &		NLmaxG,NS,NLT1,SRCGG,
     &		DISTGG,XCGG,YLGG,GGOBS,XEGG,XLSGG,XCEGG,XLLGG,XDLGG,IYRSGG,DOYSGG)
	     NS = NLT1
	     do I=1,NLT1
		XCEARTH = XMAP_SC_POS(EARTH,IYRSGG(I),DOYSGG(I),nCar,JDCar)
		XCEARGG(I) = XCEARTH		! Carrington variable of Earth at Helios obs. time
		H1LNG = (1.-(XCEGG(I)-int(XCEGG(I))))*360.0
		XCEGG(I) = XMAP_OBS_POS(XCEARTH,H1LNG)	! Earth-based Carrington variable of H1
C	Print *, 'H1 I, XCEARGG, XCEGG, Diff', I, XCEARGG(I), XCEGG(I), XCEARGG(I)-XCEGG(I)
	     end do
          end if
C-------
C If HELIOS 2 data is to be deconvolved, read the HELIOS 2 data points.
C
	  if(bh2) then
	    iSc = 2
	    LT2 = 2
	    call READ_HOSF(cWildHelios2,nCarHOS2,JDCarHOS2,XCbegHOS2,NCHOS,iYrHOS2,
     &		RRS,SPEED,0.,                       ! Radius,Speed,Eref,
     &		XElowG,XEhighG,iEorWG,XRlimG,
     &		iSc,iHOS,cHOS,iSectBeg,iSectEnd,LT2,iDayHOS,nDevHOS,dSigHOS,
     &		NLmaxG,NS,NLT2,SRCGG,
     &		DISTGG,XCGG,YLGG,GGOBS,XEGG,XLSGG,XCEGG,XLLGG,XDLGG,IYRSGG,DOYSGG)
	     do I=NLT1+1,NLT2
		XCEARTH = XMAP_SC_POS(EARTH,IYRSGG(I),DOYSGG(I),nCar,JDCar)
		XCEARGG(I) = XCEARTH		! Carrington variable of Earth at Helios obs. time
		H2LNG = (1.-(XCEGG(I)-int(XCEGG(I))))*360.0
		XCEGG(I) = XMAP_OBS_POS(XCEARTH,H2LNG)	! Earth-based Carrington variable of H2
C	Print *, 'H2 I, XCEARGG, XCEGG, Diff', I, XCEARGG(I), XCEGG(I), XCEARGG(I)-XCEGG(I)
	     end do
          end if
	  if (bh1) NLG = NLT1
	  if (bh2) NLG = NLT2
C
C XRN = heliocentric distance P (AU)
C XEN = Topocentric elongation P (deg)
C
          write (*,'(/,A,I6,A,I6,A,I6,A/)') ' There were', NLT1, ' H1 and',
     &			NLT2-NLT1,' H2 for a total of', NLG,' Helios lines of sight'
C
	  NLST = 1
	  if (bh1) then
	    call UBVConst(LT1,U,APM)
	print *, RRS, Sun__RAu, APM
	    xTemp = ThomsonlosS10Far(0.,-1.,sngl(RRS/Sun__RAu/sind(30.0)),30.0,APM,P)
	    print *, 'xTemp = ',xTemp, P
            do I=NLST,NLT1
	      if(GGOBS(I).ne.BADR4()) then 
		call UBVConst(LT1,U,APM)
		xTHOM = ThomsonlosS10Far(0.,-1.,sngl(DISTGG(I)/Sun__RAu/ABS(XEGG(I))),ABS(XEGG(I)),APM,P)
		ZZZ(I) = GGOBS(I)*xTemp/xTHOM
C	        if(I.ge.2000.and.I.le.2100) print *, I,ZZZ(I),XLSGG(I),XCEGG(I),XLLGG(I),DOYSGG(I),GGOBS(I)
C	        if(I.ge.NLT1+1800.and.I.le.NLT1+1900) print *, I,ZZZ(I),XLSGG(I),XCEGG(I),XLLGG(I),DOYSGG(I),GGOBS(I)
	      else
	        print *, 'Bad Brightness, xTemp, xTHOM', I, GGOBS(I), xTemp, xTHOM
		ZZZ(I) = BadR4()
	      end if
            end do
	    NLST = NLT1+1
	  end if
	  if(bh2) then
	    call UBVConst(LT2,U,APM)
	    xTemp = ThomsonlosS10Far(0.,-1.,sngl(RRS/Sun__RAu/sind(30.0)),30.0,APM,P)
	    print *, 'xTemp =',xTemp
            do I=NLST,NLG
	      if(GGOBS(I).ne.BADR4()) then 
	        call UBVConst(LT2,U,APM)
	        xTHOM = ThomsonlosS10Far(0.,-1.,sngl(DISTGG(I)/Sun__RAu/ABS(XEGG(I))),ABS(XEGG(I)),APM,P)
	        ZZZ(I) = GGOBS(I)*xTemp/xTHOM
C	        if(I.ge.2000.and.I.le.2100) print *, I,ZZZ(I),XLSGG(I),XCEGG(I),XLLGG(I),DOYSGG(I),GGOBS(I)
C	        if(I.ge.NLT1+1800.and.I.le.NLT1+1900) print *, I,ZZZ(I),XLSGG(I),XCEGG(I),XLLGG(I),DOYSGG(I),GGOBS(I)
	      else
	        print *, 'Bad GGOBS(I) xTemp, xTHOM', I, GGOBS(I), xTemp, xTHOM
	        ZZZ(I) = BadR4()
	      end if
	    end do
           end if
	end if
	call ArrR4Zero(NLmaxV,VSSIG)
	call ArrI4Constant(NLmaxV,1 ,NBSV )
	call ArrR4Zero(NLmaxG,GSSIG)
	call ArrI4Constant(NLmaxG,1 ,NBSG )
	do I=1,NLmaxV
C	    if (I.lt.100) print *, I, SRCV(I)
	    if (VOBS(I) .gt. VLIM) NBSV(I) = 0
	    if (VOBS(I) .lt. VLIL) NBSV(I) = 0
	end do
	ICntBadr4 = 0
	do I=1,NLmaxG
C	    if (I.lt.200) print *, I, SRCGG(I), GGOBS(I)
C	    print *, I, SRCGG(I), GGOBS(I)
C	    if (SRCGG(I).eq.' 1  2   2') print *, I, SRCGG(I), GGOBS(I)
	    if (GGOBS(I) .eq. BadR4()) then 
		NBSG(I) = 0
		ICntBadr4 = ICntBadR4 + 1
	    end if		
	    GGOBS(I) = GGOBS(I)*GLIMF
	    if (GGOBS(I) .gt. GLIMP) NBSG(I) = 0
	    if (GGOBS(I) .lt. GLIMM) NBSG(I) = 0
	end do
	if(ICntBadr4.gt.0) print *, ' The number of bad sources are = ', ICntBadr4
				
C=======
	do I=1,NLG
	  if(bGCamb) then
	    XCG  (NLG-I+1)	= XCGG(I)
	    YLG  (NLG-I+1)	= YLGG(I)
	    DISTG(NLG-I+1)	= DISTGG(I)
	    GOBS (NLG-I+1)	= GGOBS(I)
	    XEG  (NLG-I+1)	= XEGG(I)
	    XLSG (NLG-I+1)	= XLSGG(I)
	    XCEG (NLG-I+1)	= XCEGG(I)
	    XCEARG(NLG-I+1)	= XCEGG(I)
	    XLLG (NLG-I+1)	= XLLGG(I)
	    XDLG (NLG-I+1)	= XDLGG(I)
	    DOYSG(NLG-I+1)	= DOYSGG(I)
	    IYRSG(NLG-I+1)	= IYRSGG(I)
	    XLLGD5 		= XLLG(NLG-I+1)/5.0
	    XDLGD5		= XDLG(NLG-I+1)/5.0
	    if(XLLGD5.gt.0.) then
		XLLGD5 = XLLGD5 + 0.5
	    else
		XLLGD5 = XLLGD5 - 0.5
	    end if
	    if(XLLGD5.gt.0.) then
		XDLGD5 = XDLGD5 + 0.5
	    else
		XDLGD5 = XDLGD5 - 0.5
	    end if
	    IXLLGD5 = XLLGD5
	    IXDLGD5 = XDLGD5
	    write(SRCGS(1:5),'(I5)') IXLLGD5*5
	    write(SRCGS(6:9),'(I4)') IXDLGD5*5
	    SRCG(NLG-I+1) = SRCGS
	  end if
	  if(bGNago.or.bGOoty) then
	    XCG(I)	= XCGG(I)
	    YLG(I)	= YLGG(I)
	    DISTG(I)	= DISTGG(I)
	    GOBS(I)	= GGOBS(I)
	    XEG(I)	= XEGG(I)
	    XLSG(I)	= XLSGG(I)
	    XCEG(I)	= XCEGG(I)
	    XCEARG(I)	= XCEGG(I)
	    XLLG(I)	= XLLGG(I)
	    XDLG(I)	= XDLGG(I)
	    DOYSG(I)	= DOYSGG(I)
	    IYRSG(I)	= IYRSGG(I)
	  end if
	end do
C
C  Sort the Helios data into consecutive XCEARG times
C
	if(bHelio) then
C
C Find minimum and next to minimum Helios XCEARGG
 (Modified 2/26/02 BVJ)
C
	  XCeamin  = 10000.
	  XCeaminp = 10000.
	  do K=1,NLG
	      if(XCEARGG(K).lt.XCeamin) XCeamin = XCEARGG(K)
	  end do	         
	  do K=1,NLG
	      if(XCEARGG(K).lt.XCeaminp.and.XCEARGG(K).gt.XCeamin) XCeaminp = XCEARGG(K)
	  end do

	print *, 'XCeamin, XCeaminp',XCeamin, XCeaminp
	  K = 1
	  do K2=1,NLG
	      do K1=1,NLG
		   if(XCEARGG(K1) .eq. XCeamin .and. K .le. NLG) then
       			    XCG   (K)	= XCGG(K1)
     			    YLG   (K)	= YLGG(K1)
     			    DISTG (K)	= DISTGG(K1)
  		            GOBS  (K)	= GGOBS(K1)
			    ZZ    (K)	= ZZZ(K1)
    			    XEG   (K)	= XEGG(K1)
     			    XLSG  (K)	= XLSGG(K1)
			    XCEARG(K)	= XCEARGG(K1)
     		   	    XCEG  (K)	= XCEGG(K1)
     		            XLLG  (K)	= XLLGG(K1)
            		    XDLG  (K)	= XDLGG(K1)
    			    DOYSG (K)	= DOYSGG(K1)
     			    IYRSG (K)	= IYRSGG(K1)
			    SRCG  (K)	= SRCGG(K1)
C	if(K.lt.200) print *, K,K1,SRCGG(K1),SRCG(K)
,XCEARGG(K1),XCeamin,NLG
C	    if (SRCG(K).eq.' 1  1  14') print *, K, SRCG(K), GOBS(K)
C	    if (SRCG(K).eq.' 1  3 -90') print *, K, SRCG(K), GOBS(K)
C	    if (SRCG(K).eq.' 2  3 +90') print *, K, SRCG(K), GOBS(K)
C	print *, K, K1, XCeamin, XCeaminp, XCEARG(K), XEG(K)
			    K = K + 1
		  end if
	      end do
	      XCeamin = XCeaminp
	      XCeaminp = 100000.	         
	      do K3=1,NLG		! Modified 2/26/02 BVJ
	          if(XCEARGG(K3).lt.XCeaminp.and.XCEARGG(K3).gt.XCeamin) XCeaminp = XCEARGG(K3)
	      end do
	  end do
	  print *, 'The total number of Helios 1 and 2 lines of sight is',NLG
	end if								! End of helios 1 and 2 sort
       ALng = 1.*(nLng-1)/(iLng-1)
       if(bForecast) then 				! For forecast mode
	  NTV = (NTV + 1.5)/2.0 + 6
	  NTG = (NTG + 1.5)/2.0 + 6
       end if
       if(bVcon)
     &  call MkTimesvd(bDaydV,bForecast,MJDcntr,NmidHRV,aNdayV,XCtest1,XCtest2,xInc,NTV,NTmaxV,
     &  nCar,JDCar,ALng,XCintV,XCbeV,XCtbegV,XCtendV,XCtioffV,RconstV,ItoffV,IYRV,DOYV)
       if(bGcon)
     &  call MkTimesvd(bDaydG,bForecast,MJDcntr,NmidHRG,aNdayG,XCtest1,XCtest2,xInc,NTG,NTmaxG,
     &	nCar,JDCar,ALng,XCintG,XCbeGG,XCtbegG,XCtendG,XCtioffG,RconstG,ItoffG,IYRG,DOYG)
       if(.not.bGcon) then
	   nTG = nTV
	   ItoffG = ItoffV
	   XCtbegG = XCtbegV
	   XCtendG = XCtendV
	   do I=1,nTmaxG
	       XCintG(I) = XCintV(I)
	       XCbegG(1,I) = XCbeV(1,I)
	       XCbegG(2,I) = XCbeV(2,I)
	   end do
       end if
       if(.not.bVcon) then
	   nTV = nTG
	   ItoffV  = ItoffG
	   XCtbegV = XCtbegG
	   XCtendV = XCtendG
	   do I=1,nTmaxV
	       XCintV(I)  = XCintG(I)
	       XCbeV(1,I) = XCbegG(1,I)
	       XCbeV(2,I) = XCbegG(2,I)
	   end do
       end if
       if(bOnly1) then
	   nTV = nTG
	   ItoffV  = ItoffG
	   XCtbegV = XCtbegG
	   XCtendV = XCtendG
	   do I=1,nTmaxV		! The temporal spaces for XCintV, etc. better be the same as for XCintG 
	       XCintV(I)  = XCintG(I)
	       XCbeV(1,I) = XCbegG(1,I)
	       XCbeV(2,I) = XCbegG(2,I)
	   end do
	end if    
       print *, 'XCtbegV, XCtendV, New NTV', XCtbegV, XCtendV, NTV
       print *, 'XCtbegG, XCtendG, New NTG', XCtbegG, XCtendG, NTG
       print *, XCbeV, XCintV, ItoffV
       print *, XCbeGG, XCintG, ItoffG
       XCtbegF = XCtbegG				! For forecast mode
       XCtendF = XCtendG				! For forecast mode
       if(bGcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP)
       if(bVcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),VMAP)
       NLmax = max(NLmaxG,NLmaxV)
       NTVV = ItoffV
       NTGG = ItoffG
       NTVE = 0
       NTGE = 0
       K    = 0
       do KKK=1,NLmax+nTG		! Loop before deconvolution to set up VMAP, DMAP
	K = K + 1
	IVSW = 0			! Switch turn on for NTVV value of K
	if(K.le.NLV+1) then	! If K is more than the source # + 1 don't turn on switch
	  if(bVcon.and.XCEV(K).gt.XCintV(NTVV+2).and.XEV(K).le.XEhighV) then
	    NTVV = NTVV + 1
	    NTVB = NTVE + 1
	    NTVE = K - 1
	    if(NTVB.ne.NTVE+1) IVSW = 1
	    if(K.le.NLV) K = K - 1
	    NTVD = NTVE-NTVB+1
	    write (*,'(A,I4,A,2I6,2F10.3)') 'V loop ',NTVV,' - S#, K, XCEV, XCintV', NTVD,K,NCoff+XCEV(K),NCoff+XCintV(NTVV+2)
	    if(NTVV.eq.NTV) NLV = NTVE
	  end if
	end if
	if(IVSW.eq.1) then
C	print *, ' V switch ',IVSW,'  loop',NTVV
	    if (NLV .eq. 0) call Say('ipshtd','E','NoV','no velocity IPS data')
	    if (bVproxy) then			! Read proxy point-P V map from file
		I = iReadProxyMapN(1,NCoff+XCintV(NTVV),nLng,nLat,VMAP(1,1,NTVV))
		call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax)
		print *, 'The VMAP',NTVV,' min and max after read is', amin, amax
	    end if
	    call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),aminV,amaxV)
	    if (amaxV.eq.badR4()) then				! Try to calculate a proxy point-P V map
		call MapGrid(XCbeV(1,NTVV),XCbeV(2,NTVV),YLbeg,YLend,NC,
     &			NTVD,XCV(NTVB),YLV(NTVB),VOBS(NTVB),NAV,DAV,nLng,nLat,
     &			VMAP(1,1,NTVV),VDEV(1,1,NTVV),VPNT(1,1,NTVV),
     &			Vmin,Vmax,VDmin,VDmax)
C		write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintV(NTVV)
	    end if
						! Calculate weights, WTSV, along all lines of sight
	    call MkLOSWeights(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) 
	 call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAPHOLE(1,1,NTVV))
	 if(amaxV.eq.BadR4()) then
	     call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax)
	     if(amax.ne.BadR4()) then
		call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,NTVV),CONRV,3,0.0,90.0)
	    end if
	 end if
	 call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax)
	 print *, 'The VMAP' ,NTVV,' min and max is', amin, amax
	 if(NTVV.eq.ItoffV+1) then				! Fill earlier arrays
	     do N=1,ItoffV
		amaxV = BadR4()
		if(bVproxy) then
		    I = iReadProxyMapN(1,NCoff+XCintV(N),nLng,nLat,VMAP(1,1,N))
		    call arrR4getminmax(nLngLat,VMAP(1,1,N),aminV,amaxV)
		    print *, 'The VMAP' ,N,' min and max after read is', aminV, amaxV
		end if
		if(amaxV.eq.BadR4()) then
		    call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax)
		    print *, 'The VMAP' ,N,' min and max is', amin, amax
		    if(amax.ne.BadR4()) then
			call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAP(1,1,N))
			call Say('ipshtd','I','Info','Copy first source array into data file ' 
     &		//cStr(:Int2Str(N,cStr)))
		    end if
		end if
	     end do
	 end if
        end if
	IGSW = 0					! Switch turn on for NTGG value of K
	if(K.le.NLG+1) then	! If K is more than the source # +1 don't turn on switch
	 if(K.gt.0) then
	  if (.not.bHelio.and.bGcon.and.XCEG(K).gt.XCintG(NTGG+2).and.XEG(K).le.XEhighG) then
	    NTGG = NTGG + 1
	    NTGB = NTGE + 1
	    NTGE = K - 1
	    if(NTGB.ne.NTGE+1) IGSW = 1
	    if(K.le.NLG) K = K - 1
	    NTGD = NTGE-NTGB+1
	    if(.not.bTcon) write (*,'(A,I4,A,2I6,2F10.3)') 
     &		'G loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2)
	    if(bTcon) write (*,'(A,I4,A,2I6,2F10.3)') 
     &		'B loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2)
	    if(NTGG.eq.NTG) NLG = NTGE
	  end if
	  if (bHelio.and.bGcon.and.XCEARG(K).gt.XCintG(NTGG+2)) then
	    NTGG = NTGG + 1
	    NTGB = NTGE + 1
	    NTGE = K - 1
	    if(NTGB.ne.NTGE+1) IGSW = 1
	    if(K.le.NLG) K = K - 1
	    NTGD = NTGE-NTGB+1
	    write (*,'(A,I4,A,2I6,2F10.3)') 
     &		'B loop ',NTGG,' - S#, K, XCEARG, XCintG', NTGD,K,NCoff+XCEARG(K),NCoff+XCintG(NTGG+2)
	    if(NTGG.eq.NTG) NLG = NTGE
	  end if
	 end if
	end if
	if(IGSW.eq.1) then
	 if(NTGB.ne.0) then
	    if (NLG .eq. 0) call Say('ipshtd','E','NoG','no Density or G-level IPS data')
	    if (bDproxy) then			! Read proxy density map from file
		I = iReadProxyMapN(2,NCoff+XCinTG(NTGG),nLng,nLat,DMAP(1,1,NTGG))
	    end if
 	     call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),aminDMAP,amaxDMAP)
	    if (amaxDMAP.eq.badR4()) then		! Try to calculate a point-P synoptic G map
C		write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintG(NTGG)
	      if (.not.bHelio) then 
		call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC,
     &			NTGD,XCG(NTGB),YLG(NTGB),GOBS(NTGB),NAV,DAV,nLng,nLat,
     &			DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG),
     &			Gmin,Gmax,GDmin,GDmax)
	      else
		call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC,
     &			NTGD,XCG(NTGB),YLG(NTGB),ZZ(NTGB),NAV,DAV,nLng,nLat,
     &			DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG),
     &			Gmin,Gmax,GDmin,GDmax)
	      end if
	    end if

C	do NLLL=NTGB,NTGB+NTGD
C	    if (NLLL.lt.200) print *, NLLL, SRCG(NLLL), GOBS(NLLL)
C	    if (SRCG(NLLL).eq.' 1  1  14') print *, NLLL, SRCG(NLLL), GOBS(NLLL)
C	end do

C  Calculate weights, WTSG, according to NLOSWG along all lines of sight .
	    call MkLOSWeights(NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG)
	    if(.not.bTcon) then		! Deal with g-level DMAP
	      call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),PPGMAP(1,1,NTGG))
C  Convert G to density and smooth density maps if no proxy density maps are input.
	      if(amaxDMAP.eq.badR4()) then
	        call ConvertG2D(0,nLngLat,PPGMAP(1,1,NTGG),DMAP(1,1,NTGG),RRG,DEN1AU,PWG)
 	         call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax)
	        if(amax.ne.BadR4()) then
		  call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,NTGG),CONRD,3,0.0,90.0)
	        end if
	        print *, 'After G to Den convert the map min and max are', amin, amax
	      end if
	    else			! Deal with Thomson scattering DMAP
C  Smooth Thomson density maps if no proxy density maps are input.
	      if(amaxDMAP.eq.badR4()) then
	        call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax)
	        if(amax.ne.BadR4()) then
		  call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,NTGG),CONRD,3,0.0,90.0)
	        end if
	        print *, 'After Thomson to Den convert the map min and max are', amin, amax
	      end if
	    end if
	 call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax)
	 print *, 'The DMAP',NTGG,' min and max is', amin, amax
C========
	 call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAPHOLE(1,1,NTGG))	! Dummy density map for print out
	 call ArrR4TimesConstant(-nLngLat,DMAPHOLE(1,1,NTGG),RRS*RRS,DMAPHOLE(1,1,NTGG))	! n*R^2 with holes possible
	 if(NTGG.eq.ItoffG+1) then				! Fill earlier arrays
	     do N=1,ItoffG
		amaxD = BadR4()
		if(bDproxy) then
		    I = iReadProxyMapN(2,NCoff+XCintG(N),nLng,nLat,DMAP(1,1,N))
		    call arrR4getminmax(nLngLat,DMAP(1,1,N),aminD,amaxD)
		    print *, 'The DMAP',N,' min and max after read is', aminD, amaxD
		end if
		if(amaxDMAP.eq.BadR4()) then
		    call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax)
 		     print *, 'The DMAP',N,' min and max is', amin, amax
		    if(amax.ne.BadR4()) then
		        call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAP(1,1,N))
			call Say('ipshtd','I','Info','Copy first source array into data file ' 
     &		//cStr(:Int2Str(N,cStr)))
		    end if
		end if
	     end do
	 end if
	 end if
	end if
       end do			! End arrays that set up initial source values and maps
	 amaxallV = BadR4()
	 do N=1,NTV
	     call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
	     if(amax.ne.BadR4()) amaxallV = amax
	 end do
	 if(amaxallV.eq.BadR4()) then
	     do N=1,NTV
	         call ArrR4Constant(nLngLat,Speed,VMAP(1,1,N))
	     end do
	     print *, 'All velocity maps were bad and all were set to',Speed
	 end if
	 amaxallD = BadR4()
	 do N=1,NTG
	     call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
	     if(amax.ne.BadR4()) amaxallD = amax
	     amaxallD = BadR4()
	 end do
	 if(amaxallD.eq.BadR4()) then
	     TD2 = 1.0						   ! Initial total Gvalue for g-level maps
	     if(bTcon) 	TD2 = 0.66667*DEN1AU*((R1AU/RRD)**FALLOFF) ! divide total density by two-thirds for initial Thomson maps
	     do N=1,NTG
	         call ArrR4Constant(nLngLat,TD2,DMAP(1,1,N))
	     end do
	     print *, 'All density maps were bad and all were set to',TD2
	 end if
	 if(.not.bVcon) then
	    do N=1,NTV
		call ArrR4Constant(nLngLat,Speed,VMAP(1,1,N))
		call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N))
		call ArrR4Copy(nLngLat,VMAP(1,1,N),DMAPV(1,1,N))
	   end do
	   call arrR4getminmax(nLngLat,VMAP(1,1,NTV),amin,amax)
	   print *, 'All',NTV,' VMAPs beginning with .not. bVcon are', amax
	end if
 	 if(.not.bGcon) then
	  if(.not.bTcon) then
	    do N=1,NTG
		call ArrR4Constant(nLngLat,1.,DMAP(1,1,N))
		call ArrR4Copy(nLngLat,DMAP(1,1,N),PPGMAP(1,1,N))
		call ConvertG2D(0,nLngLat,PPGMAP(1,1,N),DMAP(1,1,N),RRG,DEN1AU,PWG)
		call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N))
		call ArrR4Copy(nLngLat,DMAP(1,1,N),VMAPD(1,1,N))
	    end do
	    call arrR4getminmax(nLngLat,DMAP(1,1,NTG),amin,amax)
	    print *, 'All',NTG,' beginning DMAPs for .not. bGcon are', amax
	  end if
	end if
	NLGG = 1
	do I=1,NLG
	    if(WTSG(1,I).ne.0.0) NLGG = NLGG + 1
	end do
	
	NLG = NLGG - 1
	call Say('ipshtd','I','Info',' D observations '
     &		//cStr(:Int2Str(NLG,cStr)))
	call Say('ipshtd','I','Info',' V observations '
     &		//cStr(:Int2Str(NLV,cStr)))
	call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD)
	call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV)
C Produce the initial Thomson scattering base maps if needed
	if (bTcon) then
	  DB = DEN1AU*((R1AU/RRD)**FALLOFF)
	  DB2 = DB/3.0	
	  do N=1,NTG
	      do i=1,nLng
		  do j=1,nLat
		      DBASE(i,j,N)  = DB2  		! Produce initial base map
		      DMAP(i,j,N)   = DMAP(i,j,N) + DB2	! Add initial base to DMAP
		  end do
	      end do
	  end do
	end if
C	do I=1, NLG
C	    if(I.lt.200) print *, I, SRCG(I), GOBS(I)
C	    if (SRCG(I).eq.' 1  1  14') print *, I, SRCG(I), GOBS(I)
C	end do
	call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &			nLng,nLat,NTG,NTV,ConsTG*aNdayG,DMAP,DMAPV,DDT,DDD)
	if(bVproxy) call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV)
	if(bDproxy) then		! Do this for use after velocity iteration
	    call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,
     &			nLng,nLat,nTV,nTG,ConsTV*aNdayV,VMAP,VMAPD,VVT,VVV)
	    call FillWholeT(1,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD)
	end if
       if(.not.bVproxy) then
	call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV)
	do N=1,NTV
C	  call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax)
C	  print *, 'DMAPV min max N',amin,amax,N
C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'VMAP  min max N',amin,amax,N
	  call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
	  call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
	  call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,90.0)
	  call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
	  call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
	  call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,90.0)
C	  call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax)
C	  print *, 'DMAPV min max N',amin,amax,N
C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'VMAP  min max N',amin,amax,N
	end do
	call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp)
	call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp)
	do N=1,NTV
	  call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
	  call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
	end do
       end if

C        call MkShiftd(XCbeV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,
C     &	VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
C     &	NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

        call MkShiftdn(XCbeV,XCtbegV,XCtendV,ALng,aNdayV,nLng,nLat,nMap,nTV,
     &	VMAP,DMAPV,XCshift,RconstV,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
     &	NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

	do N=1,NTV
C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'VMAP  min max N',amin,amax,N
C	  call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax)
C	  print *, 'DMAPV min max N', amin,amax,N
	  call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,1,0.0,90.0)
	end do
C
C Comment out   call MkDensity(nLng,nLat,nMap,XCshift,DDfact,VMAP,RRV,dRR)
C
	Nit  = -1			! Why three iterations counters??
	NitV = -1
	NitG = -1
	do while (Nit .lt. NiterT)
	    Nit = Nit+1
	    print *, ' '
	    call Say('ipshtd','I','Info','Iteration '//cStr(:Int2Str(Nit,cStr)))
	    NVS = 0
	    if (Nit .eq. LimitoV) NVS = 1
	    NGS = 0
	    if (Nit .eq. LimitoG) NGS = 1
	    if (bVcon) then				! V deconvolution
		NitV = NitV+1				! # V iterations
							! bGmv2=.TRUE., implies bGcon=.FALSE.
		if (bGmv2) then				! Make a density map assuming mv^2 = C
		  RRD = RRV
		  call Veleq3(nLng,nLat,NTV,VMAP,SPEEDEQ)
		  print *, 'Solar equatorial mv^2 = constant value SPEEDEQ =',SPEEDEQ 
		  do N=1,NTV
		    call Mk_V2DN(nLngLat,VMAP(1,1,N),DMAPV(1,1,N),SPEEDEQ,DEN1AU,RRD)
		  end do
		  call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,
     &				    nLng,nLat,NTV,NTG,ConsTV*aNdayV,DMAPV,VMAPD,VVT,VVV)
	    	  call ArrR4Copy(nLngLat,VMAPD(1,1,N),DMAP(1,1,N))
		  if(NitV .eq. NiterT) then
		    do N=1,NTG
	    	      call ArrR4Copy(nLngLat,VMAPD(1,1,N),DMAPHOLE(1,1,N))
		    end do
		  end if
		end if
		if(bstopw) then
		    print *, 'Before MkPostd V, V deconvolution'
		    call stopwatch(' ',' ')
		end if
		call MkPostd(XCbeV,XCtbegV,XCtendV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV,
     &			XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XCEV,XLLV,XDLV,IYRSV,DOYSV,
     &			XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM)
		  call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV,
     &			DDfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,DFAC)
		  call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV,
     &			DVfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,VFAC)
 		call MkVModeltd(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX,
     &			WTSV,DFAC,VFAC,PWV,NLV,NLOSV,NLOSVP1,VMAP,VLIM,
     &			DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij)
		call FixModeltdn(1,NLV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO)
		if(bWriteMS.and.NitV.eq.NiterT) then
		    call writeM_Scomp(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV,PWV,VRATIO,
     &			NLmaxV,NSRV,SRCV,DOYSV,XCEV,VOBS,VMOD,XLON,XCtpr,NLOSVP1,FIX,VSSIG,NBSV)
		end if
		VRATION(NitV+1) = VRATIO		! Put in 7/9/00
		if(NitV.ge.2.and.VRATION(NitV+1).gt.(VRATION(NitV)+(VRATION(NitV+1)*0.05)).and.
     &		VRATION(NitV+1).gt.(VRATION(NitV-1)+(VRATION(NitV+1)*0.10))) then
		    print *, 'THE SUMMED VELOCITY RATIO IS NO LONGER CONVERGING'
C		    NitV = NiterT
		end if
		VFACTOR = 1.0			! Use only in the following subroutine for UCSD velocities
		if(bVUCSD) VFACTOR = 0.10	! Use only in the following subroutine for UCSD velocities
		call MkVMaptdn(LON,LAT,ITIM,RPX,XEV,VWTij,FIX,NLV,NLOSV,NLOSVP1, DMAPV,XLON,XLAT,
     &			nLng,nLat,nTV,VMAP,ANLIMITV, VLIM,VSSIG,VMOD,VSIGB,NVS,NBSV,VMAPSIG,
     &			ANMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,
     &			Alng,ConstL,CONRV,CONVT,aNdayV*VFACTOR,VZtmp)
		do N=1,NTV
C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'In V determination VMAP min max',amin,amax
		   if (NitV .eq. NiterT) call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N))
		   call ArrR4SetMinMax(-nLngLat,VMAP(1,1,N),VLL,VUL)
		end do
	    end if
	    if(bstopw) then
		call stopwatch(' ',' ')
		print *, 'After V deconvolution and writing VMAPHOLE'
		print *, 'Arranging DMAP'
		call stopwatch(' ',' ')
	    end if
	    if (.not.bDproxy.or.Nit.ge.1) then
	     if (bGcon .or. bVmv2) then	! Update density and density shifts even on the last
		call FillWholeT(0,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD)
		call FillMaptN(0,XCbeGG,nLng,nLat,nTG,ConsTG,DMAP,DDD)
		do N=1,NTG
C		  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C		  print *, 'Before Grids MkShiftd DMAP  min max N',amin,amax,N
		  call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))
		  call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))
		  call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD,3,0.0,90.0)
C		  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C		  print *, 'Before MkShiftd DMAP  min max N',amin,amax,N
		end do
		call Timesmooth(nLng,nLat,nTG,DMAP, 0,1.0*CONDT/aNdayG,0.,Dtmp)
		do N=1,NTG
		  call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))
		end do
	     end if
	     if(bstopw) then
		 call stopwatch(' ',' ')
	         print *, 'After arranging DMAP'
	         print *, 'Before arranging VMAPD and D Mkshiftd'
	         call stopwatch(' ',' ')
	     end if
	     if (bGcon) then
		call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,
     &				nLng,nLat,NTV,NTG,ConsTV*aNdayV,VMAP,VMAPD,VVT,VVV)
		call FillWholeT(1,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD)
		do N=1,NTG
C		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
C		  print *, 'Before Grids MkShiftd VMAPD min max N',amin,amax,N
		  call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N))
		  call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N))
		  call GridSphere2D(ALng,nLng,nLat,1,VMAPD(1,1,N),CONRD,3,0.0,90.0)
C		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
C		  print *, 'Before MkShiftd VMAPD min max N',amin,amax,N
		end do
		call Timesmooth(nLng,nLat,nTG,VMAPD,0,1.0*CONDT/aNdayG,0.,Dtmp)
		do N=1,NTG
		  call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N))
 		 end do
	      end if
	    end if
	    if(bstopw) then
		call stopwatch(' ',' ')
		print *, 'After arranging VMAPD and before D Mkshiftd'
		call stopwatch(' ',' ')
	    end if
	    if(bGcon) then

C		call MkShiftd(XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,
C     &		VMAPD,DMAP,XCshift,RconstG,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
C     &		NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

	      if(Nit .ge. 20) then
		call MkShiftdnn(XCbeGG,XCtbegG,XCtendG,ALng,aNdayG,nLng,nLat,nMap,nTG,
     &		nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,RconstG,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL,
     &		NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)
	      else
		call MkShiftdn(XCbeGG,XCtbegG,XCtendG,ALng,aNdayG,nLng,nLat,nMap,nTG,
     &		VMAPD,DMAP,XCshift,RconstG,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
     &		NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)
	      end if

	    end if
	    if(bstopw) then
		call stopwatch(' ',' ')
		print *, 'After D Mkshiftd'
	    end if
	    if (bGcon) then				! G deconvolution
		NitG = NitG+1				! # G iterations
							! bVmv2=.TRUE., implies bVcon=.FALSE.
		if (bVmv2) then 			 ! Make a velocity map assuming mv^2 = C
		  RRV = RRD
 		   call Deneq3(nLng,nLat,NTG,DMAP,RRD,PWR,FALLOFF,DENEQ)
		   print *, 'Solar equatorial r^2 mv^2 = constant value DENEQ =',DENEQ
		  do N=1,NTG
		    call Mk_D2VN(nLngLat,DMAP(1,1,N),VMAPD(1,1,N),DENEQ,Speed,RRD)
		  end do
		  call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &				    nLng,nLat,NTG,NTV,ConsTG*aNdayG,VMAPD,DMAPV,DDT,DDD)
		  call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAP(1,1,N))
		  if(NitG .eq. NiterT) then
		    do N=1,NTV
	    	      call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAPHOLE(1,1,N))
		    end do
		  end if
		end if
		if(bstopw) then
		    print *, 'Before D MkPostd D deconvolution'
		    call stopwatch(' ',' ')
		end if
		call MkPostd(XCbeGG,XCtbegG,XCtendG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG,
     &		XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XCEARG,XLLG,XDLG,IYRSG,DOYSG,
     &		XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM)
		call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG,
     &			DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC)
C		do N=1,NTG
C	  	    call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C	  	    print *, 'Before MkGModeltdn DMAP min max N',amin,amax,N
C		end do
		if(.not.bTcon) then					! G-level model-maker from density
 		    call MkGModeltdn(XCbegG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX,
     &			WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP,
     &			nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi)
		    NGval = 2						! G-level fix
		else							! Brightness model-maker from density
		    call MkDHModeltd(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX,
     &			WTSG,DFAC,FALLOFF,NLG,NLOSG,DLOS,DISTG,NLOSGP1,DMAP,DBASE,BBASE,
     &			nLng,nLat,nTG,RRS,GMOD2,GWTij)
		    NGVal = 3		! Thomson scattering fix
		end if
		do I=1,NLG
		    GGOBS(I) = GOBS(I)  + BBASE(I)			! Each iteration, use modified observed brightnesses
		end do
		call FixModeltdn(NGVal,NLG,GGOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO)
		do I=1,NLG
		    GMOD2(I) = GMOD2(I) - BBASE(I)			! Model brightness minus base compares with GOBS in write
		end do
		if(bWriteMS.and.NitG.eq.NiterT) then 
		    call writeM_Scomp(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO,
     &			NLmaxG,NSRG,SRCG,DOYSG,XCEG,GOBS,GMOD2,XLON,XCtpr,NLOSGP1,FIX,GSSIG,NBSG)
		end if
		GRATION(NitG+1) = GRATIO		! Put in 7/9/00
		if(NitG.ge.2.and.GRATION(NitG+1).gt.(GRATION(NitG)+(GRATION(NitG+1)*0.05)).and.
     &		GRATION(NitG+1).gt.(GRATION(NitG-1)+(GRATION(NitG+1)*0.10))) then
		    print *, 'THE SUMMED DENSITY RATIO IS NO LONGER CONVERGING'
C		    Nit = NiterT
		end if
		if(.not.bTcon) then			! density map maker from g-level
 		    call MkDMaptdn(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1,
     &			nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG,
     &			ANMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,
     &			Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp)
		else					! density map maker from brightness
 		    call MkDMaptdn(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1,
     &			nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG,
     &			ANMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,
     &			Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp)
C  Calculate the new base to the density.  The base is made from the initial DEN1AU (DEN1AU/3.) and is 
C  assumed constant over the whole map.  Someday we may want to iterate this base to fit 1 AU density
C  (see Mk_Base subroutine).
		    call Mk_Base(VMAPD,DMAP,nLng,nLat,nTG,RRV,DEN1AU,FALLOFF,VEL1AU,
     &			  BadR4(),BadR4(),DBASE)
		end if
		do N=1,NTG
C	  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C	  print *, 'After MkDmaptd DMAP min max',amin,amax
		  if (NitG .eq. NiterT) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N))
		  call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG))
		end do
	    end if
	    if(bstopw) then
		call stopwatch(' ',' ')
		print *, 'After D deconvolution and writing DMAPHOLE'
		print *, 'Arranging VMAP'
		call stopwatch(' ',' ')
	    end if
	    if(Nit.ne.NiterT) then 	! Update velocity shifts and vmap except on last
 
 	       if (bVcon.or.bGmv2) then
		call FillWholeT(1,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) 
		call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV)
		do N=1,NTV
C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'After MkDmaptdn  VMAP min max N',amin,amax,N cd IPS
   		  call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
		  call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
		  call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,90.0)
		end do
		call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp)
		do N=1,NTV
 		   call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))
		end do
	      end if
	      if (.not.bOnly1) then
	       if (bVcon) then
		if(bstopw) then
		    call stopwatch(' ',' ')
		    print *, 'After arranging VMAP'
		    print *, 'Before arranging DMAPD and V Mkshiftd'
		    call stopwatch(' ',' ')
		end if
                call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &				nLng,nLat,NTG,NTV,ConsTG*aNdayG,DMAP,DMAPV,DDT,DDD)
		call FillWholeT(0,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV)
		do N=1,NTV
C		  call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax)
C	  	  print *, 'After MkDmaptdn DMAPV min max N',amin,amax,N
		  call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
		  call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
		  call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,90.0)
		end do
		call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp)
		do N=1,NTV
		  call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N))
		end do
		if(bstopw) then
		    call stopwatch(' ',' ')
		    print *, 'After arranging DMAPV and before V Mkshiftd'
		    call stopwatch(' ',' ')
		end if

C		call MkShiftd(XCbeV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,
C     &		VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
C     &		NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

		call MkShiftdn(XCbeV,XCtbegV,XCtendV,ALng,aNdayV,nLng,nLat,nMap,nTV,
     &		VMAP,DMAPV,XCshift,RconstV,DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,
     &		NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

	       end if
	      end if
	    end if
	    if(bstopw) then
		call stopwatch(' ',' ')
		print *, 'After V Mkshiftd'
	    end if
	end do						! End of iteration loop
	J = 0						! List bad IPS sources
	do I=1,NLV
	    if (NBSV(I) .eq. 0) then
		J = J+1
		write (*,'(3A,I4,A,F7.1,A,F9.3)') ' Bad V',SRCV(I),' L.O.S. removed, # =',
     &			I,', V =',VOBS(I),', elong. =',XEV(I)
	    end if
	end do
	write (*,'(A,I4,A)') ' There were',J,' bad V L.O.S. removed.' 
	J = 0	
	do I=1,NLG
	    if (NBSG(I).eq.0) then
		J = J+1
		if(.not.bTcon) write (*,'(3A,I6,A,F7.3,A,F9.3)') ' Bad',SRCG(I),' G L.O.S. removed, # =',
     &			I,', G = ',GOBS(I),', elong. =',XEG(I)
		if(bHelio) write (*,'(3A,I6,A,F7.3,A,F9.3)') ' Bad Helios',SRCG(I),' L.O.S. removed, # =',
     &			I,', G = ',GOBS(I),', elong. =',XEG(I)
	    end if
	end do
	write (*,'(A,I4,A)') ' There were',J,' bad density L.O.S. removed.' 
C Here VMAPHOLE contains the original point-P map (bVcon=.FALSE.) (   or is equal to VMAP 
	do N=1,NTV
	  call ArrR4Copy(nLngLat,VMAPHOLE(1,1,N),VMAP(1,1,N))
	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
	  print *, 'VMAP N, amin and amax = ', N, amin, amax
	  if(amax.ne.BadR4()) then
	    call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))   ! Remove 90 deg discrepancy
	    call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax)		  ! Remove small velocity holes
	    call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove 90 deg discrepancy
	    call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))   ! Remove seam
	    XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region of interest on map
	    call FillMapOLS(XCbeV(1,N),XCbeV(2,N),XCmidV,nLng,nLat,ConstL,VMAP(1,1,N))   ! Remove 180 deg bad spots
	    call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove seam
	    call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove 90 deg discrepancy
	  end if
	end do
	do N=1,NTG
	  call ArrR4Copy(nLngLat,DMAPHOLE(1,1,N),DMAP(1,1,N))
	  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
	  if(amax.ne.BadR4()) then
	    call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy
	    call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Remove small density holes
	    call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy 
	    call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
C
C  Below can't be done with Helios data since Helios is not at the location of Earth
C	    XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region of interest on Map
C	    call FillMapOLS(XCbegG(1,N),XCbegG(2,N),XCmidG,nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 180 deg bad spots
C	    call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
	    call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy
	  end if
	end do
	call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat,
     &			NTV,NTG,ConsTV*aNdayV,VMAP,VMAPD,VVT,VVV)

	do N=1,NTG
	  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
	  print *, 'VMAPD N, amin and amax = ', N, amin, amax
	end do

C	call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,nLng,nLat,
C    &			NTG,NTV,ConsTG*aNdayG,DMAP,DMAPV,DDT,DDD)
C*****  Take data out at different spacecraft (and earth) locations (bExtract tells where)
	print *, ' Before Extractd',XCbegG(1,1),XCtbegG,XCstrt
	call Extractd(bExtract,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat,
     &	      nMap,nTG,FALLOFF,VMAPD,DMAP,XCshift,RconstG,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)
	if(bWrite3D) then
C*****  Write to disk of results interpolated to the even latitudes, longitudes
C	and times over the heights (specified as nMapn) of the XCshift model.
	 RRSCON = (RRS/R1AU)**FALLOFF
	 do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	 end do
	 nTF = nTG
	 nMapn = nMap
	 Nit = Nit + 1
	 NinterC = 0
	 cPre = ' '
C The first write3D_infotd makes maps that have holes to be used to show Carrington maps 
	 print *, 'Before first write3Dinfo'
	 call write3D_infotd(0,cPre,Nit,NiterT,NinterC,NCoff,XCintG,XCintG,XCbegG,XCtbegG,XCtendG,RRS,dRR,	
     &			nLng,nLat,nMapn,nTF,nTmaxG,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D)
	do N=1,NTG
	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
	  if(amax.ne.BadR4()) then
	    call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Remove more small velocity holes
	    call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove 90 deg discrepancy
	    call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove seam
	    call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Remove more small velocity holes
	    call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove 90 deg discrepancy
	    call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove seam
	  end if
	end do
	do N=1,NTG
	  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
	  if(amax.ne.BadR4()) then
	    call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Remove more small density holes
	    call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy
	    call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
	    call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Remove more small density holes
	    call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy
	    call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
	  end if
	end do
	if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
	    call ArrR4Copy(nLngLatnTG,VMAP,VMAPD)
	else
	    call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat,
     &			NTV,NTG,ConsTV*aNdayV,VMAP,VMAPD,VVT,VVV)
	end if
	do N=1,NTG
	  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
	  print *, 'VMAPD N, amin and amax = ', N, amin, amax
	end do
C The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views
	 print *, 'Before second write3D_infotd'
	 cPre = 'nv3f'
	 NinterD = 3
	 call write3D_infotd(1,cPre,Nit,NiterT,NinterD,NCoff,XCintG,XCintG,XCbegG,XCtbegG,XCtendG,RRS,dRR,	
     &			nLng,nLat,nMapn,nTF,nTmaxG,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D)

	end if

c*TJD*added 1 to Nit to prevent _18 being added to file names
C and made NinterD = 3 to add 3 interpolations between times
	Nit = Nit + 1
	NinterD = 3
	
c*TJD* Magnetic modifications +
C	if ( bMagnetic ) then

C	if (.not.( (nTmaxV .eq. nTmaxG))) print *, 'check nTmaxV array boundaries'	

C	call MkVeltd(nLng,nLat,nMap,nTmaxV,XCshift,VMAPD, DVfact,TrueVel)

C	print *, 'nTmaxG',nTmaxG	
C	print *, 'nTV is ', nTV

C	call Write3D_bbt(1,Nit,NiterT,NinterD,NCoff,XCintG,XCbeGG,XCtbegG,
C     &XCtendG,RRS,dRR,nLng,nLat,nMap,nTmaxV,nTmaxV,XCbegMAT,XCendMAT,
C     &XCbegROI,XCendROI,XCshift,VMAPD,TrueVel,BR3D,BT3D,
C     &PWV,PWG,
C     &PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
C     &90.0,Scale,Btmp,BR2DT)

c
C	end if			     
C*TJD* Magnetic modifications -	

	call Say('IPSDtest','S','Stop','program successfully completed')
	end