C+ C NAME: C smeiipstd_p C C The latest 9/15/00 work - tried a copydtodv here and this filled C the density matrix a lot better when it was placed just before write3dinfo - C There was no apparent change in extractd from the ea_file but maybe something C was different with ulysses. At this stage, this program does not give exactly C the same thing as the current version on the vax, but it is close. C Found an error 10/15/00 (A gridsphere2d was placed outside a do loop after C the first MKshiftd. This version should now have no seam in each temporal value. C C Modifications to the current program until BVJ 03/01/02 C C To compile the program with the g77 compler type: C C ./mk_smeiipsdt C C To run the program type: C C $myfor/SMEIIPSTD/smeiipstd nagoya=nagoys,$dat/nagoya/yearly smei=smei,$dat/smei C C Program defaults should allow the November 24 1977 data files to be run to give density and C magnetic field. C C PURPOSE: C Construct deconvolved synoptic maps for the IPS observations using C a "time-dependent" technique that projects all observations to a C single lower reference surface. Note - uneven time cadences are C allowed for velocity and density data sets. C C CATEGORY: C Data processing C CALLING SEQUENCE: C run smeiipstd_p C INPUTS: C Data file $CAM:G*. C OUTPUTS: C Rectangular grid of data points C FUNCTIONS/SUBROUTINES: C Local2UT,MAP_TZERO,READ_GO, C SetGipsy C INCLUDE: C include 'sun.h' C COMMON BLOCKS: C SIDE EFFECTS: C RESTRICTIONS: C PROCEDURE: C > The modified Julian day MJD used for the IPS data is related to the C regular Julian day (as used in subroutine Julian) by C MJD = JD-2400000.5 C (i.e. the modified Julian day starts at midnight). C > Once the logical $CAM is set (during the first run of smeiipstd), the C only way to force smeiipstd 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 Dec 2003 Bernie Jackson (UCSD) (runs with SMEI data) C- program smeiipstd_p parameter (NLmaxG = 300000, ! Max # G data points (13000) (400,000) C parameter (NLmaxG = 81000, ! Max # G data points (13000) (400,000) & NLmaxV = 2000, ! Max # V data points (2000) & NGREST = 2, ! G temporal resolution factor better than 1 day (Helios=1) (2 for SMEI works) & NVREST = 2, ! V temporal resolution factor better than 1 day (on cass183) & NGRESS = 3, ! Spatial G resolution factor better than 20 degrees(Helios=2) (3 for SMEI works) & NVRESS = 3, ! Spatial V resolution factor better than 20 degrees (with temp=2 on cass183) & NF = 1, ! LOS resolution factor (Should be 1 if Spatial resolution is >= 2) & nTmaxG = 50*NGREST, ! Max # G times (50) (100) (150) & nTmaxV = 50*NVREST) ! Max # V times (50) (100) (150) parameter (NLOSG = 20*NF*NGRESS, ! # resolution elements along each G line of sight (40) (80) (125) & NLOSV = 20*NF*NVRESS, ! # resolution elements along each V line of sight (40) (80) (125) & NLOSGP1 = NLOSG+1, & NLOSVP1 = NLOSV+1, & dLOSG = 0.1/(NGRESS*NF), ! LOS G Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & dLOSV = 0.1/(NVRESS*NF), ! LOS V Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & xInc = 0.5) ! Carrington map "extra" increment parameter (nLng = 54*NGRESS+1, ! # Longitudes, resolution is 360/(nLng-1) (55) (109) (217) (325) & iLng = 18*NGRESS+1, ! # Longitudes to use in extracted map (19) (37) (73) (109) & nLat = 9*NGRESS+1, ! # Latitudes , resolution is 180/(nLat-1) (9)(10) (19) (37) (55) & nMap = 31, ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) & nLngLat = nLng*nLat) parameter (dRR = 0.1) ! Resolution (AU) in radial direction (0.1) integer iXPG (NLmaxG), ! Array indices from daily G maps & iYPG (NLmaxG), & iMJDG (NLmaxG), ! Original MJD of daily G maps & IYRFG (NLmaxG), & IRECG (NLmaxG), & IYRSG (NLmaxG), ! Year of source observation & IYRSGG (NLmaxV), ! Year of source observation & NBSG (NLmaxG), & LON (NLOSGP1,NLmaxG), & LAT (NLOSGP1,NLmaxG), & ITIM (NLOSGP1,NLmaxG) parameter (NLOC1 = 4000) real RARA(NLOC1), & DECDEC(NLOC1), & ZSMEIZSMEI(NLOC1) integer ISSISS(NLOC1) real DISTG (NLmaxG), ! G distance from Sun to Earth & DISTGG (NLmaxV), ! G distance from Sun to Earth & XLSG (NLmaxG), ! Ecliptic longitude of Sun & XLSGG (NLmaxV), ! Ecliptic longitude of Sun & XLLG (NLmaxG), ! Longitude diff. (Earth-Sun) & XLLGG (NLmaxV), ! Longitude diff. (Earth-Sun) & XDLG (NLmaxG), ! Declination of Source & XDLGG (NLmaxV), ! Declination of Source & XCEG (NLmaxG), ! Carrington variable of Earth & XCEGG (NLmaxV), ! Carrington variable of Earth & XEG (NLmaxG), ! Source Elongation & XEGG (NLmaxV), ! Source Elongation & XCG (NLmaxG), ! Carrington variable & XCGG (NLmaxV), ! Carrington variable & YLG (NLmaxG), ! Heliographic latitude (degree) & YLGG (NLmaxV), ! Heliographic latitude (degree) & GOBS (NLmaxG), ! Solar wind observed G or brightness values - final & GGOBS (NLmaxG), ! Solar wind observed G or brigntness values - initial and modified & BBASE (NLmaxG), ! Solar wind observed base brightness values & ZZ (NLmaxG), ! Thomson scattering density values for MapGRID & ZZZ (NLmaxG), ! Thomson scattering density values for MapGRID c & VOBS (NLmaxV), ! Solar wind observed velocities & VOBS (NLmaxG), ! Solar wind observed velocities & GMOD2 (NLmaxG), ! Model G values & FIX (NLmaxG), ! Fix to model to make OK & GSSIG (NLmaxG), ! G-level source sigma & WTSG (NLOSG ,NLmaxG), & GWTij (NLOSG ,NLmaxG), & XLON (NLOSGP1,NLmaxG), & XLONP (NLOSGP1,NLmaxG), & XLAT (NLOSGP1,NLmaxG), & RPX (NLOSGP1,NLmaxG), & XLproj (NLOSGP1,NLMAXG), & VFAC (NLOSGP1,NLmaxG), ! New array (3/99) & DFAC (NLOSGP1,NLmaxG) character SRCV (NLmaxV)*9, ! Velocity source name c & SRCVS*9, ! TJD, never referenced & SRCVsav(NLmaxV)*9, & SRCG (NLmaxG)*12, & SRCGG (NLmaxV)*12, & SRCGS*9,SRCGsav(NLmaxG)*9 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 & VVOBS (NLmaxG), & VMOD (NLmaxV), ! Model velocities & VSSIG (NLmaxV), ! Velocity source sigma & VWTij (NLOSV,NLmaxV), & WTSV (NLOSV,NLmaxV) real XCinTV (nTmaxV), & XCinTG (nTmaxG), & DOYV (nTmaxV), & DOYG (nTmaxG), & XCbeV (2,nTmaxV), & XCbeGG (2,nTmaxG) integer IYRV (nTmaxV), & IYRG (nTmaxG) c & NMAPD (nLng,nLat,nTmaxG), TJD, never ref. ! Scratch integer arrays c & NMAPV (nLng,nLat,nTmaxV) TJD, never ref. real VDEV (nLng,nLat,nTmaxV), & VPNT (nLng,nLat,nTmaxV), & PPGMAP (nLng,nLat,nTmaxG), ! Point-P G map (with holes) & VMAPHOLE(nLng,nLat,nTmaxV), & DM (nLng,nLat), & VM (nLng,nLat), & DMAPHOLE(nLng,nLat,nTmaxG), & DMAP (nLng,nLat,nTmaxG), ! Density map (no holes) & DBASE (nLng,nLat,nTmaxG), ! Base map to use in Thomson scattering & VMAP (nLng,nLat,nTmaxV), & VMAPD (nLng,nLat,nTmaxG), & DMAPV (nLng,nLat,nTmaxV), & GDEV (nLng,nLat,nTmaxG), & GPNT (nLng,nLat,nTmaxG), & DDfact (nLng,nLat,nMap,nTmaxG), & DVfact (nLng,nLat,nMap,nTmaxG), ! New array (3/99) & XCshift (nLng,nLat,nMap,nTmaxG,3), & V3D (nLng,nLat,nMap), & D3D (nLng,nLat,nMap), & Vtmp (nLng,nLat,nTmaxG,3), ! Scratch real arrays & VtmpV (nLng,nLat,nTmaxV,3), & VZtmp (nLng,nLat,nTmaxV), & Dtmp (nLng,nLat,nTmaxG), & XLTtmp (nLng,nTmaxG,5), & XLtmp (nLng,nTmaxG), & XLtmpt (nLng,nTmaxG), & Dtmpt (nLng,nTmaxG), & Vtmpt (nLng,nTmaxV), c & Vtp (nLng,nLat), TJD never ref c & Dtp (nLng,nLat), TJD never ref & ANMAPD (nLng,nLat,nTmaxG), & DMAD (nLng,nLat,nTmaxG), & FIXMD (nLng,nLat,nTmaxG), & WTSMD (nLng,nLat,nTmaxG), & DMEANFD (nLng,nLat,nTmaxG), & DFIXM2D (nLng,nLat,nTmaxG), & WTPD (NLOSG,NLmaxG), & GMWTi (NLmaxG), & ANMAPV (nLng,nLat,nTmaxV), & VMAV (nLng,nLat,nTmaxV), & FIXMV (nLng,nLat,nTmaxV), & WTSMV (nLng,nLat,nTmaxV), & VMEANFV (nLng,nLat,nTmaxV), & VFIXM2V (nLng,nLat,nTmaxV), & WTPV (NLOSV,NLmaxV), & VVV (nLng), & DDD (nLng), & VVT (nTmaxV), & DDT (nTmaxG), & VRATION (200), & GRATION (200), & Scale (2) byte NSIDE(9*nLng*nLat) parameter (nCar = 500) parameter (nCarHOS1 = 500) parameter (nCarHOS2 = 500) parameter (nCarSMEI = 500) real*8 JDCar(nCar), & JD, JEpoch, & dLngSun,dLatSun,dDisSun , & JDCarHOS1(nCarHOS1), & JDCarHOS2(nCarHOS2), & JDCarSMEI(nCarSMEI), & JDCX, JDCX2, JDCP, & MJD_to_JD /2400000.5d0/, ! JD = MJD+MJD_to_JD c & MJDend, ! AUTO mode: last MJD (TJD, never ref) & MJDinc /.757639d0/, ! AUTO mode: increment & DLONG8, & DLAT8, & DDIST8 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. & DOYSsav8, & dorbit real*8 XCtpr (NLOSGP1,NLMAXG), & XCtim (NLOSGP1,NLMAXG), & XCintDV (nTmaxV), & XCintDG (nTmaxG) real*8 MJDref, & MJDcntr, & XCtbegV, & XCtbegG, & XCtendV, & XCtendG C------- C String written to output files. First three chars are also used in file name construction character cStrPPV(1)*80 /'PPV: Point P velocity map with holes'/, & cStrPPG(1)*80 /'PPG: Point P G-level map with holes'/, & cStrPPD(1)*80 /'PPD: 1 AU Point P density map with holes'/, & cStrDV (1)*80 /'DVV: Deconvolved velocity map'/, & cStrDD (1)*80 /'DDD: Deconvolved density map'/, & cStrDVs(1)*80 /'DVS: Deconvolved, smoothed velocity map'/, & cStrDDs(1)*80 /'DDS: Deconvolved, smoothed density map'/, & cStrVname*14 /'vts1941_00.txt'/, & cStrGname*14 /'gts1941_00.txt'/, & cSym*17 /'LIB__SIGNAL_QUIET'/ logical bGcon, bVcon, bTcon, & bVNago, bVOoty, bVUCSD, & bGcamb, bHelio, bGOoty, bGNago, bh1, bh2, bSMEI, & bVmv2 /.FALSE./, ! Use mv^2 = constant & bGmv2 /.FALSE./, ! Use mv^2 = constant & bDaydV, bDaydG, & bOnly1 /.FALSE./, ! If .TRUE., use only one shift matrix (on density matrix) & bVproxy /.FALSE./, ! Velocity proxy map used for deconvolution & bDproxy /.FALSE./, ! Density proxy map used for deconvolution & bMagnetic /.FALSE./, ! Magnetic map available & bGWRITE /.FALSE./, ! Write out a G data file & bVWRITE /.FALSE./, ! Write out a V data file & bMirror /.FALSE./, & bjMirror /.FALSE./, & bWrite /.FALSE./, & bWrite3D /.FALSE./, & bWrite3Do /.FALSE./, & bWriteMS /.TRUE./, & bWritelos /.TRUE./, & bStopw /.FALSE./, & bForecast, & bExtract(12) /2*.FALSE.,2*.TRUE.,8*.FALSE./, & bEWRITE, ! Write out an Earth file & bMWRITE, ! Write out a Mercury file & bMAWRITE, ! Write out a Mars file & bUWRITE, ! Write out a Ulysses file & bH1WRITE, ! Write out an Helios 1 file & bH2WRITE ! Write out an Helios 2 file equivalence (bExtract( 3), bEWRITE), & (bExtract( 1), bMWRITE), & (bExtract( 4), bMAWRITE), & (bExtract(10), bUWRITE) , & (bExtract(11), bH1WRITE), & (bExtract(12), bH2WRITE) integer iEorWG /0/, & iEorWV /0/, c & MJD, c & MJDfrst /48180/, ! MJD for earliest data file c & MJDlast /49622/, ! MJD for latest data file & MJDfrst / 44100/, & MJDlast / 51800/, & LimitoV /-1/, & LimitoG /-1/, & NLOSWV /1/, & NLOSWG /1/, & NTV /44/, & NTG /44/ integer NiterT /18/, & iOffUT / 8/, & NAV / 0/, & NC / 1/, ! # rotations averaged & nDUMP / 0/, ! Threshold for #pnts/bin & nFILL / 3/, ! Min. # neighbours for fill & NmidHRV /0/, & NmidHRG /0/ real YLbeg /-90./, & YLend / 90./, & XCadj(3) /3*0./, & XCadjB(3) /3*0./, & XCmargin /.5/, & XElowG /30./, & XElowV /11.5/, & XEhighG /80./, & XEhighV /180./, & XRlimG /90./, & XRlimV /90./, & GPOWER /0./, & R1AU /1.0/, C & DEN1AU /6.8/ , & DEN1AU /5.0/, & VEL1AU /425.0/, & ANLIMITV /1.5/, & ANLIMITG /1.5/, & aNdayV /1.0/, & aNdayG /1.0/ real VMAPSIG /1.0/, & DMAPSIG /1.0/, & GSIG /1.0/, & GSIGB /3.0/, & VSIG /1.0/, & VSIGB /3.0/, & FCONSV /1.0/, & FCONSG /1.0/, & CONRV /15./, & CONRD /15./, & CONVT /0.75/, & CONDT /0.75/, & ConsTV /1.0/, & ConsTG /1.0/, & VLIM /1200./, ! Velocity upper limit for sources & VLIL /100./, ! Velocity lower limit for sources & GLIMP /4.0/, ! G-level upper limit for sources & GLIMF /1.0/, ! G-level input multiplier & GLIMM /0.2/, ! G-level lower limit for sources & HELIOmax /10000./, ! Maximum Helios brightness allowed & HELIOmin /-10000./, ! Minimum Helios brightness allowed & SMEImax /900./, ! Maximum SMEI brightness allowed & SMEImin /-900./ ! Minimum SMEI brightness allowed real RRS / -15. /, ! Reference distance for deconvolution & VUL / 1500.0 /, ! Velocity upper limit for maps & VLL / 100.0 /, ! Velocity lower limit for maps & DUL / 500.0 /, ! Density upper limit for maps & DLL / 0.0 / ! Density lower limit for maps real DAV /1./, & Speed /1000./, & PWG /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.3/ , ! Radial power of density & PWRV /0.3/, ! Radial power of velocity & FALLOFF /2.00/ ! Density falloff with solar distance (>2.00 Implies ! an acceleration of the solar wind with height.) real 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 /90./ ! Angular distance from Earth to begin D limit character cCam*4, & cMon*3, c & cStrID*11, TJD, never ref & cVar(3)*100, & cArg*100, & cStr*80, & cPre*4 character cWildNagoya*80, & cWildOoty*80, & cWildUCSD*80, & cWildHelios1*80, & cWildHelios2*580, & cWildSMEI*80 C Below specific to Helios or SMEI character cHOS(10)*80 character cSMEI(10)*80 integer iSectBeg /-15/, & iSectEnd / 16/ integer iDayHOS / 0/, ! Time filter on Helios data & nDevHOS /10/ ! Glitch removal on Helios data (doesn't matter if iDayHOS is 0) integer iDaySMEI / 0/, ! Time filter on Helios data & nDevSMEI /10/ real dSigHOS /0./ ! Offset of Helios data from mean (doesn't matter if iDayHOS is 0) real dSigSMEI /0./ character STRING*80 external iReadNagoyan, iProcessNagoyan external iReadOotyn, iProcessOotyn external ireadUCSDn, iprocessUCSDn external EARTH , HELIOS_1, HELIOS_2 external EARTH8 include 'dirspec.h' include 'sun.h' include 't3d_array.h' integer Flt2Str integer ibtt(2),iett(2) c**TJD**added for magnetic time dependent version real Btmp(nLng,nLat,3), !scratch & XCbegMAT(nTmaxG), !3-Car time frame & XCendMAT(nTmaxG), & XCbegROI(nTmaxG), !ROI time frame & XCendROI(nTmaxG), & TrueVel(nLng,nLat,nMap,nTmaxG), !real velocity array & BR3D(nLng,nLat,nMap,nTmaxG), !scratch & BR2DT(nLng,nLat,nTmaxG), !scratch & Bt3D(nLng,nLat,nMap,nTmaxG) C:\smei\exe\win\smeiipstd.exe ucsd=SANIPS,C:\dat\ucsd\ helios1=A77B4_318_(4).R08,C:\dat\Helios\ helios2=B77B4_318_(4).R08,C:\dat\Helios\ call T3D_set_mode(T3D__MODE_TIME, 1) iVar = 3 call ForeignArg(' ',iVar,cVar,cArg) NLOC = NLOC1 ! Set the maximum number of SMEI lines of sight in an orbit to the maximum possible C------- C Calculate start times for Carrington rotations between 1985 and now. C Rotation NCoff+I starts at time JDCar(I). aNdayV = aNdayV/NVREST aNdayG = aNdayG/NGREST ConsTMV = ConsTV*nTV ConsTMD = ConsTG*nTG ConstL = nLng/3 I = iDeleteSymbol(cSym,2) iYr = 1974 Doy = 330. call MAP_TZERO(EARTH,iYr,Doy,.01,nCar,JDCar) print *, nCar, JDCar(1), JDCar(500) NCoff = N_CARRINGTON(iYr,Doy) print *, 'NCoff =',NCoff NCoffHOS1 = 0 NCoffHOS2 = 0 NCoffSMEI = 0 call ArrR4Constant(NLmaxV,5000.,XCEV) ! Initialize XCEV array call ArrR4Constant(NLmaxG,5000.,XCEG) ! Initialize XCEG array call ArrR4Constant(NTmaxV,5000.,XCintV) ! Initialize XCintV array call ArrR4Constant(NTmaxG,5000.,XCintG) ! Initialize XCintG array C------- C Set range of Helios Carrington variables to work with: [XClowHOS,XChighHOS] XClowHOS1 = 1 C XClowHOS1 = 16 C XClowHOS1 = 28 C XClowHOS1 = 30 XChighHOS1 = 129 XClowHOS2 = 1 C XClowHOS2 = 16 C XClowHOS2 = 29 C XClowHOS2 = 30 XChighHOS2 = 65 XClowSMEI = 1 XChighSMEI = 300 C XCbegHOS1 = 38.5 C XCbegHOS2 = 38.5 C XCbegHOS1 = 30.0 ! Works 11/8/03 C XCbegHOS2 = 30.0 XCbegHOS1 = 38.5 XCbegHOS2 = 38.5 C------ C iOffUT is an offset time (in hours) between local time and UT. This is used only if the program C is operated in forecast mode. Note that the local time is assumed to be the same as the C computer time. The default offset is 8 hours (the difference between PST and GST). call AskYN('Forecast mode?$no',bForecast) ! Regular or forecast mode if (bForecast) call AskI4('Offset between local time and UT (hours)?',iOffUT) call AskI4('Number of iterations?',NiterT) ! # iterations call AskR4('Deconvolution distance (AU; 0=solar surface)?',RRS) if (RRS .lt. 0) RRS = -RRS*Sun__RAu ! Neg. RRS are in units of solar radii RRS = max(RRS,sngl(Sun__RAu)) ! Keep RRS > 1 solar radius ! Deconvolve velocity, G-level or both? call AskWhat('Deconvolve: velocity, density, both$1$3',I) C------- C A proxy map for density and/or velocity can be used to replace the point-P C map determined from the observations. C The map should be in files dcrot.xxx and vcrot.xxx. They must be ascii files of dimension C nLng, nLat. C C If bVmv2=TRUE then bVcon=FALSE C If bGmv2=TRUE then bGcon=FALSE bVcon = I .eq. 1 .or. I .eq. 3 bGcon = I .eq. 2 .or. I .eq. 3 if (bVcon) then ! V deconvolution: select IPS station call AskWhat('Use velocity data from: Nagoya, Ooty, UCSD$1$1',I) bVNago = I .eq. 1 bVOoty = I .eq. 2 bVUCSD = I .eq. 3 bVmv2 = .FALSE. if(bVNago) NLOSWV = 1 if(bVOoty) NLOSWV = 2 if(bVUCSD) NLOSWV = 3 else ! No V deconvolution bVNago = .FALSE. bVOoty = .FALSE. bVUCSD = .FALSE. call AskYN('Do you want to make a velocity map using mv^2 = con$no',bVmv2) end if if (bGcon) then ! G or density deconvolution: select IPS station or Thomson scattering source call AskWhat('Use density or G-level data from: SMEI, Helios, Camb, Ooty, Nagoya$1$1',I) bHelio = I .eq. 2 if(bHelio) bh1 = .TRUE. if(bHelio) bh2 = .TRUE. bSMEI = I .eq. 1 bTcon = .FALSE. if(bHelio) call ASKYN('HELIOS 1',bh1) if(bHelio) call ASKYN('HELIOS 2',bh2) if(bHelio) GLIMP = HELIOmax if(bHelio) GLIMM = HELIOmin if (bh1) call MAP_TZERO(HELIOS_1,iYr,Doy,.01,nCarHOS1,JDCarHOS1) if (bh2) call MAP_TZERO(HELIOS_2,iYr,Doy,.01,nCarHOS2,JDCarHOS2) if (bSMEI) call MAP_TZERO(EARTH,iYr,Doy,.01,nCarSMEI,JDCarSMEI) if(bSMEI) GLIMP = SMEImax if(bSMEI) GLIMM = SMEImin if(bHelio) bTcon = .TRUE. if(bSMEI) bTcon = .TRUE. bGCamb = I .eq. 3 bGOoty = I .eq. 4 bGNago = I .eq. 5 bGmv2 = .FALSE. cCam = cEnvi//'CAM' if(bh1) NLOSWG = 5 ! B light used in Helios 1 analysis if(bh2) NLOSWG = 5 ! B light used in Helios 2 analysis if(bSMEI) NLOSWG = 6 ! V light Thomson scattering weights for SMEI if(bGCamb) NLOSWG = 2 if(bGOoty) NLOSWG = 1 if(bGNago) NLOSWG = 1 else ! No g deconvolution bHelio = .FALSE. bGCamb = .FALSE. bGOoty = .FALSE. bGNago = .FALSE. bSMEI = .FALSE. call AskYN('Do you want to make a density map using mv^2 = con',bGmv2) end if if (bGCamb) call SetGipsy(cCam,iOffUT,iEdt,MJDfrst,MJDlast) if (bVNago .or. bGNago ) call ForeignFile(iVar,cVar,'nagoya' ,cWildNagoya) if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty) if (bVUCSD) call ForeignFile(iVar,cVar,'ucsd' ,cWildUCSD) if (bh1) call ForeignFile(iVar,cVar,'helios1',cWildHelios1) if (bh2) call ForeignFile(iVar,cVar,'helios2',cWildHelios2) if (bSMEI) call ForeignFile(iVar,cVar,'smei' ,cWildSMEI) C if (bSMEI) cWildSMEI='/home/bjackson/dat/smei/jan_2005/points_200??????????????????????' C if (bSMEI) cWildSMEI='/mnt/storage/smei/data1/base/points_2005_????????????????????' C if (bSMEI) cWildSMEI='/mnt/storage2/smei/data/base/points_2003_????????????????????' print *, cWildSMEI if(bVcon) then bDaydV = .TRUE. call AskYN('Deconvolve V from one day to the next$yes',bDaydV) if(bDaydV) then if(bVUCSD) aNdayV = 12. call AskR4('Number of days to combine for velocity?',aNdayV) if(bVNago) NmidHRV = -10 call AskI4('Number of hours from V UT midnight?',NmidHRV) end if call AskI4('Number of velocity time values?',NTV) ! # NTV time itervals if(NTV.gt.NTmaxV-2) then NTV = NTmaxV-2 call AskI4('Number of velocity time values?',NTV) ! Ask again end if end if if(bGcon) then bDaydG = .TRUE. call AskYN('Deconvolve density or G from one day to the next$yes',bDaydG) if(bDaydG) then call AskR4('Number of days to combine for density or G-level?',aNdayG) if(bGNago) NmidHRG = -10 call AskI4('Number of hours from density or G UT midnight?',NmidHRG) end if call AskI4('Number of density or G-level time values?',NTG) ! # NTG time itervals if(NTG.gt.NTmaxG-2) then NTG = NTmaxG-2 call AskI4('Number of density or G-level time values?',NTG) ! Ask again end if end if C I = iFilePath('/home/bjackson/dat/nagoya/',1,'yearly','nagoya.'//cWildChar(:4),cWildNagoya) C I = iFilePath(cEnvi//'DAT',1,'Ooty' ,'Ooty.' //cWildChar(:4),cWildOoty ) C------- C Select the rotation for Helios 1. C if(bHelio) then if(bh1) then XCbegHOS1 = max(XClowHOS1,NCoffHOS1+XCbegHOS1) write (STRING,'(2(A,F7.2),A)') '$',XClowHOS1,'$',XChighHOS1,'$0' call ASKR4('Helios 1 rotation (0 = STOP)'//STRING,XCbegHOS1) if (XCbegHOS1 .le. 0) STOP 'Stop in IPSH12Dtd' XCbegHOS1 = min(max(XCbegHOS1,XClowHOS1),XChighHOS1) XCbegHOS1 = XCbegHOS1-NCoffHOS1 C------ C For the moment use only a single Helios rotation c call ASKI4('# Rotations to be averaged',NCHOS) c NCHOS = min(NCHOS,1+int(XChighHIGH-(NCoffHOS1+XCbegHOS1)) NCHOS = 1 C C------- Info only: start time of Carrington rotation C I = int(XCbegHOS1) JDCX = JDCarHOS1(I) JDCP = JDCarHOS1(I+NCHOS) NCoffP = NCoffHOS1 + I C call JULIAN(1,IYR,DOY,JDCX,JEPOCH) IDOY = DOY call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY) JDCX2 = JDCarHOS1(I+1) call JULIAN(1,IYR2,DOY2,JDCX2,JEPOCH) IDOY2 = DOY2 call DATE_DOY(1,IYR2,CMON2,IMON2,IDAY2,IDOY2) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/,30X,A,I2,3A,I4,A,I3,A)') 'Helios 1 rotation ',NCoffP, & ' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )', & ' and ends on date ',IDAY2,'-',CMON2,'-',IYR2,' ( DOY ',IDOY2,' )' iYrHOS1 = IYR end if C------- C Select the rotation for Helios 2. C if(bh2) then XCbegHOS2 = max(XClowHOS2,NCoffHOS2+XCbegHOS2) write (STRING,'(2(A,F7.2),A)') '$',XClowHOS2,'$',XChighHOS2,'$0' call ASKR4('Helios 2 rotation (0 = STOP)'//STRING,XCbegHOS2) if (XCbegHOS2 .le. 0) stop 'Stop in IPSH12Dtd' XCbegHOS2 = min(max(XCbegHOS2,XClowHOS2),XChighHOS2) XCbegHOS2 = XCbegHOS2-NCoffHOS2 C------ C For the moment use only a single Helios rotation c call ASKI4('# Rotations to be averaged',NCHOS) c NCHOS = min(NCHOS,1+int(XChighHIGH-(NCoffHOS2+XCbegHOS2)) NCHOS = 1 C C------- Info only: start time of Carrington rotation C I = int(XCbegHOS2) JDCX = JDCarHOS2(I) JDCP = JDCarHOS2(I+NCHOS) NCoffP = NCoffHOS2 + I call JULIAN(1,IYR,DOY,JDCX,JEPOCH) IDOY = DOY call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY) JDCX2 = JDCarHOS2(I+1) call JULIAN(1,IYR2,DOY2,JDCX2,JEPOCH) IDOY2 = DOY2 call DATE_DOY(1,IYR2,CMON2,IMON2,IDAY2,IDOY2) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/,30X,A,I2,3A,I4,A,I3,A)') 'Helios 2 rotation ',NCoffP, & ' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )', & ' and ends on date ',IDAY2,'-',CMON2,'-',IYR2,' ( DOY ',IDOY2,' )' iYrHOS2 = IYR end if C------- C NEarth is the Earth-based Carrington rotation closest to the Helios rotation C XCbeg. JD = 0.5d0*(JDCP+JDCX) ! Halfway Helios data range call JULIAN(1,IYR,DOY,JD,JEPOCH) NEarth = N_CARRINGTON(IYR,DOY) ! Nearest Carrington rotation IDOY = DOY call DATE_DOY(1,IYR,CMON,IMON,IDAY,IDOY) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A)') 'Nearest IPS rotation ',NEarth, & ' starts on date ',IDAY,'-',CMON,'-',IYR,' ( DOY ',IDOY,' )' XCbeg = NEarth-NCoff-NC/2 ! Center IPS range on Helios range write (*,'(/,10X,2(A,F6.1))') 'Carrington range for IPS data ',NCoff+XCbeg,'-',NCoff+XCbeg+NC C C Set up file parameters for different data sets output C IROT = NEarth + 0.00001 RDUM = IROT IHGT = RRS*10.0 + 0.00001 NCOFH1 = int(XCbegHOS1) NCOFH2 = int(XCbegHOS2) NCOFI = int(XCbeg) PLUS = 0.0 C C------- Get modified Carrington variable for center of plot C Get heliographic coordinates for center of plot (Earth position) call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) XCEarth = XCbeg+.5 ! Center of plot I = XCEarth ! Rotation containing plot center MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I)) ! Julian day for plot center call JULIAN(1,IYR,DOY,MJDcntr,JEPOCH) ! IYR,DOY for plot center MJDcntr = MJDcntr-MJD_to_JD ! Modified Julian day for plot center call SunNewcomb(0,IYR,DOY,DLONG8,DLAT8,DDIST8) EARTH_LNG = mod(sngl(DLONG8)+180.,360.) EARTH_LAT = sngl(DLAT8) call ECLIPTIC_HELIOGRAPHIC(0,IYR,DOY,EARTH_LNG,EARTH_LAT) ! Heliographic coord. for plot center XCtest1 = XCbeg-XCmargin ! Search limits XCtest2 = XCbeg+NC+XCmargin I = XCtest2 MJDref = JDCar(I)-MJD_to_JD+(XCtest2-I)*(JDCar(I+1)-JDCar(I)) ! Start MJD for data search else ! Non-Helios mode C Calculate range of IPS Carrington variables to work with: [XClo,XChi] C XClo is based on the first data file available (MJDfrst). C XChi is based on the current time. if(bVUCSD.and..not.bGcon) then Doy = Doy + 30. else call Julian(11,iYr,Doy,dble(MJDfrst),JEpoch) end if print *, 'before JD', iYr,Doy,XClo,xchi XClo = NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XClo = XClo-.5 call Local2UT(iOffUT,iYr,Doy) XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo) print *, 'before JD', XClo,xchi end if RRV = RRS RRG = RRS C------- C See if SMEI data exists C C------- C Deal with IPS velocity C if (bVcon) then if(NiterT.gt.8) LimitoV = NiterT/2 call AskI4('Iteration at which to limit V-source deviations, -1=none',LimitoV) C------- C Fix the reference distance used for the V map formed. C Get E or W or EW; Set the elongation angles call AskWhat('All velocity data, East data, West data ?$0',iEorWV) if (iEorWV .eq. 2) iEorWV = -1 call AskR4('Minimum elongation for V-data$0$180$',XElowV) call AskR4('Maximum elongation for V-data$'//cStr(:Flt2Str(XElowV,1,cStr))//'$180$',XEhighV) call AskR4('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 limit for V deconvolution?',VLIM) end if if (bGcon) then if(NiterT.gt.8) LimitoG = NiterT/2 call AskI4('Iteration at which to limit D or G-source deviations, -1=none',LimitoG) C Get E or W or EW; Set the elongation angles call AskWhat('All D or G data, East data, West data ?$0',iEorWG) if (iEorWG .eq. 2) iEorWG = -1 if(bHelio) XElowG = 15. if(bSMEI) XElowG = 15. C if(bSMEI) XElowG = 70. if(bGNago) XElowG = 5. call AskR4('Minimum elongation for D or G-data$0$180$',XElowG) if(bHelio) XEhighG = 180. if(bSMEI) XEhighG = 180. if(bGNago) XEhighG = 180. call AskR4('Maximum elongation for D or 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) if(bHelio) XRlimG = 180. if(bSMEI) XRlimG = 180. call AskR4('Maximum RA relative to Sun for D or G-data$0$180$',XRlimG) RRG = RRS call AskR4('D or G-map reference distance in AU (0=solar surface)?',RRG) if (RRG .lt. 0) RRG = -RRG*Sun__RAu RRG = max(RRG,sngl(Sun__RAu)) call AskYN('Do you want to use a proxy map for density?',bDproxy) call AskR4('Maximum limit for D or G deconvolution?',GLIMP) call AskR4('Minimum limit for D or G deconvolution?',GLIMM) call AskR4('Factor by which to modify brightness?',GLIMF) end if CONRD = CONRD/NGRESS call AskR4('Density spatial filter (degrees)?',CONRD) CONDT = CONDT*aNdayG call AskR4('Density temporal filter (days)?',CONDT) call AskR4('Density time hole filter?',ConsTG) if (bGcon) then ANLIMITG = ANLIMITG call AskR4('# of D or G-points per bin that constitute a deconvolution',ANLIMITG) end if if(bVNago.and.bGNago) CONRV = CONRD call AskR4('Velocity spatial filter (degrees)?',CONRV) if(bVNago.and.bGNago) CONVT = CONDT call AskR4('Velocity temporal filter (days)?',CONVT) call AskR4('Velocity time hole filter?',ConsTV) if (bVcon) then ANLIMITV = ANLIMITV*aNdayV*CONRV/20.0/NVRESS call AskR4('# of V-points per bin that constitute a deconvolution',ANLIMITV) end if if(bGCamb) PWG = 0.15 if(.not.bTcon) call AskR4(' G-level power to fit density?',PWG) if(bTcon) call AskR4('1 AU average density to fit base?',DEN1AU) if(bVNago.and.bGNago) PWV = PWG call AskR4('G-level power to fit velocity?',PWV) if(.not.bTcon) call AskR4('Radial G-level falloff to fit G density?',PWRG) if(bTcon) call AskR4('Radial Density falloff?',FALLOFF) if(bVNago.and.bGNago) PWRV = PWRG call AskR4('Radial G-level falloff to fit V density?',PWRV) call AskR4('Solar wind speed for G-level traceback?',Speed) RRD = RRG C------- Prompts for time if (bForecast) then ! Forecast mode call Local2UT(iOffUT,iYr,Doy) ! Doy = universal time call Julian(10,iYr,Doy,JD,JEpoch) ! Current time in MJD write (*,'(/,10X,A,I6,A,I6,A,F11.3)') 'First MJD =',MJDfrst, & ' Last MJD =',MJDlast,' Current MJD =',JD MJDref = JD ! MJD at center of plot (Earth position) write (cStr,'(A,I5,A,F7.1,A)') '$',MJDfrst,'$',JD,'$0$' call AskR8('Modified Julian Day (0 = STOP)'//cStr,MJDref) if (MJDref .le. 0) call Say('smeiipstd','I','StopPlay','StopPlay') C------- Make sure iYr,Doy match MJDref C Get heliographic coordinates for center of plot (Earth position) call Julian(11,iYr,Doy,MJDref,JEpoch) call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) C------- Get modified Carrington variable for center of plot C [Note that EarthLng=(1-frac(XCEarth))*360] XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XCbeg = XCEarth-.5 call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) XCtest1 = XCbeg-XCmargin ! Search limits XCtest2 = XCbeg+NC+XCmargin XCend= XCbeg + NC XCstrt = XCbeg XCbeg= XCbeg - xInc XCend= XCend + xInc MJDcntr = MJDref else ! Normal mode if(bVUCSD.and..not.bHelio) XCbeg = 1661.6 C print *, 'Here XCbeg, NCoff =',XCbeg, NCoff C if(bVNago.and.bSMEI) XCbeg = 2008.8 - NCoff C if(bSMEI) XCbeg = 2025.2 - NCoff ! January2005 period if(bSMEI) XCbeg = 2015.8 - NCoff ! Apr_May_2004 period if(bHelio) then XClo = NCoff+XCbeg-1 XChi = NCoff+XCbeg+1 end if if(bHelio) write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$' if(bHelio) XCbeg = 1661.6 if(bGCamb) XCbeg = 1884. if(bVNago.and.bSMEI) XCbeg = XCbeg + NCoff if(.not.bVNago.and.bSMEI) XCbeg = XCbeg + NCoff call AskR4('Carrington rotation (0 = STOP)'//cStr,XCbeg) if (XCbeg .le. 0) call Say('smeiipstd','I','StopPlay','StopPlay') C NORBITST = 10400 C call AskI4('Beginning orbit number',NORBITST) C MX = 600 C call AskI4('Total number of orbits',MX) NintF = 1 if(bForecast) then nTF = nTG - 1 NintF = 0 end if call AskI4('# Intermediate rotations interpolated',NintF) call AskI4('# Rotations averaged',NC) if(bHelio) NC = min(NC,1+int(XChi-XCbeg) ) ! Stay within available data range XCbeg = XCbeg-NCoff C------- Info only: start time of Carrington rotation call Julian(1,iYr,Doy,JDCar(int(XCbeg)),JEpoch) iDoy = Doy call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+int(XCbeg), & ' starts on date ',iDay,'-',cMon,'-',iYr,' ( Doy ',iDoy,' )' C------- Get modified Carrington variable for center of plot C Get heliographic coordinates for center of plot (Earth position) if(bVNago.and.bSMEI) call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) if(.not.bVNago.and.bSMEI) call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) 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 print *, 'XCtest1, XCtest2, NC, XCmargin', XCtest1, XCtest2, 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 if(bSMEI) then XCORBEG = 1957.9786 ! SMEI beginning orbit Carrington rotation, doy = 0.0, year 2000, doy 1, 0 UT ibtt(1) = ((XCbeg + NCoff) - XCORBEG)*27.275 ! Approx. days since the beginning orbit C ibtt(1) = (2015.2259 - XCORBEG)*27.275 ! about orbit 6600 (actually 6591) print *, XCbeg + NCoff,ibtt(1) ibtt(2) = 0 iett(1) = ((XCbeg + NCoff + 1.2 + xInc) - XCORBEG)*27.275 ! Approx. days since the ending orbit iett(2) = 0 call smei_orbit2(ibtt,NORBITST,dorbit) C print *, dorbit call smei_orbit2(iett,iorbite,dorbit) MX = iorbite - NORBITST + 1 print *, ' Beginning SMEI orbit # =', NORBITST,' Number of orbits =',MX end if C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCbeg call AskYN('Do you want to output model-source files$no',bWriteMS) call AskYN('Do you want 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 a Three D file$yes',bWrite3D) ! 3D file write call AskYN('Output completely filled Three D files$no',bWrite3Do) ! Completely filled 3D file write call AskYN('Output an Earth parameter file$yes',bEWRITE) ! Earth extraction call AskYN('Output a Mercury parameter file$no',bMWRITE) ! Mercury extraction call AskYN('Output a Mars parameter file$no',bMAWRITE) ! Mars extraction call AskYN('Output a Ulysses parameter file$no',bUWRITE) ! Ulysses extraction call AskYN('Output a Helios 1 parameter file$no',bH1WRITE) ! Helios 1 extraction call AskYN('Output a Helios 2 parameter file$no',bH2WRITE) ! Helios 2 extraction C Do you want to include a magnetic map at the reference surface height? C If yes, then read the source surface map (at 2.5Rs) at C the resolution of the reference arrays. call AskYN('Do you want to make a magnetic map?$no',bMagnetic) C------- C The range of the Carrington variable plotted is [XCbeg,XCend]. The target range C [XCtest1,XCtest2] includes an extra safety margin to make sure that all relevant C data are read from file. if (bVNago) call ReadVIPSn8(iReadNagoyan,iProcessNagoyan,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,SRCV,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(iReadOotyn, iProcessOotyn,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,SRCV,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,XEVsav,XCVsav,YLVsav,VVsav,GGsav) if (bGOoty) call ReadVIPSn8(iReadOotyn, iProcessOotyn,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, c & RRV,XElowG,XEhighG,iEorWG,XRlimG, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,SRCGG,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, c & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTV,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(iReadNagoyan,iProcessNagoyan,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,SRCGG,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) C C Set up the next two limits for the deconvolution just in case they are used C C print *, 'Before read_smei8' if(bSMEI) then NS = 0 call READ_SMEI8(cWildSMEI,nCar,JDCar,XCtest1,nint(XCtest2-XCtest1),iYrSMEI,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & iSMEI,cSMEI,iDaySMEI,nDevSMEI,dSigSMEI, & RARA,DECDEC,ZSMEIZSMEI,ISSISS, & NORBITST,MX,NLOC,NLmaxG,NS,NLT1, & SRCG,DISTG,XCG,YLG,GOBS,XEG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8) NLG = NLT1 C print *, 'Out of SMEIREAD' do I=1,NLG if(GOBS(I).ne.BADR4()) then C call UBVConst(LT2,U,APM) C xTHOM = ThomsonlosS10Far(0.,-1.,sngl(DISTG(I)/Sun__RAu/ABS(XEG(I))),ABS(XEG(I)),APM,P) C ZZ(I) = GOBS(I)*xTemp/xTHOM C if(I.ge.2000.and.I.le.2100) print *, I,ZZ(I),XLSG(I),XCEG(I),XLLG(I),DOYSG8(I),GOBS(I) C if(I.ge.NLT1+1800.and.I.le.NLT1+1900) print *, I,ZZ(I),XLSG(I),XCEG(I),XLLG(I),DOYSG8(I),GOBS(I) ZZ(I) = GOBS(I)*DEN1AU/DISTG(I) else print *, 'Bad GOBS(I) xTemp, xTHOM', I, GOBS(I), xTemp, xTHOM ZZ(I) = BadR4() end if end do end if C print *, 'After read_smei8', NLG C call ArrR4Zero(NLmaxV,VSSIG) call ArrI4Constant(NLmaxV,1 ,NBSV ) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLmaxG,1 ,NBSG ) do I=1,NLmaxV C if (I.lt.100) print *, I, SRCV(I) if (VOBS(I) .gt. VLIM) NBSV(I) = 0 if (VOBS(I) .lt. VLIL) NBSV(I) = 0 end do ICntBadr4 = 0 do I=1,NLmaxG if (GOBS(I) .eq. BadR4()) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if GOBS(I) = GOBS(I)*GLIMF if (GOBS(I) .gt. GLIMP) NBSG(I) = 0 if (GOBS(I) .lt. GLIMM) NBSG(I) = 0 end do if(ICntBadr4.gt.0) print *, ' The number of bad sources are = ', ICntBadr4 C======= do I=1,NLG if(bGNago.or.bGOoty) then XCG(I) = XCGG(I) YLG(I) = YLGG(I) DISTG(I) = DISTGG(I) GOBS(I) = GGOBS(I) XEG(I) = XEGG(I) XLSG(I) = XLSGG(I) XCEG(I) = XCEGG(I) XLLG(I) = XLLGG(I) XDLG(I) = XDLGG(I) DOYSG8(I) = DOYSGG8(I) IYRSG(I) = IYRSGG(I) C SRCG(I) = SRCGG(I) end if end do nTminTS = 1 nTmaxTS = nTG ALng = 1.*(nLng-1)/(iLng-1) if(bForecast) then ! For forecast mode NTV = (NTV + 1.5)/2.0 + 6 NTG = (NTG + 1.5)/2.0 + 6 end if ItoffG = 0 ItoffV = 0 if(NLV.eq.0) bVcon = .FALSE. print *, 'Before MktimesV', nCar, JDCar(1), JDCar(500) if(bVcon) & call MkTimes(bForecast,MJDcntr,NmidHRV,aNdayV,NTmaxV,ALng,nLng,nCar,JDCar, & NTV,IYRBV,IdaybegV,XCintDV,XCintV,XCbeV,XCtbegV,XCtendV,IYRV,DOYV) C & call MkTimesvd(bDaydV,bForecast,MJDcntr,NmidHRV,aNdayV,XCtest1,XCtest2,xInc,NTV,NTmaxV, C & nCar,JDCar,ALng,XCintV,XCbeV,XCtbegV,XCtendV,XCtioffV,RconstV,ItoffV,IYRV,DOYV) if(NLG.eq.0) bGcon = .FALSE. if(bGcon) & call MkTimes(bForecast,MJDcntr,NmidHRG,aNdayG,NTmaxG,ALng,nLng,nCar,JDCar, & NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG) C print *, 'After MktimesG' C & call MkTimesvd(bDaydG,bForecast,MJDcntr,NmidHRG,aNdayG,XCtest1,XCtest2,xInc,NTG,NTmaxG, C & nCar,JDCar,ALng,XCintG,XCbeGG,XCtbegG,XCtendG,XCtioffG,RconstG,ItoffG,IYRG,DOYG) C print *, '7000th source elongation, day and brightness = ', XEG(7000),DOYSG8(7000),GOBS(7000) if(.not.bGcon) then nTG = nTV 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 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 ! One shift matrix per iteration 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 print *, 'XCtbegV, XCtendV, New NTV', XCtbegV, XCtendV, NTV print *, 'XCtbegG, XCtendG, New NTG', XCtbegG, XCtendG, NTG print *, 'XCbeV =' , XCbeV print *, 'XCintV =', XCintV print *, 'XCbeGG =', XCbeGG print *, 'XCintG =', XCintG print *, 'XCintDV =', XCintDV print *, 'XCintDG =', XCintDG C C Find the minimum and maximum time index value in the tiseries C C nTmin = 0 ! Old version C nTmax = NTG C do I=1,NTG C if(XCintDG(I).ge.DOYSG8(1).and.nTmin.eq.0) then C nTminTS = I C nTmin = I C end if C if(XCintDG(I).ge.DOYSG8(NLG).and.nTmax.eq.NTG) then C nTmaxTS = I-1 C nTmax = I C end if C end do print *, 'IYRBG',IYRBG nTmin = 0 nTmax = nTG nTminTS = 1 nTmaxTS = nTG-1 Adaybeg = dble(IdaybegG) tXC1 = sngl(DOYSG8(1)) if(XCintDG(nTG).gt.Adaybeg.and.DOYSG8(1).lt.dble(180.)) tXC1 = sngl(DOYSG8(1) + Adaybeg) tXCN = sngl(DOYSG8(NLG)) if(XCintDG(nTG).gt.Adaybeg.and.DOYSG8(NLG).lt.dble(180.)) tXCN = sngl(DOYSG8(NLG) + Adaybeg) C print *, nTG,XCintDG(nTG),Adaybeg,NLG,DOYSG8(NLG),tXCN do I=1,nTG if(sngl(XCintDG(I)).ge.tXC1.and.nTmin.eq.0) then nTminTS = I nTmin = I end if if(sngl(XCintDG(I)).ge.tXCN.and.nTmax.eq.nTG) then nTmaxTS = I-1 nTmax = I end if end do print *, 'nTminTS, nTmaxTS',nTminTS, nTmaxTS XCtbegF = XCtbegG ! For forecast mode XCtendF = XCtendG ! For forecast mode if(bGcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP) if(bVcon) call ArrR4Constant(nLng*nLat*nTV,BadR4(),VMAP) NLmax = max(NLmaxG,NLmaxV) NTVV = ItoffV NTGG = ItoffG NTVE = 0 NTGE = 0 K = 0 do KKK=1,NLmax+nTG ! Loop before deconvolution to set up VMAP, DMAP K = K + 1 IVSW = 0 ! Switch turn on for NTVV value of K if(K.le.NLV+1) then ! If K is more than the source # + 1 don't turn on switch if(bVcon.and.XCEV(K).gt.XCintV(NTVV+2).and.XEV(K).le.XEhighV) then NTVV = NTVV + 1 NTVB = NTVE + 1 NTVE = K - 1 if(NTVB.ne.NTVE+1) IVSW = 1 if(K.le.NLV) K = K - 1 NTVD = NTVE-NTVB+1 write (*,'(A,I4,A,2I6,2F10.3)') 'V loop ',NTVV,' - S#, K, XCEV, XCintV', NTVD,K,NCoff+XCEV(K),NCoff+XCintV(NTVV+2) if(NTVV.eq.NTV) NLV = NTVE end if end if if(IVSW.eq.1) then C print *, ' V switch ',IVSW,' loop',NTVV if (NLV .eq. 0) call Say('smeiipstd','E','NoV','no velocity IPS data') if (bVproxy) then ! Read proxy point-P V map from file I = iReadProxyMapN(1,NCoff+XCintV(NTVV),nLng,nLat,VMAP(1,1,NTVV)) call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP',NTVV,' min and max after read is', amin, amax end if call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),aminV,amaxV) if (amaxV.eq.badR4()) then ! Try to calculate a proxy point-P V map call MapGrid(XCbeV(1,NTVV),XCbeV(2,NTVV),YLbeg,YLend,NC, & NTVD,XCV(NTVB),YLV(NTVB),VOBS(NTVB),NAV,DAV,nLng,nLat, & VMAP(1,1,NTVV),VDEV(1,1,NTVV),VPNT(1,1,NTVV), & Vmin,Vmax,VDmin,VDmax) C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintV(NTVV) end if ! Calculate weights, WTSV, along all lines of sight call MkLOSWeights(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAPHOLE(1,1,NTVV)) if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,NTVV),CONRV,3,0.0,90.0) end if end if call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,NTVV,' min and max is', amin, amax if(NTVV.eq.ItoffV+1) then ! Fill earlier arrays do N=1,ItoffV amaxV = BadR4() if(bVproxy) then I = iReadProxyMapN(1,NCoff+XCintV(N),nLng,nLat,VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),aminV,amaxV) print *, 'The VMAP' ,N,' min and max after read is', aminV, amaxV end if if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,N,' min and max is', amin, amax if(amax.ne.BadR4()) then call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAP(1,1,N)) call Say('smeiipstd','I','Info','Copy first source array into data file ' & //cStr(:Int2Str(N,cStr))) end if end if end do end if end if IGSW = 0 ! Switch turn on for NTGG value of K if(K.le.NLG+1) then ! If K is more than the source # +1 don't turn on switch if(K.gt.0) then if (.not.bHelio.and.bGcon.and.XCEG(K).gt.XCintG(NTGG+2).and.XEG(K).le.XEhighG) then NTGG = NTGG + 1 NTGB = NTGE + 1 NTGE = K - 1 if(NTGB.ne.NTGE+1) IGSW = 1 if(K.le.NLG) K = K - 1 NTGD = NTGE-NTGB+1 if(.not.bTcon) write (*,'(A,I4,A,2I6,2F10.3)') & 'G loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2) if(bTcon) write (*,'(A,I4,A,2I6,2F10.3)') & 'B loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2) if(NTGG.eq.NTG) NLG = NTGE end if end if end if if(IGSW.eq.1) then if(NTGB.ne.0) then if (NLG .eq. 0) call Say('smeiipstd','E','NoG','no Density or G-level IPS data') if (bDproxy) then ! Read proxy density map from file I = iReadProxyMapN(2,NCoff+XCinTG(NTGG),nLng,nLat,DMAP(1,1,NTGG)) end if call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),aminDMAP,amaxDMAP) if (amaxDMAP.eq.badR4()) then ! Try to calculate a point-P synoptic G map C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintG(NTGG) C if (.not.bHelio) then if (.not.bSMEI) then call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC, & NTGD,XCG(NTGB),YLG(NTGB),GOBS(NTGB),NAV,DAV,nLng,nLat, & DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG), & Gmin,Gmax,GDmin,GDmax) else call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC, & NTGD,XCG(NTGB),YLG(NTGB),ZZ(NTGB),NAV,DAV,nLng,nLat, & DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG), & Gmin,Gmax,GDmin,GDmax) end if end if C do NLLL=NTGB,NTGB+NTGD C if (NLLL.lt.200) print *, NLLL, SRCG(NLLL), GOBS(NLLL) C if (SRCG(NLLL).eq.' 1 1 14') print *, NLLL, SRCG(NLLL), GOBS(NLLL) C end do C Calculate weights, WTSG, according to NLOSWG along all lines of sight . call MkLOSWeights(NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) print *, 'Sample L.O.S. wt. ',NTGB,NTGD,5*NGRESS,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),WTSG(5*NGRESS,NTGB) if(.not.bTcon) then ! Deal with g-level DMAP call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),PPGMAP(1,1,NTGG)) C Convert G to density and smooth density maps if no proxy density maps are input. if(amaxDMAP.eq.badR4()) then call ConvertG2D(0,nLngLat,PPGMAP(1,1,NTGG),DMAP(1,1,NTGG),RRG,DEN1AU,PWG) call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,NTGG),CONRD,3,0.0,90.0) end if print *, 'After G to Den convert the map min and max are', amin, amax end if else ! Deal with Thomson scattering DMAP C Smooth Thomson density maps if no proxy density maps are input. if(amaxDMAP.eq.badR4()) then call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,NTGG),CONRD,3,0.0,90.0) end if print *, 'After Thomson to Den convert the map min and max are', amin, amax end if end if call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) print *, 'The DMAP',NTGG,' min and max is', amin, amax C======== call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAPHOLE(1,1,NTGG)) ! Dummy density map for print out call ArrR4TimesConstant(-nLngLat,DMAPHOLE(1,1,NTGG),RRS*RRS,DMAPHOLE(1,1,NTGG)) ! n*R^2 with holes possible if(NTGG.eq.ItoffG+1) then ! Fill earlier arrays do N=1,ItoffG amaxD = BadR4() if(bDproxy) then I = iReadProxyMapN(2,NCoff+XCintG(N),nLng,nLat,DMAP(1,1,N)) call arrR4getminmax(nLngLat,DMAP(1,1,N),aminD,amaxD) print *, 'The DMAP',N,' min and max after read is', aminD, amaxD end if if(amaxDMAP.eq.BadR4()) then call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) print *, 'The DMAP',N,' min and max is', amin, amax if(amax.ne.BadR4()) then call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAP(1,1,N)) call Say('smeiipstd','I','Info','Copy first source array into data file ' & //cStr(:Int2Str(N,cStr))) end if end if end do end if end if end if end do ! End arrays that set up initial source values and maps amaxallV = BadR4() do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) amaxallV = amax end do if(amaxallV.eq.BadR4()) then do N=1,NTV call ArrR4Constant(nLngLat,Speed,VMAP(1,1,N)) end do print *, 'All velocity maps were bad and all were set to',Speed end if amaxallD = BadR4() do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) amaxallD = amax amaxallD = BadR4() end do if(amaxallD.eq.BadR4()) then TD2 = 1.0 ! Initial total Gvalue for g-level maps if(bTcon) TD2 = 0.66667*DEN1AU*((R1AU/RRD)**FALLOFF) ! divide total density by two-thirds for initial Thomson maps do N=1,NTG call ArrR4Constant(nLngLat,TD2,DMAP(1,1,N)) end do print *, 'All density maps were bad and all were set to',TD2 end if if(.not.bVcon) then do N=1,NTV call ArrR4Constant(nLngLat,Speed,VMAP(1,1,N)) call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N)) end do call arrR4getminmax(nLngLat,VMAP(1,1,NTV),amin,amax) print *, 'All',NTV,' VMAPs beginning with .not. bVcon are', amax end if if(.not.bGcon) then if(.not.bTcon) then do N=1,NTG call ArrR4Constant(nLngLat,1.,DMAP(1,1,N)) call ArrR4Copy(nLngLat,DMAP(1,1,N),PPGMAP(1,1,N)) call ConvertG2D(0,nLngLat,PPGMAP(1,1,N),DMAP(1,1,N),RRG,DEN1AU,PWG) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) call ArrR4Copy(nLngLat,DMAP(1,1,N),VMAPD(1,1,N)) end do call arrR4getminmax(nLngLat,DMAP(1,1,NTG),amin,amax) print *, 'All',NTG,' beginning DMAPs for .not. bGcon are', amax end if end if C FINAL CUT or COUNT NLGG = 1 do I=1,NLG C if(DOYSG8(I).le.149.25) then if(WTSG(1,I).ne.0.0) NLGG = NLGG + 1 C end if end do NLG = NLGG - 1 call Say('smeiipstd','I','Info',' D observations ' & //cStr(:Int2Str(NLG,cStr))) call Say('smeiipstd','I','Info',' V observations ' & //cStr(:Int2Str(NLV,cStr))) ffact = 15./48. ! about fifteen of 48 days total are being reconstructed redundant = float(NLG)/((nLng/Alng)*(nLat)*(nTG*ffact)) print *, 'The D observations are',redundant,' more than the digital resolution requires.' call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) C Produce the initial Thomson scattering base maps if needed if (bTcon) then DB = DEN1AU*((R1AU/RRD)**FALLOFF) DB2 = DB/3.0 do N=1,NTG do i=1,nLng do j=1,nLat DBASE(i,j,N) = DB2 ! Produce initial base map DMAP(i,j,N) = DMAP(i,j,N) + DB2 ! Add initial base to DMAP end do end do end do end if C do I=1, NLG C if(I.lt.200) print *, I, SRCG(I), GOBS(I) C if (SRCG(I).eq.' 1 1 14') print *, I, SRCG(I), GOBS(I) C end do call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTV*aNdayV,DMAP,DMAPV,Vtmpt,DDT,DDD) if(bVproxy) call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV) if(bDproxy) then ! Do this for use after velocity iteration call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,nTV,nTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV) call FillWholeT(1,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD) end if if(.not.bVproxy) then call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV) do N=1,nTV call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) print *, 'DMAPV min max N',amin,amax,N call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP min max N',amin,amax,N call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,90.0) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,90.0) C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N',amin,amax,N C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,nTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP min max N',amin,amax,N end do end if call MkShiftdn(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N', amin,amax,N call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,1,0.0,90.0) end do C C Comment out call MkDensity(nLng,nLat,nMap,XCshift,DDfact,VMAP,RRV,dRR) C Nit = -1 ! Why three iterations counters?? NitV = -1 NitG = -1 do while (Nit .lt. NiterT) Nit = Nit+1 print *, ' ' call Say('smeiipstd','I','Info','Iteration '//cStr(:Int2Str(Nit,cStr))) NVS = 0 if (Nit .eq. LimitoV) NVS = 1 NGS = 0 if (Nit .eq. LimitoG) NGS = 1 if (bVcon) then ! V deconvolution NitV = NitV+1 ! # V iterations ! bGmv2=.TRUE., implies bGcon=.FALSE. if (bGmv2) then ! Make a density map assuming mv^2 = C RRD = RRV call Veleq3(nLng,nLat,NTV,VMAP,SPEEDEQ) print *, 'Solar equatorial mv^2 = constant value SPEEDEQ =',SPEEDEQ do N=1,NTV call Mk_V2DN(nLngLat,VMAP(1,1,N),DMAPV(1,1,N),SPEEDEQ,DEN1AU,FALLOFF,RRD) end do call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,NTV,NTG,ConsTG*aNdayG,DMAPV,VMAPD,Dtmpt,VVT,VVV) do N=1,NTG call ArrR4Copy(nLngLat,VMAPD(1,1,N),DMAP(1,1,N)) end do if(NitV .eq. NiterT) then do N=1,NTG call ArrR4Copy(nLngLat,VMAPD(1,1,N),DMAPHOLE(1,1,N)) end do end if end if if(bstopw) then print *, 'Before MkPostdn V, V deconvolution' call stopwatch(' ',' ') end if call MkPostdn(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) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DVfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,VFAC) call MkVModeltd(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX, & WTSV,DFAC,VFAC,PWV,NLV,NLOSV,NLOSVP1,VMAP,VLIM, & DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij) call FixModeltdn(1,NLV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) if(bWriteMS.and.NitV.eq.NiterT) then call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV,PWV,VRATIO, & NLmaxV,NSRV,SRCV,DOYSV8,XCEV,VOBS,VMOD,XLON,XCtpr,NLOSVP1,FIX,VSSIG,NBSV) end if VRATION(NitV+1) = VRATIO ! Put in 7/9/00 if(NitV.ge.2.and.VRATION(NitV+1).gt.(VRATION(NitV)+(VRATION(NitV+1)*0.05)).and. & VRATION(NitV+1).gt.(VRATION(NitV-1)+(VRATION(NitV+1)*0.10))) then print *, 'THE SUMMED VELOCITY RATIO IS NO LONGER CONVERGING' C NitV = NiterT end if VFACTOR = 1.0 ! Use only in the following subroutine for UCSD velocities if(bVUCSD) VFACTOR = 0.10 ! Use only in the following subroutine for UCSD velocities call MkVMaptdn(LON,LAT,ITIM,RPX,XEV,VWTij,FIX,NLV,NLOSV,NLOSVP1,XLON,XLAT, & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VMOD,VSIGB,NVS,NBSV,VMAPSIG, & ANMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV, & Alng,ConstL,CONRV,CONVT,aNdayV*VFACTOR,VZtmp) do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'In V determination VMAP min max',amin,amax if (NitV .eq. NiterT) then call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N)) call arrR4getminmax(nLngLat,VMAPHOLE(1,1,N),amin,amax) if (amax.ne.BadR4().and.CONRV.gt.ALng*360./float(nLng)) then call GridSphere2D(ALng,nLng,nLat,1,VMAPHOLE(1,1,N),CONRV,0,0.0,0.0) end if end if call ArrR4SetMinMax(-nLngLat,VMAP(1,1,N),VLL,VUL) end do end if if(bstopw) then call stopwatch(' ',' ') print *, 'After V deconvolution and writing VMAPHOLE' print *, 'Arranging DMAP' call stopwatch(' ',' ') end if C C Condition the next DMAP density set to change and to use in mkshiftd. C if (.not.bDproxy.or.Nit.ge.1) then if (bGcon .or. bVmv2) then ! Update density and density shifts even on the last C print *, 'Before first density FillWholeT' call FillWholeT(0,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTG,ConsTG,DMAP,DDD) do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before Grids MkShiftd DMAP min max N',amin,amax,N 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,90.0) C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkShiftd DMAP min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTG,DMAP, 0,1.0*CONDT/aNdayG,0.,Dtmp) do N=1,NTG call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) end do end if if(bstopw) then call stopwatch(' ',' ') print *, 'After arranging DMAP' print *, 'Before arranging VMAPD and D Mkshiftd' call stopwatch(' ',' ') end if if (bGcon) then C C Condition the next velocity map set for density (VMAPD) to use in mkshiftd. C call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,NTV,NTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV) C print *, 'Before second density FillWholeT' call FillWholeT(1,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD) do N=1,NTG C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before Grids MkShiftd VMAPD min max N',amin,amax,N call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAPD(1,1,N),CONRD,3,0.0,90.0) C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before MkShiftd VMAPD min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTG,VMAPD,0,1.0*CONDT/aNdayG,0.,Dtmp) do N=1,NTG call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) end do end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After arranging VMAPD and before D Mkshiftd' call stopwatch(' ',' ') end if if(bGcon) then if(Nit .ge. 20) then C call MkShiftdnn(XCbeGG,XCtbegG,XCtendG,ALng,aNdayG,nLng,nLat,nMap,nTG, C & nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,RconstG,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, C & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) else call MkShiftdn(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG, & nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After D Mkshiftd' end if if (bGcon) then ! G deconvolution NitG = NitG+1 ! # G iterations ! bVmv2=.TRUE., implies bVcon=.FALSE. if (bVmv2) then ! Make a velocity map assuming mv^2 = C RRV = RRD call Deneq3(nLng,nLat,NTG,DMAP,RRD,FALLOFF,DENEQ) print *, 'Solar equatorial 1 AU mv^2 = constant value DENEQ =',DENEQ do N=1,NTG call Mk_D2VN(nLngLat,DMAP(1,1,N),VMAPD(1,1,N),DENEQ,Speed,FALLOFF,RRD) end do call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTV*aNdayV,VMAPD,DMAPV,Vtmpt,DDT,DDD) do N=1,NTV call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAP(1,1,N)) end do if(NitG .eq. NiterT) then do N=1,NTV call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAPHOLE(1,1,N)) end do end if end if if(bstopw) then print *, 'Before D MkPostdn D deconvolution' call stopwatch(' ',' ') end if call MkPostdn(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) if(NitG.eq.NiterT.and.bWritelos) call writelosmap(NLOSGP1,NLG,nTG,NBSG,ITIM,XLON,XLAT) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC) C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkGModeltdn DMAP min max N',amin,amax,N C end do if(.not.bTcon) then ! G-level model-maker from density call MkGModeltdn(XCbegG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP, & nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi) NGval = 2 ! G-level fix else ! Brightness model-maker from density call MkDHModeltd(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,FALLOFF,NLG,NLOSG,DISTG,NLOSGP1,DMAP,DBASE,BBASE, & nLng,nLat,nTG,RRS,GMOD2,GWTij) NGVal = 3 ! Thomson scattering fix end if C C print *, 'XCINTG(78), NLOSG, dLOSG, DMAP(115,18,70) + ', XCINTG(78), NLOSG, dLOSG, DMAP(115,18,70), C & DMAP(115,18,69),DMAP(115,18,71),DMAP(100,18,69),DMAP(100,18,70) do I=1,NLG C GGOBS(I) = GOBS(I) + BBASE(I) ! Each iteration, use modified observed brightnesses end do call FixModeltdn(NGVal,NLG,GGOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) do I=1,NLG GMOD2(I) = GMOD2(I) - BBASE(I) ! Model brightness minus base compares with GOBS in write end do if(bWriteMS.and.NitG.eq.NiterT) then call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO, & NLmaxG,NSRG,SRCG,DOYSG8,XCEG,GOBS,GMOD2,XLON,XCtpr,NLOSGP1,FIX,GSSIG,NBSG) end if GRATION(NitG+1) = GRATIO ! Put in 7/9/00 if(NitG.ge.2.and.GRATION(NitG+1).gt.(GRATION(NitG)+(GRATION(NitG+1)*0.05)).and. & GRATION(NitG+1).gt.(GRATION(NitG-1)+(GRATION(NitG+1)*0.10))) then print *, 'THE SUMMED DENSITY RATIO IS NO LONGER CONVERGING' C Nit = NiterT end if C C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkDmaptdn DMAP min max',amin,amax C end do C if(.not.bTcon) then ! density map maker from g-level call MkDMaptdn(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) else ! density map maker from brightness call MkDMaptdn(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) C Calculate the new base to the density. The base is made from the initial DEN1AU (DEN1AU/3.) and is C assumed constant over the whole map. Some day we may want to iterate this base to fit 1 AU density C (see Mk_Base subroutine). C C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'After MkDmaptdn DMAP min max',amin,amax C end do call Mk_Base(VMAPD,DMAP,nLng,nLat,nTG,nTminTS,nTmaxTS,RRV,DEN1AU,FALLOFF,VEL1AU, & BadR4(),BadR4(),DBASE) end if do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'After Mk_Base DMAP min max',amin,amax if (NitG .eq. NiterT) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG)) end do end if if(bstopw) then call stopwatch(' ',' ') print *, 'After D deconvolution and writing DMAPHOLE' print *, 'Arranging VMAP' call stopwatch(' ',' ') end if if(Nit.ne.NiterT) then ! Update velocity shifts and vmap except on last C C Condition the next VMAP velocity set to change and to use in mkshiftd. C if (bVcon.or.bGmv2) then C print *, 'Before first velocity FillWholeT' call FillWholeT(1,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV) do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'After MkDmaptdn VMAP min max N',amin,amax,N cd IPS call 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,90.0) end do call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) end do end if if (bVcon) then if(bstopw) then call stopwatch(' ',' ') print *, 'After arranging VMAP' print *, 'Before arranging DMAPD and V Mkshiftd' call stopwatch(' ',' ') end if C C Condition the next density map set for velocity (DMAPV) to use in mkshiftd and for use in the C next velocity model. C C print *, 'DMAP(50,9,76),DMAP(50,9,77),DMAPV(50,9,39)',DMAP(50,9,76),DMAP(50,9,77),DMAPV(50,9,39) call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTV*aNdayV,DMAP,DMAPV,Vtmpt,DDT,DDD) C print *, 'DMAP(50,9,76),DMAP(50,9,77),DMAPV(50,9,39)',DMAP(50,9,76),DMAP(50,9,77),DMAPV(50,9,39) C print *, 'Before second velocity FillWholeT' call FillWholeT(0,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV) do N=1,NTV C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'After MkDmaptdn DMAPV min max N',amin,amax,N call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,90.0) end do call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) end do if(bstopw) then call stopwatch(' ',' ') print *, 'After arranging DMAPV and before V Mkshiftd' call stopwatch(' ',' ') end if if (.not.bOnly1) then ! limit analysis to a density iteration only call MkShiftdn(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) end if end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After V Mkshiftd' end if end do ! End of iteration loop J = 0 ! List bad IPS sources do I=1,NLV if (NBSV(I) .eq. 0) then J = J+1 write (*,'(3A,I4,A,F7.1,A,F9.3)') ' Bad V',SRCV(I),' L.O.S. removed, # =', & I,', V =',VOBS(I),', elong. =',XEV(I) end if end do write (*,'(A,I4,A)') ' There were',J,' bad V L.O.S. removed.' J = 0 do I=1,NLG if (NBSG(I).eq.0) then J = J+1 if(.not.bTcon) write (*,'(3A,I6,A,F7.3,A,F9.3)') ' Bad',SRCG(I),' G L.O.S. removed, # =', & I,', G = ',GOBS(I),', elong. =',XEG(I) if(bHelio) write (*,'(3A,I6,A,F7.3,A,F9.3)') ' Bad Helios',SRCG(I),' L.O.S. removed, # =', & I,', B (S10) = ',GOBS(I),', elong. =',XEG(I) if(bSMEI) write (*,'(A,1X,2A,I7,A,F7.3,A,F9.3)') ' Bad SMEI',SRCG(I),' L.O.S. removed, # =', & I,', B (S10) = ',GOBS(I),', elong. =',XEG(I) end if end do write (*,'(A,I5,A)') ' There were',J,' bad density L.O.S. removed.' C Here VMAPHOLE contains the original point-P map (bVcon=.FALSE.) ( or is equal to VMAP do N=1,NTV call ArrR4Copy(nLngLat,VMAPHOLE(1,1,N),VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'Before VMAP N, amin and amax = ', N, amin, amax if(amax.ne.BadR4()) then call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy C print *, 'Before Gridfill Zmin,Zmax', Zmin,Zmax call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove small velocity holes C print *, 'Gridfill Zmin,Zmax', Zmin,Zmax 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 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 C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'After VMAP N, amin and amax = ', N, amin, amax end if end do do N=1,NTG call ArrR4Copy(nLngLat,DMAPHOLE(1,1,N),DMAP(1,1,N)) call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam C Below can't be done with Helios data since Helios is not at the location of Earth C XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region of interest on Map 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) C call FillMapOLS(XCbegG(1,N),XCbegG(2,N),XCmidG,nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 180 deg bad spots 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 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) 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 C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,nLng,nLat, C & NTG,NTV,ConsTV*aNdayV,DMAP,DMAPV,Vtmpt,DDT,DDD) C***** Take data out at different spacecraft (and earth) locations (bExtract tells where) print *, ' Before the first Extractd',XCbegG(1,1),XCtbegG,XCstrt call Extractd(bExtract,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C if(bWrite3D) then C***** Write to disk of results interpolated to the even latitudes, longitudes C and times over the heights (specified as nMapn) of the XCshift model. RRSCON = (RRS/R1AU)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do nTF = nTG Nit = Nit + 1 if(bWrite3D) then NinterC = 0 cPre = ' ' C The first write3D_infotd makes maps that have holes to be used to show Carrington maps print *, 'Before first write3Dinfo' call write3D_infotd3DM(0,cPre,Nit,NiterT,NinterC,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) end if print *, 'Before first 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(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 end if end do print *, 'Before second 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) ! Changed to DMAP by BVJ - 4/20/07 if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes AEmax = AEmaxD/(RRS*RRS) AEmin = AEminD/(RRS*RRS) SpeedD = 0.0 call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if end if end do print *, 'After second filling' if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1).and.XCtbegG.eq.XCtbegV) then 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 do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do do N=1,NTG call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) print *, 'VMAPD N, amin and amax = ', N, amin, amax end do C The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views if(bWrite3D) then print *, 'Before second write3D_infotd' cPre = 'nv3f' NinterD = 3 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM(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) 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) 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 = AEminD/(RRS*RRS) 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' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do do N=1,NTG call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) print *, 'VMAPD N, amin and amax = ', N, amin, amax end do 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(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 end if ************************************************************************************************************************************* C Change the density back since the falloff at distance is produced in xtractdn RRSCON = (R1AU/RRS)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do print *, ' Before 2nd Extractd',XCbegG(1,1),XCtbegG,XCstrt call Extractdn(bExtract,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) if ( bMagnetic ) then end if call Say('smeiipstd','S','Stop','program successfully completed') end