C+
C NAME:
C	ipstd_20n_inp_mag3_v23
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. This analysis now includes
C     a way to input DSCOVR in-situ data.
C     This now incorporates a analysis to provide pseudo forecasts from the last input ISEE data.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	run ipstd_20n_inp_mag3_v23
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	 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 extractdv 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.  Prior I tried
C	using. Found an error 10/15/00 (A gridsphere2d was placed outside a do loop after 
C	the first MKshiftdn0n. This version should now have no seam in each temporal value.
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 IPSD2020), the
C	only way to force IPSD2020 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	Jan 2020 Bernie Jackson (UCSD) v20
C     Feb 2022 Bernie Jackson (UCSD) v22
C     Jan 2023 Bernie Jackson (UCSD) v23
C-
	program ipstd_20n_inp_mag3_v23

	parameter         (NLmaxG	=  30000,	! Max # G data points  (13000) (400,000)
     &			NLmaxV	=  30000,	! Max # V data points  (30000 should be more than adequate for any IPS difficulties)
     &			NGREST	=  1,		! G temporal resolution factor better than 1 day
     &			NVREST	=  1,		! V temporal resolution factor better than 1 day
     &			NGRESS	=  1,		! Spatial G resolution factor better than 20 degrees
     &			NVRESS	=  1,		! Spatial V resolution factor better than 20 degrees
     &			NF	      =  2,		! LOS resolution factor (Should be 1 if Spatial resolution is >= 2)
     &			nTmaxG	=  55*NGREST,	! Max # G times	(50) (100) (150)
     &			nTmaxV	=  55*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,
     &			 nLngLatnTG =nLng*nLat*nTMaxG,
     &			 nLngLatnTV =nLng*nLat*nTMaxV)

	parameter	      (NintLn	= 3,		! High resolution output parameters
     &			 NintLa = 3,
     &			 NintH  = 3)

	parameter	      (dRR	= 0.1)		! Resolution (AU) in radial direction 		(0.1)
	parameter	      (NTP	= 19)		      ! number of body locations possible to extract.

	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),
     &			LON_in	(NLmaxG),
     &			LAT_in	(NLmaxG),
     &			ITIM_in	(NLmaxG)
		
	real		      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
     &			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)
     &			GGGOBS(NLmaxG),	! Dummy Solar wind observed G values
     &			GOBS	(NLmaxG),	! Solar wind observed G values
     &			DGOBS (NLmaxG),	! Solar wind observed G value errors
     &			GGOBS	(NLmaxG),	! Dummy Solar wind observed G values
     &			DGGOBS(NLmaxG),	! Dummy Solar wind observed G value errors
     &			GOBSS	(NLmaxG),	! Solar wind observed G values     
     &			GSSIZE(NLmaxG),	! Dummy Source g-level size
     &			GSIZE	(NLmaxG),	! Source g-level size
     &			GFFREQ(NLmaxG),	! Dummy source g-level frequency
     &			GFREQ	(NLmaxG),	! Source g-level frequency
     &			VVOBS	(NLmaxV),	! Dummy Solar wind observed velocities
     &			VOBS	(NLmaxV),	! True Solar wind observed velocities
     &			DVOBS	(NLmaxV),	! True Solar wind observed velocity errors
     &			VSIZE	(NLmaxV),	! Source velocity size
     &			VFREQ	(NLmaxV),	! Source velocity frequency
     &			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),
     &			XLproj(NLOSGP1,NLMAXG),
     &			VFAC	(NLOSGP1,NLmaxG),	! New array (3/99)
     &			DFAC	(NLOSGP1,NLmaxG),

     &			XLON_in     (NLmaxG),
     &			XLONP_in    (NLmaxG),
     &			XLAT_in     (NLmaxG),
     &			RPX_in      (NLmaxG),
     &			XLproj_in   (NLmaxG),
     &			VFAC_in	(1,NLmaxG),
     &			DFAC_in	(1,NLmaxG)

	character	      SRCV	(NLmaxV)*9,		! Velocity source name
     &			SRCVLAST*9,
     &			SRCVG	(NLmaxV)*115,		! Full source info and g-value and velocity file
     &                  cSRCVG*115,
     &			SRCGV	(NLmaxG)*115,		! Full source info and g-value and velocity file
     &                  cSRCGV*115,
     &			SRCVGsav(NLmaxV)*115,
     &			SRCGVsav(NLmaxG)*115,
     &			SRCVsav (NLmaxV)*9,
     &			SRCGsav (NLmaxG)*9,
     &			SRCGG	(NLmaxG)*9,		! G-level source name
     &			SRCGLAST*9,
     &			SRCALAST*9,
     &			cvgfiles(NLmaxG)*11,
     &			cggfiles(NLmaxG)*11,
     &			SRCGS*9,
     &                  cSite*2

	integer		IYRFV	(NLmaxV),
     &			IRECV	(NLmaxV),
     &			NSRV	(NLmaxV),
     &			NSRG	(NLmaxG),
     &			IYRSV	(NLmaxV),		! Year of source observation
     &			NBSV	(NLmaxV)

	real		      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
     &			VMOD	(NLmaxV),	! Model velocities
     &			VSSIG	(NLmaxV),	! Velocity source sigma	
     &			VWTij	(NLOSV,NLmaxV),
     &			WTSV	(NLOSV,NLmaxV),
     &			VWTij_in(NLmaxV),
     &			GWTij_in(NLmaxG),     
     &			WTSV_in	(NLmaxV),
     &			WTSD_in	(NLmaxG),     
     &			WTPV_in	(NLmaxV),
     &			WTPD_in	(NLmaxG)

	real		      XCinTV(nTmaxV),
     &			XCinTG(nTmaxG),
     &			DOYV	(nTmaxV),
     &			DOYG	(nTmaxG),
     &			XCbeV	(2,nTmaxV),
     &			XCbeGG(2,nTmaxG),
     &			XCbe	(2,4*nTmaxG)

	integer		IYRV	(nTmaxV),
     &			IYRG	(nTmaxG),
     &			NMAPD	(nLng,nLat,nTmaxG),	! Scratch integer arrays
     &			NMAPV	(nLng,nLat,nTmaxV)

	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),
     &			DMHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1),
     &			VMHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1),
     &			DMAPHOLE(nLng,nLat,nTmaxG),
     &			DMAP	(nLng,nLat,nTmaxG),	! Density map (no holes)
     &			VMAP	(nLng,nLat,nTmaxV),
     &			VMAPD	(nLng,nLat,nTmaxG),
     &			VMAPF	(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),
     &			V3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &			D3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),     
     &			BR3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &			BT3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &			BN3DHR	(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1),
     &			V3DT	(nLng,nLat,nMap,4*nTmaxG),
     &			D3DT	(nLng,nLat,nMap,4*nTmaxG),
     &			Vtmp 	(nLng,nLat,nTmaxG,3),		! Scratch real arrays
     &			VtmpV	(nLng,nLat,nTmaxV,3),
     &			VZtmp	(nLng,nLat,nTmaxV),
     &			Dtmp 	(nLng,nLat,nTmaxG),
     &			XLTtmp	(nLng,nTmaxG,5),
     &			XLtmp	(nLng,nTmaxG),
     &			XLtmpt	(nLng,nTmaxG),
     &			Dtmpt	(nLng,nTmaxG),
     &			Vtmpt	(nLng,nTmaxG),
     &			Vtp	(nLng,nLat),
     &			Dtp	(nLng,nLat),
     &			ANMAPD	(nLng,nLat,nTmaxG),
     &			AWTMAPD	(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),
     &			AWTMAPV	(nLng,nLat,nTmaxV),
     &			ANDMAPV	(nLng,nLat,nTmaxG),
     &			AWTDMAPV	(nLng,nLat,nTmaxG),
     &			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),
C     &			Scale	(2)
     &			Scale	(4)

	byte		      NSIDE(9*nLng*nLat)

	parameter	      (nCar	= 2000)
	real*8		JDCar(nCar),
     &			JD, JEpoch,
     &			dLngSun,dLatSun,dDisSun
     
	real*8 		DOYSG8	(NLmaxG),	! DOY (and fraction) of source obs.
     &			DOYSGG8	(NLmaxV),	! DOY (and fraction) of source obs.
     &			DOYSV8	(NLmaxV),	! DOY (and fraction) of source obs.
     &			RADUMG8	(NLMAXG),	! Dummy array for glevel sort
     &			RADUMV8	(NLMAXV),	! Dummy array for velocity sort
     &			DOYSsav8,
     &			dorbit
          
	real*8		XCtpr	(NLOSGP1,NLMAXG),
     &			XCtim	(NLOSGP1,NLMAXG),
     &			XCtpr_in(NLMAXG),
     &			XCtim_in(NLMAXG),
     &			XCintDV	(nTmaxV),
     &			XCintDG	(nTmaxG)
	
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,
     &			bVNago, bVOoty, bVUCSD, bVEISC, bVGen,
     &			bGUCSD, bGCamb, bGOoty, bGNago, bGGen,
     &			bVACE,  bDACE,  bFACE,  bGACE,
     &			bVFACE, bVACE2, bDACE2,
     &			bVWind, bDWind,
     &			bVCELIAS,bDCELIAS,
     &			bVDSCOVR,bDDSCOVR,
     &                  bsinglesitev, bmultisitev, bsinglesiteg, bmultisiteg,
     &                  bsinglesitevo, bsinglesitevn, bsinglesitevk, bsinglesitevm, bsinglesitevb, bsinglesitevl,
     &                  bsinglesitego, bsinglesitegn, bsinglesitegk, bsinglesitegm, bsinglesitegb, bsinglesitegl,
     &                  bmultisitevo,  bmultisitevn,  bmultisitevk,  bmultisitevm,  bmultisitevb,  bmultisitevl, 
     &                  bmultisitego,  bmultisitegn,  bmultisitegk,  bmultisitegm,  bmultisitegb,  bmultisitegl,
     &                  bVERR0OK, bGERR0OK,
     &                  bVERR0OKO, bVERR0OKN, bVERR0OKK, bVERR0OKM, bVERR0OKB, bVERR0OKL,
     &                  bGERR0OKO, bGERR0OKN, bGERR0OKK, bGERR0OKM, bGERR0OKB, bGERR0OKL,
     &			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	      /.TRUE./,	! Magnetic map available
     &			bMagnet		/.FALSE./,	! Magnetic map check
     &			bmagexx		/.TRUE./,	! Produce 3D magnetic field data files
     &			bmagex		/.TRUE./,	! Extract magnetic field data as well as velocity and density data
     &			bAdjust		/.TRUE./,	! Adjust the magnetic field BF factors for solar cycle and instrument types
     &			bGWRITE		/.FALSE./,	! Write out a G data file
     &			bVWRITE		/.FALSE./,	! Write out a V data file
     &			bMirror		/.FALSE./,
     &			bjMirror	      /.FALSE./,
     &			bWrite		/.FALSE./,
     &			bWrite3D1	      /.FALSE./,
     &			bWrite3D2	      /.FALSE./,
     &			bWrerr		/.FALSE./,
     &			bDdener		/.TRUE./,
     &			bVdener		/.FALSE./,
C     &			bVdener		/.TRUE./,
C     &			bVverer		/.FALSE./,
     &			bVverer		/.TRUE./,
     &			bDverer		/.FALSE./,    
     &			bDdenerb	      /.TRUE./,
     &			bVdenerb	      /.FALSE./,
C     &			bVdenerb	      /.TRUE./,
C     &			bVvererb	      /.FALSE./,
     &			bVvererb	      /.TRUE./,
     &			bDvererb	      /.FALSE./,
     &			bWrite3D_HR	      /.TRUE./,
     &			bWrite3D_HRb	/.FALSE./,
     &			bWrite3Do	      /.TRUE./,
     &			bWriteMS	      /.FALSE./,
     &			bWritelos	      /.FALSE./,
     &			bOneOr1		/.TRUE./,
     &			bOneOr2		/.TRUE./,
     &			bFillmag	      /.TRUE./,
     &			bForecast,
     &			bSource	      /.TRUE./,	! Write out good sources
     &			bExtract(NTP)	/2*.FALSE.,.TRUE.,9*.FALSE.,7*.FALSE./,
     &			bEWRITE,			! Write out an Earth file
     &			bMWRITE,			! Write out a Mercury file
     &			bVEWRITE,			! Write out a Venus file
     &			bMAWRITE,			! Write out a Mars file     
     &			bUWRITE,			! Write out a Ulysses file
     &			bSAWRITE,			! Write out a STEREO A file
     &			bSBWRITE,			! Write out a STEREO B file
     &			bMSWRITE,			! Write out a MESSENGER file
     &                  bPSWRITE,               ! Write out a Parker Solar Probe file
     &                  bBCWRITE,               ! Write out a BepiColumbo file
     &                  bSOWRITE,               ! Write out a Solar Orbiter file
     &			b67WRITEF,			! Write out a comet 67P/CG file
     &			bExtractf(NTP)	/2*.FALSE.,.TRUE.,9*.FALSE.,7*.FALSE./,
     &			bEWRITEF,			! Write out an Earth filled file
     &			bMWRITEF,			! Write out a Mercury filled file
     &			bV2WRITE,			! Write out a Venus filled file
     &			bMAWRITEF,			! Write out a Mars filled file     
     &			bUWRITEF,			! Write out a Ulysses filled file
     &			bSAWRITEF,			! Write out a STEREO A filled file
     &			bSBWRITEF,			! Write out a STEREO B filled file
     &			bMXWRITEF,			! Write out a MESSENGER filled file
     &                  bPSWRITEF,              ! Write out a Parker Solar Probe filled file
     &                  bBCWRITEF,              ! Write out a BepiColumbo filled file
     &                  bSOWRITEF,              ! Write out a Solar Orbiter filled file
     &			bC67WRITEF			! Write out a comet 67P/CG filled file
	
	equivalence	      (bExtract( 3), bEWRITE),
     &			(bExtract( 1), bMWRITE),
     &			(bExtract( 2), bVEWRITE),
     &			(bExtract( 4), bMAWRITE),
     &			(bExtract(10), bUWRITE),
     &			(bExtract(13), bSAWRITE),
     &			(bExtract(14), bSBWRITE),
     &			(bExtract(15), bMSWRITE),
     &			(bExtract(16), b67WRITE),
     &			(bExtract(17), bPSWRITE),
     &			(bExtract(18), bBCWRITE),
     &			(bExtract(19), bSOWRITE),        
     &			(bExtractf( 3), bEWRITEF),
     &			(bExtractf( 1), bMWRITEF),
     &			(bExtractf( 2), bV2WRITEF),
     &			(bExtractf( 4), bMAWRITEF),
     &			(bExtractf(10), bUWRITEF),     
     &			(bExtractf(13), bSAWRITEF),
     &			(bExtractf(14), bSBWRITEF),
     &			(bExtractf(15), bMXWRITEF),
     &			(bExtractf(16), bC67WRITEF),
     &			(bExtractf(17), bPSWRITEF),
     &			(bExtractf(18), bBCWRITEF),
     &			(bExtractf(19), bSOWRITEF)      
     
	character 	cPrefix (NTP)*2	/'Me','Ve','Ea','Ma','Ju','Sa','Ur','Ne','Pl','Ul','H1','H2','Sa','Sb','MS',
     &			'Cg','Ps','Bc','So'/
	character 	cPrefixf(NTP)*2	/'m1','v2','e3','m4','j5','s6','u7','n8','p9','uy','ha','hb','s1','s2','mx',
     &			'67','p1','b1','o1'/ ! p1 = Parker, solar, b1 = Bepi, o1 = Solar Orbiter, 
 
	integer		iEorWG	/0/,
     &			iEorWV	/0/,
c     &			MJD,
c     &			MJDfrst		/48180/, ! MJD for earliest data file
c    &			MJDlast		/49622/, ! MJD for latest data file
     &			MJDfrst	/ 46100/,
     &			MJDlast	/ 60000/,
     &			LimitoV	/-1/,
     &			LimitoG	/-1/,
     &			NLOSWV	/1/,
     &			NLOSWG	/1/,
     &			NTV		/44/,
     &			NTG		/44/

	parameter	      iVarp	=	27	! Number of potential data files accessed via a parameter list (6 IPS + 5 in-situ + 16 mag)
	parameter	      iSyst	= 	8	! total number of IPS systems accessed at one time
	integer		iSys        (iSyst)   ! System number 1-UCSD, 2-Eiscat, 3-Ooty, 4-Nagoya, 5-KSWC, 6-MEXART  7-BSA3, 8-LOFAR
      real              VERRLS      (iSyst)   ! General system velocity error ratio allowed from 0.0 - 1.0: 0.0 - nothing, 1.0 - all    
      real              GERRLS      (iSyst)   ! General system g-level error ratio allowed from 0.0 - 1.0: 0.0 - nothing, 1.0 - all

	integer		NiterT	/18/,
C     &			iOffUT	/ 0/,
     &			iOffUT	/ 7/,		! Changed 5/2/11 BVJ
     &			NAV		/ 0/,
     &			NC		/ 1/,		! # rotations averaged
     &			nDUMP		/ 0/,		! Threshold for #pnts/bin
     &			nFILL		/ 3/,		! Min. # neighbours for fill
     &			NmidHRV	/0/,
     &			NmidHRG	/0/,
     &			NLVmin	/50/,		! Default minimum number of velociy lines of sight
     &			NLGmin	/50/		! Default minimum number of g-level lines of sight

	real		      YLbeg		/-90./,
     &			YLend		/ 90./,
     &			XCadj(3)	/3*0./,
     &			XCadjB(3)	/3*0./,
     &			XCmargin 	/.5/,
     &			XElowG	/11.5/,
     &			XElowV	/11.5/,
     &			XElowNG	/11.5/,
     &			XElowNV	/11.5/,
     &			XElowOG	/11.5/,
     &			XElowOV	/11.5/,
     &			XElowKG	/11.5/,
     &			XElowKV	/11.5/,
     &			XElowMG	/21.0/,
     &			XElowMV	/21.0/,
     &			XElowPG	/26.1/,
     &			XElowPV	/26.1/,
     &			XElowLG	/20.0/,
     &			XElowLV	/20.0/,
     &			XEhighG	/80./,
     &			XEhighV	/180./,
     &			XRlimG	/90./,
     &			XRlimV	/90./,
     &                  GPOWER      /0./,
     &			R1AU		/1.0/,
     &			DEN1AU	/5.0/,
     &			ANLIMITV	/1.2/,
     &			ANLIMITG	/1.2/,
     &			ERDLOS	/1.2/,
     &			ERVLOS	/1.2/,    
     &			ERDLOSB	/1.2/,
     &			ERVLOSB	/1.2/,
     &			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		/14./,
     &			CONRD		/14./,
     &			CONRO		/8.0/,
     &			CONVT		/0.65/,
     &			CONDT		/0.65/,
     &			CONOT		/0.45/,
     &			ConsTV	/1.0/,
     &			ConsTG	/1.0/,
     &			VLIMM		/1200./,	! Velocity upper limit for sources
     &			VLIML		/201./,		! Velocity lower limit for sources
     &			VLIM		/1200./,	! Velocity upper limit for tomography
     &			VLIL		/100./,		! Velocity lower limit for tomography
     &			GLIMP		/4.0/,		! G-level upper limit for sources
     &			GLIMOF	/0.5/,		! Ooty g-level input excursion multiplier
     &			GLIMNF	/1.0/,		! Nagoya g-level input excursion multiplier
     &			GLIMKF	/1.0/,		! KSWC g-level input excursion multiplier
     &			GLIMMF	/0.2/,		! MEXART g-level excursion input multiplier
     &			GLIMLF	/1.0/,		! LOFAR g-level excursion input multiplier
     &			GLIMPF	/1.0/,		! BSA (Pushchino) g-level input excursion multiplier
     &			GLIMF		/1.0/,		! General g-level input excursion multiplier
     &			GLIMM		/0.2/,		! G-level lower limit for sources
     &			DMAX		/200./,		! Density upper limit for in-situ inputs
     &			DMIN		/0.1/,		! Density lower limit for in-situ inputs
     &			VVMAX		/1500./,	      ! Velocity upper limit for in-situ inputs	
     &			VVMIN		/100./		! Velocity lower limit for in-situ inputs
			

 	real		      RRS		/-15./,	      ! Reference distance for deconvolution
     &			RRX		/1.0/,	      ! Reference distance for Map extraction	
     &			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./,
     &			FALLOFFD	/2.0/,
     &			FALLOFFV	/0.0/,
     &			FALLOFFBR	/2.0/,
     &			FALLOFFBT	/1.0/,
     &			FALLOFFBN	/1.34/,
     &			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.20/,			 	! Radial power of density
     &			PWRV		/0.20/,				! Radial power of velocity
     &			AEmaxV	/800.0/,			      ! Anti-Earth maximum velocity
     &			AEmaxD	/10.0/,			      ! Anti-Earth maximum density at 1 AU
     &			AEminV	/200.0/,			      ! Anti-Earth minimum velocity
     &			AEminD	/2.0/,			      ! Anti-Earth minimum density at 1 AU
     &			AEangV	/135./,			      ! Angular distance from Earth to begin V limit
     &			AEangD	/135./			      ! Angular distance from Earth to begin D limit

	real*8		MJDref,
     &			MJDcntr,
     &			XCtbegV,
     &			XCtbegG,
     &			XCtendV,
     &			XCtendG,
     &                  DIDOYS,
     &                  AIDOYS

	character	      cCam*4,
     &			cMon*3,
     &			cStrID*11,
     &			cVar(27)*256,
     &			ccVar*256,
     &			cArg*100,
     &			cStr*80,
     &			cPre*4

	character	      cWildNagoya*100,
     &			cWildOoty*100,
     &			cWildGen*80,
     &			cWildGens(8)*100,			! No of IPS systems general files
     &                  cWildACE*80,
     &                  cWildACE2*80,
     &                  cWildWind*80,
     &                  cWildVCELIAS*80,
     &                  cWildDCELIAS*80,
     &                  cWildDSCOVR*80,
     &			cWild*80

	external	      iReadNagoyan8, iProcessNagoyan8
	external	      iReadOotyn8,   iProcessOotyn8
	external	      iProcessGen8ss_n, iProcessGen8ms_n
	external          iReadGen8ss_n, iReadGen8ms_n     
 	external 	      EARTH

	include		'dirspec.h'
	include		'sun.h'
	include		't3d_array_3.h'
	include		'b3d_param_3.h'
C	include		't3d_array.h'
	character	      cWildMAG(B3D__NN)*256
	logical		bWildMAG(B3D__NN)
      logical           bForinput    /.FALSE./
	integer		Flt2Str

c**TJD**added for magnetic time dependant version
	real 		      Btmp(nLng,nLat,3),		!scratch
c    &			XCbegMAT(nTmaxG),		!3-Car time frame
c    &			XCendMAT(nTmaxG),
     &			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),		!scratch
     & 			BR2DT(nLng,nLat,nTmaxG),	!scratch
     &			BT3D(nLng,nLat,nMap),		!scratch
     &			BN3D(nLng,nLat,nMap),		!scratch
     & 			BRR2DT(nLng,nLat,nTmaxG),	!scratch
     & 			BRC2DT(nLng,nLat,nTmaxG),	!scratch	
     & 			BTC2DT(nLng,nLat,nTmaxG),	!scratch
     & 			BNC2DT(nLng,nLat,nTmaxG),	!scratch
     &			TT3D	(4*nTmaxG),		      !holds XC corresponding to 6 hour int. steps output by write3dinfo_td3d
     &			XC3DT	(nLng,nLat,nMap,4*nTmaxG,3),							
     &			BRcheck (nLng,nLat),
     &			Vcheck5 (nLng,nLat),
     &			Vcheck7 (nLng,nLat)
     

C           real         aimagess(2048/2,1024/2,1)
            real         aimagess(1536,768)
            print*, 'Before gridsphere test'
            icirx = 1536
            iciry = 768
            CONRV = 2.0
            do I=1,icirx
              do J=1,iciry
                aimagess(I,J) = 1.0        
              end do
            end do

            print*, 'Got to here 1',icirx,iciry,CONRV
            call GridSphere2D(1.0,icirx,iciry,1,aimagess,CONRV,0,0.0,0.0)               ! Smooth the initial image by Gridsphere2D
            print*, 'Got to here 2',icirx,iciry,CONRV
      STOP

	call T3D_set_mode(T3D__MODE_TIME, 1)
C     ALSTGFREQ = 0.0
	iVar = iVarp				! Up to 6 IPS pointers, 5 in-situ pointers, 16 magnetic field input = 27 total
	call ForeignArg(' ',iVar,cVar,cArg)

	call ArrI4Constant(iSyst,0,iSys)	! Initialize iSys array to 0 (the total possible number of IPS systems, iSyst = 7)

	do I=1,iVar
	write(*,'(A,I3,1X,A64)') 'I, cVar(I)', I, cVar(I)
	  ccVar  = cVar(I)
	  if(ccVar(1:9) .eq.'UCSD=UCSD')	iSys(1) = 1
	  if(ccVar(1:8) .eq.'gen=UCSD') then
		call ForeignFile(1,ccVar,'gen',cwildGens(1))
		iSys(1) = 1
	  end if
	  if(ccVar(1:13).eq.'EISCAT=EISCAT')	iSys(2) = 2
	  if(ccVar(1:10).eq.'gen=EISCAT') then
		call ForeignFile(1,ccVar,'gen',cwildGens(2))
		iSys(2) = 2
	  end if
	  if(ccVar(1:9) .eq.'ooty=ooty')	iSys(3) = 3
	  if(ccVar(1:8) .eq.'gen=ooty') then
		call ForeignFile(1,ccVar,'gen',cwildGens(3))
		iSys(3) = 3
	  end if
	  if(ccVar(1:13).eq.'nagoya=nagoya')	iSys(4) = 4
	  if(ccVar(1:10).eq.'gen=nagoya') then
		call ForeignFile(1,ccVar,'gen',cwildGens(4))
		iSys(4) = 4
	  end if
	  if(ccVar(1:8) .eq.'gen=KSWC') then
		call ForeignFile(1,ccVar,'gen',cwildGens(5))
  		iSys(5) = 5
	  end if
	  if(ccVar(1:10).eq.'gen=MEXART') then
		call ForeignFile(1,ccVar,'gen',cwildGens(6))
	  	iSys(6) = 6
	  end if
	  if(ccVar(1:8).eq.'gen=BSA3') then
		call ForeignFile(1,ccVar,'gen',cwildGens(7))
	  	iSys(7) = 7
	  end if
	  if(ccVar(1:9).eq.'gen=LOFAR') then
		call ForeignFile(1,ccVar,'gen',cwildGens(8))
	  	iSys(8) = 8
	  end if
	end do

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

	DGOBSS  = 1.0

	Scale(1) = 0.0
	Scale(2) = 1.0
	Scale(3) = 0.0
	Scale(4) = 1.0

	I = iDeleteSymbol(cSym,2)
C	iYr = 1997
C	iYr = 1985
	iYr = 1985
	Doy = 1.
	call MAP_TZERO(EARTH,iYr,Doy,.01,nCar,JDCar)
	NCoff = N_CARRINGTON(iYr,Doy)

	call ArrR4Constant(NLmaxV,327.0,VFREQ)	! Initialize VFREQ array	! Nagoya or Ooty default
	call ArrR4Constant(NLmaxG,327.0,GFREQ)	! Initialize GFREQ array 	! Nagoya or Ooty defaut

	call ArrR4Constant(NLmaxV,0.1,VSIZE)	! Initialize VSIZE array 	! Nagoya or Ooty default
	call ArrR4Constant(NLmaxG,0.1,GSIZE)	! Initialize GSIZE array 	! Nagoya or Ooty default

	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 iOffUT is an offset time (in hours) between the local computer 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 7 hours (the difference between PDT and GST).

	call AskYN('Forecast mode?$no',bForecast)		! Regular or forecast mode
	if (bForecast) then
	  call AskI4('Offset between local computer time and UT (hours)?',iOffUT)
	  NLVmin = NLVmin/2
	  NLGmin = NLGmin/2
	end if

C	niterT = 2
	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, G-level, both,$1$3',I)
C-------
C A proxy map for g-values and/or velocity can be used to replace the point-P
C map determined from the observations.
C The map should be in files DMAP.DAT and VMAP.DAT. 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: UCSD, EISCAT, Ooty, Nagoya, Gen,?$1$5',I)
	    bVUCSD = I .eq. 1
	    bVEISC = I .eq. 2
	    bVOoty = I .eq. 3
	    bVNago = I .eq. 4
	    bVGen  = I .eq. 5

	    if(bVUCSD) XElowV = 35.0
	    if(bVEISC) XElowV = 3.0
	    if(bVOoty) XElowV = 11.5
	    if(bVNago) XElowV = 11.5

	    if(bVUCSD) then 
	      Print *, 'Sorry: The UCSD Velocity reconstructions are not yet supported. Try another.'
	      stop
	    end if

	    noSys  = 0
	    do ii=1, iSyst
	      if(iSys(ii).ne.0) noSys = noSys + 1             ! If noSys is zero, Gen input is not used
	    end do		

	    if(noSys.eq.0) then                               ! If noSys is zero, Gen input is not used
	      Print*, 'You asked for a General input velocity IPS file but did not specify its path. Please specify the path.'
	      stop
	    end if

          bVERR0OK = .TRUE.
	    if(noSys.eq.0) call ASKYN('Is use of IPS velocity data with NO error restriction OK?',bVERR0OK)  ! Gen is not specified
          if(bVERR0OK) VERRL = 1.0                          ! An error limit of 1.0 means all valid Velocities values are accepted'
          if(.not.bVERR0OK) then
            VERRL = 0.2                                     ! An error limit of 1.0 means all valid Velocities values are accepted'
	      call ASKR4('IPS velocity fractional error limit restriction',VERRL)
          end if

          if(bVGen) then
            do ii=1,iSyst
	        if(iSys(ii).ne.0) then

                if(iSys(ii).eq.3) then                      ! Ooty General Format IPS data
                  bVERR0OKO = .TRUE.
                  call ASKYN('Is use of Ooty IPS velocity data with NO error restriction OK?',bVERR0OKO)
                  if(bVERR0OKO) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKO) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('Ooty IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  bsinglesitevo = .TRUE.
                  bmultisitevo  = .FALSE.
                  print*, 'Ooty V',VERRLS(ii), bVERR0OKO, bsinglesitevo, bmultisitevo 
                end if

                if(iSys(ii).eq.4) then                      ! Nagoya General Format IPS data
                  bsinglesitevn = .FALSE.
                  bmultisitevn  = .FALSE.
	            call AskWhat('Nagoya single-site, multi-site velocity data, or both: SS, MS, BS$1$2',ISMBVN)
	            bsinglesitevn = ISMBVN .eq. 1
                  bmultisitevn  = ISMBVN .eq. 2
                  if(ISMBVN.eq.3) then
                    bsinglesitevn = .TRUE.
                    bmultisitevn  = .TRUE.
                  end if
                  bVERR0OKN = .TRUE.
                  call ASKYN('Is use of Nagoya IPS velocity data with NO error restriction OK?',bVERR0OKN)
                  if(bVERR0OKN) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKN) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('Nagoya IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  print*, 'ISEE V',VERRLS(ii), bVERR0OKN, bsinglesitevn, bmultisitevn 
                end if

                if(iSys(ii).eq.5) then                      ! MEXART General Format IPS data
                  bVERR0OKK = .TRUE.
                  call ASKYN('Is use of KSWC IPS velocity data with NO error restriction OK?',bVERR0OKK)
                  if(bVERR0OKK) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKK) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('KSWC IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  bsinglesitevk = .TRUE.
                  bmultisitevk  = .FALSE.
                  print*, 'KSWC V',VERRLS(ii), bVERR0OKK, bsinglesitevk, bmultisitevk 
                end if

               if(iSys(ii).eq.6) then                      ! MEXART General Format IPS data
                  bVERR0OKM = .TRUE.
                  call ASKYN('Is use of MEXART IPS velocity data with NO error restriction OK?',bVERR0OKM)
                  if(bVERR0OKM) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKM) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('MEXART IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  bsinglesitevm = .TRUE.
                  bmultisitevm  = .FALSE.
                  print*, 'MEXART V',VERRLS(ii), bVERR0OKM, bsinglesitevm, bmultisitevm 
                end if

                if(iSys(ii).eq.7) then                      ! BSA3 General Format IPS data
                  bVERR0OKOB = .TRUE.
                  call ASKYN('Is use of BSA3 IPS velocity data with NO error restriction OK?',bVERR0OKB)
                  if(bVERR0OKB) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKB) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('BSA3 IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  bsinglesitevb = .FALSE.
                  bmultisitevb  = .FALSE.                  
                  print*, 'BSA3 V',VERRLS(ii), bVERR0OKB, bsinglesitevb, bmultisitevb 
                end if

               if(iSys(ii).eq.8) then                     ! LOFAR General Format IPS data
                  bsinglesitevl = .FALSE.
                  bmultisitevl  = .FALSE.
	            call AskWhat('LOFAR single-site, multi-site velocity data, or both: SS, MS, BS$1$2',ISMBVL)
	            bsinglesitevl = ISMBVL .eq. 1
                  bmultisitevl  = ISMBVL .eq. 2
                  if(ISMBVL.eq.3) then
                    bsinglesitevl = .TRUE.
                    bmultisitevl  = .TRUE.
                  end if

                  bVERR0OKL = .TRUE.
                  call ASKYN('Is use of LOFAR IPS velocity data with NO error restriction OK?',bVERR0OKL)
                  if(bVERR0OKL) VERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid Velocities values are accepted'
	            if(.not.bVERR0OKL) then
                    VERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid Velocities values are accepted'
	              call ASKR4('LOFAR IPS velocity fractional error limit restriction',VERRLS(ii))
                  end if
                  print*, 'LOFAR V',VERRLS(ii), bVERR0OKL, bsinglesitevl, bmultisitevl 
                end if
              end if
            end do
          end if

	    bVmv2  = .FALSE.

	    if(bVEISC) NLOSWV = 4
	    if(bVNago) NLOSWV = 1
	    if(bVOoty) NLOSWV = 1
	    if(bVGen)  NLOSWV = 1

	    call AskI4('Minimum # of velocity lines of sight to use in the reconstruction?',NLVmin)

	    call AskYN('Use in-situ velocity in the 3D reconstructions?$yes',bVFACE)

	    bVACE    = .FALSE.
	    bVACE2   = .FALSE.
	    bVWind   = .FALSE.
	    bVCELIAS = .FALSE.
	    bVDSCOVR = .FALSE.
	    if(bVFACE) then
		call AskWhat('Use in-situ velocity from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$4',J)
	    	bVACE    = J .eq. 1
	    	bVACE2   = J .eq. 2
	    	bVWind   = J .eq. 3
	    	bVCELIAS = J .eq. 4
	    	bVDSCOVR = J .eq. 5
	    end if

	    if(bVFACE) call AskR4('Minimum velocity for in-situ velocity values?',VVMIN)
	    if(bVFACE) call AskR4('Maximum velocity for in-situ velocity values?',VVMAX)

	else				! No V deconvolution
	    bVEISC = .FALSE.	
	    bVNago = .FALSE.
	    bVOoty = .FALSE.
	    call AskYN('Do you want to make a velocity map using mv^2 = con?$no',bVmv2)
	end if
	if(.not.bVFACE) bVACE  = .FALSE.

	if (bGcon) then			! G deconvolution: select IPS station
	    if(bVOoty) call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$3',I)
	    if(bVNago) call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$4',I)
	    if(bVGen)  call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$5',I)
	    if(.not.bVUCSD.and..not.bVNago.and..not.bVOoty.and..not.bVGen) 
     &		       call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$5',I)
	    bGUCSD = I .eq. 1
	    bGCamb = I .eq. 2
	    bGOoty = I .eq. 3
	    bGNago = I .eq. 4
	    bGGen  = I .eq. 5

	    if(bGUCSD) then
	      Print *, 'Sorry: The UCSD g-level reconstructions are not supported. Try another.'
	      stop
	    end if

	    noSys  = 0
	    do ii=1, iSyst
	      if(iSys(ii).ne.0) noSys = noSys + 1
	    end do		
	    if(noSys.eq.0) then
	      Print*, 'You asked for a General input g-level IPS file but did not specify its path. Please specify the path.'
	      stop
	    end if

          bGERR0OK = .TRUE.
	    if(noSys.eq.0) call ASKYN('Is use of IPS g-level data with NO error restriction OK?',bGERR0OK)
          if(bGERR0OK)  GERRL = 1.0                         ! An error limit of 1.0 means all valid g-level values are accepted.'
          if(.not.bGERR0OK) then
            GERRL = 0.2                                     ! An error limit of 1.0 means all valid g-level values are accepted'
	      call ASKR4('Ooty IPS g-level fractional error limit restriction',GERRL)
          end if

          if(bGGen) then
            do ii=1, iSyst
	        if(iSys(ii).ne.0) then

                if(iSys(ii).eq.3) then                     ! Ooty General Format IPS data
                  call ASKYN('Is use of Ooty IPS g-level data with NO error restriction OK?',bGERR0OKO)
                  if(bVERR0OKO) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKO) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('Ooty IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  ISMBGO = 1
                  bsinglesitego = .TRUE.
                  bmultisitego  = .FALSE.
                  print*, 'Ooty G',GERRLS(ii), bGERR0OKO, bsinglesitego, bmultisitego 
                end if

                if(iSys(ii).eq.4) then                      ! Nagoya General Format IPS data
                  bsinglesitegn = .FALSE.
                  bmultisitegn  = .FALSE.
	            call AskWhat('Nagoya single-site, multi-site g-level data, or both: SS, MS, BS$1$1',ISMBGN)
	            bsinglesitegn = ISMBGN .eq. 1
                  bmultisitegn  = ISMBGN .eq. 2
                  if(ISMBGN.eq.3) then
                    bsinglesitegn = .TRUE.
                    bmultisitegn  = .TRUE.
                  end if
                  bGERR0OKN = .TRUE.
                  call ASKYN('Is use of Nagoya IPS g-level data with NO error restriction OK?',bGERR0OKN)
                  if(bGERR0OKN) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKN) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('Nagoya IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  print*, 'ISEE G',GERRLS(ii), bGERR0OKN, bsinglesitegn, bmultisitegn 
                end if

                if(iSys(ii).eq.5) then                      ! KSWC General Format IPS data
                  bGERR0OKK = .TRUE.
                  call ASKYN('Is use of KSWC IPS g-level data with NO error restriction OK?',bGERR0OKK)
                  if(bGERR0OKK) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKK) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('KSWC IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  bsinglesitegk = .TRUE.
                  bmultisitegk  = .FALSE.
                  print*, 'KSWC G',GERRLS(ii), bGERR0OKK, bsinglesitegk, bmultisitegk 
                end if

                if(iSys(ii).eq.6) then                      ! MEXART General Format IPS data
                  bGERR0OKM = .TRUE.
                  call ASKYN('Is use of MEXART IPS g-level data with NO error restriction OK?',bGERR0OKM)
                  if(bGERR0OKM) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKM) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('MEXART IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  bsinglesitegm = .TRUE.
                  bmultisitegm  = .FALSE.
                  print*, 'MEXART G',GERRLS(ii), bGERR0OKM, bsinglesitegm, bmultisitegm 
                end if

                if(iSys(ii).eq.7) then                      ! BSA3 General Format IPS data
                  bGERR0OKB = .TRUE.
                  call ASKYN('Is use of BSA3 IPS g-level data with NO error restriction OK?',bGERR0OKB)
                  if(bGERR0OKB) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKB) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('BSA3 IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  bsinglesitegb = .TRUE.
                  bmultisitegb  = .FALSE.
                  print*, 'BSA3 G',GERRLS(ii), bGERR0OKB, bsinglesitegb, bmultisitegb 
                end if

                if(iSys(ii).eq.8) then                     ! LOFAR General Format IPS data
                  bsinglesitegl = .FALSE.
                  bmultisitegl  = .FALSE.
	            call AskWhat('LOFAR single-site, multi-site g-level data, or both: SS, MS, BS$1$1',ISMBGL)
	            bsinglesitegl = ISMBGL .eq. 1
                  bmultisitegl  = ISMBGL .eq. 2
                  if(ISMBGL.eq.3) then
                    bsinglesitegl = .TRUE.
                    bmultisitegl  = .TRUE.
                  end if
                  bGERR0OKL = .TRUE.
                  call ASKYN('Is use of LOFAR IPS g-level data with NO error restriction OK?',bGERR0OKL)
                  if(bGERR0OKL) GERRLS(ii) = 1.0           ! An error limit of 1.0 means all valid g-level values are accepted'
	            if(.not.bGERR0OKL) then
                    GERRLS(ii) = 0.2                       ! An error limit of 1.0 means all valid g-level values are accepted'
	              call ASKR4('LOFAR IPS g-level fractional error limit restriction',GERRLS(ii))
                  end if
                  print*, 'LOFAR G',GERRLS(ii), bGERR0OKL, bsinglesitegl, bmultisitegl 
                end if
              end if
            end do
          end if

	    bGmv2  = .FALSE.
	    cCam = cEnvi//'CAM'
	    if(bGGen ) NLOSWG = 1
	    if(bGNago) NLOSWG = 1
	    if(bGOoty) NLOSWG = 1
	    if(bGcamb) NLOSWG = 2
	    if(bGUCSD) NLOSWG = 1

	    call AskI4('Minimum # of g-level lines of sight to use in the reconstruction?',NLGmin)
	    
	else				! No g deconvolution
	    bGUCSD = .FALSE.
	    bGCamb = .FALSE.
	    bGOoty = .FALSE.
	    bGNago = .FALSE.
	    bGgen  = .FALSE.
	    call AskYN('Do you want to make a density map using mv^2 = con?',bGmv2)
	end if

	    if(bGcon) call AskYN('Use in-situ density in the 3D reconstructions?$yes',bFACE)

	    bDACE    = .FALSE.
	    bDACE2   = .FALSE.
	    bDWind   = .FALSE.
	    bDCELIAS = .FALSE.

	    if(bGcon.and.bFACE) then
	      I = 0
	      if(bVACE)    call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$1',I)
	      if(bVACE2)   call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$2',I)
	      if(bVWind)   call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$3',I)
	      if(bVCELIAS) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$4',I)
	      if(bVDSCOVR) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$5',I)
C
	      if(I.eq.0)   call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$1',I)
	      bDACE    = I .eq. 1
	      bDACE2   = I .eq. 2
	      bDWind   = I .eq. 3
	      bDCELIAS = I .eq. 4 
	      bDDSCOVR = I .eq. 5 
	    end if
	    if(.not.bGcon.and.bFACE) then
	      if(bVACE)    call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$1',I)
	      if(bVACE2)   call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$2',I)
	      if(bVWind)   call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$3',I)
	      if(bVCELIAS) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$4',I)
	      if(bVDSCOVR) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS, DSCOVR$1$5',I)
	      bDACE    = I .eq. 1
	      bDACE2   = I .eq. 2
	      bDWind   = I .eq. 3
	      bDCELIAS = I .eq. 4
	      bDDSCOVR = I .eq. 5
	    end if
	    if(bDACE.and.bFACE)    call AskYN('Fit in-situ density data to ACE level 0 data?$yes',bFACE)
	    if(bDACE2.and.bFACE)   call AskYN('Fit in-situ density data to ACE level 2 data?$yes',bFACE)
	    if(bDWind.and.bFACE)   call AskYN('Fit in-situ density data to Wind?$yes',bFACE)
	    if(bDCELIAS.and.bFACE) call AskYN('Fit in-situ density data to CELIAS?$yes',bFACE)
	    if(bDDSCOVR.and.bFACE) call AskYN('Fit in-situ density data to DSCOVR?$yes',bFACE)

	    if(bFACE) call AskR4('Minimum density for fit in-situ density values?',DMIN)
	    if(bFACE) call AskR4('Maximum density for fit in-situ density values?',DMAX)
	    DLIM = DMAX/(RRS*RRS)

	    if(bDACE)  call AskYN('Fit mean g-level to ACE_L0?$yes',bGACE)
	    if(bDACE2) call AskYN('Fit mean g-level to ACE_L2?$yes',bGACE)
	    if(bDWind) call AskYN('Fit mean g-level to Wind?$yes',bGACE)
	    if(bDCELIAS) call AskYN('Fit mean g-level to CELIAS?$yes',bGACE)
	    if(bDDSCOVR) call AskYN('Fit mean g-level to DSCOVR?$yes',bGACE)

	if (bGCamb) call SetGipsy(cCam,0,iEdt,MJDfrst,MJDlast)
	if (bVNago .or. bGNago .or. bVEISC) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya)
	if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty'  ,cWildOoty  )
C	if (bVGen  .or. bGGen  ) call ForeignFile(iVar,cVar,'gen'   ,cWildGen   )
      if (bVACE .or .bDACE) then
	  call ForeignFile(-iVar-1,cVar,'ace',cWildACE)
           if (cWildACE .eq. ' ') then
C                cWildACE = '/zshare/dat/insitu/acesw_????.hravg'
                cWildACE = '/home/soft/dat/insitu/acesw_????.hravg'
           end if
           write(*,'(A75)') cWildACE
	end if
      if (bVACE2 .or .bDACE2) then
	  call ForeignFile(-iVar-1,cVar,'swace',cWildACE2)
           if (cWildACE2 .eq. ' ') then
C                cWildACE2 = '/zshare/dat/insitu/swace_????.hravg'
                cWildACE2 = '/home/soft/insitu/swace_????.hravg'
           end if
           write(*,'(A75)') cWildACE2
	end if
      if (bVWind .or .bDWind) then
	  call ForeignFile(-iVar-1,cVar,'wind',cWildWind)
           if (cWildWind .eq. ' ') then
C                cWildWind = '/zshare/dat/insitu/wind_swe/????_WIND_hourly_averages'
                cWildWind = '/home/soft/dat/insitu/wind_swe/????_WIND_hourly_averages'
           end if
           write(*,'(A75)') cWildWind
	end if
      if (bVCELIAS) then
	  call ForeignFile(-iVar-1,cVar,'celias',cWildVCELIAS)
           if (cWildVCELIAS .eq. ' '.and..not.bForecast) then
C                cWildVCELIAS = '/zshare/dat/insitu/celias_????.hravg'
                cWildVCELIAS = '/home/soft/dat/insitu/celias_????.hravg'
           end if
           if (cWildVCELIAS .eq. ' '.and.bForecast) then
C                cWildVCELIAS = '/zshare/dat/insitu/realcelias/celias_realtime.hravg*'
                cWildVCELIAS = '/home/soft/dat/insitu/realcelias/celias_realtime.hravg*'
           end if
           write(*,'(A75)') cWildVCELIAS
	end if
      if (bDCELIAS) then
	  call ForeignFile(-iVar-1,cVar,'celias',cWildDCELIAS)
           if (cWildDCELIAS .eq. ' '.and..not.bForecast) then
C                cWildDCELIAS = '/zshare/dat/insitu/celias_????.hravg'
                cWildDCELIAS = '/home/soft/dat/insitu/celias_????.hravg'
           end if
           if (cWildDCELIAS .eq. ' '.and.bForecast) then
C                cWildDCELIAS = '/zshare/dat/insitu/realcelias/celias_realtime.hravg*'
                cWildDCELIAS = '/home/soft/dat/insitu/realcelias/celias_realtime.hravg*'
           end if
           write(*,'(A75)') cWildDCELIAS
	end if
      if (bVDSCOVR .or .bDDSCOVR) then
	  call ForeignFile(-iVar-1,cVar,'DSCOVR',cWildDSCOVR)
        if (cWildDSCOVR .eq. ' ') then
          cWildDSCOVR = '/home/soft/dat/insitu/dscovr_????.hravg'
        end if
        write(*,'(A75)') cWildDSCOVR
	end if

	if(bVcon) then
	  bDaydV = .TRUE.
	  call AskYN('Deconvolve V from one day to the next?$yes',bDaydV)
	  if(bDaydV) then
	      call AskR4('Number of days to combine for velocity?',aNdayV)
	      if(bVEISC) NmidHRV = -1
	      if(bVNago) NmidHRV = -9
	      if(bVGen)  NmidHRV = -9
	      if(bVOoty) NmidHRV = -6
	      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 G from one day to the next$yes',bDaydG)
	  if(bDaydG) then
	      call AskR4('Number of days to combine for g-level?',aNdayG)
	      if(bGNago) NmidHRG = -9
	      if(bGGen)  NmidHRG = -9
	      if(bGOoty) NmidHRG = -6
	      call AskI4('Number of hours from g UT midnight?',NmidHRG)
	  end if
	  call AskI4('Number of g-level time values?',NTG) ! # NTG time itervals
	  if(NTG.gt.NTmaxG-2) then
	      NTG = NTmaxG-2
	      call AskI4('Number of 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 Calculate range of 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.

	call Julian(11,iYr,Doy,dble(MJDfrst),JEpoch)
	XClo = NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)
	XClo = XClo-.5
C	call Local2UT(iOffUT,iYr,Doy)
	call Local2UT(NmidHRV,iYr,Doy)					! Changed 5/2/11 BVJ
	XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo)

	RRV = RRS
	RRG = RRS

	if (bVcon) then

	    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

	    if(bVGen) then
	      if(iSys(3).eq.3) call AskR4('Minimum elongation for Ooty V-data?$0$11.5$',XElowOV)
	      if(iSys(4).eq.4) call AskR4('Minimum elongation for Nagoya V-data?$11.5$',XElowNV)
	      if(iSys(5).eq.5) call AskR4('Minimum elongation for KSWC V-data?$11.5$',XElowKV)
	      if(iSys(6).eq.6) call AskR4('Minimum elongation for MEXART V-data?$21.0$',XElowMV)	
	      if(iSys(7).eq.7) call AskR4('Minimum elongation for BSA (Pushchino) V-data?$21.0$',XElowPV)	
	      if(iSys(8).eq.8) call AskR4('Minimum elongation for LOFAR V-data?$20.0$',XElowLV)
	    else
	      call AskR4('Minimum elongation for V-data??$0$180$',XElowV)
	    end if

	    call AskR4('Maximum elongation for V-data?$'//cStr(:Flt2Str(XElowV,1,cStr))//'$180$',XEhighV)

	    call AskR4('Anti-Earth maximum velocity at 1 AU?$'//cStr(:Flt2Str(600.0,1,cStr))//'$1000.0$',AEmaxV)
	    call AskR4('Anti-Earth minimum velocity at 1 AU?$'//cStr(:Flt2Str(0.,1,cStr))//'$400.0$',AEminV)
	    call AskR4('Angular distance from Earth to begin V limit?$'//cStr(:Flt2Str(90.,1,cStr))//'$180.0$',AEangV)

	    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 upper limit for tomography velocity?',VLIM)   

	    call AskR4('Maximum lower limit for tomography velocity?',VLIL)
	    
	    call AskR4('Maximum upper limit for source velocity?',VLIMM)   

	    if(bVNago) VLIML = VLIL
	    if(bVGen)  VLIML = VLIL
	    call AskR4('Maximum lower limit for source velocity?',VLIML)   
 
	end if 
 
	if (bGcon) then 
 
	    call AskR4('# of G-points per bin that constitute a deconvolution?',ANLIMITG)

	    if(NiterT.gt.8) LimitoG = NiterT/2
	    call AskI4('Iteration at which to limit G-source deviations, -1=none',LimitoG) 
 
C  Get E or W or EW; Set the elongation angles

	    call AskWhat('All G-level data, East data, West data ?$0',iEorWG) 
	    if (iEorWG .eq. 2) iEorWG = -1  
	    if(bGNago) XElowG = XElowV
	    if(bGGen)  XElowG = XElowV
	    if(bGOoty) XElowG = XElowV
	    if(bGCAMB) XElowG = 30.0
	    if(bGUCSD) XElowG = XElowG

	    if(bGGen) then
	      if(iSys(3).eq.3) call AskR4('Minimum elongation for Ooty G-data?$0$11.5$',XElowOG)
	      if(iSys(4).eq.4) call AskR4('Minimum elongation for Nagoya G-data?$11.5$',XElowNG)
	      if(iSys(5).eq.5) call AskR4('Minimum elongation for KSWC G-data?$11.5$',XElowKG)
	      if(iSys(6).eq.6) call AskR4('Minimum elongation for MEXART G-data?$21.0$',XElowMG)	
	      if(iSys(7).eq.7) call AskR4('Minimum elongation for BSA (Pushchino) G-data?$21.0$',XElowPG)	
	      if(iSys(8).eq.8) call AskR4('Minimum elongation for LOFAR G-data?$20.0$',XElowLG)
	    else
	      call AskR4('Minimum elongation for G-data?$0$180$',XElowG)
	    end if

	    if(bGNago) XEhighG = 180. 
	    if(bGGen)  XEhighG = 180. 
	    if(bGOoty) XEhighG = 180.
	    if(bGCAMB) XEhighG = 80.0
	    if(bGUCSD) XEhighG = 80.0
	    call AskR4('Maximum elongation for G-data?$'//cStr(:Flt2Str(XElowG,1,cStr))//'$180$',XEhighG) 

	    call AskR4('Anti-Earth maximum density at 1 AU?$'//cStr(:Flt2Str(5.,1,cStr))//'$30.$',AEmaxD)
	    call AskR4('Anti-Earth minimum density at 1 AU?$'//cStr(:Flt2Str(0.,1,cStr))//'$4.0$',AEminD)
	    call AskR4('Angular distance from Earth to begin D limit?$'//cStr(:Flt2Str(90.,1,cStr))//'$180.0$',AEangD)
 
	    call AskR4('Maximum RA relative to Sun for G-data?$0$180$',XRlimG) 

	    RRG = RRS
	    call AskR4('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 G deconvolution?',GLIMP)
	    call AskR4('Minimum limit for G deconvolution?',GLIMM)

	    NDUP = 0
	    if(bGGen) then
	      if(iSys(3).eq.3) call AskR4('Factor by which to modify Ooty g-level excursions?',GLIMOF)
	      if(iSys(4).eq.4) call AskR4('Factor by which to modify Nagoya g-level excursions?',GLIMNF)
	      if(iSys(5).eq.5) call AskR4('Factor by which to modify KSWC g-level excursions?',GLIMKF)
	      if(iSys(6).eq.6) call AskR4('Factor by which to modify MEXART g-level excursions?',GLIMMF)	
	      if(iSys(6).eq.6) call AskI4('Duplicate MEXART g-level excursions?',NDUP)
	      if(iSys(7).eq.7) call AskR4('Factor by which to modify BSA (Pushchin g-level excursions?',GLIMPF)
	      if(iSys(8).eq.8) call AskR4('Factor by which to modify LOFAR g-level excursions?',GLIMLF)	
	    else
	      call AskR4('Factor by which to modify g-level excursions?',GLIMF)
	    end if

	end if

	call AskR4('G-level time hole filter?',ConsTG)
	call AskR4('Velocity time hole filter?',ConsTV)
  	CONRD = CONRD/float(NGRESS)
	call AskR4(' G-level spatial filter?',CONRD)
	if(bVEISC) CONRV = CONRD
	if(bVNago) CONRV = CONRD 
	if(bVGen)  CONRV = CONRD 
  	if(bVOoty) CONRV = CONRD
  	if(CONRD.ne.CONRV) CONRV = CONRD
	call AskR4('Velocity spatial filter?',CONRV)
  	CONDT = CONDT/float(NGREST)	
	call AskR4(' G-level temporal filter?',CONDT)
	if(bVEISC) CONVT = CONDT
	if(bVNago) CONVT = CONDT
	if(bVGen)  CONVT = CONDT
 	if(bVOoty) CONVT = CONDT
 	if(CONVT.ne.CONDT) CONVT = CONDT
	call AskR4('Velocity temporal filter?',CONVT)
	if(bGNago)	PWG = 0.35				! Returned back to fit data 5/30/2011 BVJ
C	if(bGNago)	PWG = 0.20				! Changed 4/28/2011 to fit data BVJ
	if(bGGen)	PWG = 0.35
	if(bGOoty) 	PWG = 0.35
	if(bGCamb) 	PWG = 0.15
	if(bGUCSD) 	PWG = 0.15
	call AskR4(' G-level power to fit density?',PWG)
 	if(bGcon) call AskR4('1 AU average density to fit base?',DEN1AU)
	if(bVEISC)	PWV = 0.35	
	if(bVNago)	PWV = 0.35
	if(bVGen) 	PWV = 0.35
	if(bVOoty)	PWV = 0.35
	if(PWV.ne.PWG)  PWV = PWG
	call AskR4('G-level power to fit velocity?',PWV)
	call AskR4('Radial G-level falloff to fit G density?',PWRG)
	if(bVEISC)            PWRV = PWRG
	if(bVNago.and.bGNago) PWRV = PWRG
	if(bVGen.and.bGGen)   PWRV = PWRG
	if(PWRV.ne.PWRG)      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

	call AskR4('Map Extraction distance (AU; 0=solar surface)?',RRX)
	if (RRX .lt. 0) RRX = -RRX*SUN__RAU			! Neg. RRX are in units of solar radii
	RRX = max(RRX,sngl(SUN__RAU))				! Keep RRX > 1 solar radius

C-------	Prompts for time

	if (bForecast) then			! Forecast mode

	    call Local2UT(iOffUT,iYr,Doy)		! calls the local computer time clock. With the offset calculates the Doy = universal time
	    call Julian(10,iYr,Doy,JD,JEpoch)		! Current time in MJD

	    write (*,'(/,10X,A,I6,A,I6,A,F10.3)') 'First MJD =',MJDfrst,
     &			'   Last MJD =',MJDlast,'   Current MJD =',JD

	    MJDref = JD					! MJD at center of plot (Earth position)
	    write (cStr,'(A,I5,A,F9.3,A)') '$',MJDfrst,'$',JD,'$0$'
          print*,' '
	    call AskR8('Modified Julian Day (0 = STOP)'//cStr,MJDref)
	    if (MJDref .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay')

C         print*,'MJDref, JD, JD- =',MJDref, JD, JD - 0.01d0
          if (MJDref .lt. (JD - 0.02d0)) then 
            if(iSys(4).eq.4.or.bVNago.or.bGnago) then         
              print *,' '
              if(.not.bForinput) call AskYN('Limit in-situ measurements to < DOY 0.79861?$yes',bForinput)
              if(bForinput) then
                print *,' '
                print *, 'bForinput is set .true. ISEE IPS inputs after DOY 0.79861 are removed.'
                print *,' '
                cWildVCELIAS = '/home/soft/dat/insitu/celias_????.hravg'            ! In the case CELIAS in-situ V is used
                cWildDCELIAS = '/home/soft/dat/insitu/celias_????.hravg'            ! In the case CELIAS in-situ D is used
              end if
            end if
          end if

C          if (cWildVCELIAS .eq. ' '.and.bForecast.and.bForinput) then
C             cWildVCELIAS = '/home/soft/dat/insitu/celias_????.hravg'
C          end if
C          if (cWildDCELIAS .eq. ' '.and.bForecast.and.bForinput) then
C             cWildDCELIAS = '/home/soft/dat/insitu/celias_????.hravg'
C          end if

C-------	Make sure iYr,Doy match MJDref
C		Get heliographic coordinates for center of plot (Earth position)
C	MJDref = MJDref + 7
	    call Julian(11,iYr,Doy,MJDref,JEpoch)
          DOYS = Doy                                  ! DOYS is used in the Forecast restriction if necessary
	    Idoy = int(Doy)
	    HR   = (Doy-float(Idoy))*24.0
	    IH   = int(HR)
	    IMN   = nint((HR-float(IH))*60.0)
	    call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy)
	    write(*,'(A,I2.2,A,I2.2,A,I2.2,A,A,I5,A,F8.3,A)') 
     &		'      (This MJD is ',IH,':',IMN,' UT ',iDay,' ',cMon,iYr,' -- DOY',Doy,')' 
	    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

	    call AdjustJDCar(XCbeg, nCar, JDCar, NCoff)

	    call T3D_set(T3D__TNOW, 0, XCbeg+0.5)

	    XCtest1 = XCbeg-XCmargin		! Search limits
	    XCtest2 = XCbeg+NC+XCmargin
	    XCend= XCbeg + NC
	    XCstrt = XCbeg
	    XCbeg= XCbeg - xInc
	    XCend= XCend + xInc
	    MJDcntr = MJDref
	    NC = 0
	    NCC = NC

	else						! Normal mode
	    write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$'
	    XCbeg = max(XClo,NCoff+XCbeg)
C	    XCbeg = 1964.6		! Bastille-day  event
C	    XCbeg = 2055.85 ! Used for EISCAT data run
C	    XCbeg = 2061.
C	    XCbeg = 2147.
C	    XCbeg = 2153.
C	    XCbeg = 2156.7
C	    XCbeg = 2056.
C	    XCbeg = 2209.8
C	    XCbeg = 2226.3
	    XCbeg = 2182.4

	    if(bGCamb) XCbeg = 1884
	    call AskR4('Carrington rotation (0 = STOP)?'//cStr,XCbeg)
	    if (XCbeg .le. 0) call Say('IPSD2020','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
	    NCC = NC

	    XCbeg = XCbeg-NCoff

	    call AdjustJDCar(XCbeg, nCar, JDCar, 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

      print*, ' '
	call AskYN('Do you want to output model-source files?$no',bWriteMS)
	call AskYN('Do you want to output l.o.s. source surface maps?$no',bWritelos)
	call AskYN('Do you want one shift matrix per iteration?$no',bOnly1)
	call AskYN('Output good source files?$yes',bSource)
      call AskYN('Output three D error files?$no',bWrerr)
	NOneOr1 = -2
	NOneOr2 = -2
	NinterC = 3										! Default first extract matrix
	NinterD = 3										! Default forecast extract matrix
        call AskYN('Output 3D files with holes?$no',bWrite3D1) 	  	! 3D 1 file write
	    if(bWrite3D1) then
		call AskI4('Number of intermediate rotation densities with holes?',NinterC)	! Affects first extract
		call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr1)
		NOneOr1 = 0
		if(bOneOr1) NOneOr1 = 1
	    end if
        call AskYN('Output 3D mostly filled files?$no',bWrite3D2) 	  			! 3D "f" file write
	    if(bWrite3D2) then
		call AskI4('Number of intermediate filled rotation densities?',NinterD)		! Affects forecast extract
		call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr2)
		NOneOr2 = 0
		if(bOneOr2) NOneOr2 = 1
	    end if

	nMapb   = 1
	nMapHRb = 1
	nMapHRMb = nMap

        call AskYN('Output 3D high resolution files?$yes',bWrite3D_HR)! 3D "h" file write
        if(bWrite3D_HR) then

	  RRRS = RRS
	  call AskR4('Minimum high-resolution height in AU?',RRRS)
	  if(RRRS.le.RRS) then
	    nMapb = 1
	  else
	    nMapb = nint((RRRS-RRS)/dRR) + 1
	  end if
	  if(nMapb.gt.nMap) nMapb = nMap				! Limit minimum height to maximum available 
	  RRRSM = 3.0
	  call AskR4('Maximum high-resolution height in AU?',RRRSM)
	  nMapHR = nint((RRRSM-RRS)/dRR) + 1
	  if(nMapHR.gt.nMap) nMapHR = nMap				! Limit maximum height to maximum available
	  if(nMapHR.lt.nMapb) nMapHR = nMapb				! Don't be stupid!

          call AskYN('Do you want the density error to limit HR densities?$no',bDdener)
          if(bDdener) then
            call ASKR4('What density L.O.S. crossing limit?',ERDLOS)
            call AskYN('Do you want the density error to limit HR velocities?$no',bVdener)
          end if
          call AskYN('Do you want the velocity error to limit HR velocities?$no',bVverer)
          if(bVverer) then
            bVdener = .FALSE. 
            call ASKR4('What velocity L.O.S. crossing limit?',ERVLOS)
            call AskYN('Do you want the velocity error to limit HR densities?$no',bDverer)
            if(bDverer) bDdener = .FALSE.
          end if
        end if

        call AskYN('Output completely filled 3D files?$no',bWrite3Do) ! Completely filled 3D "o" file write

        call AskYN('Output a 3D high resolution base?$yes',bWrite3D_HRb)! Completely filled 3D base "b" file write
        if(bWrite3D_HRb) then

	  RRRRS = RRS
	  call AskR4('Minimum high-resolution base height in AU?',RRRRS)
	  if(RRRRS.le.RRS) then
	    nMapHRb = 1
	  else
	    nMapHRb = nint((RRRRS-RRS)/dRR) + 1
	  end if
          if(nMapHRb.gt.nMap) nMapHRb=nMap				! Limit maximum height to maximum available

	  RRRRSM = 3.0
	  call AskR4('Maximum high-resolution base height in AU?',RRRRSM)
	  nMapHRMb = nint((RRRRSM-RRS)/dRR) + 1
	  if(nMapHRMb.gt.nMap) nMapHRMb = nMap				! Limit maximum height to maximum available
	  if(nMapHRMb.lt.nMapHRb) nMapHRMb = nMapHRb			! Don't be stupid!
          call AskYN('Do you want the den error to limit HR base densities?$no',bDdenerb)
          if(bDdenerb) then
            call ASKR4('What base density L.O.S. crossing limit?',ERDLOSB)
            call AskYN('Do you want the den error to limit HR base velocities?$no',bVdenerb)
          end if
          call AskYN('Do you want the vel error to limit HR base velocities?$no',bVvererb)
          if(bVvererb) then
            bVdenerb = .FALSE. 
            call ASKR4('What base velocity L.O.S. crossing limit?',ERVLOSB)
            call AskYN('Do you want the vel error to limit HR base densities?$no',bDvererb)
            if(bDvererb) bDdenerb = .FALSE.
          end if
        end if

C
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.

	bMagnetic = Bfield_choose_3(iVarp,cVar)
        call AskYN('Do you want to make magnetic source surface maps?$yes',bMagnetic)
	if (bMagnetic) then
	  bMagnet = .TRUE.
	  bMagnetic = .FALSE.
	  do J=1, B3D__NN
	    bWildMag(J) = .FALSE.
	    call ForeignFile(-iVarp-1,cVar,B3D__STRCMD(J),cWildMag(J))
	    if (itrim(cWildMag(J)) .gt. 0) then
	      bMagnetic   = .TRUE.
	      bWildMag(J) = .TRUE.
	    end if
	  end do
	  if (bMagnet.and..not.bMagnetic) then
	    print *, "You requested magnetic maps but didn't specify one correct path to find them."
	    stop
	  end if
	  if (bMagnetic) call AskYN('Do you want to produce magnetic field 3D files?$yes',bmagexx)
	  if (bMagnetic) call AskYN('Do you want to extract magnetic fields?$yes',bmagex)

	  BFACTORRR = 0.0
	  if(bWildMag(1).or.bWildMag(5).or.bWildMag(9).or.bWildMag(13))  BFACTORRR = 1.0	! radial open component factor
	  BFACTORRC = 0.0
	  if(bWildMag(2).or.bWildMag(6).or.bWildMag(10).or.bWildMag(14)) BFACTORRC = 1.0	! radial closed component factor
	  BFACTORTC = 0.0
	  if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.1.0) then
	    BFACTORRC = 2.0
	    if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0	! tangential closed component factor
	  end if	  
	  if(BFACTORRR.EQ.1.0.and.BFACTORRC.EQ.0.0) then
	    BFACTORRR = 2.0
	    if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0
	  end if

	  if((BFACTORRR+BFACTORRC).GE.1.0) then
	    if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0
	  end if
	  if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.0.0) then
	    if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 2.0
	  end if
	    BFACTORNC = 0.0
	  if(bWildMag(4).or.bWildMag(8).or.bWildMag(12).or.bWildMag(16)) BFACTORNC = 2.0

	  if (bMagnetic.and.BFACTORRR.ne.0.0) call AskR4('Multiplicative factor for CSSS fields?',BFACTORRR)
	  if (bMagnetic.and.BFACTORRC.ne.0.0) call AskR4('Multiplicative factor for radial closed fields?',BFACTORRC)
	  if (bMagnetic.and.BFACTORTC.ne.0.0) call AskR4('Multiplicative factor for tangential closed fields?',BFACTORTC)
	  if (bMagnetic.and.BFACTORNC.ne.0.0) call AskR4('Multiplicative factor for normal closed fields?',BFACTORNC)

	  if (bMagnetic.and.(BFACTORRR.ne.0.0.or.BFACTORRC.ne.0.0.or.BFACTORTC.ne.0.0.or.BFACTORNC.ne.0.0))
     &	  call AskYN('Adjust field factors for possible known effects?$no',bAdjust)
	  if(bAdjust) call Adjust_BF(bWildMag,B3D__NN,iYr,Doy,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC)
	end if

C     Regular Output files
      call AskYN('Output an Earth parameter file?$yes',bEWRITE)		      ! Earth extraction
      call AskYN('Output a Mercury parameter file?$yes',bMWRITE)		      ! Mercury extraction
      call AskYN('Output a Venus parameter file?$yes',bVEWRITE)		      ! Venus extraction
      call AskYN('Output a Mars parameter file?$yes',bMAWRITE)		      ! Mars extraction
	if(iYr.lt.2009) then
          call AskYN('Output a Ulysses parameter file?$no',bUWRITE)		! Ulysses extraction
	end if
        if(iYr.ge.2007) then
          call AskYN('Output a STEREO A parameter file?$yes',bSAWRITE)		! STEREO A extraction
          call AskYN('Output a STEREO B parameter file?$no',bSBWRITE)		! STEREO B extraction
        end if

	if(iYr.ge.2007) then
          call AskYN('Output a MESSENGER parameter file?$no',bMSWRITE)		! MESSENGER extraction
	end if
          if (iYr .eq. 2018) then
            if (iDoy .ge. 225) then
              call AskYN('Output a Parker Solar parameter file$no',bPSWRITE)
            else if (iDoy .ge. 294) then
              call AskYN('Output a Parker Solar parameter file$no',bPSWRITE)
              call AskYN('Output a BepiColumbo parameter file$no',bBCWRITE)
            end if
          end if
C
          if (iYr .gt. 2018) then
            call AskYN('output a Parker Solar parameter file$yes',bPSWRITE)
            call AskYN('Output a BepiColumbo parameter file$yes',bBCWRITE)
          end if
C
          if (iYr .eq. 2020) then
            if (iDoy .ge. 42) then
              call AskYN('Output a Solar Orbiter parameter file$yes',bSOWRITE)
            end if
          else if (iYr .gt. 2020) then
            call AskYN('Output a Solar Orbiter paramter file$yes',bSOWRITE)
          end if

	call AskYN('Output a Comet 67P parameter file?$no',b67WRITE)		! Comet 67P/GC extraction

C     Filled Files

      call AskYN('Output an Earth parameter filled file?$yes',bEWRITEF)	! Earth filled extraction
      call AskYN('Output a Mercury parameter filled file?$yes',bMWRITEF)	! Mercury filled extraction
      call AskYN('Output a Venus parameter filled file?$yes',bV2WRITEF)	! Venus extraction
      call AskYN('Output a Mars parameter filled file?$yes',bMAWRITEF)	! Mars filled extraction
	if(iYr.lt.2009) then
          call AskYN('Output a Ulysses parameter filled file?$no',bUWRITEF)	! Ulysses filled extraction
	end if
        if(iYr.ge.2007) then
          call AskYN('Output a STEREO A parameter filled file?$yes',bSAWRITEF)	! STEREO A filled extraction
          call AskYN('Output a STEREO B parameter filled file?$no',bSBWRITEF)	! STEREO B filled extraction
        end if
	if(iYr.ge.2007) then
          call AskYN('Output a MESSENGER parameter filled file?$no',bMXWRITEF)	! MESSENGER filled extraction
	end if

          if (iYr .eq. 2018) then
            if (iDoy .ge. 225) then
              call AskYN('Output a Parker Solar parameter filled file$no',bPSWRITEF)
            else if (iDoy .ge. 294) then
              call AskYN('Output a Parker Solar parameter filled file$no',bPSWRITEF)
              call AskYN('Output a BepiColumbo parameter filled file$no',bBCWRITEF)
            end if
          end if
C
          if (iYr .gt. 2018) then
            call AskYN('output a Parker Solar parameter filled file$yes',bPSWRITEF)
            call AskYN('Output a BepiColumbo parameter filled file$yes',bBCWRITEF)
          end if
C
          if (iYr .eq. 2020) then
            if (iDoy .ge. 42) then
              call AskYN('Output a Solar Orbiter parameter filled file$yes',bSOWRITEF)
            end if
          else if (iYr .gt. 2020) then
            call AskYN('Output a Solar Orbiter paramter filled file$yes',bSOWRITEF)
          end if

	call AskYN('Output a Comet 67P filled parameter file?$no',bC67WRITEF)	! Comet 67P/GC filled extraction

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.

	NLV  = 0
	NLG  = 0	
	do I=1,NLmaxV
	  IYRSV(I) = 0
	end do
	do I=1,NLmaxG
	  IYRSGG(I) = 0
	end do

	if (bVNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bVEISC) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bVOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEVsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bVGen) then

	  do ii=1, iSyst

	    if(iSys(ii).ne.0) then

	      if(ii.eq.3) then                          ! Ooty IPS velocities
	        XElowV       = XElowOV
	        bsinglesitev = bsinglesitevo
              bmultisitev  = bmultisitevo
              bVERR0OK     = bVERR0OKO
              VERRL        = VERRLS(ii)
	      end if
	      if(ii.eq.4) then                          ! ISEE IPS velocities
	        XElowV       = XElowNV
	        bsinglesitev = bsinglesitevn
              bmultisitev  = bmultisitevn
              bVERR0OK     = bVERR0OKN
              VERRL        = VERRLS(ii)
	      end if
	      if(ii.eq.5) then                          ! KSWC IPS velocities
	        XElowV       = XElowKV
	        bsinglesitev = bsinglesitevk
              bmultisitev  = bmultisitevk
              bVERR0OK     = bVERR0OKK
              VERRL        = VERRLS(ii)
	      end if
	      if(ii.eq.6) then                          ! MEXART IPS velocities
	        XElowV       = XElowMV
	        bsinglesitev = bsinglesitevm
              bmultisitev  = bmultisitevm
              bVERR0OK     = bVERR0OKM
              VERRL        = VERRLS(ii)
	      end if
	      if(ii.eq.7) then                          ! BSA3 IPS velocities
	        XElowV       = XElowPV
	        bsinglesitev = bsinglesitevb
              bmultisitev  = bmultisitevb
              bVERR0OK     = bVERR0OKB
              VERRL        = VERRLS(ii)
	      end if
	      if(ii.eq.8) then                          ! LOFAR IPS velocities
	        XElowV       = XElowLV
	        bsinglesitev = bsinglesitevl
              bmultisitev  = bmultisitevl
              bVERR0OK     = bVERR0OKL
              VERRL        = VERRLS(ii)
	      end if
	      
            NLVSS = 0
            NLVMS = 0
	      NLVB = NLV + 1

            if(bsinglesitev) then
	        call ReadVIPSn8ss_n(iSys(ii),iReadGen8ss_n,iProcessGen8ss_n,cWildGens(ii),
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,DVOBS,GGOBS,DGGOBS,VFREQ,VSIZE,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,VFsav,VSsav)
            end if
C
            NLVBS = NLV                         ! Intermediate velocity reads
            NLVSS = NLVBS - NLVB  + 1           ! Number of single-site speed reads

            if(bmultisitev) then
	        call ReadVIPSn8ms_n(iSys(ii),iReadGen8ms_n,iProcessGen8ms_n,cWildGens(ii),
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRV,XElowV,XEhighV,iEorWV,XRlimV,
     &		NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav,
     &		IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,DVOBS,GGOBS,DGGOBS,VFREQ,VSIZE,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,VFsav,VSsav)
            end if

   	      NLVE = NLV
   	      
            NLVMS = NLVE - NLVBS               !  Number of multi-site velocity reads

	      do III=NLVB,NLVE
C             
	        if(ii.eq.4) then					! This is ISEE velocity DATA
                if(III.eq.NLVB) then
                  print*,' '
                  write(*,'(A,I6)') 'ISEE single-site speed data read  ',NLVSS
                  write(*,'(A,I6)') 'ISEE multi-site velocity data read',NLVMS
		    end if
		    VNOerr = DVOBS(III)/VOBS(III)
                if(VNOerr.lt.0.0.and..not.bVERR0OK) VOBS(III) = -99.999
		    if(VNOerr.gt.VERRL) VOBS(III) = -99.999                       ! if bVERR0OK .TRUE. VERRL should be set to 1.0	
		    if(III.eq.NLVB) then                        ! First velocity value read in allows a Nagoya count and info print.
		      IVNAGOYA = 0
                  write(*,'(A,I4,A,I4)') 'ISEE velocity or speed data begin at ',NLVB,' and end at ',NLVE
		      if(VERRL.gt.0.9999) write(*,'(A)') 'ISEE velocity or speed data are input. No V or speed error is imposed.'
		      if(VERRL.le.0.9999) write(*,'(A,F10.7)') 'ISEE velocity or speed data are input. VERRL =', VERRL
		    end if
		    if(VOBS(III).gt.VLIML.and.VOBS(III).lt.VLIMM) then
C	            write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)')
C     &		    'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III)		
                  IVNAGOYA  = IVNAGOYA + 1		
                end if  
		    if(III.eq.NLVE) then
		      write(*,'(A,I4,A)') 'There are ',IVNAGOYA,' valid General Nagoya IPS velocity or speed sources for use.'
		    end if	
	        end if
C
	        if(ii.eq.6) then					! This is MEXART Speed DATA
                if(III.eq.NLVB) then
                  print*,' '
                  write(*,'(A,I6)')   'MEXART single-site speed data read  ',NLVSS
                  write(*,'(A,I6,A)') 'MEXART multi-site velocity data read',NLVMS,' (none are available).'
		    end if
		    VMEXerr = DVOBS(III)/VOBS(III)
                if(VMEXrr.lt.0.0.and..not.bVERR0OK) VOBS(III) = -99.999
		    if(VMEXerr.gt.VERRL) VOBS(III) = -99.999                       ! if bVERR0OK .TRUE. VERRL should be set to 1.0
		    if(III.eq.NLVB) then
		      IVMEX = 0
                  write(*,'(A,I4,A,I4)') 'MEXART IPS speed data begin at ',NLVB,' and end at ',NLVE
		      if(VERRL.gt.0.9999) write(*,'(A)') 'MEXART IPS speed data are input. No speed speed error is imposed.'
		      if(VERRL.le.0.9999) write(*,'(A,F10.7)') 'MEXART IPS speed data are input. VERRL =', VERRL
		    end if
	          write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)') 
     &		    'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III)
		    if(VOBS(III).gt.VLIML.and.VOBS(III).lt.VLIMM) IVMEX  = IVMEX + 1
		    if(III.eq.NLVE) then
		      write(*,'(A,I4,A)') 'There are ',IVMEX,' valid General MEXART IPS speed sources for use.'
		    end if
	        end if
C
	        if(ii.eq.7) then		! This is BSA (Pushchino) speed data
                if(III.eq.NLVB) then
                  print*,' '
                  write(*,'(A,I6)')   'Pushchino single-site speed data read  ',NLVSS
                  write(*,'(A,I6,A)') 'Pushchino multi-site velocity data read',NLVMS,' (none are available).'
		    end if
		    VPUSHerr = DVOBS(III)/VOBS(III)
                if(VPUSHerr.lt.0.0.and..not.bVERR0OK) VOBS(III) = -99.999
		    if(VPUSHerr.gt.VERRL) VOBS(III) = -99.999                       ! if bVERR0OK .TRUE. VERRL should be set to 1.0
		    if(III.eq.NLVB) then
		      IVPUSH = 0
                  write(*,'(A,I4,A,I4)') 'Pushchino IPS speed data begin at ',NLVB,' and end at ',NLVE
		      if(VERRL.gt.0.9999) write(*,'(A)') 'Pushchino IPS speed data are input. No speed speed error is imposed.'
		      if(VERRL.le.0.9999) write(*,'(A,F10.7)') 'Puschino IPS speed data are input. VERRL =', VERRL
		    end if
		    if(VOBS(III).le.VLIMM.and.VOBS(III).ge.VLIML) then
	            write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)') 
     &		    'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III)
		      IVPUSH  = IVPUSH + 1
		    end if	  		
		    if(III.eq.NLVE) then
		      write(*,'(A,I4,A)') 'There are ',IVPUSH,'  valid General Pushchino IPS speed sources for use.'
		    end if
	        end if
C
	        if(ii.eq.8) then		! This is LOFAR Velocity DATA
                if(III.eq.NLVB) then
                  print*,' '
                  write(*,'(A,I6)') 'LOFAR single-site speed data read  ',NLVSS
                  write(*,'(A,I6)') 'LOFAR multi-site velocity data read',NLVMS
		    end if
		    VLOerr = DVOBS(III)/VOBS(III)
                if(VLOerr.lt.0.0.and..not.bVERR0OK) VOBS(III) = -99.999
		    if(VLOerr.gt.VERRL) VOBS(III) = -99.999                       ! if bVERR0OK .TRUE. VERRL should be set to 1.0
		    if(III.eq.NLVB) then
		      IVLOFAR = 0
                  write(*,'(A,I4,A,I4)') 'LOFAR IPS velocity or speed data begin at ',NLVB,' and end at ',NLVE
		      if(VERRL.gt.0.9999) write(*,'(A)') 'LOFAR IPS velocity or speed data are input. No V or speed error is imposed.'
		      if(VERRL.le.0.9999) write(*,'(A,F10.7)') 'LOFAR IPS velocity or speed data are input. VERRL =', VERRL
		    end if
		    if(VOBS(III).le.VLIMM.and.VOBS(III).ge.VLIML) then
C	            write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)') 
C     &		    'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III)
		      IVLOFAR = IVLOFAR + 1
		    end if
		    if(III.eq.NLVE) then
		      write(*,'(A,I4,A)') 'There are ',IVLOFAR,' valid General LOFAR IPS velocity or speed sources for use.'
		    end if
	        end if
C
	      end do

	    end if

	  end do	
        print*,' '
        print*,'End of velocity reads and winnowing.'
C	  Sort the IPS velocity array input parameters
	  
	  call sort_IPS_data3(NLV,NLMAXV,RADUMV8,cvgfiles,SRCVG,SRCV,IYRFV,IRECV,
     &			     XCV,YLV,DISTV,VOBS,DVOBS,VFREQ,VSIZE,XEV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8)

	end if

	print *, ' '

	if (bGCamb) call ReadGIPS8('$CAM',nCar,JDCar,
     &		iEdt,.FALSE.,bForecast,
     &		XCtest1,XCtest2,MJDref,MJDfrst,MJDlast,
     &		RRG,Speed,GPOWER,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxG,NLG,
     &		iXPG,iYPG,iMJDG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,GGOBS,
     &		1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,
     &		XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bGOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bGNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya,
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	if (bGGen) then

	  do ii=1, iSyst
	    if(iSys(ii).ne.0) then

	      if(ii.eq.3) then                          ! Ooty IPS g-levels
	        XElowG       = XElowOG
	        bsinglesiteg = bsinglesitego
              bmultisiteg  = bmultisitego
              bGERR0OK     = bGERR0OKO
              GERRL        = GERRLS(ii)
	      end if
	      if(ii.eq.4) then                          ! Nagoya IPS g-levels
	        XElowG       = XElowNG
	        bsinglesiteg = bsinglesitegn
              bmultisiteg  = bmultisitegn
              bGERR0OK     = bGERR0OKN
              GERRL        = GERRLS(ii)
	      end if
	      if(ii.eq.5) then                          ! KSWC IPS g-levels
	        XElowG       = XElowKG
	        bsinglesiteg = bsinglesitegk
              bmultisiteg  = bmultisitegk
              bGERR0OK     = bGERR0OKK
              GERRL        = GERRLS(ii)
	      end if
	      if(ii.eq.6) then                          ! MEXART IPS g-levels
	        XElowG       = XElowMG
	        bsinglesiteg = bsinglesitegm
              bmultisiteg  = bmultisitegm
              bGERR0OK     = bGERR0OKM
              GERRL        = GERRLS(ii)
	      end if
	      if(ii.eq.7) then                          ! Pushchino IPS g-levels
	        XElowG       = XElowPG
	        bsinglesiteg = bsinglesitegb
              bmultisiteg  = bmultisitegb
              bGERR0OK     = bGERR0OKB
              GERRL        = GERRLS(ii)
	      end if
	      if(ii.eq.8) then                          ! LOFAR IPS g-levels
	        XElowG       = XElowLG
	        bsinglesiteg = bsinglesitegl
              bmultisiteg  = bmultisitegl
              bGERR0OK     = bGERR0OKL
              GERRL        = GERRLS(ii)
	      end if
	      
            NLGSS = 0
            NLGMS = 0
	      NLGB = NLG + 1

            if(bsinglesiteg) then
	        call ReadVIPSn8ss_n(-iSys(ii),iReadGen8ss_n,iProcessGen8ss_n,cWildGens(ii),
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,DVOBS,GGOBS,DGGOBS,GFFREQ,GSSIZE,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,GFsav,GSSsav)
            end if

            NLGBS = NLG                               !  Intermediate g-level read value
            NLGSS = NLGBS - NLGB  + 1                 !  Number of single-site g-level reads
C
            if(bmultisiteg) then
	        call ReadVIPSn8ms_n(-iSys(ii),iReadGen8ms_n,iProcessGen8ms_n,cWildGens(ii),
     &		nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref,
     &		RRG,XElowG,XEhighG,iEorWG,XRlimG,
     &		NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav,
     &		IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,DVOBS,GGOBS,DGGOBS,GFFREQ,GSSIZE,
     &		1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,GFsav,GSSsav)            
            end if           
	      NLGE = NLG  	      
            NLGMS = NLGE - NLGBS     	            !  Number of multi-site g-level reads
C
	      do III=NLGB,NLGE
	        if(ii.eq.3) then						            ! Ooty IPS g-level data
                if(III.eq.NLGB) then
                  print*,' '
                  write(*,'(A,I6)')   'Ooty single-site g-level data read',NLGSS
                  write(*,'(A,I6,A)') 'Ooty multi-site g-level data read ',NLGMS,' (none are available).'
                end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMOF + 1.0
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMOF
		    GOotyerr = DGGOBS(III)/GGOBS(III)
                if(GOotyerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GOotyerr.gt.GERRL) GGOBS(III) = -99.999                     ! if bGERR0OK .TRUE. GERRL should be set to 1.0
		    if(III.eq.NLGB) then
		      IGOoty = 0
                  write(*,'(A,I4,A,I4)') 'Ooty IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'Ooty g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'Ooty g-level data are input. GERRL Fraction =', GERRL
		    end if
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') 
     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		      IGOoty  = IGOoty + 1		  		
		    end if
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGOoty,' valid General Ooty IPS g-level sources for use.'
		    end if
		  end if
C
	        if(ii.eq.4) then						            ! ISEE IPS g-level data
                if(III.eq.NLGB) then
                  print*,' ' 
                  write(*,'(A,I6)') 'ISEE single-site g-level data read',NLGSS
                  write(*,'(A,I6)') 'ISEE multi-site g-level data read ',NLGMS
                end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMNF + 1.0	      
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMNF
		    GNOerr = DGGOBS(III)/GGOBS(III)
                if(GNOerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GNOerr.gt.GERRL) GGOBS(III) = -99.999                       ! if bGERR0OK .TRUE. GERRL should be set to 1.0	
		    if(III.eq.NLGB) then
		      IGNAGOYA  = 0
                  write(*,'(A,I4,A,I4)') 'Nagoya IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'Nagoya g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'Nagoya g-level data are input. GERRL Fraction =', GERRL
		    end if 		
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
		      IGNAGOYA = IGNAGOYA + 1
C	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') 
C     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		    end if
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGNAGOYA,' valid General Nagoya IPS g-level sources for use.'
		    end if
		  end if
C
		  if(ii.eq.5) then						             ! KSWC IPS g-level data
                if(III.eq.NLGB) then
                  print*,' '
                  write(*,'(A,I6)')   'KSWC single-site g-level data read',NLGSS
                  write(*,'(A,I6,A)') 'KSWC multi-site g-level data read ',NLGMS,' (none are available).'
		    end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMKF + 1.0	      
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMKF
		    GKSWCerr = DGGOBS(III)/GGOBS(III)
                if(GKSWCerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GKSWCerr.gt.GERRL) GGOBS(III) = -99.999                     ! if bGERR0OK .TRUE. GERRL should be set to 1.0
		    if(III.eq.NLGB) then
		      IGKSWC = 0
                  write(*,'(A,I4,A,I4)') 'KSWC IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'KSWC g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'KSWC g-level data are input. GERRL Fraction =', GERRL
		    end if
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') 
     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		      IGMEX  = IGMEX + 1		  		
		    end if
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGMEX,' valid General KSWC IPS g-level sources for use.'
		    end if
		  end if
C
	        if(ii.eq.6) then						            ! MEXART IPS g-level DATA
                if(III.eq.NLGB) then
                  print*,' '
                  write(*,'(A,I6)')   'MEXART single-site g-level data read',NLGSS
                  write(*,'(A,I6,A)') 'MEXART multi-site g-level data read ',NLGMS,' (none are available).'
		    end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMMF + 1.0	      
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMMF
		    GMEXerr = DGGOBS(III)/GGOBS(III)
                if(GMEXerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GMEXerr.gt.GERRL) GGOBS(III) = -99.999                     ! if bGERR0OK .TRUE. GERRL should be set to 1.0
		    if(III.eq.NLGB) then
		      IGMEX = 0
                  write(*,'(A,I4,A,I4)') 'MEXART IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'MEXART g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'MEXART g-level data are input. GERRL Fraction =', GERRL
		    end if
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') 
     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		      IGMEX  = IGMEX + 1		  		
		    end if
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGMEX,' valid General MEXART IPS g-level sources for use.'
		    end if
	        end if
C
	        if(ii.eq.7) then							        ! BSA (Pushchino) IPS g-level DATA
                if(III.eq.NLGB) then
                  print*,' '
                  write(*,'(A,I6)')   'Pushchino single-site g-level data read',NLGSS
                  write(*,'(A,I6,A)') 'Pushchino multi-site g-level data read ',NLGMS,' (none are available).'
		    end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMPF + 1.0
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMPF	      
		    GBSAerr = DGGOBS(III)/GGOBS(III)
                if(GBSAerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GBSAerr.gt.GERRL) GGOBS(III) = -99.999                       ! if bGERR0OK .TRUE. GERRL should be set to 1.0
		    if(III.eq.NLGB) then
		      IGPUSH = 0
                  write(*,'(A,I4,A,I4)') 'BSA IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'BSA g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'BSA g-level data are input. GERRL Fraction =', GERRL
		    end if
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)')  
     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		      IGPUSH  = IGPUSH + 1
		    end if
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGPUSH,' valid General BSA IPS g-level sources for use.'
		    end if
	        end if
C
	        if(ii.eq.8) then							        ! LOFAR IPS g-level DATA
                if(III.eq.NLGB) then       
                  print*,' '
                  write(*,'(A,I6)') 'LOFAR single-site g-level data read',NLGSS
                  write(*,'(A,I6)') 'LOFAR multi-site g-level data read ',NLGMS
		    end if
		    GGOBS(III)  = (GGOBS(III)-1.0)*GLIMLF + 1.0	      
		    if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMLF
		    GLOerr = DGGOBS(III)/GGOBS(III)
                if(GLOerr.lt.0.0.and..not.bGERR0OK) GGOBS(III) = -99.999
		    if(GLOerr.gt.GERRL) GGOBS(III) = -99.999                        ! if bGERR0OK .TRUE. GERRL should be set to 1.0	
		    if(III.eq.NLGB) then
		      IGLOFAR = 0
                  write(*,'(A,I4,A,I4)') 'LOFAR IPS g-level data begin at ',NLGB,' and end at ',NLGE
		      if(GERRL.gt.0.9999) write(*,'(A)') 'LOFAR g-level data are input. No g-level error is imposed.'
		      if(GERRL.le.0.9999) write(*,'(A,F10.7)') 'LOFAR g-level data are input. GERRL Fraction =', GERRL
		    end if
		    if(GGOBS(III).le.GLIMP.and.GGOBS(III).ge.GLIMM) then
C	            write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') 
C     &		    'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III)
		      IGLOFAR  = IGLOFAR + 1
		    end if		  		
		    if(III.eq.NLGE) then
		      write(*,'(A,I4,A)') 'There are ',IGLOFAR,' valid General LOFAR IPS g-level sources for use.'
	          end if
	        end if
	      end do
	    end if
	  end do
        print*,' '
        print*,'End of g-level reads and winnowing.'
C	  Sort the IPS g-level array input parameters

	  call sort_IPS_data3(NLG,NLMAXG,RADUMG8,cggfiles,SRCGV,SRCGG,IYRFG,IRECG,
     &				XCGG,YLGG,DISTGG,GGOBS,DGGOBS,GFFREQ,GSSIZE,XEGG,XLSGG,XCEGG,XLLGG,XDLGG,IYRSGG,DOYSGG8)
	end if

	print *, ' '
	print *, 'Before the in-situ reads bVACE,bVACE2,bVWind,bVCELIAS,bVDSCOVR', bVACE, bVACE2, bVWind, bVCELIAS, bVDSCOVR
	print *, 'Before the in-situ reads bDACE,bDACE2,bDWind,bDCELIAS,bDDSCOVR', bDACE, bDACE2, bDWind, bDCELIAS, bDDSCOVR
	print *, ' '

	N_INV = 0

	if (bVACE) call ReadACE8_v(cWildACE,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV,
     &		NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)
     
	if (bVACE2) call ReadACE28_v(cWildACE2,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV,
     &		NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)

	if (bVWind) call ReadWind8_v(cWildWind,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV,
     &		NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)

	if (bVCELIAS.and..not.bForecast) call ReadCELIAS8_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,
     &		RRV,NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)

	if (bVCELIAS.and.bForecast.and..not.bForinput) call ReadCELIAS8R_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,
     &		RRV,NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)

	if (bVCELIAS.and.bForecast.and.bForinput) call ReadCELIAS8_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,
     &		RRV,NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)

	if (bVDSCOVR) call ReadDSCOVR8_v(cWildDSCOVR,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV,
     &		NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV,
     &		IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS)
	N_ING = 0

	if (bDACE) call ReadACE8_d(cWildACE,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG,
     &		NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	if (bDACE2) call ReadACE28_d(cWildACE2,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG,
     &		NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	if (bDWind) call ReadWind8_d(cWildWind,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG,
     &		NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	if (bDCELIAS.and..not.bForecast) call ReadCELIAS8_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,
     &		RRG,NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	if (bDCELIAS.and.bForecast.and..not.bForinput) call ReadCELIAS8R_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,
     &		RRG,NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	if (bDCELIAS.and.bForecast.and.bForinput) call ReadCELIAS8_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,
     &		RRG,NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)
     
	if (bDDSCOVR) call ReadDSCOVR8_d(cWildDSCOVR,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG,
     &		NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG,
     &		IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS)

	call ArrR4Zero(NLmaxV,VSSIG)
	call ArrI4Constant(NLmaxV,1 ,NBSV )
	call ArrR4Zero(NLmaxG,GSSIG)
	call ArrI4Constant(NLmaxG,1 ,NBSG )

	if(bforecast) write(*,'(A,F10.4,A,I5)') 'The Forecast cut-off time is DOY',Doy,' of',iYr

	write(*,'(A,4I10)') 'NLG, NLV, N_ING, N_INV',NLG,NLV,N_ING,N_INV

	LASTVVALUE = 0
	ALASTVDOY  = 0.0
	LASTViYr   = 0
	do I=1,NLV
C	    if (I.lt.100) print *, I, SRCV(I), SRCGG(I)
	    if (VOBS(I) .gt. VLIMM) then	
		NBSV(I) = 0
		VOBS(I) = BadR4()
	    end if
	    if (VOBS(I) .lt. VLIML) then
		NBSV(I) = 0
		VOBS(I) = BadR4()
	    end if
	    if(bforecast) then						! Provides a sharp cut-off in the velocity forecast lines of sight.
	      if(IYRSV(I).ge.iYr.and.DOYSV8(I).gt.dble(Doy).or.IYRSV(I).eq.0) then
	        if(LASTVVALUE.eq.0) then
	          write(*,'(A,A,A,F10.4,A,I5)') 'The last valid IPS velocity LOS was ',SRCVLAST,' on',ALASTVDOY,' of',LASTViYr
	          LASTVVALUE = 1
	        end if
	        NBSV(I) = 0
		VOBS(I) = BadR4()
	      end if
	      if(I.le.NLV.and.NBSV(I).ne.0) then
	        ALASTVDOY = sngl(DOYSV8(I))
	        LASTViYr  = IYRSV(I)
	        SRCVLAST  = SRCV(I)
	      end if
	    end if
	end do

	LASTGVALUE = 0
	ALASTGDOY  = 0.0
	LASTGiYr   = 0
	do I=1,NLG
	    if(GLIMF.ne.1.0) then						! g-level site specific multiplicative factor not specified
		  GGOBS(I)  = (GGOBS(I)-1.0)*GLIMF + 1.0	      
		  if(DGGOBS(I).gt.0.0) DGGOBS(I) = DGGOBS(I)*GLIMF
	    end if	
	    if (GGOBS(I) .gt. GLIMP) then
		NBSG(I) = 0
		GGOBS(I) = BadR4()
	    end if
	    if (GGOBS(I) .lt. GLIMM) then
		NBSG(I) = 0
		GGOBS(I) = BadR4()
	    end if
	    if(bforecast) then						! Provides a sharp cut-off in the g-level forecast lines of sight.
C	    print *, I, NLmaxG, IYRSGG(I), iYr, DOYSGG8(I), dble(Doy), LASTGVALUE 
	      if(IYRSGG(I).ge.iYr.and.DOYSGG8(I).gt.dble(Doy).or.IYRSGG(I).eq.0) then
	        if(LASTGVALUE.eq.0) then
	          write(*,'(A,A,A,F10.4,A,I5)') 'The last valid IPS g-level LOS was  ',SRCGLAST,' on',ALASTGDOY,' of',LASTGiYr
	          LASTGVALUE = 1
	        end if
	        NBSG(I) = 0
		  GGOBS(I) = BadR4()
	      end if      
	      if(I.le.NLG.and.NBSG(I).ne.0) then
	        ALASTGDOY = sngl(DOYSGG8(I))
	        LASTGiYr  = IYRSGG(I)
	        SRCGLAST  = SRCGG(I)
	      end if
	    end if
	end do

	LASTAVALUE = 0
	ALASTADOY  = 0.0
	LASTAiYr   = 0
	do I=NLV+1,NLmaxV
	    if(bforecast) then						! Provides a sharp cut-off in the in-situ V forecast measurements.
	      if(IYRSV(I).ge.iYr.and.DOYSV8(I).gt.dble(Doy)) then
	        if(LASTAVALUE.eq.0) then
	          write(*,'(A,A,A,F10.4,A,I5)'), 'The last valid in-situ velocity was ',SRCALAST,' on',ALASTADOY,' of',LASTAiYr
	          LASTAVALUE = 1
	        end if
	        NBSV(I) = 0
		VOBS(I) = BadR4()
	      end if	    
	      if(I.le.NLV+N_INV.and.NBSV(I).ne.0) then
	        ALASTADOY = sngl(DOYSV8(I))
	        LASTAiYr  = IYRSV(I)
	        SRCALAST  = SRCV(I)
	      end if
	    end if
	end do
	
	LASTAVALUE = 0
	ALASTADOY  = 0.0
	LASTAiYr   = 0
	NINSITUBEND = 0
	do I=NLG+1,NLmaxG
	    if(bforecast) then						! Provides a sharp cut-off in the in-situ D forecast measurements.
	      if(IYRSGG(I).ge.iYr.and.DOYSGG8(I).gt.dble(Doy)) then
	        if(LASTAVALUE.eq.0) then
	          write(*,'(A,A,A,F10.4,A,I5)') 'The last valid in-situ density was ',SRCALAST,' on',ALASTADOY,' of',LASTAiYr
	          LASTAVALUE = 1
	        end if
	        NBSG(I) = 0
		NINSITUBEND = NINSITUBEND + 1
C               write(*,*) 'GGOBS(I)', GGOBS(I), ' IYRSGG(I)',IYRSGG(I), ' DOYSGG8(I)', DOYSGG8(I)
		GGOBS(I) = BadR4()
	      end if	    
	      if(I.le.NLG+N_ING.and.NBSG(I).ne.0) then
	        ALASTADOY = sngl(DOYSGG8(I))
	        LASTAiYr  = IYRSGG(I)
	        SRCALAST  = SRCV(I)
	      end if
	    end if
	    if (GGOBS(I) .gt. DMAX) then
		NBSG(I) = 0
		GGOBS(I) = BadR4()
	    end if
	    if (GGOBS(I) .lt. DMIN) then
		NBSG(I) = 0
		GGOBS(I) = BadR4()
	    end if
	end do

	if(NINSITUBEND.gt.0) write(*,'(A,I5)') 'The number of in-situ density points beyond the cut-off =',NINSITUBEND

	if(bFACE) write(*,'(A,F8.3)') 'The 1 AU upper limit of the in-situ density is:',DMAX
	if(bFACE) write(*,'(A,F8.3)') 'The 1 AU lower limit of the in-situ density is:',DMIN

C=======
     
	II   = 0
	III  = 0
	do I=1,NLG+N_ING
	  II = II + 1
	  if(bGCamb) then
	    XCG  (NLG-II+1)	= XCGG(I)
	    YLG  (NLG-II+1)	= YLGG(I)
	    DISTG(NLG-II+1)	= DISTGG(I)
	    GOBS (NLG-II+1)	= GGOBS(I)
	    DGOBS (NLG-II+1)	= DGGOBS(I)
	    GFREQ (NLG-II+1)	= GFFREQ(I)
	    GSIZE (NLG-II+1)	= GSSIZE(I)
	    XEG  (NLG-II+1)	= XEGG(I)
	    XLSG (NLG-II+1)	= XLSGG(I)
	    XCEG (NLG-II+1)	= XCEGG(I)
	    XLLG (NLG-II+1)	= XLLGG(I)
	    XDLG (NLG-II+1)	= XDLGG(I)
	    DOYSG8(NLG-II+1)	= DOYSGG8(I)
	    IYRSG(NLG-II+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
	    SRCGG(NLG-II+1) = SRCGS

	  end if
	  if(bGNago.or.bGOoty.or.bGGen) then
	    XCG(II)	   = XCGG(I)
	    YLG(II)	   = YLGG(I)
	    DISTG(II)  = DISTGG(I)
	    GOBS(II)   = GGOBS(I)
	    DGOBS(II)  = DGGOBS(I)
	    GFREQ(II)  = GFFREQ(I)
	    GSIZE(II)  = GSSIZE(I)
	    XEG(II)	   = XEGG(I)
	    XLSG(II)   = XLSGG(I)
	    XCEG(II)   = XCEGG(I)
	    XLLG(II)   = XLLGG(I)
	    XDLG(II)   = XDLGG(I)
	    DOYSG8(II) = DOYSGG8(I)
	    IYRSG(II)  = IYRSGG(I)

	    if(bGGen.and.GFFREQ(I).eq.140.0.and.NDUP.ne.0) then		! MEXART and other data
	      do IDUP=1,NDUP-1
		III = III + 1
		II  = II  + 1
		XCG(II)		= XCGG(I)
		YLG(II)		= YLGG(I)
		DISTG(II)	= DISTGG(I)
		GOBS(II)	= GGOBS(I)
		DGOBS(II)	= DGGOBS(I)
		GFREQ(II)	= GFFREQ(I)
		GSIZE(II)	= GSSIZE(I)
		XEG(II)		= XEGG(I)
		XLSG(II)	= XLSGG(I)
		XCEG(II)	= XCEGG(I)
		XLLG(II)	= XLLGG(I)
		XDLG(II)	= XDLGG(I)
		DOYSG8(II)	= DOYSGG8(I)
		IYRSG(II)	= IYRSGG(I)
	      end do
	    end if

	  end if
	end do
	NLG = NLG + III

       ALng = 1.*(nLng-1)/(iLng-1)
       
C       print *, 'Spot check of IPS data - 30', IYRSG(30),IYRSV(30),DOYSG8(30),DOYSV8(30),XEG(30),XEV(30),GOBS(30),VOBS(30)

        if(bForecast) then				! For forecast mode
	  NTV = nint((float(NTV) + 3.0)/2.0) + 7
	  NTG = nint((float(NTG) + 3.0)/2.0) + 7
        end if

	print *, ' '
        write(*,'(A,2I5)') ' Before MkTimes NTV, NTG ',NTV, NTG
	write(*,'(2(A,I9))') ' NLV =',NLV,'  NLG =',NLG
	print *, ' '

	NLVVV = 0
	do I=1,NLV
          if (VOBS(I) .ge. VLIML .and. VOBS(I) .le. VLIMM) NLVVV = NLVVV + 1
	end do
	do I=NLV+1,NLV+N_INV
          if (VOBS(I) .ge. VLIML .and. VOBS(I) .le. VLIMM) NLVVV = NLVVV + 1
	end do
	if (NLVVV .eq. 0) then
	  print *, 'There are no valid IPS or in-situ velocity measurements. Set bVcon to false.'
	  bVcon = .FALSE.
	end if

	NLGGG = 0
	do I=1,NLG
          if (GOBS(I) .ge. GLIMM .and. GOBS(I) .le. GLIMP) NLGGG = NLGGG + 1
	end do
	do I=NLG+1,NLG+N_ING
          if (GOBS(I) .ge. DMIN .and. GOBS(I) .le. DMAX) NLGGG = NLGGG + 1
	end do
	if (NLGGG .eq. 0) then
	  print *, 'There are no valid g-level or in-situ density measurements. Set bGcon to false.'
	  bGcon = .FALSE.
	end if

        if(.not.bVcon.and..not.bGcon) then
	  print *, ' STOP!!!!!!!'
          print *, ' There are no sources for either a velocity' 
	  print *, ' or a g-level deconvolution during this interval'
	  stop
        end if

	print *, 'Before MkTimes MJDcntr =', MJDcntr
	print *, ' '
	
	ItoffG = 0
	ItoffV = 0
	IdaybegV  = 365
	IdaybegG  = 365
	IYRBV  = IYRSV(1)
	IYRBG  = IYRSG(1)

        if(NLV.eq.0) bVcon = .FALSE.
	if(bVcon)
     &	    call MkTimes(bForecast,MJDcntr,NmidHRV,aNdayV,NTmaxV,ALng,nLng,nCar,JDCar,
     &		0,NTV,IYRBV,IdaybegV,XCintDV,XCintV,XCbeV,XCtbegV,XCtendV,IYRV,DOYV)
     
        if(NLG.eq.0) bGcon = .FALSE.
	if(bGcon)
     &	    call MkTimes(bForecast,MJDcntr,NmidHRG,aNdayG,NTmaxG,ALng,nLng,nCar,JDCar,
     &		0,NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG)

C       print *, 'SRCVG(30)',SRCVG(30)
C       print *, 'Spot check of IPS data - 30', IYRSG(30),IYRSV(30),DOYSG8(30),DOYSV8(30),XEG(30),XEV(30)

C       print *, 'SRCVG(409)',SRCVG(409)
C       print *, 'Spot check of IPS data - 409', IYRSG(409),IYRSV(409),DOYSG8(409),DOYSV8(409),XEG(409),XEV(409)

       write (*,'(A,4I5)') ' After MkTimes IYRBV, IYRBG, IdaybegV, IdaybegG',IYRBV,IYRBG,IdaybegV,IdaybegG

       if(.not.bGcon) then
	   nTG     = nTV
	   IYRBG   = IYRBV
	   ItoffG  = ItoffV
	   XCtbegG = XCtbegV
	   XCtendG = XCtendV
	   do I=1,nTmaxG
	       XCintG(I)   = XCintV(I)
	       XCintDG(I)  = XCintDV(I)
	       XCbegG(1,I) = XCbeV(1,I)
	       XCbegG(2,I) = XCbeV(2,I)
	   end do
       end if

       if(.not.bVcon) then
	   nTV     = nTG
	   IYRBV   = IYRBG
	   ItoffV  = ItoffG
	   XCtbegV = XCtbegG
	   XCtendV = XCtendG
	   do I=1,nTmaxV
	       XCintV(I)  = XCintG(I)
	       XCintDV(I) = XCintDG(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 
	       XCintDV(I) = XCintDG(I)
	       XCintV(I)  = XCintG(I)
	       XCbeV(1,I) = XCbegG(1,I)
	       XCbeV(2,I) = XCbegG(2,I)
	   end do
       end if

       write(*,'(A,2F12.5,I12)') ' XCtbegV, XCtendV, NTV', XCtbegV, XCtendV, NTV
       write(*,'(A,2F12.5,I12)') ' XCtbegG, XCtendG, NTG', XCtbegG, XCtendG, NTG

       write(*,'(/A)') 	 ' XCbeV ='
       write(*,'(8F9.6)')  XCbeV       
       write(*,'(/A)')	 ' XCbeGG ='
       write(*,'(8F9.6)')  XCbeGG
       write(*,'(/A)') 	 ' XCintV ='
       write(*,'(8F9.6)')  XCintV
       write(*,'(/A)')	 ' XCintG ='
       write(*,'(8F9.6)')  XCintG
       write(*,'(/A)')	 ' XCintDV ='
       write(*,'(8F9.4)')  XCintDV
       write(*,'(/A)')	 ' XCintDG ='
       write(*,'(8F9.4)')  XCintDG

       XCtbegF = XCtbegG				! For forecast mode
       XCtendF = XCtendG				! For forecast mode

       if(NLVVV.lt.NLVmin.and.N_INV.lt.NLVmin) then
         bVcon = .FALSE.
         NLV = 0
	 write(*,'(A,I3,A)') 'There are too few (',NLVVV+N_INV,') sources for a velocity deconvolution. NLV set to zero.'
	 print *, ' '
       end if
       
       if(NLGGG.lt.NLGmin.and.N_ING.lt.NLGmin) then
	 bGcon = .FALSE.
         NLG = 0
	 write(*,'(A,I3,A)') 'There are too few (',NLGGG+N_ING,') sources for a g-level deconvolution. NLG set to zero.'
	 print *, ' '
       end if

       if(.not.bVcon.and..not.bGcon) then
	 print *, ' STOP!!!!!!!'
	 write(*,'(2(A,I9))') 'NLV =',NLV,' NLG =', NLG
         print *, ' There are too few sources for either a velocity' 
	 print *, ' or a g-level deconvolution'
	 print *, ' You will need to set the lower limit of either the velocity'
	 print *, ' or g-level LOS (or both) to a smaller number.'
	 stop
       end if

       if(bGcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP)
       if(bVcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),VMAP)

       NLmax = max(NLmaxG,NLmaxV)

       NLVBEG = 0
       NLGBEG = 0    
       NLVEND = NLV
       NLGEND = NLG    

	print *, 'NLVEND =', NLVEND, 'NLGEND =',NLGEND

       NTVV = ItoffV
       NTGG = ItoffG
       NTVE = 0
       NTGE = 0
       K    = 0
       Nbadvel = 0
       Nbadden = 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.NLVEND+1) then	! If K is more than the source # + 1 don't turn on switch
C	print *, ' Before V switch loop K, NTVV, NLmax', K,NTVV,NLmax,XCEV(K),IYRSV(NTVV+2),DOYV(NTVV+2),XCintV(NTVV+2)
	  if(bVcon.and.(XCEV(K).gt.XCintV(NTVV+2).and.XEV(K).le.XEhighV).or.(XEV(K).le.XEhighV.and.K.ge.NLV)) then
	    NTVV = NTVV + 1
	    NTVB = NTVE + 1
	    NTVE = K - 1
	    if(NTVB.ne.NTVE+1) IVSW = 1
	    if(K.le.NLVEND) K = K - 1
	    NTVD = NTVE-NTVB+1
	    if(NTTV.eq.1) NLVBEG = NTVD
	    if(NTVD.ne.0)
     &	    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) NLVEND = NTVE
	  end if
	end if

	if(IVSW.eq.1) then

C	print *, ' V switch ',IVSW,'  loop',NTVV

	    if (NLVEND .eq. 0) call Say('IPSD2020','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
C	    call MkLOSWeightsm(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV)
	    call MkLOSWeightsx(NLOSWV,dLOSV,NLOSV,NTVD,VFREQ(NTVB),VSIZE(NTVB),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,0.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(amax.ne.BadR4()) Nbadvel = Nbadvel + 1

	 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('IPSD2020','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.NLGEND+1) then	! If K is more than the source # +1 don't turn on switch
	 if(K.gt.0) then
C	print *, ' Before G switch loop K, NTGG, NLmax', K,NTGG,NLmax,XCEG(K),IYRSG(NTGG+2),DOYG(NTGG+2),XCintG(NTGG+2)
	  if(bGcon.and.(XCEG(K).gt.XCintG(NTGG+2).and.XEG(K).le.XEhighG).or.(XEG(K).le.XEhighG.and.K.ge.NLG)) then
	    NTGG = NTGG + 1
	    NTGB = NTGE + 1
	    NTGE = K - 1
	    if(NTGB.ne.NTGE+1) IGSW = 1
	    if(K.le.NLGEND) K = K - 1
	    NTGD = NTGE-NTGB+1
	    if(NTGG.eq.1) NLVBEG = NTGD
	    if(NTGD.ne.0)
     &	    write (*,'(A,I4,A,2I6,2F10.3)') 'G loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2)
	    if(NTGG.eq.NTG) NLGEND = NTGE
	  end if
	 end if
	end if

	if(IGSW.eq.1) then
	 if(NTGB.ne.0) then
	    if (NLGEND .eq. 0) call Say('IPSD2020','E','NoG','no 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
		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)

C		write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintG(NTGG)

	    end if
						! Calculate weights, WTSG, along all lines of sight

C	    call MkLOSWeightsm(NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG)
	    call MkLOSWeightsx(NLOSWG,dLOSG,NLOSG,NTGD,GFREQ(NTGB),GSIZE(NTGB),XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG)

C          if(GFREQ(NTGB).ne.ALSTGFREQ) then
C           write (*,'(A,I5,2F12.4)') 'NTGB, GFREQ(NTGB), GSIZE(NTGB)', NTGB, GFREQ(NTGB), GSIZE(NTGB)
C           ZEROVALUEWT = 0.0
C           write (*,'(F12.5)')      ZEROVALUEWT
C           do NNNWTS=1,40
C             write (*,'(F12.5)')     WTSG(NNNWTS,NTGB)/WTSG(20,NTGB)
C           end do          
C           ALSTGFREQ = GFREQ(NTGB)
C         end if


	 call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),PPGMAP(1,1,NTGG))
						! Convert G to density
	 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,0.0)
	     end if
	 end if

	 call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax)
	 print *, 'The DMAP',NTGG,' min and max is', amin, amax
	 if(amax.ne.BadR4()) Nbadden = Nbadden + 1

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(amaxD.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('IPSD2020','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

       if(Nbadvel .eq. 0 .and. .not. bVproxy) then	! If all else fails, fill the velocity arrays with a valid values
	 do N=1,nTV
           do I=1,Nlng
             do J=1,Nlat
               VMAP(I,J,N) = 400.0
	     end do
	   end do
	 end do
	 print *, ' '
	 print *, 'The velocity arrays all needed to be filled with valid values.'
	 print *, ' '
       end if

       if(Nbadden .eq. 0 .and. .not. bDproxy) then	! If all else fails, fill the density arrays with a valid values
	 do N=1,nTG
           do I=1,Nlng
             do J=1,Nlat
               DMAP(I,J,N) = DEN1AU*R1AU/RRG
	     end do
	   end do
	 end do
	 print *, ' '
	 print *, 'The density arrays all needed to be filled with valid values.'
	 print *, ' '
       end if

       if(bVFACE) then
	 do I=1,N_INV
	   WTSV_in(I) = 10.0
	 end do
       end if
	
       if(bFACE) then
	 do I=1,N_ING
	   WTSD_in(I) = 10.0
	 end do
        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
	    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,' DMAPs beginning with .not. bGcon are', amax
	end if

	NLGG = 0
	NLGGL = 0
	do L=1,NLG
	   if(L.lt.NLGBEG.or.L.gt.NLGEND) then
	     GOBS(L) = BadR4()
	     NBSG(L) = 0
	   end if 
         if(bForinput.and.L.eq.NLGEND.and.NLGG.ne.0) then         ! if this option is chosen and L is equal to the last data point then check.

           IDOYS = int(DOYS)                                      ! DOYS is the day of year of the run time

C          If the run time is greater than  0.79861 on the run time day day then the last time data came in was at 0.79861 on the run time day
           if(DOYS.gt.(IDOYS + 0.79861))  DIDOYS = dfloat(IDOYS) + 0.79861d0

C          If the run time is less than 0.79861d0 on the run time day then the last time data came in was was at 0.79861 on the day before the run time day
           if(DOYS.lt.(IDOYS + 0.79861))  DIDOYS = dfloat(IDOYS-1) + 0.79861d0
C           write(*,'(A,3I5,3F8.3,I3)')'Internal end loop D0', L, NLGEND, IDOYS, DOYS, DIDOYS, DOYSG8(NLGEND),NBSG(NLGEND)
           if(DOYSG8(NLGEND).le.DIDOYS) then                       ! Is the last data point less than the last run time of the day?
             go to 1120                                            ! If yes go on and let the normal process treat the last data point.          
           else
             do LL= 1,300                                          ! Step back one source at a time to remove sources prior to the run time.             
               LLL = NLGEND-LL+1
	         if(DOYSG8(LLL).gt.DIDOYS) then         
C                 write(*,'(A,3I5,2F8.3,I3)')'Internal end loop D1', L, LL, LLL, DIDOYS, DOYSG8(LLL), NBSG(LLL)
	           GOBS(LLL) = BadR4()
	           NBSG(LLL) = 0
                 if((NLGG - LL + 1).gt.0) then                     ! This works to NOT make a mistake if there are too few good images needed removing
                   NLGG  = NLGG - LL + 1
                   NLGGL = LLL - 1
                 else
                   NLGGL = 0
                 end if
C                 write(*,'(A,3I5,2F8.3,I3)')'Finished end loop D3', L, NLGG, NLGGL, DIDOYS, DOYSG8(LLL), NBSG(LLL)
               end if
             end do
             go to 1121
           end if
         end if
1120     continue
	   if(GOBS(L) .ne. BadR4() .and. GOBS(L) .ne. 0.0) then
	     NLGG = NLGG + 1
	     NLGGL = L
	   end if
1121     continue
	end do

	NLVV = 0
	NLVVL = 0
	do L=1,NLV
	   if(L.lt.NLVBEG.or.L.gt.NLVEND) then
	     VOBS(L) = BadR4()
	     NBSV(L) = 0
	   end if
         if(bForinput.and.L.eq.NLVEND.and.NLVV.ne.0) then          ! if this option is chosen and L is equal to the last data point then check.

           IDOYS = int(DOYS)                                       ! DOYS is the day of year of the run time

C          If the run time is greater than 0.79861 on the run time day day then the last time data came in was at 0.79861 on the run time day
           if(DOYS.gt.(IDOYS + 0.79861))  DIDOYS = dfloat(IDOYS) + 0.79861d0

C          If the run time is less than 0.79861d0 on the run time day then the last time data came in was was at 0.79861 on the day before the run time day
           if(DOYS.lt.(IDOYS + 0.79861))  DIDOYS = dfloat(IDOYS-1) + 0.79861d0
           if(DOYSV8(NLVEND).le.DIDOYS) then                       ! Is the last data point less than the last run time of the day?
             go to 1130                                            ! If yes go on and let the normal process treat the last data point.          
           else
             do LL= 1,300                                          ! Step back one source at a time to remove sources prior to the run time.             
               LLL = NLVEND-LL+1
	         if(DOYSV8(LLL).gt.DIDOYS) then
	           VOBS(LLL) = BadR4()
	           NBSV(LLL) = 0
                 if((NLVV - LL + 1).gt.0) then                     ! This hopefully works to NOT make a mistake if there are too few good images needed removing
                   NLVV  = NLVV - LL + 1
                   NLVVL = LLL - 1
                 else
                   NLVVL = 0
                 end if
               end if
             end do
             go to 1131
           end if
         end if
1130     continue
	   if(VOBS(L) .ne. BadR4()) then
	      NLVV  = NLVV + 1
	      NLVVL = L
	   end if
1131     continue
	end do

	if(NLGG.lt.NLGmin) then
	  print *, ' '
	  write(*,'(A,I3,A)') 'There are too few g-level (',NLGG,') remote sources for a good deconvolution.'
	  write(*,'(A)') 'All remote sources set bad.'
	  do L=1,NLG
	    NBSG(L) = 0
	  end do
	  NLGG = 0
        end if

	NLVV = 0
	NLVVL = 0
	do L=1,NLV
	   if(L.lt.NLVBEG.or.L.gt.NLVEND) then
	     VOBS(L) = BadR4()
	     NBSV(L) = 0
	   end if
	   if(VOBS(L) .ne. BadR4() .and. VOBS(L) .ne. 0.0) then
	      NLVV  = NLVV + 1
	      NLVVL = L
	   end if
	end do

	if(NLVV.lt.NLVmin) then
	  print *, ' '
	  write(*,'(A,I3,A)') 'There are too few velocity (',NLVV,') remote sources for a good deconvolution.'
	  write(*,'(A)') 'All remote sources set bad.'
	  do L=1,NLV
	    NBSV(L) = 0
	  end do
	  NLVV  = 0
        end if

	NLGGG  = 0
	NLGGGL = 0
	do L=NLG+1,NLG+N_ING
	  if(GOBS(L) .ne. BadR4() .and. GOBS(L) .ne. 0.0) then
	    NLGGG =  NLGGG + 1
	    NLGGGL = L
	  end if
	end do

	NLVVV  = 0
	NLVVVL = 0
	do L=NLV+1,NLV+N_INV
	   if(VOBS(L) .ne. BadR4() .and. VOBS(L) .ne. 0.0) then
	     NLVVV  = NLVVV + 1
	     NLVVVL = L
	   end if
	end do

	if(bGcon.and.NLGGL.ne.0) then
	    print *, ' '
	    write (*,'(1X,A)') 'The last IPS g-value was:'
	    write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') 
     &	         'I =',NLGGL,' Year =',IYRSG(NLGGL),' DOY =',DOYSG8(NLGGL),' IPS g-level =',GOBS(NLGGL)
	end if

	if(bFACE.and.NLGGGL.ne.0) then
	    print *, ' '
	    write (*,'(1X,A)') 'The last good in-situ density value was:'
	    write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') 
     &	         'I =',NLGGGL,' Year =',IYRSG(NLGGGL),' DOY =',DOYSG8(NLGGGL),' In-situ density =',GOBS(NLGGGL)
	end if

	if(bVcon.and.NLVVL.ne.0) then
	  print *, ' '
	  write (*,'(1X,A)') 'The last good IPS velocity value was:'
	  write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') 
     &	         'I =',NLVVL,' Year =',IYRSV(NLVVL),' DOY =',DOYSV8(NLVVL),' IPS velocity =',VOBS(NLVVL)
	end if

	if(bVFACE.and.NLVVVL.ne.0) then
	  print *, ' '
	  write (*,'(1X,A)') 'The last good in-situ velocity value was:'
	  write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') 
     &	         'I =',NLVVVL,' Year =',IYRSV(NLVVVL),' DOY =',DOYSV8(NLVVVL),' In-situ velocity =',VOBS(NLVVVL)
	end if

	do N=1,NTV
          do I=1,Nlng
            do J=1,Nlat
              VMAP(I,J,N) = 400.0
	    end do
	  end do
	end do

	print *, ' '
	call Say('ipsdt_20n_inp','I','Info',' # good IPS V observations         '
     &		//cStr(:Int2Str(NLVV,cStr)))

	call Say('ipsdt_20n_inp','I','Info',' # good IPS g-level observations   '
     &		//cStr(:Int2Str(NLGG,cStr)))
     
	call Say('ipsdt_20n_inp','I','Info',' # good in-situ V observations '
     &		//cStr(:Int2Str(NLVVV,cStr)))
     
	call Say('ipsdt_20n_inp','I','Info',' # good in-situ D observations '
     &		//cStr(:Int2Str(NLGGG,cStr)))
	print *, ' '
	
C      stop

	call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD)

	call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV)

	if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
	    call ArrR4Copy(nLngLatnTG,DMAP,DMAPV)
	else
	    call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &			nLng,nLat,NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD)
	end if

	if(bVproxy) call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV)

	if(bDproxy) then		! Do this for use after velocity iteration

	    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,VMAP,VMAPD,Dtmpt,VVT,VVV)
	    end if
	    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,0.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,0.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

        call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed,
     &	nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,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,0.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('IPSD2020','I','Info','Iteration '//cStr(:Int2Str(Nit,cStr)))

C	    write(*,'(A,I3,4F7.3)') 'nit, CONRV, CONVT, CONSTV CONSTMV',nit, CONRV, CONVT, CONSTV, CONSTMV

	    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
		  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,RRV)
		  end do

		  if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
		      call ArrR4Copy(nLngLatnTG,DMAPV,DMAP)
		  else
		      call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,
     &				    nLng,nLat,NTV,NTG,ConsTV,DMAPV,VMAPD,Dtmpt,VVT,VVV)

		      call ArrR4Copy(nLngLatnTG,VMAPD,DMAP)
		  end if

		  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

		call MkPostd_in(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV,
     &			XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8,
     &			XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM,
     &			N_INV,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in)

C		print *, 'Spot check after Mkpostd_V - 1,30', ITIM(1,1),ITIM(1,30),IYRSV(1),IYRSV(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30)

		  if(NitV.eq.NiterT.and.bWritelos) 
     & 			call writelosmapvg_in(0,NLOSVP1,NLV,N_INV,nTV,NBSV,NCOFF,ITIM,ITIM_in,XLON,XLON_in,XLAT,XLAT_in)

		  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 Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV,
     &			DDfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in)
		  call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV,
     &			DVfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,VFAC_in)

 		call MkVModeltd_in(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX,
     &			WTSV,DFAC,VFAC,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSV_in,DFAC_in,VFAC_in,PWV,
     &			NLV,N_INV,NLOSV,NLOSVP1,VMAP,VLIM,
     &			DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij,VWTij_in)

		do I=1,NLV+N_INV
C		    if(I.le.NLV) print *, 'Before  NLV',I,VMOD(I),VOBS(I),WTSV(20,I),NBSV(I),VWTij(20,I)
C		    if(I.gt.NLV) print *, 'Through NLV',I,VMOD(I),VOBS(I),WTSV_in(I-NLV),NBSV(I),VWTij_in(I-NLV)
 		    if(WTSV(1,I).eq.0.and.I.le.NLV) then
		      NBSV(I) = 0
 		      VMOD(I) = Speed
 		    end if
 		    if(WTSV_in(I-NLV).eq.0.and.I.gt.NLV) then
		      NBSV(I) = 0
 		      VMOD(I) = Speed
 		    end if
		end do    

C		if(nit.eq.18) stop

C		  if(VMOD(I).eq.BadR4().or.VOBS(I).eq.BadR4().and.I.gt.NLV) then
C		    NBSV(I) =  0
C		  else
C		    if(I.le.NLV) print *, 'Before  NLV',I,VMOD(I),VOBS(I),WTSV(20,I),NBSV(I),VWTij(20,I)
C		    if(VMOD(I).gt.2000..or.VMOD(I).lt.100.) print *, I,VMOD(I),VOBS(I)
C		    if(I.gt.NLV) print *, 'Through NLV',I,VMOD(I),VOBS(I),WTSV_in(I-NLV),NBSV(I),VWTij_in(I-NLV)
C		    NBSV(I) =  1
C		  end if
C		end do

C		  if(I.gt.367.and.I.le.NLV) then
C		if(WTSV(1,I).eq.0.and.I.le.NLV) then
C		    print *, 'Before the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I)
C		    NBSV(I) = 0
C		  end if
C		if(WTSV_in(I-NLV).eq.0.and.I.gt.NLV) then
C		    print *, 'Before the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I)
C		    NBSV(I) = 0
C		  end if
C		end do

		call FixModeltdn(1,NLV+N_INV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO)

C		do I=1,NLV+N_INV
C		  if(I.gt.367.and.I.le.NLV) then
C		    print *, 'After the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I)		    
C		  end if
C		end do

		if(NitV .eq. 0 .or. NitV .eq. NiterT) then

		  do I=1,NLV+N_INV,24
		    if(I.le.NLV) then		
C		      print *, 'Observed and model V:', I,' ',SRCV(I),VOBS(I),VMOD(I),FIX(I),NBSV(I)	
		    else
C		      print *, 'Observed and model V:', I,' ',SRCV(I),VOBS(I),VMOD(I),FIX(I),NBSV(I),XLON_in(I-NLV),XLAT_in(I-NLV),XCtim_in(I-NLV)
		    end if		
		  end do

		end if

		if(bWriteMS.and.NitV.eq.NiterT) then
		    call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV+N_INV,PWV,VRATIO,
     &			NLmaxV,NSRV,SRCV,DOYSV8,XCEV,VOBS,VMOD,XLON,XLproj,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

		if(NVREST.eq.1) then

		  call MkVMaptdn0n_in(LON,LAT,ITIM,RPX,LON_in,LAT_in,ITIM_in,RPX_in,XEV,VWTij,VWTij_in,
     &			FIX,NLV,N_INV,NLOSV,NLOSVP1,
     &			nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VSIGB,NVS,NBSV,VMAPSIG,
     &			ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,WTPV_in,
     &			Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp)

		end if
		if(NVREST.gt.1) then

		  call MkVMaptdn0n_inm(LON,LAT,ITIM,RPX,LON_in,LAT_in,ITIM_in,RPX_in,XEV,VWTij,VWTij_in,
     &			FIX,NLV,N_INV,NLOSV,NLOSVP1,
     &			nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VSIGB,NVS,NBSV,VMAPSIG,
     &			ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,WTPV_in,
     &			Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp)

		end if

		do N=1,NTV

C			call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	     		if(amax.eq.BadR4()) then
C		 		print *, 'VMAP is bad before map determination N =', N,' Iter =', Nit
C	    		end if 

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 (.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 MkShiftdn0n DMAP  min max N',amin,amax,N

		  AEmax = AEmaxD/(RRS*RRS)
		  AEmin = AEminD/(RRS*RRS)
		  call FillMapOLS(2,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,Speed,ConstL,AEmax,AEmin,AEangD,DMAP)

		  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,0.0)

C		  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C		  print *, 'Before MkShiftdn0n 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 (bGcon) then

C		if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
C		    call ArrR4Copy(nLngLatnTG,VMAP,VMAPD)
C		else
		    call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,
     &			nLng,nLat,NTV,NTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV)
C		end if

		do N=1,NTG

C	  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
C	  print *, 'Before FillWholeT VMAPD min max',amin,amax

		end do

		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 MkShiftdn0n 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,0.0)

C		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
C		  print *, 'Before Timesmooth and MkShiftdn0n 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))
C		  call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax)
C		  print *, 'Before MkShiftdn0n VMAPD min max N',amin,amax,N

		 end do

	      end if
	    end if

	    if(bGcon) then

        	call MkShiftdn0n(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG,Speed,
     &		nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL,
     &		NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

C	print *, VMAPD  (30,5,20),VMAPD  (31,5,20),VMAPD  (30,6,20),VMAPD  (31,6,20),VMAPD  (30,5,21),VMAPD  (31,5,21),
C     &		 DMAP   (30,5,20),DMAP   (31,5,20),DMAP   (30,6,20),DMAP   (31,6,20),DMAP   (30,5,21),DMAP   (31,5,21),
C     &		 XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11),
C     &		 DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11),
C     &		 DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11)

	    end if

	    if (bGcon) then				! G deconvolution
		NitG = NitG+1				! # G iterations
							! bVmv2=.TRUE., implies bVcon=.FALSE.
		if (bVmv2) then
		  RRV = RRD
		   call Deneq3(nLng,nLat,NTG,DMAP,RRD,DENEQ)
		  print *, 'Solar equatorial 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

		  if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
		      call ArrR4Copy(nLngLatnTG,VMAPD,VMAP)
		  else
		      call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &				    nLng,nLat,NTG,NTV,ConsTG,VMAPD,DMAPV,Vtmpt,DDT,DDD)

		      call ArrR4Copy(nLngLatnTV,DMAPV,VMAP)

		  end if

		  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
		call MkPostd_in(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG,
     &			XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8,
     &			XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM,
     &			N_ING,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in)

C		call MkPostd_in(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV,
C     &			XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8,
C     &			XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM,
C     &			N_INV,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in)

C		call MkPostd(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG,
C     &			XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8,
C     &			XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM)

C		print *, 'Spot check after Mkpostd_G - 1,30', ITIM(1,1),ITIM(1,30),IYRSG(1),IYRSG(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30)

		if(NitG.eq.NiterT.and.bWritelos) 
     &			call writelosmapvg_in(1,NLOSGP1,NLG,N_ING,nTG,NBSG,NCOFF,ITIM,ITIM_in,XLON,XLON_in,XLAT,XLAT_in)

C		if(NitG.eq.NiterT.and.bWritelos) call writelosmapvg(1,NLOSGP1,NLG,nTG,NBSG,NCOFF,ITIM,XLON,XLAT)

		call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG,
     &			DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC)
     
		call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG,
     &			DDfact,N_ING,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in)

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
 		call MkGModeltd_in(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX,
     &			WTSG,DFAC,PWG,DEN1AU,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSD_in,DFAC_in,
     &			NLG,N_ING,NLOSG,NLOSGP1,DMAP,
     &			nLng,nLat,nTG,RRS,GMOD2,GWTij,GWTij_in,WTPD,GMWTi)

C 		call MkGModeltdn(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX,
C     &			WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP,
C     &			nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi)

		do I=1,NLG+N_ING
		    if(I.le.NLG) then
		      if(GOBS(I).eq.badR4()) then
		        GOBSS(I) = badR4()
		        NBSG(I) = 0
		      else
		        GOBSS(I) = GOBS(I)

			if(bFACE) then
			  if(DGOBSS.le.1.0) then
		            if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)*DGOBSS
		            end if
		            if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)-(1.0-DGOBSS)
		            end if
			  else
		            if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)-(1.0-DGOBSS)
		            end if
		            if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)*DGOBSS
		            end if
			  end if			
			end if

		      end if
		    else
		      if(GOBS(I).eq.badR4()) then
		        GOBSS(I) = badR4()
		        NBSG(I) = 0
		      else

		        GOBSS(I) = (GOBS(I)/(DEN1AU*((1.0/RPX_in(I-NLG))**2)))**(2.0*PWG)

			if(bFACE) then
			  if(DGOBSS.le.1.0) then
		            if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)*DGOBSS
		            end if
		            if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)-(1.0-DGOBSS)
		            end if
			  else
		            if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)-(1.0-DGOBSS)
		            end if
		            if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then
		              GOBSS(I) = GOBSS(I)*DGOBSS
		            end if
			  end if
			end if

		      end if
C		      print *, GOBSS(I),GOBS(I),DEN1AU,RPX_in(I-NLG),PWG
		    end if
 		    if(WTSG(1,I).eq.0.and.I.le.NLG) then
		      NBSG(I) = 0
 		      GMOD2(I) = 1.0
 		    end if
 		    if(WTSD_in(I-NLG).eq.0.and.I.gt.NLG) then
		      NBSG(I) = 0
 		      GMOD2(I) = 1.0
 		    end if
 		end do

		call FixModeltdn(2,NLG+N_ING,GOBSS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO)

C		call FixModeltdn(2,NLG,GOBSS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO)

		if(NitG .eq. 0 .or. NitG .eq. NiterT) then

		  do I=1,NLG+N_ING,24
		    if(I.le.NLG) then	
C		      print *, 'Observed and model g:', I,' ',SRCGG(I),GOBSS(I),GMOD2(I),FIX(I),NBSG(I)
		    else
C		      print *, 'Observed and model g:', I,' ',SRCGG(I),GOBSS(I),GMOD2(I),FIX(I),NBSG(I),XLON_in(I-NLG),XLAT_in(I-NLG),XCtim_in(I-NLG)
		    end if		
		  end do

		end if

		if(bWriteMS.and.NitG.eq.NiterT) then
		    call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG+N_ING,PWG,GRATIO,
     &			NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GOBSS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG)
		end if

C		if(bWriteMS.and.NitG.eq.NiterT) then 
C		    call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO,
C     &			NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GOBSS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG)
C		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

C	print *, 'Into MkDmaptdn0n_in',nLng,nLat,nTG,NLOSG,NLG

		if(NGREST.eq.1) then

		  call MkDMaptdn0n_in(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in,
     &			PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1,
     &			nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG,
     &			ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in,
     &			Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp)

		end if
		if(NGREST.gt.1) then

		  call MkDMaptdn0n_inm(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in,
     &			PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1,
     &			nLng,nLat,nTG,DMAP,ANLIMITG,DLIM,GSSIG,GSIGB,NGS,NBSG,DMAPSIG,
     &			ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in,
     &			Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp)

		end if

C	print *, 'Out of MkDmaptdn0n_in'

C 		call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1,
C     &			nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG,
C     &			ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,
C     &			Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp)

		do N=1,NTG

C	  call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
C	  print *, 'After MkDmaptd0 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(Nit.ne.NiterT) then 	! Update velocity shifts and vmap except on last

C
C 	    If in-situ density is available define DEN1AU except on the last if bFACE is yes
C

	       DEN1  = 0.0
	       ADEN1 = 0.0
	       if(bFACE) then
	         do I=NLG+1,NLG+N_ING
 		    if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then
		      DEN1  = DEN1 + GOBS(I)
		      ADEN1 = ADEN1 + 1.0
 		    end if
 	         end do
 	         DEN1AUO = DEN1AU
 	         DEN1AU = DEN1/ADEN1
 	         if(abs(DEN1AU-DEN1AUO).gt.0.000001) print *, 'The old 1 AU density',DEN1AUO,' has been changed to',DEN1AU
 	       end if

	       if(bFACE.and.Nit.eq.NiterT) print *, 'The final 1 AU density has been set to',DEN1AU

	        if(bGACE) then
 	          DGOBSSO = DGOBSS

	          DGOBSS  = 0.0
	          ADGOBSS = 0.0
C	          do I=1,NLG
C 		    if(GOBS(I).le.GLIMP.and.GOBS(I).ge.GLIMM.and.NBSG(I).ne.0) then
C		      DGOBSS  = DGOBSS + GOBS(I)
C		      ADGOBSS = ADGOBSS + 1.0
C 		    end if
C 	          end do
	          do I=NLG+1,NLG+N_ING
 		    if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then
		      DGOBSS  = DGOBSS + (GOBS(I)/(DEN1AU*((1.0/RPX_in(I-NLG))**2)))**(2.0*PWG)
		      ADGOBSS = ADGOBSS + 1.0
		    end if
 	          end do
 	          DGOBSS = DGOBSS/ADGOBSS


 	          if(abs(DGOBSS-DGOBSSO).gt.0.000001) print *, 'The old G-level mean',DGOBSSO,' has been changed to',DGOBSS
 	        end if

	       if(bGACE.and.Nit.eq.NiterT) print *, 'The final G-level mean of 1.0 has been set to',DGOBSS

	       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)

C	if(niTV.eq.2) then
C	  print *, 'Line 2828 iteration 2 stop'
C	  stop
C	end if

		do N=1,NTV

C	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'After MkDmaptdn0  VMAP min max N',amin,amax,Ncd IPS

		  call FillMapOLS(2,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,Speed,ConstL,AEmaxV,AEminV,AEangV,VMAP)

   		  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,0.0)

		end do

		call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) 	! prior to 2016 BVJ
C		call Timesmooth(nLng,nLat,nTV,VMAP, 2,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

C		if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then
C		    call ArrR4Copy(nLngLatnTG,DMAP,DMAPV)
C		else
                call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,
     &				nLng,nLat,NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD)
C		end if

		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 MkDmaptdn0 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,0.0)

		end do

		call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) 	! prior to 2016 BVJ
C		call Timesmooth(nLng,nLat,nTV,DMAPV,2,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

        	call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed,
     &		nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL,
     &		NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt)

C	print *, VMAP   (30,5,20),VMAP   (31,5,20),VMAP   (30,6,20),VMAP   (31,6,20),VMAP   (30,5,21),VMAP   (31,5,21),
C     &		 DMAPV  (30,5,20),DMAPV  (31,5,20),DMAPV  (30,6,20),DMAPV  (31,6,20),DMAPV  (30,5,21),DMAPV  (31,5,21),
C     &		 XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11),
C     &		 DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11),
C     &		 DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11)

	       end if
	      end if
	    end if

	end do						! End of iteration loop

C *********************************************************************
C	The below is needed at least for header writes in the NV3 files
C
C	    call T3D_set (T3D__CLIPLNG		, 0, 90.0	)
C	    call T3D_iset(T3D__ITERATIONS	, 0, NiterT	)
C	    call T3D_iset(T3D__NCOFF		, 0, NCoff	)

C	    call T3D_set (T3D__RR		, 0, RRS	)
C	    call T3D_set (T3D__DRR		, 0, dRR	)

C	    call T3D_set (T3D__TT		, 0, TT		)
C	    call T3D_set (T3D__DTT		, 0, dTT	)

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

C	    call T3D_set (T3D__PWN_V		, 0, PWV	)
C	    call T3D_set (T3D__PWN_G		, 0, PWG	)
C	    call T3D_set (T3D__PWR_V		, 0, PWRV	)
C	    call T3D_set (T3D__PWR_G		, 0, PWRG	)

C	    call T3D_set (T3D__D1AU		, 0, DEN1AU	)
C	    call T3D_set (T3D__SMOOTH_V		, 0, CONRV	)
C	    call T3D_set (T3D__SMOOTH_D		, 0, CONRD	)
C	    call T3D_set (T3D__FILL_V		, 0, CONSTV	)
C	    call T3D_set (T3D__FILL_D		, 0, CONSTG	)
C	    call T3D_set (T3D__SMOOTH_TIME_V	, 0, CONVT	)
C	    call T3D_set (T3D__SMOOTH_TIME_D	, 0, CONDT	)

	    call T3D_iset(T3D__NL_V		, 0, NLV	)
	    call T3D_iset(T3D__NL_G		, 0, NLG	)
	    call T3D_iset(T3D__NLOS_V		, 0, NLOSV	)
	    call T3D_iset(T3D__NLOS_G		, 0, NLOSG	)
	    call T3D_set (T3D__DLOS_V		, 0, dLOSV	)
	    call T3D_set (T3D__DLOS_G		, 0, dLOSG	)

	    call T3D_set (T3D__BINX_V		, 0, ANLIMITV	)
	    call T3D_set (T3D__BINX_D		, 0, ANLIMITG	)

C *********************************************************************

	if(bSource) then				! Write out good sources used	

	  ididg = 0
	  ididv = 0
	  print *, ' '
	  write(*,'(A,L3,4I8,4L3)') 'Before call writeallgoodsources', bSource, NLG, NLGG, NLV, NLVV, bGcon, bVcon, bGgen, bVgen

	  if(bVgen.or.bGgen) then
	      write(*,'(A,I5,A,I5)') 'Write all sources (good or bad) from the general format. NLV =', NLV,' NLG =',NLG
	      itype=1						!goodsources.txt file output
 	      call writeallgoodsources1(itype,cggfiles,cvgfiles,NLMAXV,NLMAXG,SRCGV,SRCVG,NBSV,NBSG,NLV,NLG,N_INV,N_ING)
	  else
	      write(*,'(A,I5,A,I5)') 'Write all sources (good or bad) in the old nagoya fmt. NLV =', NLV,' NLG =',NLG
	      itype=0						!goodsources.txt file output
 	      call writeallgoodsources(itype,cggfiles,cvgfiles,NLMAXV,NLMAXG,SRCGV,SRCVG,NBSV,NBSG,NLV,NLG,N_INV,N_ING)
	  end if

C	  if(NLGG.eq.0.and.NLVV.ne.0.and.bVgen) then
C	    write(*,'(A,I5)') 'NLGG equals 0. writegoodsources. NLV =', NLV
C 	    call writegoodsources(NLMAXV,SRCVG,NBSV,NBSG,NLV,N_INV)
C	    ididg = 1
C	    ididv = 1
C	  end if

C	  if(NLGG.ne.0.and.NLVV.eq.0.and.bGgen) then
C	    write(*,'(A,I5)') 'NLVV equals 0. writegoodsources. NLG =',NLG
C	    call writegoodsources(NLMAXG,SRCGV,NBSV,NBSG,NLG,N_ING)
C	    ididg = 1
C	    ididv = 1
C	  end if

C	  if(NLG.eq.NLV) then
C	    if(ididg.eq.0.and.ididv.eq.0) print *, 'NLG equals NLV for writegoodsource. Both =',NLG
C	    istation = 1
C	    if(bVUCSD) istation = 0
C	    if(bGNago) istation = 1
C	    if(bGOoty) istation = 2
C	    call writegoodsourcegv(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NBSG,NLV,N_INV)

C	    if(bGgen.and.ididg.eq.0.and.ididv.eq.0) call writegoodsources(NLMAXV,SRCVG,NBSV,NBSG,NLV,N_INV)
	    
C	    if(.not.bGgen.and.ididg.eq.0.and.ididv.eq.0) call writegoodsourcegv(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NBSG,NLV,N_INV)
C	    ididv = 1
C	    ididg = 1
C	  else

C	    istation = 1
C	    if(bGcon.and.ididg.eq.0) then
C	      if(bVUCSD) istation = 0
C	      if(bGNago) istation = 1
C	      if(bGOoty) istation = 2
C	      print *, 'NLG does not equal NLV. Writegoodsourceg. NLG =',NLG,' NLV =',NLV
C	      call writegoodsourceg(istation,NLMAXG,cggfiles,SRCGV,NBSG,NLG,N_ING)
C	      ididg = 1
C	    end if
C	    if(bVcon.and.ididv.eq.0) then
C	      if(bVUCSD) istation = 0
C	      if(bGNago) istation = 1
C	      if(bGOoty) istation = 2
C	      print *, 'NLG does not equal NLV. Writegoodsourcev. NLG =',NLG,' NLV =',NLV
C	      call writegoodsourcev(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NLV,N_INV)
C	      ididv = 1
C	    end if
C	  end if

C	  if(NLGG.ne.0.and.NLVV.ne.0.and.bVgen.and.bGgen) then

C	    istation = 1
C	    if(ididv.eq.0) then
C	      write(*,'(A,I5,A,I5)') 'Make writegoodsourcev. NLG =',NLG, ' NLV =', NLV
C	      call writegoodsourcev(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NLV,N_INV)
C	      ididv = 1
C	    end if
C	    if(ididv.eq.0) then
C	      write(*,'(A,I5,A,I5)') 'Make writegoodsourceg. NLG =',NLG, ' NLV =', NLV
C	      call writegoodsourceg(istation,NLMAXG,cggfiles,SRCGV,NBSG,NLG,N_ING)
C	      ididg = 1
C	    end if
C	  end if

	end if

	J = 0						! List bad IPS sources
	NGOODSDVAL = 0
      NGOODVISEE = 0 
      NGOODVLOFA = 0 
      NGOODVBSA3 = 0
      NGOODVOOTY = 0
      NGOODVMEXA = 0

	do I=1,NLV+N_INV
          if (NBSV(I) .eq. 1 .and. I .le. NLV .and. bVGen) then
            cSRCVG = SRCVG(I)
C            write(*,'(A,A)') 'cSRCVG(43:44) = ', cSRCVG(43:44)            
            read(cSRCVG(43:44),'(A)') cSITE
C            write(*,'(A,A)') 'Site = ', cSITE
            if(cSITE.eq.'S ') NGOODVISEE = NGOODVISEE + 1      
            if(cSITE.eq.'L ') NGOODVLOFA = NGOODVLOFA + 1           
            if(cSITE.eq.'P ') NGOODVBSA3 = NGOODVBSA3 + 1      
            if(cSITE.eq.'O ') NGOODVOOTY = NGOODVOOTY + 1 
            if(cSITE.eq.'M ') NGOODVMEXA = NGOODVMEXA + 1 
          end if
	    if (NBSV(I) .eq. 0) then
		J = J+1
		if(VOBS(I).eq.BadR4()) VOBS(I) = -9999.9
		write (*,'(3A,I5,A,F7.1)') ' Bad V ', SRCV(I),' source removed, # =',
     &			I,', elongation =',XEV(I)
	    end if
	    if(I.gt.NLV.and.NBSV(I).eq.1) then
	      if(NGOODSDVAL.eq.0) FIRSTDOY = sngl(DOYSV8(I))
	      NGOODSDVAL = NGOODSDVAL + 1
	      FINALDOY = sngl(DOYSV8(I))
	    end if
	end do
	write (*,'(A,I8,A)') ' There were',J,' bad V sources removed.' 
      if(NLV.ne.0) then
            if(isys(3).eq.3) write (*,'(A,I7,A)') ' There were',NGOODVOOTY,' good Ooty IPS velocities used.'
            if(isys(4).eq.4) write (*,'(A,I7,A)') ' There were',NGOODVISEE,' good ISEE IPS velocities used.'
            if(isys(6).eq.6) write (*,'(A,I7,A)') ' There were',NGOODVMEXA,' good MEXART IPS velocities used.'
            if(isys(7).eq.7) write (*,'(A,I7,A)') ' There were',NGOODVBSA3,' good BSA3 IPS velocities used.'
            if(isys(8).eq.8) write (*,'(A,I7,A)') ' There were',NGOODVLOFA,' good LOFAR IPS velocities used.'
      end if
	if(N_INV.ne.0) then
	  write (*,'(A,I7,2(A,F6.1))') 
     &		' There were',NGOODSDVAL," good in-situ velocities. 1'st DOY =",FIRSTDOY ,' End DOY =',FINALDOY
	end if

	J = 0	
	NGOODSDVAL = 0
      NGOODGISEE = 0 
      NGOODGLOFA = 0 
      NGOODGBSA3 = 0
      NGOODGOOTY = 0
      NGOODGMEXA = 0
	do I=1,NLG+N_ING
          if (NBSG(I) .eq. 1 .and. I .le. NLG .and. bGGen) then
            cSRCGV = SRCGV(I)
C            write(*,'(A,A)') 'cSRCGV(43:44) = ', cSRCGV(43:44)            
            read(cSRCGV(43:44),'(A)') cSITE
C            write(*,'(A,A)') 'Site = ', cSITE
            if(cSITE.eq.'S ') NGOODGISEE = NGOODGISEE + 1      
            if(cSITE.eq.'L ') NGOODGLOFA = NGOODGLOFA + 1           
            if(cSITE.eq.'P ') NGOODGBSA3 = NGOODGBSA3 + 1      
            if(cSITE.eq.'O ') NGOODGOOTY = NGOODGOOTY + 1 
            if(cSITE.eq.'M ') NGOODGMEXA = NGOODGMEXA + 1 
          end if
	    if (NBSG(I).eq.0) then
		J = J+1
		if(GOBSS(I).eq.BadR4()) GOBSS(I) = -99.999
		write (*,'(3A,I5,A,F7.3)') ' Bad G ',SRCGG(I),' source removed, # =',
     &			I,', elongation =',XEG(I)
	    end if
	    if(I.gt.NLG.and.NBSG(I).eq.1) then
	      if(NGOODSDVAL.eq.0) FIRSTDOY = sngl(DOYSG8(I))
	      NGOODSDVAL = NGOODSDVAL + 1
	      FINALDOY = sngl(DOYSG8(I))
	    end if
	end do
	write (*,'(A,I8,A)') ' There were',J,' bad g-level sources removed.' 
      if(NLG.ne.0) then
            if(isys(3).eq.3) write (*,'(A,I7,A)') ' There were',NGOODGOOTY,' good Ooty IPS g-values used.'
            if(isys(4).eq.4) write (*,'(A,I7,A)') ' There were',NGOODGISEE,' good ISEE IPS g-values used.'
            if(isys(6).eq.6) write (*,'(A,I7,A)') ' There were',NGOODGMEXA,' good MEXART IPS g-values used.'
            if(isys(7).eq.7) write (*,'(A,I7,A)') ' There were',NGOODGBSA3,' good BSA3 IPS g-values used.'
            if(isys(8).eq.8) write (*,'(A,I7,A)') ' There were',NGOODGLOFA,' good LOFAR IPS g-values used.'
      end if
	if(N_ING.ne.0) then
	  write (*,'(A,I7,2(A,F6.1))') 
     &		' There were',NGOODSDVAL," good in-situ densities. 1'st DOY =",FIRSTDOY ,' End DOY =',FINALDOY
	end if

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),DMAPV(1,1,N))			    ! Place into a dummy array
	end do

C	Changed 9/21/01 BVJ
C	if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then			   ! Do this if smoother V maps needed
C	    call ArrR4Copy(nLngLatnTV,DMAPV,VMAP)				   ! In any case, replace dummy array
C	else				   					   
C	    call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat,
C     &			nTV,nTV,ConsTV,DMAPV,VMAP,Vtmpt,VVT,VVV)
C	end if
	    call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat,
     &			nTV,nTV,ConsTV,DMAPV,VMAP,Vtmpt,VVT,VVV)
	do N=1,NTV
	  call arrR4getminmax(nLngLat,VMAP(1,1,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
C	    XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region (Change 1/18 BVJ)
C	    XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/27.275) ! Mid region of interest on map

	    call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,Speed,ConstL,AEmaxV,AEminV,AEangV,VMAP)

C	    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),VMAPD(1,1,N))			    ! Place into a dummy array
	end do

C	Changed below 9/21/01 BVJ
C	if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then			   ! Do this if smoother D maps needed
C	    call ArrR4Copy(nLngLatnTG,VMAPD,DMAP)				   ! Replace dummy array
C	else
C	    call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat,
C     &			nTG,nTG,ConsTG,VMAPD,DMAP,Vtmpt,DDT,DDD)
C	end if
	    call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat,
     &			nTG,nTG,ConsTG,VMAPD,DMAP,Vtmpt,DDT,DDD)
	do N=1,NTG
	  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	    XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region (Changed 1/18/02)
C	    XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/27.275) ! Mid region of interest on Map

	    AEmax = AEmaxD/(RRS*RRS)
	    AEmin = AEminD/(RRS*RRS)
	    call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,Speed,ConstL,AEmax,AEmin,AEangD,DMAP)

	    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

	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,VMAP,VMAPD,Dtmpt,VVT,VVV)
	end if

C	call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,nLng,nLat,
C    &			NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD)

	print *, ' Before the first Extractdvd',XCbegG(1,1),XCtbegG,XCstrt

C*****  Take data out at different spacecraft (and earth) locations (bExtract tells where)

	NinterD = 3

C	call Extractdvd(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat,
C     &	      nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)

	call Extractdvd(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat,
     &	      nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)

C	call Extractdnn(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat,nMap,nTG,NinterD,FALLOFFD,
C     &	      VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)


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.

	nTF = nTG
	nMapn = nMap

	Nit = Nit + 1

	cPre = 'nv3d'

C The first write3D_infotd makes maps that have holes to be used to show Carrington maps

	call write3D_infotd3D_3(NOneOr1,bForecast,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,	
     &			nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT)

	do I=1,4*nTmaxG						! Will not work for Extractd3D if XCbegG's are different
	    XCbe(1,I) = XCbegG(1,1)
	    XCbe(2,I) = XCbegG(2,1)
	end do

C	if(bWrerr) then							! Changed BVJ 18/5/11
C These (mode 2) write3D_infotd calls make 3D arrays with time error maps to show reconstruction reliability

	 do N=1,nTG
	     do I=1,nLng
	     	 do J=1,nLat
		     if(ANMAPD(I,J,N).eq.0.0) ANMAPD(I,J,N)   = badR4()
		     VDEV(I,J,N) = ANMAPD(I,J,N)                            ! VDEV will go out of bounds if NTG not = to NTV BVJ 2020
		     if(AWTMAPD(I,J,N).eq.0.0) AWTMAPD(I,J,N) = badR4()
		     GDEV(I,J,N) = AWTMAPD(I,J,N)
		 end do
	     end do
	 end do
	 call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat,
     &			nTG,nTG,ConsTG,VDEV,ANMAPD,Vtmpt,DDT,DDD)
	 call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat,
     &			nTG,nTG,ConsTG,GDEV,AWTMAPD,Vtmpt,DDT,DDD)
	 PIRSQD = (CONDT/aNdayG)*3.14159*((CONRD/(20.0/float(NGRESS)))**2)	! normalize weights and lines of sight to spatial resolution
	 do N=1,nTG
	   call arrR4getminmax(nLngLat,ANMAPD(1,1,N),amin,amax)
	   if(amax.ne.BadR4()) then
	     call GridFill(2,nLng,nLat,ANMAPD(1,1,N),NSIDE,Zmin,Zmax)		! Remove small ANMAPD holes
	     call ArrR4TimesConstant(-nLngLat,ANMAPD(1,1,N),PIRSQD,ANMAPD(1,1,N))
	   end if	     
	   call arrR4getminmax(nLngLat,AWTMAPD(1,1,N),amin,amax)
	   if(amax.ne.BadR4()) then	   
	     call GridFill(2,nLng,nLat,AWTMAPD(1,1,N),NSIDE,Zmin,Zmax)		! Remove small AWTMAPD holes)	     	
	     call ArrR4TimesConstant(-nLngLat,AWTMAPD(1,1,N),PIRSQD,AWTMAPD(1,1,N))
	   end if
	 end do

	 do N=1,nTV
	     do I=1,nLng
	     	 do J=1,nLat
		     if(ANMAPV(I,J,N).eq.0.0) ANMAPV(I,J,N)   = badR4()
		     VDEV(I,J,N) = ANMAPV(I,J,N)
		     if(AWTMAPV(I,J,N).eq.0.0) AWTMAPV(I,J,N) = badR4()
		     GDEV(I,J,N) = AWTMAPV(I,J,N)		     
		 end do
	     end do
	 end do    	     	
	 call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat,     ! I think this is wrong (should be NTV,NTG) BVJ 2020
     &			nTV,nTV,ConsTV,VDEV,ANMAPV,Vtmpt,VVT,VVV)             ! But nTV, nTV OK doesn't change much (ANMAPV not used)
	 call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat,     ! Above should be ANDMAPV that is used in 3-D writes
     &			nTV,nTV,ConsTV,GDEV,AWTMAPV,Vtmpt,VVT,VVV)

	 PIRSQV = (CONVT/aNdayV)*3.14159*((CONRV/(20.0/float(NVRESS)))**2)	! normalize weights and lines of sight to spatial resolution

	 do N=1,nTV
	   call arrR4getminmax(nLngLat,ANMAPV(1,1,N),amin,amax)
	   if(amax.ne.BadR4()) then
	     call GridFill(2,nLng,nLat,ANMAPV(1,1,N),NSIDE,Zmin,Zmax)		! Remove small ANMAPV holes
	     call ArrR4TimesConstant(-nLngLat,ANMAPV(1,1,N),PIRSQV,ANMAPV(1,1,N))
	   end if
	   call arrR4getminmax(nLngLat,AWTMAPV(1,1,N),amin,amax)
	   if(amax.ne.BadR4()) then
	     call GridFill(2,nLng,nLat,AWTMAPV(1,1,N),NSIDE,Zmin,Zmax)		! Remove small AWTMAPV holes
	     call ArrR4TimesConstant(-nLngLat,AWTMAPV(1,1,N),PIRSQV,AWTMAPV(1,1,N))
	   end if
	 end do


	if(bWrerr) then								! Changed BVJ 18/5/11

	 NinterCC = 0
	 cPre = 'er1_'
	 print *, 'Before the error1 (L.O.S. crossings) write'

	 call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,	
     &			nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,ANMAPV,ANMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D)	

	 print *, 'Before the error2 (projected weights) write'

	 NinterCC = 0
	 cPre = 'er2_'
	 call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,	
     &			nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,AWTMAPV,AWTMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D)

	end if

	if(.not. bForecast) then

	 print *, ' Before Extractd3D',XCbegG(1,1),XCtbegG,XCstrt

	 call Extractd3D(bExtract,RRS,dRR,XCbe,XCtbegG,XCtendG,nLng,nLat,
     &	      nMap,nTG,NinterC+1,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar)

	end if

	call Timesmooth(nLng,nLat,nTV,VMAP,2,CONVT/aNdayV,0.75*CONVT/aNdayV,VZtmp) ! Fill in times one time away

	do N=1,NTV
	  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
	    SpeedV = 0.0
	    call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP)
	    call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))    ! Remove seam
	  end if

	  call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV/2.0,0,0.0,0.0)	   ! Change 5/9/08 to smooth polar holes

	end do

	call Timesmooth(nLng,nLat,nTG,DMAP,2,CONDT/aNdayG,0.75*CONDT/aNdayG,Dtmp)  ! Fill in times one time away

	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
	    AEmax = AEmaxD/(RRS*RRS)
	    AEmin = AEmin/(RRS*RRS)
	    SpeedD = 0.0
	    call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP)
	    call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
	  end if

	  call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD/2.0,0,0.0,0.0)	   ! Change 5/9/08 to smooth polar holes

	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,VMAP,VMAPD,Dtmpt,VVT,VVV)
	end if

C The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views

C The first write3D_infotd makes maps that have holes to be used to show Carrington maps 


	cPre = 'nv3f'

	call write3D_infotd3D_3(NOneOr2,bForecast,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,	
     &		nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar,
     &		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &		90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT)

C Must change to low value in current version

	 RRSCON = (RRS/R1AU)**FALLOFFD
	 do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	 end do

C End change

	if(bWrite3D_HR) then
	  print *, 'Before second write3D_infotd at High Resolution'
	  cPre = 'nv3h'
	  NintLng = NintLn
	  NintLat = NintLa
	  NintHt  = NintH
	  NinterD = 3							! In original map coordinates (so nMapHR = 1.5 AU).
	  if(nMapHR*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHR = nMap	! Make the setting of nMapHR idiot proof   
	  if(nMapb.gt.nMapHR) nMapb = nMapHR
	print *, RRRS, RRRSM, nMapb, nMapHR	 
C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU.

	  call write3D_infotd3DM_HR_3(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,
     &			nLng,nLat,nMapb,nMapHR,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			bDdener,bVdener,bVverer,bDverer,ERDLOS,ERVLOS,ANMAPD,ANDMAPV,
     &			90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR)
	 end if

C Must change back to regular value in current version

	 RRSCON = (R1AU/RRS)**FALLOFFD
	 do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	 end do

C End change

	if(bForecast.and..not.bmagex) then

C	 print *, 'DATA EXTRACTED FROM FINAL MATRIX FOR FORECAST'

	  print *, ' Before 2nd Extractdvd',XCbegG(1,1),XCtbegG,XCstrt

	  NinterD = 3

	  call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat,
     &	        nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)

	end if

*************************************************************************************************************************************

C Produce a set of nv3o files that are completely filled as well as we can make them.

C The third write3D_infotd makes maps that have no holes to be used to provide the best possible boundary for a 3D MHD model

        if(bWrite3Do.or.bWrite3D_HRb.or..not.bForecast) then
C	ixxxxxxx = 0
C	if(ixxxxxxx.eq.1) then 			! ixxxxxxx loop
     	 print *, 'Before third filling'
	 do N=1,NTV
	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	  print *, 'VMAP N, amin and amax = ', 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 arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
C	    print *, 'VMAP N, amin and amax = ', N, amin, amax
	    if(amax.ne.BadR4()) then
	      call GridFill(0,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax)		   ! Fill ALL velocity holes in any map with a hole
C	      call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))  ! Remove 90 deg discrepancy
C
C	      SpeedV = 0.0
C	      call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP)
      
C	      call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N))  ! Remove seam
	    end if
	  end if
	 end do
	 call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV)		     ! Fill all the velocity times
	 do N=1,NTV
	  call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax)
	  print *, 'VMAP N, amin and amax = ', N, amin, amax
	 end do
	 print *, 'Before forth filling'
	 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 arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax)
	    if(amax.ne.BadR4()) then
	      call GridFill(0,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax)		     ! Fill ALL density holes in any map with a hole

C	      AEmax = AEmaxD/(RRS*RRS)
C	      AEmin = AEmin/(RRS*RRS)	(needs fixing)
C	      AEmin = AEminD/(RRS*RRS)	(fixed 2/29/2008 BVJ)
C	      SpeedD = 0.0
C	      call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP)
	      
C	      call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove 90 deg discrepancy
C	      call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N))  ! Remove seam
	    end if
	  end if
	 end do
	 call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD)		     ! Fill all the density times
	 print *, 'After forth filling'
	 if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1).and.XCtbegG.eq.XCtbegV) then
	    print *, 'Copy array from VMAP to VMAPD directly'
	    call ArrR4Copy(nLngLatnTG,VMAP,VMAPD)
	 else
C	 print *, 'VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78)
	    call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat,
     &			NTV,NTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV)
C	 print *, 'VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78)
	 end if
C	end if 			!ixxxxxxx loop
	 print *, 'Before third write3D_infotd'

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

         if(bWrite3Do) then

C Must change to low value in current version

	   RRSCON = (RRS/R1AU)**FALLOFFD
	   do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	   end do

C End change

	   cPre = 'nv3o'
	   NinterD = 3

C	   if(ixxxxxxx.eq.1) then 			! ixxxxxxx loop	 
C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU.

	   call write3D_infotd3DM_3(0,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,	
     &			nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D)
C	   end if 					!ixxxxxxx loop

C Change the density back since the falloff at distance is produced in Extractdvd
	
	   RRSCON = (R1AU/RRS)**FALLOFFD
	   do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	   end do

C End change

	 end if

	 if(bWrite3D_HRb) then

C Must change to low value in current version

	   RRSCON = (RRS/R1AU)**FALLOFFD
	   do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	   end do

C End change

	   print *, ' '
	   print *, 'Before third write3D_infotd of High Resolution base'
	   cPre = 'nv3b'
	   NintLng = NintLn
	   NintLat = NintLa
	   NintHt  = NintH
	   NinterD = 3
	   nMapb = nMapHRb							! In original map coordinates (so nMapHRMb = 0.25 AU).
	   if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb	! Make the setting of nMapHR idiot proof   
	   if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb

	print *, RRRRS, RRRRSM, nMapb, nMapHRMb
	 
C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU.

	   call write3D_infotd3DM_HR_3(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,
     &			nLng,nLat,nMapb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC,
     &			PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANDMAPV,
     &			90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR)

C Change the density back since the falloff at distance is produced elsewhere
	
	   RRSCON = (R1AU/RRS)**FALLOFFD
	   do N=1,nTG
	     call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N))
	   end do

C End change

	 end if
	end if

C*************************************************************************************************************************************

	if ( bMagnetic ) then

C	    inexc = 0
C	    if(inexc.eq.0) go to 4567
	    dTT = (XCtendG-XCtbegG)/(sngl((nTF-1)*(NinterD+1)))		! Time increment
	    call T3D_set (T3D__TT		, 0, sngl(XCtbegG))	! Needed for header write in Write3D_bbtm_3.
	    call T3D_set (T3D__DTT		, 0, dTT 	)

	    call T3D_iset(T3D__NCOFF		, 0, NCoff	)   ! Needed to pass these values, and for the header write in Write3D_bbtm_3
	    call T3D_set (T3D__RR		, 0, RRS	)
	    call T3D_set (T3D__DRR		, 0, dRR	)
	    call T3D_iset(T3D__NLNG		, 0, nLng	)
	    call T3D_iset(T3D__NLAT		, 0, nLat	)

	    call T3D_iset(T3D__LAT		, 0, -90	)
	    call T3D_iset(T3D__LAT+1		, 0, +90	)

	    call T3D_iset(T3D__NRAD		, 0, nMap	)
	    call T3D_iset(T3D__NTIM		, 0, nTim	)

C 4567	continue

	    call T3D_iset(T3D__NCOFF		, 0, NCoff	)   ! Must pass these values to Bfield_Get_3 in Get_bbtm_3
	    call T3D_set (T3D__RR		, 0, RRS	)
	    call T3D_iset(T3D__NLNG		, 0, nLng	)
	    call T3D_iset(T3D__NLAT		, 0, nLat	)
	    call T3D_iset(T3D__NRAD		, 0, nMap	)
	    call T3D_iset(T3D__NTIM		, 0, nTim	)
	    call T3D_iset(T3D__LAT		, 0, -90	)
	    call T3D_iset(T3D__LAT+1		, 0, +90	)

	write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC

	    call get_bbtm_3(XCbeGG,nLng,nLat,nTG,RRS,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,
     &						XCINTG,BRR2DT,BRC2DT,BTC2DT,BNC2DT,iBflag)

	    if(bmagexx) then

	      print *, ' '
	      print *, 'Before write3D_bbtm_HR_3 after HR base'

	      NintLng = NintLn
	      NintLat = NintLa
	      NintHt  = NintH
	      Ninter = 3
	      ClipLng = 90.0
	      nMapbb = nMapHRb							! In original map coordinates (so nMapHRMb = 0.25 AU).
	      if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb	! Make the setting of nMapHR idiot proof   
	      if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb

	      IradialBRR = 0
	      IradialBRC = 0
	      IradialBTC = 0
	      IradialBNC = 0
	      do N=1,nTF+1
		call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax)
		if(amax.ne.BadR4()) then
	          call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD/2.0,2,0.0,0.0)	   ! Change 6/27/16 to smooth polar holes
		  IradialBRR = IradialBRR + 1						   ! Case 1
		end if
		call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax)
		if(amax.ne.BadR4()) then
	          call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD/2.0,2,0.0,0.0)	   ! Change 6/27/16 to smooth polar holes
		  IradialBRC = IradialBRC + 1						   ! Case 2
		end if
		call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax)
		if(amax.ne.BadR4()) then
	          call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD/2.0,2,0.0,0.0)	   ! Change 6/27/16 to smooth polar holes
		  IradialBTC = IradialBTC + 1						   ! Case 3
		end if
		call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax)
		if(amax.ne.BadR4()) then
	          call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD/2.0,2,0.0,0.0)	   ! Change 6/27/16 to smooth polar holes
		  IradialBNC = IradialBNC + 1						   ! Case 4
		end if
	      end do

	      if(bforecast) then
	        if(IradialBRR.ne.0.and.IradialBRR.ne.(nTF+1)) then		   ! Case 1

		    write(*,'(A,I4,A)') ' There are no radial magnetic fields for',(nTF+1)-IradialBRR,' BRR source surfaces.'
		    print *, 'These are filled in the following code.'
		    print *, ' '

	          call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRR2DT,DDD)
		    call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRR2DT,DDD)

	          Iradial = 0
	          do N=1,nTF+1
		    call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax)
		    if(amax.ne.BadR4().and.amin.eq.BADR4()) then
		      call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD,4,0.0,0.0)
		      write(*,'(A,I4,A)') 'BRR source surface',N,' was completely filled by gridsphere'
		    else
		      if(amin.eq.BadR4().and.amax.eq.BadR4()) then
		        write(*,'(A,I4,A)') 'BRR source surface',N,' was completely empty'
		        Iradial = Iradial + 1
		      end if
		    end if
		  end do
		  if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRR source surfaces could not be filled with gridsphere'
	        end if

	        if(IradialBRC.ne.0.and.IradialBRC.ne.(nTF+1)) then		   ! Case 2

		    write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBRC,' BRC source surfaces.'
		    print *, 'These are filled in the following code.'
		    print *, ' '

	          call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRC2DT,DDD)
		    call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRC2DT,DDD)

	          Iradial = 0
	          do N=1,nTF+1
		    call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax)
		    if(amax.ne.BadR4().and.amin.eq.BADR4()) then
		      call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD,4,0.0,0.0)
		      write(*,'(A,I4,A)') 'BRC source surface',N,' was completely filled by gridsphere'
		    else
		      if(amax.ne.BadR4()) then
		        write(*,'(A,I4,A)') 'BRC source surface',N,' was completely empty'
		        Iradial = Iradial + 1
		      end if
		    end if
		  end do
		  if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRC source surfaces were not filled with gridsphere'
	        end if

	        if(IradialBTC.ne.0.and.IradialBTC.ne.(nTF+1)) then		   ! Case 3

		  write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBTC,' BTC source surfaces.'
		  print *, 'These are filled in the following code.'
		  print *, ' '

	          call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BTC2DT,DDD)
		  call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BTC2DT,DDD)

	          Iradial = 0
	          do N=1,nTF+1
		    call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax)
		    if(amax.ne.BadR4().and.amin.eq.BADR4()) then
		      call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD,4,0.0,0.0)
		      write(*,'(A,I4,A)') 'BTC source surface',N,' was completely filled by gridsphere'
		    else
		      if(amax.ne.BadR4()) then
		        write(*,'(A,I4,A)') 'BTC source surface',N,' was completely empty'
		        Iradial = Iradial + 1
		      end if
		    end if
		  end do
		  if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BTC source surfaces were not filled with gridsphere'
	        end if

	        if(IradialBNC.ne.0.and.IradialBNC.ne.(nTF+1)) then		   ! Case 4

		  write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBNC,' BNC source surfaces.'
		  print *, 'These are filled in the following code.'
		  print *, ' '

	          call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BNC2DT,DDD)
		  call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BNC2DT,DDD)

	          Iradial = 0
	          do N=1,nTF+1
		    call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax)
		    if(amax.ne.BadR4().and.amin.eq.BADR4()) then
		      call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD,4,0.0,0.0)
		      write(*,'(A,I4,A)') 'BNC source surface',N,' was completely filled by gridsphere'
		    else
		      if(amax.ne.BadR4()) then
		        write(*,'(A,I4,A)') 'BNC source surface',N,' was completely empty'
		        Iradial = Iradial + 1
		      end if
		    end if
		  end do
		  if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BNC source surfaces were not filled with gridsphere'
	        end if
	      end if

	      write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC

C	BFACTORRR = 1.0
C	BFACTORRC = 1.0
C	BFACTORTC = 1.0
C	BFACTORNC = 0.0

C	      do N=1,nTF+1
C		call ArrR4Constant(nLngLat,BFACTORRR,BRR2DT(1,1,N))
C		call ArrR4Constant(nLngLat,BFACTORRC,BRC2DT(1,1,N))
C		call ArrR4Constant(nLngLat,BFACTORTC,BTC2DT(1,1,N))
C		call ArrR4Constant(nLngLat,BFACTORNC,BNC2DT(1,1,N))
C	      end do

C	write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC

	      do N=1,nTF
	        if(bFillmag) then
		  call ArrR4Constant(nLngLat,VMAPD(1,1,N),VMAPF(1,1,N))
		else
		  call ArrR4Constant(nLngLat,VMAP(1,1,N),VMAPF(1,1,N))
	        end if
	      end do

	      call write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG,iYrBG,
     &			RRS,dRR,nLng,nLat,nMapbb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC,
     &	      		PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT,
     &			ClipLng,VMAPF,XCshift,DVfact,
     &			BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR)

	    end if

	    if(bmagex) then

	      write(*,'(A,F10.6,F12.7,F10.6)') ' Before Final Extractdvdm_3',XCbegG(1,1),XCtbegG,XCstrt

C	      ifieldw2 = 1
	      ifieldw2 = 0
	      if(ifieldw2.ne.0) then
		print *, ' '
		write (*,'(A,100(I4,E12.5))') 'BRR2DT', (L,BRR2DT(10,4,L),L=1,ntG)
		print *, ' '
		print *, ' '
		write (*,'(A,100(I4,E12.5))') 'BRC2DT', (L,BRC2DT(10,4,L),L=1,ntG)
		print *, ' '
		print *, ' '
		write (*,'(A,100(I4,E12.5))') 'BTC2DT', (L,BTC2DT(10,4,L),L=1,ntG)
		print *, ' '
		print *, ' '
		write (*,'(A,100(I4,E12.5))') 'BNC2DT', (L,BNC2DT(10,4,L),L=1,ntG)
	      end if

	      print *, ' '

	      NinterD = 3

	      call Extractdvdm_3(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat,
     &	        nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift,
     &		DVfact,DDfact,TT3D,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar)

	     end if

	end if			     

	if(.not. bMagnetic) then

	  print *, ' Before Final Extractdvd',XCbegG(1,1),XCtbegG,XCstrt

	  NinterD = 3
	  call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat,
     &	        nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar)
	end if

	call Say('IPSD2020','S','Stop','program successfully completed')

	stop
	end