C+ C NAME: C smeiipstd0_intel 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 MKshiftdn0. 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 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 smeiipstd0_intel parameter (NLmaxG = 385000, ! Max # G data points (13000) (400,000) C parameter (NLmaxG = 81000, ! Max # G data points (13000) (400,000) & NLmaxV = 8000, ! Max # V data points (10000) & NGREST = 2, ! G temporal resolution factor better than 1 day (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 (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) & nTmaxx = 55, ! was 48 & nTmaxG = nTmaxx*NGREST+1, ! Max # G times (50) (100) (150) & nTmaxV = nTmaxx*NVREST+1) ! 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 (NintLn = 3, ! High resolution output parameters & NintLa = 3, & NintH = 3) parameter (dRR = 0.1) ! Resolution (AU) in radial direction (0.1) parameter (NTP = 14) ! 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 (NLmaxV), ! Year of source observation & NBSG (NLmaxG), & LON (NLOSGP1,NLmaxG), & LAT (NLOSGP1,NLmaxG), & ITIM (NLOSGP1,NLmaxG) parameter (NLOC1 = 4000) real RARA(NLOC1), & DECDEC(NLOC1), & ORBITVV(NLOC1), & ZSMEIZSMEI(NLOC1) integer ISSISS(NLOC1), & ICIC(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 & VOBS (NLmaxV), ! 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 & SRCVG (NLmaxV)*115, ! Full source info and g-value and velocity file & SRCVGsav(NLmaxV)*115, & SRCVsav (NLmaxV)*9, & SRCB (NLmaxG)*34, & SRCGG (NLmaxV)*12, & cvgfiles(NLmaxG)*11, & cgbfiles(NLmaxG)*32, & 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) 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) & 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), & 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), & 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), & 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,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), & Scale (2) byte NSIDE(9*nLng*nLat) parameter (nCar = 500) parameter (nCarSMEI = 500) real*8 JDCar(nCar), & JD, JEpoch, & dLngSun,dLatSun,dDisSun, & 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, & bGOoty, bGNago, 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./, & bWrerr /.FALSE./, & bDdener /.TRUE./, & bVdener /.TRUE./, & bVverer /.FALSE./, & bDverer /.FALSE./, & bDdenerb /.TRUE./, & bVdenerb /.TRUE./, & bVvererb /.FALSE./, & bDvererb /.FALSE./, & bWrite3D_HR /.TRUE./, & bWrite3D_HRb /.FALSE./, & bWrite3Do /.TRUE./, & bWriteMS /.TRUE./, & bWritelos /.TRUE./, & bStopw /.FALSE./, & bForecast, & bGSource /.TRUE./, ! Write out good sources & bExtract(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,2*.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 & bSAWRITE, ! Write out a STEREO A file & bSBWRITE, ! Write out a STEREO B file & bExtractf(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,2*.FALSE./, & bEWRITEF, ! Write out an Earth filled file & bMWRITEF, ! Write out a Mercury 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 equivalence (bExtract( 3), bEWRITE), & (bExtract( 1), bMWRITE), & (bExtract( 4), bMAWRITE), & (bExtract(10), bUWRITE), & (bExtract(13), bSAWRITE), & (bExtract(14), bSBWRITE), & (bExtractf( 3), bEWRITEF), & (bExtractf( 1), bMWRITEF), & (bExtractf( 4), bMAWRITEF), & (bExtractf(10), bUWRITEF), & (bExtractf(13), bSAWRITEF), & (bExtractf(14), bSBWRITEF) character cPrefix (NTP)*2 /'ME','VE','EA','MA','JU','SA','UR','NE','PL','UL','H1','H2','SA','SB'/ character cPrefixf(NTP)*2 /'m1','v2','e3','m4','j5','s6','u7','n8','p9','u1','hA','hb','s1','s2'/ 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/ parameter (NTVmaxNV = NTmaxx*NVREST-1, & NTGmaxNG = NTmaxx*NGREST-1) integer NTV /NTVmaxNV/, & NTG /NTGmaxNG/ integer NiterT /18/, & iOffUT / 0/, & NAV / 0/, & NC / 1/, ! # rotations averaged & nDUMP / 0/, ! Threshold for #pnts/bin & nFILL / 3/, ! Min. # neighbours for fill & NmidHRV /0/, & NmidHRG /0/, & NLV /0/, & NLVmin /200/, ! Default minimum number of velociy lines of sight & NLG /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/, & ERDLOS /10.0/, & ERVLOS /1.5/, & ERDLOSB /10.0/, & ERVLOSB /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/, & CONRD /14./, & CONRV /14./, & CONRVN /14./, & CONRDN /14./, & CONRVO /8.0/, & CONRDO /8.0/, & CONVT /0.75/, & CONDT /0.75/, & CONDOT /0.45/, & CONVOT /0.45/, & CONVNT /0.75/, & CONDNT /0.75/, & 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 & GLIMF /1.0/, ! G-level input multiplier & GLIMM /0.2/, ! G-level lower limit for sources & SMEImax /5./, ! Maximum SMEI brightness in S10 allowed & SMEImin /-2./ ! Minimum SMEI brightness in S10 allowed real RRS / -15. /, ! Reference distance for deconvolution & VUL / 1500.0 /, ! Velocity upper limit for maps & VLL / 100.0 /, ! Velocity lower limit for maps & DUL / 500.0 /, ! Density upper limit for maps & DLL / 0.0 / ! Density lower limit for maps real DAV /1./, & Speed /400./, & PWG /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.3/ , ! Radial power of density & PWRV /0.3/, ! Radial power of velocity & FALLOFF /2.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, & cWildSMEI*80 C Below specific to SMEI character cSMEI(10)*80 integer iSectBeg /-15/, & iSectEnd / 16/ integer iDaySMEI / 0/, & nDevSMEI /10/ real dSigSMEI /0./ character STRING*80 external iReadNagoyan8, iProcessNagoyan8 external iReadOotyn8, iProcessOotyn8 external EARTH 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) 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 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 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$1$1',I) bVNago = I .eq. 1 bVOoty = I .eq. 2 if(bVOoty) NLVmin = 500 call AskI4('Minimum number of velocity lines of sight?',NLVmin) bVmv2 = .FALSE. if(bVNago) NLOSWV = 1 if(bVOoty) NLOSWV = 2 else ! No V deconvolution bVNago = .FALSE. bVOoty = .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, Nagoya, Ooty$1$1',I) bSMEI = I .eq. 1 bTcon = .FALSE. if (bSMEI) call MAP_TZERO(EARTH,iYr,Doy,.01,nCarSMEI,JDCarSMEI) if(bSMEI) GLIMP = SMEImax if(bSMEI) GLIMM = SMEImin if(bSMEI) bTcon = .TRUE. bGNago = I .eq. 2 bGOoty = I .eq. 3 bGmv2 = .FALSE. cCam = cEnvi//'CAM' if(bSMEI) NLOSWG = 6 ! V light Thomson scattering weights for SMEI if(bGOoty) NLOSWG = 1 if(bGNago) NLOSWG = 1 else ! No g deconvolution bGOoty = .FALSE. bGNago = .FALSE. bSMEI = .FALSE. call AskYN('Do you want to make a density map using mv^2 = con',bGmv2) end if if (bVNago .or. bGNago ) call ForeignFile(iVar,cVar,'nagoya' ,cWildNagoya) if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty) 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 call AskR4('Number of days to combine for velocity?',aNdayV) if(bVNago) 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 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 = -9 if(bGOoty) NmidHRG = -6 if(bSMEI) NmidHRG = 0 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 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. call Julian(11,iYr,Doy,dble(MJDfrst),JEpoch) 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 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 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 call AskR4('Maximum lower limit for source velocity?',VLIML) 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(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(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(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) CONRV = CONRVN if(bVOoty) CONRV = CONRVO call AskR4('Velocity spatial filter (degrees)?',CONRV) if(bVNago) CONVT = CONDNT if(bVOoty) CONVT = CONVOT 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(bVOoty.or.bGOoty) PWG = 0.35 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) PWV = PWG if(bVOoty) 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 C print *, 'XCtest1, XCtest2, NC, XCmargin, Sun__MJDtoJD, MJDref', XCtest1, XCtest2, NC, XCmargin,Sun__MJDtoJD,MJDref XCend= XCbeg + NC XCstrt = XCbeg XCbeg= XCbeg - xInc XCend= XCend + xInc MJDcntr = MJDref C print *, 'MJDref',MJDref else ! Normal mode C if(bSMEI) XCbeg = 2025.2 - NCoff ! January2005 period C if(bSMEI) XCbeg = 2015.8 - NCoff ! Apr_May_2004 period if(bSMEI) XCbeg = 2063.0 - NCoff ! Nov_2007 period. 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') call AskI4('# Rotations averaged',NC) 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 C print *, 'XCtest1, XCtest2, NC, XCmargin, Sun__MJDtoJD, MJDref, MJDcntr', XCtest1, XCtest2, NC, XCmargin,Sun__MJDtoJD,MJDref,MJDcntr 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 C print *, 'MJDref',MJDref end if if(bSMEI) then XCORBEG = 1957.9786 ! SMEI beginning orbit Carrington rotation, doy = 0.0, year 2000, doy 1, 0 UT C ibtt(1) = ((XCbeg + NCoff) - XCORBEG)*27.2753 ! Approx. days since the beginning orbit (old <11/11/08 BVJ) C ibtt(1) = ((XCbeg+.2 + NCoff) - XCORBEG)*27.2753 ! Approx. days since the beginning orbit (new >11/11/08 BVJ) ibtt(1) = ((XCbeg+.15 + NCoff) - XCORBEG)*27.2753 ! Approx. days since the beginning orbit (intel >03/05/09 BVJ) C ibtt(1) = (2015.2259 - XCORBEG)*27.2753 ! about orbit 6600 (actually 6591) print *, XCbeg + NCoff,ibtt(1) ibtt(2) = 0 C iett(1) = ((XCbeg + NCoff + 1.2 + xInc) - XCORBEG)*27.2753 ! Approx. days since the ending orbit (11/11/08 BVJ) iett(1) = ((XCbeg + NCoff + 1.3 + xInc) - XCORBEG)*27.2753 ! Approx. days since the ending orbit (intel >03/05/09 BVJ) 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 good source files$yes',bGSource) call AskYN('Output three D error files$yes',bWrerr) call AskYN('Output a three D file$yes',bWrite3D) ! 3D "d" and "f" file write call AskYN('Output a three D high resolution file$yes',bWrite3D_HR)! 3D "h" file write if(bWrite3D_HR) then nMapb = 1 nMapHR = 15 call AskI4('Minimum high resolution height number - 1-31?',nMapb) call AskI4('Maximum high resolution height number - 1-31?',nMapHR) if(nMapHR.gt.nMap) then nMapHR=nMap nMapb=nMapHR end if call AskYN('Do you want the density error to limit HR densities$yes',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$yes',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 three D files$yes',bWrite3Do) ! Completely filled 3D "o" file write call AskYN('Output a Three D high resolution base$no',bWrite3D_HRb)! Completely filled 3D base "b" file write if(bWrite3D_HRb) then nMapHRb = 5 call AskI4('Minimum height number - 1-31?',nMapHRb) if(nMapHRb.gt.nMap) then nMapHRb=nMaps nMapb=nMapHRb end if call AskYN('Do you want the density error to limit HR densities$yes',bDdenerb) if(bDdenerb) then call ASKR4('What density L.O.S. crossing limit?',ERDLOSB) call AskYN('Do you want the density error to limit HR velocities$yes',bVdenerb) end if call AskYN('Do you want the velocity error to limit HR velocities$no',bVvererb) if(bVvererb) then bVdenerb = .FALSE. call ASKR4('What velocity L.O.S. crossing limit?',ERVLOSB) call AskYN('Do you want the velocity error to limit HR densities$no',bDvererb) if(bDvererb) bDdenerb = .FALSE. end if end if 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 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$yes',bSBWRITE) ! STEREO B extraction end if call AskYN('Output an Earth parameter filled file$yes',bEWRITEF) ! Earth filled extraction call AskYN('Output a Mercury parameter filled file$no',bMWRITEF) ! Mercury filled extraction call AskYN('Output a Mars parameter filled file$no',bMAWRITEF) ! Mars filled extraction call AskYN('Output a Ulysses parameter filled file$no',bUWRITEF) ! Ulysses filled extraction 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$yes',bSBWRITEF) ! STEREO B filled extraction end if 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(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,XEVsav,XCVsav,YLVsav,VVsav,GGsav) if (bGOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG,SRCVGsav,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, & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG,SRCVGsav,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, & ICIC,RARA,DECDEC,ORBITVV,ZSMEIZSMEI,ISSISS, & NORBITST,MX,NLOC,NLmaxG,NS,NLT1, & cgbfiles,SRCB,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. VLIMM) NBSV(I) = 0 if (VOBS(I) .lt. VLIML) 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)*((abs(XEG(I))**2.5)/(90.0**2.5)) .gt. GLIMP) NBSG(I) = 0 if (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) .lt. GLIMM) NBSG(I) = 0 C if (I .lt. 500) print *, 'Sample value of XEG =', XEG(I), GOBS(I), GOBS(I)*((XEG(I)**2.5)/(90.0**2.5)), GLIMP, GLIMM 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) 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.lt.NLVmin) then bVcon = .FALSE. NLV = 0 print *, 'There are too few velocity sources for a velocity deconvolution. NLV set to zero.' end if print *, 'Before MktimesV', nCar, JDCar(1), JDCar(500) if(bVcon) & call MkTimes(bForecast,MJDcntr,NmidHRV,aNdayV,NTmaxV,ALng,nLng,nCar,JDCar, & NTmaxV,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, & NTmaxG,NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG) 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 timeseries C print *, 'IYRBG',IYRBG,'IdaybegV',IdaybegV,'IdaybegG',IdaybegG 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 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,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(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 (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) 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 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,0.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,0.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 if(WTSG(1,I).ne.0.0) NLGG = NLGG + 1 end do NLG = NLGG - 1 call Say('smeiipstd','I','Info',' D observations ' & //cStr(:Int2Str(NLG,cStr))) C **** write (*,'(1X,A)') & 'The last 10 SMEI brightness values were:' NLGM10 = NLG-10 do I=NLGM10,NLG write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3,A,F8.3)') & 'I =',I,' Year =',IYRSG(I),' DOY =',DOYSG8(I),' Brightness =',GOBS(I),' Elong =',XEG(I) end do C **** 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 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,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)) call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP min max N',amin,amax,N end do end if call MkShiftdn0(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,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('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 MkPostd V, V deconvolution' call stopwatch(' ',' ') end if call MkPostd(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) if(NitV.eq.NiterT.and.bWritelos) call writelosmapvg(0,NLOSVP1,NLV,nTV,NBSV,NCOFF,ITIM,XLON,XLAT) 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,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 VFACTOR = 1.0 ! Use only in the following subroutine for UCSD velocities call MkVMaptdn0n(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,AWTMAPV,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 mkshiftdn0. 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 MkShiftdn0 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 MkShiftdn0 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 Mkshiftdn0' call stopwatch(' ',' ') end if if (bGcon) then C C Condition the next velocity map set for density (VMAPD) to use in mkshiftdn0. C do N=1,nTV vmaptot = 0. vmapsum = 0. do I=1,nlng do J=1,nlat if(VMAP(I,J,N).ne.badR4()) then vmaptot = vmaptot + VMAP(I,J,N) vmapsum = vmapsum + 1. end if end do end do vmapave = 0. if(vmapsum.ne.0) vmapave = vmaptot/vmapsum C if(N.eq.30) print *, 'N and Vmap average',N,vmapave 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) do N=1,nTG vmapave1 = 0. if((N+1)/2.lt.nTV) then vmaptot1 = 0. vmapsum1 = 0. do I=1,nlng do J=1,nlat if(VMAP(I,J,(N+1)/2).ne.badR4()) then vmaptot1 = vmaptot1 + VMAP(I,J,(N+1)/2) vmapsum1 = vmapsum1 + 1. end if end do end do if(vmapsum1.ne.0) vmapave1 = vmaptot1/vmapsum1 end if vmapave2 = 0. if((N+1)/2+1.lt.nTV) then vmaptot2 = 0. vmapsum2 = 0. do I=1,nlng do J=1,nlat if(VMAP(I,J,(N+1)/2+1).ne.badR4()) then vmaptot2 = vmaptot2 + VMAP(I,J,(N+1)/2+1) vmapsum2 = vmapsum2 + 1. end if end do end do if(vmapsum2.ne.0) vmapave2 = vmaptot2/vmapsum2 end if vmaptot = 0. vmapsum = 0. do I=1,nlng do J=1,nlat if(VMAPD(I,J,N).ne.badR4()) then vmaptot = vmaptot + VMAPD(I,J,N) vmapsum = vmapsum + 1. end if end do end do vmapave = 0. if(vmapsum.ne.0) vmapave = vmaptot/vmapsum C if(N.eq.60) print *, 'N and Vmap average',N,vmapave,vmapave1,vmapave2 end do 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 MkShiftdn0 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 MkShiftdn0 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 Mkshiftdn0' 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 MkShiftdn0(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 Mkshiftdn0' 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 MkPostd D deconvolution' call stopwatch(' ',' ') end if call MkPostd(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 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) 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,SRCB,DOYSG8,XCEG,GOBS,GMOD2,XLON,XLproj,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 MkDmaptdn0 DMAP min max',amin,amax C end do C if(.not.bTcon) then ! density map maker from g-level call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) else ! density map maker from brightness call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,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 MkDmaptdn0 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 mkshiftdn0. 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 MkDmaptdn0 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,0.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 Mkshiftdn0' call stopwatch(' ',' ') end if C C Condition the next density map set for velocity (DMAPV) to use in mkshiftdn0 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 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) 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 Mkshiftdn0' call stopwatch(' ',' ') end if if (.not.bOnly1) then ! limit analysis to a density iteration only call MkShiftdn0(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 Mkshiftdn0' end if end do ! End of iteration loop if(bGSource) then ! Write out good sources used call writegoodsourceb(cgbfiles,SRCB,GOBS,NBSG,NLG) istation = 1 if(bVNago) istation = 1 if(bVOoty) istation = 2 call writegoodsourcev(istation,cvgfiles,SRCVG,VOBS,NBSV,NLV) end if 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',SRCGG(I),' G L.O.S. removed, # =', & I,', G = ',GOBS(I),', elong. =',XEG(I) if(bSMEI) write (*,'(A,1X,2A,F7.3)') ' Bad SMEI L.O.S. at',SRCB(I),' removed, B(S10) = ',GOBS(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 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 Extractdn_intel',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdn_intel(bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,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. 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 = 'nv3d' 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 if(bWrerr) then 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() if(AWTMAPD(I,J,N).eq.0.0) AWTMAPD(I,J,N) = badR4() end do end do end do PIRSQD = (CONDT/aNdayG)*3.14159*((CONRD/(20.0/float(NGRESS)))**2) ! normalize weights and lines of sight to temporat and 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 CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,VDEV,ANMAPV,Dtmpt,VVT,VVV) call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,GDEV,AWTMAPV,Dtmpt,VVT,VVV) PIRSQV = (CONVT/aNdayV)*3.14159*((CONRV/(20.0/float(NVRESS)))**2) ! normalize weights and lines of sight to temporal and 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 NinterC = 0 cPre = 'er1_' print *, 'Before the error1 (L.O.S. crossings) write' call write3D_infotd3DM(2,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,ANMAPV,ANMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) print *, 'Before the error2 (projected weights) write' NinterC = 0 cPre = 'er2_' call write3D_infotd3DM(2,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,AWTMAPV,AWTMAPD,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) C AEmin = AEmin/(RRS*RRS) (needs fixing) AEmin = AEminD/(RRS*RRS) ! (fixed 2/29/2008 BVJ) 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 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 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 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 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) C The second high resolution write3D_infotd makes maps that have fewer holes to be used to show smoothly as C remote observer views and planar cuts. end if 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 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR(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, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdener,bVdener,bVverer,bDverer,ERDLOS,ERVLOS,ANMAPD,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if if(bforecast) then 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 Extractdn_intel',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdn_intel(bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) RRSCON = (RRS/R1AU)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do 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) 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 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 if(bWrite3D_HRb) then 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 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 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR(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, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if end if C************************************************************************************************************************************* C Change the density back since the falloff at distance is produced in Extractdn if(.not.bforecast) then 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 Extractdn_intel',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdn_intel(bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if if ( bMagnetic ) then end if call Say('smeiipstd','S','Stop','program successfully completed') end