C+ C NAME: C smeiipstd0nvhrpg_inpv24 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 MKshiftdn0n. 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 ./smeiipstd0npg_inpv24_intel 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 smeiipstd0npg_inpv24 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 April 2020 Bernie Jackson (UCSD) (runs on 64-node Leela machine) C- program smeiipstd0nvhrpg_inpv24 parameter (NLmaxG = 5000000, ! 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 = 8, ! 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 = 8, ! 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, & 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) & nLngLat = nLng*nLat) parameter (NintLn = 3, ! High resolution output parameters & NintLa = 3, & NintH = 3) parameter (dRR = 0.02) ! Resolution (AU) in radial direction (0.1) parameter (nMap = 30*(0.1/dRR)+1) ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) 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), & LON_in(NLmaxG), & LAT_in(NLmaxG), & ITIM_in(NLmaxG) parameter (NLOC1 = 20000) 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 & GOBSS (NLmaxG), ! Solar wind observed G values & 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 & GMODD_in(NLmaxV), ! Model density values & FIX (NLmaxG), ! Fix to model to make OK & GSSIG (NLmaxG), ! G-level source sigma & WTSG (NLOSG ,NLmaxG), & GWTij (NLOSG ,NLmaxG), & XLON (NLOSGP1,NLmaxG), & XLONP (NLOSGP1,NLmaxG), & XLAT (NLOSGP1,NLmaxG), & RPX (NLOSGP1,NLmaxG), & XLproj(NLOSGP1,NLMAXG), & VFAC (NLOSGP1,NLmaxG), ! New array (3/99) & DFAC (NLOSGP1,NLmaxG), & XLON_in (NLmaxG), & XLONP_in (NLmaxG), & XLAT_in (NLmaxG), & RPX_in (NLmaxG), & XLproj_in(NLmaxG), & VFAC_in (1,NLmaxG), & DFAC_in (1,NLmaxG) character SRCV (NLmaxV)*9, ! Velocity source name & SRCVLAST*9, & SRCVG (NLmaxV)*115, ! Full source info and g-value and velocity file & SRCVGsav(NLmaxV)*115, & SRCVsav (NLmaxV)*9, & SRCB (NLmaxG)*34, & SRCGG (NLmaxV)*12, & SRCBA (NLmaxV)*115, & SRCGA (NLmaxV)*9, & SRCGLAST*12, & SRCBLAST*34, & SRCALAST, & cvgfiles(NLmaxG)*11, & cggfiles(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 & VMOD (NLmaxV), ! Model velocities & VSSIG (NLmaxV), ! Velocity source sigma & VWTij (NLOSV,NLmaxV), & WTSV (NLOSV,NLmaxV), & VWTij_in(NLmaxV), & GWTij_in(NLmaxV), & WTSV_in (NLmaxV), & WTSD_in (NLmaxV), & WTPV_in (NLmaxV), & WTPD_in (NLmaxV) real XCinTV(nTmaxV), & XCinTG(nTmaxG), & DOYV (nTmaxV), & DOYG (nTmaxG), & XCbeV (2,nTmaxV), & XCbeGG(2,nTmaxG), & XCbe (2,4*nTmaxG) integer IYRV (nTmaxV), & IYRG (nTmaxG) real VDEV (nLng,nLat,nTmaxV), & VDEVD (nLng,nLat,nTmaxG), & 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), & VMAPF (nLng,nLat,nTmaxG), & DMAPV (nLng,nLat,nTmaxV), & GDEV (nLng,nLat,nTmaxG), & GPNT (nLng,nLat,nTmaxG), & DDfact(nLng,nLat,nMap,nTmaxG), & DVfact(nLng,nLat,nMap,nTmaxG), ! New array (3/99) & XCshift(nLng,nLat,nMap,nTmaxG,3), & V3D (nLng,nLat,nMap), & D3D (nLng,nLat,nMap), & V3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & D3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BR3DHR(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BT3DHR(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BN3DHR(nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & V3DT (nLng,nLat,nMap,4*nTmaxG), & D3DT (nLng,nLat,nMap,4*nTmaxG), & Vtmp (nLng,nLat,nTmaxG,3), ! Scratch real arrays & VtmpV (nLng,nLat,nTmaxV,3), & VZtmp (nLng,nLat,nTmaxV), & Dtmp (nLng,nLat,nTmaxG), & XLTtmp(nLng,nTmaxG,5), & XLtmp (nLng,nTmaxG), & XLtmpt(nLng,nTmaxG), & Dtmpt (nLng,nTmaxG), & Vtmpt (nLng,nTmaxG), 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,nTmaxV), & ANDMAPV(nLng,nLat,nTmaxG), & AWTDMAPV(nLng,nLat,nTmaxG), & VMAV (nLng,nLat,nTmaxV), & FIXMV (nLng,nLat,nTmaxV), & WTSMV (nLng,nLat,nTmaxV), & VMEANFV(nLng,nLat,nTmaxV), & VFIXM2V(nLng,nLat,nTmaxV), & WTPV (NLOSV,NLmaxV), & VVV (nLng), & DDD (nLng), & VVT (nTmaxV), & DDT (nTmaxG), & VRATION (200), & GRATION (200), & Scale (2) byte NSIDE(9*nLng*nLat) parameter (nCar = 700) parameter (nCarSMEI = 700) 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), & XCtpr_in(NLMAXG), & XCtim_in(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, & bVACE, bDACE, bFACE, bDACEn, bBFIX, & bVFACE, bVACE2, bDACE2, & bVWind, bDWind, & bVCELIAS,bDCELIAS, & 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 & bGWRITE /.FALSE./, ! Write out a G data file & bVWRITE /.FALSE./, ! Write out a V data file & bMirror /.FALSE./, & bjMirror /.FALSE./, & bWrite3D1 /.FALSE./, & bWrite3D2 /.FALSE./, & bWrerr /.FALSE./, & bDdener /.TRUE./, & bVdener /.FALSE./, & bVverer /.TRUE./, & bDverer /.FALSE./, & bDdenerb /.TRUE./, & bVdenerb /.FALSE./, & bVvererb /.TRUE./, & bDvererb /.FALSE./, & bWrite3D_HR /.TRUE./, & bWrite3D_HRb /.FALSE./, & bWrite3Do /.TRUE./, & bWriteMS /.TRUE./, & bWritelos /.TRUE./, & bOneOr1 /.TRUE./, & bOneOr2 /.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/ parameter iVarp = 27 ! Number of potential data files accessed via a parameter list (6 IPS + 5 in-situ + 16 mag) integer NiterT /18/, & iOffUT / 7/, & 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 /100/, ! 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 & DMAX /100./, ! Density upper limit for in-situ inputs - set low for spikes & DMIN /0.01/, ! Density lower limit for in-situ inputs & VVMAX /1500./, ! Velocity upper limit for in-situ inputs & VVMIN /100./ ! Velocity lower limit for in-situ inputs real RRS / -15. /, ! Reference distance for deconvolution & VUL / 1500.0 /, ! Velocity upper limit for maps & VLL / 100.0 /, ! Velocity lower limit for maps & DUL / 500.0 /, ! Density upper limit for maps & DLL / 0.0 / ! Density lower limit for maps real DAV /1./, & Speed /400./, & FALLOFFD /2.0/, ! Density falloff with solar distance (>2.00 Implies & FALLOFFV /0.0/, ! an acceleration of the solar wind with height.) & FALLOFFBR /2.0/, & FALLOFFBT /1.0/, & FALLOFFBN /1.34/, & PWG /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.3/ , ! Radial power of density & PWRV /0.3/ ! Radial power of velocity 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 limit character cCam*4, & cMon*3, c & cStrID*11, TJD, never ref & cVar(27)*256, & cArg*100, & cStr*80, & cPre*4 character cWildNagoya*80, & cWildOoty*80, & cWildSMEI*80, & cWildACE*80, & cWildACE2*80, & cWildWind*80, & cWildVCELIAS*80, & cWildDCELIAS*80, & cWild*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 real Bbegin3(3), ! Interval select for nv3h files & Eend3(3) real Bbeginb3(3), ! Interval select for nv3b files & Eendb3(3) real Bbeginm3(3), ! Interval select for magnetic files & Eendm3(3) 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 include 'b3d_param_3.h' character cWildMAG(B3D__NN)*256 logical bWildMAG(B3D__NN) logical bMagnetic /.TRUE./, ! Magnetic map available & bMagnet /.FALSE./, ! Magnetic map check & bmagexx /.TRUE./, ! Produce 3D magnetic field data files & bmagbase /.TRUE./, ! Produce the same intervals as the base files & bAdjust /.TRUE./, ! Adjust the magnetic field BF factors for solar cycle and instrument types & bmagex /.TRUE./, ! Extract magnetic field data as well as velocity and density data & bFillmag /.TRUE./ 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), !scratch & BN3D(nLng,nLat,nMap), !scratch & BRR2DT(nLng,nLat,nTmaxG), !scratch & BRC2DT(nLng,nLat,nTmaxG), !scratch & BTC2DT(nLng,nLat,nTmaxG), !scratch & BNC2DT(nLng,nLat,nTmaxG), !scratch & TT3D (4*nTmaxG), !holds XC corresponding to 6 hour int. time steps output by write3dinfo_td3d & XC3DT (nLng,nLat,nMap,4*nTmaxG,3), & BRcheck (nLng,nLat), & Vcheck5 (nLng,nLat), & Vcheck7 (nLng,nLat) save SRCVGsav,SRCsav,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav save XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav call T3D_set_mode(T3D__MODE_TIME, 1) iVar = iVarp 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) write(*,'(A,I5,2F17.8)') 'The max CR # and lower and upper JD',nCar, JDCar(1), JDCar(nCar) 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 7 hours (the difference between PDST and GST). call AskYN('Forecast mode?$no',bForecast) ! Regular or forecast mode if (bForecast) then call AskI4('Offset between local time and UT (hours)?',iOffUT) NLVmin = NLVmin/2 end if 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 call AskYN('Use in-situ velocity data in 3-D reconstructions?$yes',bVFACE) bVACE = .FALSE. bVACE2 = .FALSE. bVWind = .FALSE. bVCELIAS = .FALSE. if(bVFACE) then call AskWhat('Use in-situ velocity data from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) bVACE = I .eq. 1 bVACE2 = I .eq. 2 bVWind = I .eq. 3 bVCELIAS = I .eq. 4 end if if(bVFACE) call AskR4('Minimum velocity for in-situ velocity values',VVMIN) if(bVFACE) call AskR4('Maximum velocity for in-situ velocity values',VVMAX) else ! No V deconvolution bVNago = .FALSE. bVOoty = .FALSE. call AskYN('Do you want to make a velocity map using mv^2 = con$no',bVmv2) end if if(.not.bVFACE) bVACE = .FALSE. if (bGcon) then ! G 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. if(bSMEI) bGcon = .FALSE. 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 call AskYN('Use in-situ density data in 3-D reconstructions?$yes',bFACE) if(bGcon.and.bFACE) then call AskWhat('Use in-situ density data from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) bDACE = I .eq. 1 bDACE2 = I .eq. 2 bDWind = I .eq. 3 bDCELIAS = I .eq. 4 end if if(bTcon.and.bFACE) then call AskWhat('Use in-situ density data from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) bDACE = I .eq. 1 bDACE2 = I .eq. 2 bDWind = I .eq. 3 bDCELIAS = I .eq. 4 end if if(.not.bTcon.and.bFACE) then call AskWhat('Use in-situ density data from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) bDACE = I .eq. 1 bDACE2 = I .eq. 2 bDWind = I .eq. 3 bDCELIAS = I .eq. 4 end if if(bDACE.and.bFACE) call AskYN('Fit in-situ density data to ACE level 0 data$yes',bFACE) if(bDACE2.and.bFACE) call AskYN('Fit in-situ density data to ACE level 2 data$yes',bFACE) if(bDWind.and.bFACE) call AskYN('Fit in-situ density data to Wind$yes',bFACE) if(bDCELIAS.and.bFACE) call AskYN('Fit in-situ density data to CELIAS$yes',bFACE) if(.not.bFACE) bDCELIAS = .FALSE. if(bFACE) call AskR4('Minimum density for fit in-situ density values',DMIN) if(bFACE) call AskR4('Maximum density for fit in-situ density values',DMAX) DLIM = DMAX/(RRS*RRS) if(bFACE) call AskYN('Fit average in-situ density data to brigntness?$yes',bBFIX) if(bFACE) call ASkYN('Use in-situ density in the 3-D reconstructions?$yes',bDACEn) 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) if(bTcon) then print *, ' ' print *, 'The directory used for the SMEI tomography analysis follows:' print *, cWildSMEI end if if (bVACE .or .bDACE) then call ForeignFile(-iVar-1,cVar,'ace',cWildACE) if (cWildACE .eq. ' ') then cWildACE = '/home/soft/dat/insitu/acesw_????.hravg' end if write(*,'(A75)') cWildACE end if if (bVACE2 .or .bDACE2) then call ForeignFile(-iVar-1,cVar,'swace',cWildACE2) if (cWildACE2 .eq. ' ') then cWildACE2 = '/home/soft/insitu/swace_????.hravg' end if write(*,'(A75)') cWildACE2 end if if (bVWind .or .bDWind) then call ForeignFile(-iVar-1,cVar,'wind',cWildWind) if (cWildWind .eq. ' ') then cWildWind = '/home/soft/dat/insitu/wind_swe/????_WIND_hourly_averages' end if write(*,'(A75)') cWildWind end if if (bVCELIAS) then call ForeignFile(-iVar-1,cVar,'celias',cWildVCELIAS) if (cWildVCELIAS .eq. ' '.and..not.bForecast) then cWildVCELIAS = '/home/soft/dat/insitu/celias_????.hravg' end if if (cWildVCELIAS .eq. ' '.and.bForecast) then cWildVCELIAS = '/home/soft/dat/insitu/realcelias/celias_realtime.hravg*' end if write(*,'(A75)') cWildVCELIAS end if if (bDCELIAS) then call ForeignFile(-iVar-1,cVar,'celias',cWildDCELIAS) if (cWildDCELIAS .eq. ' '.and..not.bForecast) then cWildDCELIAS = '/home/soft/dat/insitu/celias_????.hravg' end if if (cWildDCELIAS .eq. ' '.and.bForecast) then cWildDCELIAS = '/home/soft/dat/insitu/realcelias/celias_realtime.hravg*' end if write(*,'(A75)') cWildDCELIAS end if if(bVcon) then bDaydV = .TRUE. call AskYN('Deconvolve V from one day to the next?$yes',bDaydV) if(bDaydV) then call AskR4('Number of days to combine for velocity?',aNdayV) if(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(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.or.bTcon) 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) C print *, 'before JD', iYr,Doy,XClo,xchi XClo = NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XClo = XClo-.5 C call Local2UT(iOffUT,iYr,Doy) call Local2UT(NmidHRV,iYr,Doy) ! Changed 5/2/11 BVJ XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo) C 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.or.bTcon) 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.or.bTcon) 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?',FALLOFFD) 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) Idoy = int(Doy) HR = (Doy-float(Idoy))*24.0 IH = int(HR) IMN = nint((HR-float(IH))*60.0) call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy) write(*,'(A,I2.2,A,I2.2,A,I2.2,A,A,I5,A,F8.3,A)') & ' (This MJD is ',IH,':',IMN,' UT ',iDay,' ',cMon,iYr,' -- DOY',Doy,')' call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) C------- Get modified Carrington variable for center of plot C [Note that EarthLng=(1-frac(XCEarth))*360] XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XCbeg = XCEarth-.5 call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) 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 NC = 0 NCC = NC 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 C if(bSMEI) XCbeg = 2063.0 - NCoff ! Nov_2007 period C if(bSMEI) XCbeg = 2003.3 - NCoff ! May_2003 period C if(bSMEI) XCbeg = 2061.0 - NCoff ! Eiscat paper period if(bSMEI) XCbeg = 2003.0 - NCoff ! May_2003 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) ibtt(1) = ((XCbeg+.2 + NCoff) - XCORBEG)*27.2753 ! Approx. days since the beginning orbit (new >11/11/08 BVJ) C ibtt(1) = (2015.2259 - XCORBEG)*27.2753 ! about orbit 6600 (actually 6591) C print *, XCbeg + NCoff,ibtt(1) ibtt(2) = 0 iett(1) = ((XCend + NCoff) - XCORBEG)*27.2753 ! Approx. days since the ending orbit (intel >02/02/11 BVJ) iett(2) = 0 call smei_orbit2(ibtt,NORBITST,dorbit) C print *, dorbit call smei_orbit2(iett,iorbite,dorbit) MX = iorbite - NORBITST + 1 write(*,'(A,I8,A,I5/)') ' 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$no',bGSource) call AskYN('Output three D error files$no',bWrerr) NinterC = 3 ! Default first extract matrix NinterD = 3 ! Default forecast extract matrix C The following parameter tells how many intermediate low resolution error file times to write NinterCC = 0 NOneOr1 = -2 ! tells d 3D write not to write call AskYN('Output 3D files with holes$no',bWrite3D1) ! 3D "d" file write if(bWrite3D1) then call AskI4('Number of intermediate rotation densities with holes?',NinterC) ! Affects first extract call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr1) NOneOr1 = 0 ! tells d 3D write to write if(bOneOr1) NOneOr1 = 1 end if NOneOr2 = -2 ! tells f 3D write not to write call AskYN('Output 3D mostly filled files?$no',bWrite3D2) ! 3D "f" file write if(bWrite3D2) then call AskI4('Number of intermediate filled rotation densities?',NinterD) ! Affects forecast extract call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr2) NOneOr2 = 0 ! tells f 3D write to write if(bOneOr2) NOneOr2 = 1 end if call AskYN('Output a three D high resolution file$yes',bWrite3D_HR) ! 3D "h" file write if(bWrite3D_HR) then NBE3S = 1 do I=1,3 BBegin3(I) = -1.0 EEnd3(I) = -1.0 end do call AskI4('How many intervals do you want to specify (up to 3)?',NBE3S) if(NBE3S.ne.0) then do I=1,NBE3S BBegin3(I) = XCbeg + float(NCOFF) + xInc EEnd3(I) = XCend + float(NCOFF) - xInc if(I.eq.1) then BBegin3(I) = 2003.5990 EEnd3(I) = 2003.6100 call ASKR4(' What is the beginning first limit?',BBegin3(I)) call ASKR4(' What is the ending first limit?',EEnd3(I)) BBegin3(I) = BBegin3(I) - float(NCOFF) EEnd3(I) = EEnd3(I) - float(NCOFF) BBegin3(2) = BBegin3(I) EEnd3(2) = EEnd3(I) BBegin3(3) = BBegin3(I) EEnd3(3) = EEnd3(I) end if if(I.eq.2) then call ASKR4('What is the beginning second limit?',BBegin3(I)) call ASKR4(' What is the ending second limit?',EEnd3(I)) BBegin3(I) = BBegin3(I) - float(NCOFF) EEnd3(I) = EEnd3(I) - float(NCOFF) BBegin3(3) = BBegin3(I) EEnd3(3) = EEnd3(I) end if if(I.eq.3) then call ASKR4(' What is the beginning third limit?',BBegin3(I)) call ASKR4(' What is the ending third limit?',EEnd3(I)) BBegin3(I) = BBegin3(I) - float(NCOFF) BBegin3(I) = EEnd3(I) - float(NCOFF) end if end do end if C RRRS = RRS call AskR4('Minimum high-resolution height in AU?',RRRS) if(RRRS.le.RRS) then nMapb = 1 else nMapb = nint((RRRS-RRS)/dRR) + 1 end if if(nMapb.gt.nMap) nMapb = nMap ! Limit minimum height to maximum available RRRSM = 1.5 call AskR4('Maximum high-resolution height in AU?',RRRSM) nMapHR = nint((RRRSM-RRS)/dRR) + 1 if(nMapHR.gt.nMap) nMapHR = nMap ! Limit maximum height to maximum available if(nMapHR.lt.nMapb) nMapHR = nMapb ! Don't be stupid! if(bForecast) bDdener = .FALSE. call AskYN('Do you want the density error to limit HR densities?',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 if(bForecast) bVdener = .FALSE. 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$no',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 NBE3S = 1 do I=1,3 BBeginb3(I) = -1.0 EEndb3(I) = -1.0 end do call AskI4('How many intervals do you want to specify (up to 3)?',NBE3S) if(NBE3S.ne.0) then do I=1,NBE3S BBeginb3(I) = XCbeg + float(NCOFF) + xInc EEndb3(I) = XCend + float(NCOFF) - xInc if(I.eq.1) then BBeginb3(I) = 2003.57 EEndb3(I) = 2003.64 call ASKR4(' What is the beginning first limit?',BBeginb3(I)) call ASKR4(' What is the ending first limit?',EEndb3(I)) BBeginb3(I) = BBeginb3(I) - float(NCOFF) EEndb3(I) = EEndb3(I) - float(NCOFF) BBeginb3(2) = BBeginb3(I) EEndb3(2) = EEndb3(I) BBeginb3(3) = BBeginb3(I) EEndb3(3) = EEndb3(I) end if if(I.eq.2) then call ASKR4('What is the beginning second limit?',BBeginb3(I)) call ASKR4(' What is the ending second limit?',EEndb3(I)) BBeginb3(I) = BBeginb3(I) - float(NCOFF) EEndb3(I) = EEndb3(I) - float(NCOFF) BBeginb3(3) = BBeginb3(I) EEndb3(3) = EEndb3(I) end if if(I.eq.3) then call ASKR4(' What is the beginning third limit?',BBeginb3(I)) call ASKR4(' What is the ending third limit?',EEndb3(I)) BBeginb3(I) = BBeginb3(I) - float(NCOFF) BBeginb3(I) = EEndb3(I) - float(NCOFF) end if end do end if C RRRRS = RRS call AskR4('Minimum high-resolution base height in AU?',RRRRS) if(RRRRS.le.RRS) then nMapbb = 1 else nMapbb = nint((RRRRS-RRS)/dRR) + 1 end if if(nMapbb.gt.nMap) nMapbb=nMap ! Limit maximum height to maximum available RRRRSM = 1.5 call AskR4('Maximum high-resolution base height in AU?',RRRRSM) nMapHRb = nint((RRRRSM-RRS)/dRR) + 1 if(nMapHRb.gt.nMap) nMapHRb = nMap ! Limit maximum height to maximum available if(nMapHRb.lt.nMapbb) nMapHRb = nMapbb ! Don't be stupid! if(bForecast) bDdener = .FALSE. 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.S2/7/2011. 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 C C The following parameter tells how many intermediate high resolution times to write, however ask C Since main densities are produced every 1.0/NGREST days, and this parameter means the density is C written out every 1.0/(NGREST*(NinterDD+1)) days, Thus for C NGREST = 2 and NinterDD = 3, a file is written every 3 hours. (preferred) C NGREST = 3 and NinterDD = 3, a file is written every 2 hours (OK) C NGREST = 4 and NinterDD = 2, a file is written every 2 hours (preferred) C NGREST = 4 and NinterDD = 3, a file is written every 1.5 hours (every 1 hr and 40 min) C NGREST = 5 and NinterDD = 2, a file is written every 1.66667 hours (every 1 hr and 20 min) C NGREST = 6 and NinterDD = 2, a file is written every 1.3333 hours (every 1 hr and 20 min) C NGREST = 6 and NinterDD = 3, a file is written every 1 hour (preferred) C NGREST = 8 and NinterDD = 2, a file is written every 1 hour (preferred) C NGREST= 12 and NinterDD = 1, a file is written every 1 hour (preferred) C NinterDD = 3 if(NGREST.gt.2) then ! ask for number of intermediate steps if NGREST > 2, but give preferred if(NGREST.eq.3) NinterDD = 3 if(NGREST.eq.4) NinterDD = 2 if(NGREST.eq.6) NinterDD = 3 if(NGREST.eq.8) NinterDD = 2 if(NGREST.eq.12) NinterDD = 1 if(bWrite3D_HR.or.bWrite3D_HRb) call AskI4('Number of temporal intermediate steps?',NinterDD) end if C C Do you want to include a magnetic map at the reference surface height? C If yes, then read the source surface map (at 2.5Rs) at C the resolution of the reference arrays. C bMagnetic = Bfield_choose_3(iVarp,cVar) call AskYN('Do you want to make magnetic source surface maps?$no',bMagnetic) if (bMagnetic) then bMagnet = .TRUE. bMagnetic = .FALSE. do J=1, B3D__NN bWildMag(J) = .FALSE. call ForeignFile(-iVarp-1,cVar,B3D__STRCMD(J),cWildMag(J)) if (itrim(cWildMag(J)) .gt. 0) then bMagnetic = .TRUE. bWildMag(J) = .TRUE. end if call ForeignFile(-iVar-1,cVar,'nso_sorr',cWildDCELIAS) if (cWildDCELIAS .eq. ' '.and..not.bForecast) then C cWildDCELIAS = '/zshare/dat/insitu/celias_????.hravg' cWildDCELIAS = '/home/soft/dat/insitu/celias_????.hravg' end if end do if (bMagnet.and..not.bMagnetic) then print *, "You requested magnetic maps but didn't specify one correct path to find them." stop end if if (bMagnetic) call AskYN('Do you want to produce magnetic field 3D files?$yes',bmagexx) if (bmagexx.and.bWrite3D_HRb) call AskYN('Do you want to duplicate the base file intervals?$yes',bmagbase) if(bmagbase) then do I=1,3 BBeginm3(I) = BBeginm3(I) EEndm3(I) = EEndm3(I) end do else if(bmagexx) then NBE3S = 1 do I=1,3 BBeginm3(I) = -1.0 EEndm3(I) = -1.0 end do call AskI4('How many intervals do you want to specify (up to 3)?',NBE3S) if(NBE3S.ne.0) then do I=1,NBE3S BBeginm3(I) = XCbeg + float(NCOFF) + xInc EEndm3(I) = XCend + float(NCOFF) - xInc if(I.eq.1) then BBeginm3(I) = 2003.57 EEndm3(I) = 2003.64 call ASKR4(' What is the beginning first limit?',BBeginm3(I)) call ASKR4(' What is the ending first limit?',EEndm3(I)) BBeginm3(I) = BBeginm3(I) - float(NCOFF) EEndm3(I) = EEndm3(I) - float(NCOFF) BBeginm3(2) = BBeginm3(I) EEndm3(2) = EEndm3(I) BBeginm3(3) = BBeginm3(I) EEndm3(3) = EEndm3(I) end if if(I.eq.2) then call ASKR4('What is the beginning second limit?',BBeginm3(I)) call ASKR4(' What is the ending second limit?',EEndm3(I)) BBeginm3(I) = BBeginm3(I) - float(NCOFF) EEndm3(I) = EEndm3(I) - float(NCOFF) BBeginm3(3) = BBeginm3(I) EEndm3(3) = EEndm3(I) end if if(I.eq.3) then call ASKR4(' What is the beginning third limit?',BBeginm3(I)) call ASKR4(' What is the ending third limit?',EEndm3(I)) BBeginm3(I) = BBeginm3(I) - float(NCOFF) BBeginm3(I) = EEndm3(I) - float(NCOFF) end if end do end if end if end if if (bMagnetic) call AskYN('Do you want to extract magnetic fields?$no',bmagex) BFACTORRR = 0.0 if(bWildMag(1).or.bWildMag(5).or.bWildMag(9).or.bWildMag(13)) BFACTORRR = 1.0 ! radial open component factor BFACTORRC = 0.0 if(bWildMag(2).or.bWildMag(6).or.bWildMag(10).or.bWildMag(14)) BFACTORRC = 1.0 ! radial closed component factor BFACTORTC = 0.0 if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.1.0) then BFACTORRC = 2.0 if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 ! tangential closed component factor end if if(BFACTORRR.EQ.1.0.and.BFACTORRC.EQ.0.0) then BFACTORRR = 2.0 if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 end if if((BFACTORRR+BFACTORRC).GE.1.0) then if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 end if if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.0.0) then if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 2.0 end if BFACTORNC = 0.0 if(bWildMag(4).or.bWildMag(8).or.bWildMag(12).or.bWildMag(16)) BFACTORNC = 2.0 if (bMagnetic.and.BFACTORRR.ne.0.0) call AskR4('Multiplicative factor for CSSS fields?',BFACTORRR) if (bMagnetic.and.BFACTORRC.ne.0.0) call AskR4('Multiplicative factor for radial closed fields?',BFACTORRC) if (bMagnetic.and.BFACTORTC.ne.0.0) call AskR4('Multiplicative factor for tangential closed fields?',BFACTORTC) if (bMagnetic.and.BFACTORNC.ne.0.0) call AskR4('Multiplicative factor for normal closed fields?',BFACTORNC) if (bMagnetic.and.(BFACTORRR.ne.0.0.or.BFACTORRC.ne.0.0.or.BFACTORTC.ne.0.0.or.BFACTORNC.ne.0.0)) & call AskYN('Adjust field factors for possible known effects?$no',bAdjust) if(bAdjust) call Adjust_BF(bWildMag,B3D__NN,iYr,Doy,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC) end if C C End of magnetic field questions C 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------- 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. print*,' ' print*,'Before ReadVIPSm8 (For either the old ifort or new PG Fortran compiler)' if (bVNago) call ReadVIPSm8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS) if (bGNago) call ReadVIPSm8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS) if (bVOoty) call ReadVIPSm8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS) if (bGOoty) call ReadVIPSm8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS) C For only old ifort compiler C if (bVNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, C & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, C & RRV,XElowV,XEhighV,iEorWV,XRlimV, C & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, C & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, C & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) C if (bVOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, C & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, C & RRV,XElowV,XEhighV,iEorWV,XRlimV, C & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, C & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, C & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav) C if (bGOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, C & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, C & RRG,XElowG,XEhighG,iEorWG,XRlimG, C & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG,SRCVGsav,SRCGsav, C & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, C & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) C if (bGNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, C & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, C & RRG,XElowG,XEhighG,iEorWG,XRlimG, C & NLmaxV,NLG,cvgfiles,SRCVG,SRCGG,SRCVGsav,SRCGsav, C & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, C & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) 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' end if print *, 'Before the in-situ reads bVACE, bVACE2, bVWind, bVCELIAS', bVACE, bVACE2, bVWind, bVCELIAS print *, 'Before the in-situ reads bDACE, bDACE2, bDWind, bDCELIAS', bDACE, bDACE2, bDWind, bDCELIAS N_INV = 0 if (bVACE) call ReadACE8_v(cWildACE,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVACE2) call ReadACE28_v(cWildACE2,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVWind) call ReadWind8_v(cWildWind,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVCELIAS.and..not.bForecast) call ReadCELIAS8_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVCELIAS.and.bForecast) call ReadCELIAS8R_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) N_ING = 0 if (bDACE) call ReadACE8_d(cWildACE,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCBA,SRCGA, & IYRSG,DOYSG8,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS) if (bDACE2) call ReadACE28_d(cWildACE2,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCBA,SRCGA, & IYRSG,DOYSG8,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS) if (bDWind) call ReadWind8_d(cWildWind,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCBA,SRCGA, & IYRSG,DOYSG8,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS) if (bDCELIAS.and..not.bForecast) call ReadCELIAS8_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCBA,SRCGA, & IYRSG,DOYSG8,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS) if (bDCELIAS.and.bForecast) call ReadCELIAS8R_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCBA,SRCGA, & IYRSG,DOYSG8,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS) call ArrR4Zero(NLmaxV,VSSIG) call ArrI4Constant(NLmaxV,1 ,NBSV ) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLmaxG,1 ,NBSG ) if(bforecast) write(*,'(A,F10.4,A,I5)') 'The Forecast cut-off time is DOY',Doy,' of',iYr print *, 'After the in-situ reads bVACE, bVACE2, bVWind, bVCELIAS', bVACE, bVACE2, bVWind, bVCELIAS print *, 'After the in-situ reads bDACE, bDACE2, bDWind, bDCELIAS', bDACE, bDACE2, bDWind, bDCELIAS Write (*,'(A,4I8)') ' After the in-situ reads NLG, NLV, N_ING, N_INV',NLG,NLV,N_ING,N_INV C call ArrR4Zero(NLmaxV,VSSIG) call ArrI4Constant(NLmaxV,1 ,NBSV ) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLmaxG,1 ,NBSG ) LASTVVALUE = 0 ALASTVDOY = 0.0 LASTViYr = 0 do I=1,NLV C if (I.lt.100) print *, I, SRCV(I), SRCGG(I) if (VOBS(I) .gt. VLIMM) NBSV(I) = 0 if (VOBS(I) .lt. VLIML) NBSV(I) = 0 if(bforecast) then ! Provides a sharp cut-off in the velocity forecast lines of sight. if(IYRSV(I).gt.iYr.or.DOYSV8(I).gt.dble(Doy).or.IYRSV(I).eq.0) then NBSV(I) = 0 end if if(NBSV(I).ne.0) then ALASTVDOY = sngl(DOYSV8(I)) LASTViYr = IYRSV(I) SRCVLAST = SRCV(I) LASTVVALUE = LASTVVALUE + 1 end if end if end do if(LASTVVALUE.ne.0) then write(*,'(A,I5,A,A,A,F9.4,A,I5)') ' The last',LASTVVALUE,' valid IPS velocity LOS ',SRCVLAST,' was',ALASTVDOY,' of',LASTViYr else print *, 'No IPS velocity LOS values were used' end if ICntBadr4 = 0 LASTGVALUE = 0 ALASTGDOY = 0.0 LASTGiYr = 0 do I=1,NLG if(.not.bForecast) then if(.not.bSMEI) then if(GGOBS(I) .eq. BadR4()) NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 if(NBSG(I).ne.0) LASTGVALUE = LASTGVALUE + 1 end if if(bSMEI) then if (GOBS(I) .ne. BadR4()) GOBS(I) = GOBS(I)*GLIMF if (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) .gt. GLIMP) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if if (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) .lt. GLIMM) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if if(NBSG(I).ne.0) LASTGVALUE = LASTGVALUE + 1 end if end if 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 if(bForecast) then if(.not.bSMEI) then ! Provides a sharp cut-off in the g-level forecast lines of sight. C print *, I, NLmaxG, IYRSGG(I), iYr, DOYSGG8(I), dble(Doy), LASTGVALUE if(IYRSGG(I).gt.iYr.or.DOYSGG8(I).gt.dble(Doy).or.IYRSGG(I).eq.0) then NBSG(I) = 0 end if if(NBSG(I).ne.0) then ALASTGDOY = sngl(DOYSGG8(I)) LASTGiYr = IYRSGG(I) SRCGLAST = SRCGG(I) LASTGVALUE = LASTGVALUE + 1 end if end if if(bSMEI) then if(IYRSG(I).gt.iYr.or.DOYSG8(I).ge.dble(Doy).or.IYRSG(I).eq.0) then NBSG(I) = 0 end if if(NBSG(I).ne.0) then ! Limit bad sources count to those within the forecast interval if (GOBS(I) .ne. BadR4()) GOBS(I) = GOBS(I)*GLIMF if (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) .gt. GLIMP) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if if (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) .lt. GLIMM) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if end if if(NBSG(I).ne.0) then ALASTGDOY = sngl(DOYSG8(I)) LASTGiYr = IYRSG(I) SRCBLAST = SRCB(I) LASTGVALUE = LASTGVALUE + 1 end if end if end if end do if(.not.bSMEI) then if(LASTGVALUE.ne.0) then write(*,'(A,I5,A,A,A,F9.4,A,I5)') ' The last',LASTGVALUE,' valid IPS g-levely LOS ',SRCGLAST,' was',ALASTGDOY,' of',LASTGiYr else print *, 'No IPS g-level LOS values were used' end if end if if(bSMEI) then if(LASTGVALUE.ne.0) then write(*,'(A,I7,A,A,A,F9.4,I5)') ' The last',LASTGVALUE,' SMEI LOS ',SRCBLAST,' was on',ALASTGDOY,LASTGiYr else print *, 'No SMEI LOS values were used' end if end if if(ICntBadr4.gt.0) write(*,'(A,I7)') ' The number of bad sources are = ', ICntBadr4 if(bVFACE) then print *,' ' print *, 'The 1 AU upper limit of in-situ velocity is:',VVMAX print *, 'The 1 AU lower limit of in-situ velocity is:',VVMIN end if print *, 'Doy',Doy LASTAVALUE = 0 ALASTADOY = 0.0 LASTAiYr = 0 do I=NLV+1,N_INV if(bforecast) then ! Provides a sharp cut-off in the ACE forecast measurements. if(IYRSV(I).gt.iYr.or.DOYSV8(I).gt.dble(Doy)) then NBSV(I) = 0 C VOBS(I) = BadR4() end if if (VOBS(I) .gt. VVMAX) NBSV(I) = 0 if (VOBS(I) .lt. VVMIN) NBSV(I) = 0 if(NBSV(I).ne.0) then ALASTVDOY = sngl(DOYSV8(I)) LASTViYr = IYRSV(I) SRCALAST = SRCV(I) LASTAVALUE = LASTAVALUE + 1 end if end if end do if(LASTAVALUE.ne.0) then write(*,'(A,I5,A,A,A,F9.4,A,I5)') & ' The last',LASTAVALUE,' valid in-situ forecast velocity value ',SRCALAST,' was',ALASTVDOY,' of',LASTViYr end if if(bFACE) then print *,' ' print *, 'The 1 AU upper limit of in-situ density is:',DMAX print *, 'The 1 AU lower limit of in-situ density is:',DMIN end if LASTAVALUE = 0 ALASTADOY = 0.0 LASTAiYr = 0 do I=NLG+1,N_ING if(bforecast) then ! Provides a sharp cut-off in the ACE forecast measurements. if(IYRSG(I).gt.iYr.or.DOYSG8(I).gt.dble(Doy)) then NBSG(I) = 0 C GOBS(I) = BadR4() end if if (GOBS(I) .lt. DMIN) NBSG(I) = 0 if (GOBS(I) .gt. DMAX) NBSG(I) = 0 if(NBSG(I).ne.0) then ALASTADOY = sngl(DOYSG8(I)) LASTAiYr = IYRSG(I) SRCALAST = SRCGA(I-NLG) LASTAVALUE = LASTAVALUE + 1 end if end if end do if(LASTAVALUE.ne.0) then write(*,'(A,I5,A,A,A,F9.4,A,I5)') ' The last',LASTAVALUE,' valid forecast in-situ density value ',SRCALAST,' was',ALASTADOY,' of',LASTAiYr else print *, 'No in-situ density values were used' end if C do I=307700,307710 C print *, 'In smeiipstd0n_inp', I, GOBS(I),NBSG(I) C end do 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 C print *, 'Bad GOBS(I) xTemp, xTHOM', I, GOBS(I), xTemp, xTHOM ZZ(I) = BadR4() end if end do C end if C print *, 'After read_smei8', NLG C call ArrR4Zero(NLmaxV,VSSIG) call ArrR4Zero(NLmaxG,GSSIG) do I=1,NLmaxV C if (I.lt.100) print *, I, SRCV(I) if (VOBS(I) .gt. VVMAX) NBSV(I) = 0 if (VOBS(I) .lt. VVMIN) NBSV(I) = 0 end do ICntBadr4 = 0 do I=1,NLG if(.not.bSMEI) then if (GGOBS(I) .eq. BadR4()) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if end if if(bSMEI) then 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 if (GOBS(I) .eq. BadR4().or.NBSG(I).eq.0) then NBSG(I) = 0 ICntBadr4 = ICntBadR4 + 1 end if end if end do if(ICntBadr4.gt.0.and.bsmei) print *, ' The number of bad SMEI brightness sources are = ', ICntBadr4 if(ICntBadr4.gt.0.and..not.bsmei) print *, ' The number of bad g-level 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 = nint((float(NTV) + (2.0*float(NVREST)+1.0))/2.0) + 5*NVREST + 1 NTG = nint((float(NTG) + (2.0*float(NGREST)+1.0))/2.0) + 5*NGREST + 1 end if ItoffG = 0 ItoffV = 0 print *, ' ' print *, 'Before Mktimes', nCar, JDCar(1), JDCar(500) print *, 'Before MkTimes MJDcntr =', MJDcntr print *, ' ' ItoffG = 0 ItoffV = 0 IdaybegV = 365 IdaybegG = 365 IYRBV = IYRSV(1) IYRBG = IYRSG(1) NTmaxVmkt = NTmaxV NTmaxGmkt = NTmaxG if(bForecast) then NTmaxVmkt = 0 NTmaxGmkt = 0 end if print *, ' ' print *, 'Before Mktimes', nCar, JDCar(1), JDCar(500) print *, 'Before MkTimes MJDcntr =', NTV, NTG, MJDcntr print *, 'Before MkTimes NTmaxG =', NTmaxG, NTmaxGmkt print *, ' ' if(bVcon) & call MkTimes(bForecast,MJDcntr,NmidHRV,aNdayV,NTmaxV,ALng,nLng,nCar,JDCar, & NTmaxVmkt,NTV,IYRBV,IdaybegV,XCintDV,XCintV,XCbeV,XCtbegV,XCtendV,IYRV,DOYV) if(NLG.eq.0) bGcon = .FALSE. if(bGcon.or.bTcon) & call MkTimes(bForecast,MJDcntr,NmidHRG,aNdayG,NTmaxG,ALng,nLng,nCar,JDCar, & NTmaxGmkt,NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG) print *, 'After MkTimes NTmaxG =', NTV, NTG, NTmaxG, NTmaxGmkt write (*,'(A,4I5)') ' IYRBV, IYRBG, IdaybegV, IdaybegG',IYRBV,IYRBG,IdaybegV,IdaybegG if(.not.bGcon.and..not.bTcon) then nTG = nTV if(nTG.lt.nTV) then print *, ' nTG (nTmaxG) must be larger than nTV (nTmaxV)' print *, ' if density is not deconvolved.' stop end if 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 print*, 'bVcon = ', bVcon if(.not.bVcon) then C nTV = nTG ! set here in case only density is deconvolved. Otherwise comment out. BVJ 4/30/2024 if(nTV.ne.nTG) then nTGdnTV = nTG/nTV if(nTG.ne.(nTV*nTGdnTV - 1)) then print *, ' The current density programming does not support non-multiple' print *, ' temporal resolution differences between density and velocity' print *, ' with density the higher resolution, if velocity is not deconvolved.' stop end if end if 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 C nTV = nTG if(nTV.ne.nTG) then nTGdnTV = nTG/nTV if(nTG.ne.(nTV*nTGdnTV - 1)) then print *, ' The current single mkshift programming does not support non-multiple' print *, ' temporal resolution differences between density and velocity' print *, ' with density the higher resolution, if there is only one' print *, ' shift matrix per iteration used.' stop end if end if ItoffV = ItoffG XCtbegV = XCtbegG XCtendV = XCtendG do I=1,nTmaxV ! The temporal spaces for XCintV, etc. better be the same as for XCintG XCintDV(I) = XCintDG(I) XCintV(I) = XCintG(I) XCbeV(1,I) = XCbegG(1,I) XCbeV(2,I) = XCbegG(2,I) end do end if write(*,'(A,2F12.5,I12)') ' XCtbegV, XCtendV, NTV', XCtbegV, XCtendV, NTV write(*,'(A,2F12.5,I12)') ' XCtbegG, XCtendG, NTG', XCtbegG, XCtendG, NTG write(*,'(/A)') ' XCbeV =' write(*,'(8F9.6)') XCbeV write(*,'(/A)') ' XCbeGG =' write(*,'(8F9.6)') XCbeGG write(*,'(/A)') ' XCintV =' write(*,'(8F9.6)') XCintV write(*,'(/A)') ' XCintG =' write(*,'(8F9.6)') XCintG write(*,'(/A)') ' XCintDV =' write(*,'(8F9.4)') XCintDV write(*,'(/A)') ' XCintDG =' write(*,'(8F9.4)') XCintDG C C Find the minimum and maximum time index value in the timeseries C 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 write(*,'(/A,2I5/)') ' For the density base subroutine nTminTS, nTmaxTS',nTminTS, nTmaxTS NLVVV = 0 do I=1,NLV if (VOBS(I) .ge. VLIML .and. VOBS(I) .le. VLIMM) NLVVV = NLVVV + 1 end do if(NLVVV.lt.NLVmin) then bVcon = .FALSE. write(*,'(A,I3,A)') ' There are too few (',NLVVV,') velocity measurements. NLV set to zero.' print *, ' ' NLV = 0 NLVVV = 0 end if if(NLVVV.lt.200.and.NLVVV.ne.0) then ! If the velocity source numbers are few make this constant bigger CONSTV = CONSTV*float(200)/float(NLVVV) ! Hopefully this fills in velocities from other times if the gap print *, 'NLVVV, CONSTV', NLVVV, CONSTV ! is not to big! BVJ 5/3/11 end if if(bGcon.or.bTcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP) if(bVcon) call ArrR4Constant(nLng*nLat*nTV,BadR4(),VMAP) NLmax = max(NLmaxG,NLmaxV) NLVBEG = 0 NLGBEG = 0 NLVEND = NLV NLGEND = NLG 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.NLVEND+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.NLVEND) K = K - 1 NTVD = NTVE-NTVB+1 if(NTTV.eq.1) NLVBEG = NTVD write (*,'(A,I4,A,2I6,2F10.3)') 'V loop ',NTVV,' - S#, K, XCEV, XCintV', NTVD,K,NCoff+XCEV(K),NCoff+XCintV(NTVV+2) if(NTVV.eq.NTV) NLVEND = NTVE end if end if if(IVSW.eq.1) then C print *, ' V switch ',IVSW,' loop',NTVV if (NLVEND .eq. 0) call Say('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 ArrR4Constant(nLngLat,Speed,VMAP(1,1,NTVV)) C call MapGrid(XCbeV(1,NTVV),XCbeV(2,NTVV),YLbeg,YLend,NC, C & NTVD,XCV(NTVB),YLV(NTVB),VOBS(NTVB),NAV,DAV,nLng,nLat, C & VMAP(1,1,NTVV),VDEV(1,1,NTVV),VPNT(1,1,NTVV), C & 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 MkLOSWeightsm(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) C 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.NLGEND+1) then ! If K is more than the source # +1 don't turn on switch if(K.gt.0) then if (bGcon.or.bTcon.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.NLGEND) K = K - 1 NTGD = NTGE-NTGB+1 if(NTGG.eq.1) NLGBEG = NTGD 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,2I8,2F10.3)') & 'B loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2) if(NTGG.eq.NTG) NLGEND = NTGE end if end if end if if(IGSW.eq.1) then if(NTGB.ne.0) then if (NLGEND .eq. 0) call Say('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.bSMEI) then C call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC, C & NTGD,XCG(NTGB),YLG(NTGB),GOBS(NTGB),NAV,DAV,nLng,nLat, C & DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG), C & Gmin,Gmax,GDmin,GDmax) C else C call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC, C & NTGD,XCG(NTGB),YLG(NTGB),ZZ(NTGB),NAV,DAV,nLng,nLat, C & DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG), C & Gmin,Gmax,GDmin,GDmax) C end if TD2 = 1.0 ! Initial total Gvalue for g-level maps TD2 = 0.66667*DEN1AU*((R1AU/RRD)**FALLOFFD) ! divide total density by two-thirds for initial Thomson maps call ArrR4Constant(nLngLat,TD2,DMAP(1,1,NTGG)) end if C Calculate weights, WTSG, according to NLOSWG along all lines of sight. call MkLOSWeightsm(NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) C call MkLOSWeights (NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) C 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 C 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() aminallV = BadR4() do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) if(amax.le.VLIMM.or.amax.ne.BadR4()) amaxallV = amax if(amin.ge.VLIML.or.amin.ne.BadR4()) aminallV = amin end do if(amaxallV.eq.BadR4().or.aminallV.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)**FALLOFFD) ! 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 = 0 NLGGL = 0 do L=1,NLG if(L.lt.NLGBEG.or.L.gt.NLGEND) then GOBS(L) = BadR4() NBSG(L) = 0 end if if(GOBS(L) .ne. BadR4()) then NLGG = NLGG + 1 NLGGL = L end if end do NLVV = 0 NLVVL = 0 do L=1,NLV if(L.lt.NLVBEG.or.L.gt.NLVEND) then VOBS(L) = BadR4() NBSV(L) = 0 end if if(VOBS(L) .ge. VLIML .and. VOBS(L) .le. VLIMM .and. WTSV(1,L).ne.0.0.and.NBSV(L).ne.0) then NLVV = NLVV + 1 NLVVL = L else NBSV(L) = 0 end if end do NLGGG = 0 NLGGGL = 0 do L=NLG+1,NLG+N_ING if(GOBS(L) .ne. BadR4().and.NBSG(L).ne.0) then NLGGG = NLGGG + 1 NLGGGL = L end if end do NLVVV = 0 NLVVVL = 0 do L=NLV+1,NLV+N_INV if(VOBS(L) .ne. BadR4().and.NBSV(L).ne.0) then NLVVV = NLVVV + 1 NLVVVL = L end if end do if(NLVV.lt.NLVMIN.and.NLVVV.lt.NLVMIN) then bVCon = .FALSE. else bVCon = .TRUE. end if do I=1,N_INV if(bVACE) then C WTSV_in(I) = float(1*NLVV) + 10.0 WTSV_in(I) = 10.0 else NBSV(I+NLV) = 0 end if end do C print *, 'bVCon =',bVCon if(bFACE) then do I=1,N_ING C WTSD_in(I) = float(1*NLGG) + 100.0 WTSD_in(I) = 10000.0 if(.not.bDACEn) then NBSG(I+NLG) = 0 end if end do end if C **** if(bTcon) then print *, ' ' write (*,'(1X,A)') 'The last good of 10 SMEI brightness values were:' NLGM10 = NLGGL-10 do I=NLGM10,NLGGL if(NBSG(I).ne.0) 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 if(bFACE) then print *, ' ' write (*,'(1X,A)') 'The last good in-situ density value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3,A,F8.3)') & 'I =',NLGGGL,' Year =',IYRSG(NLGGGL),' DOY =',DOYSG8(NLGGGL),' In-situ density =',GOBS(NLGGGL) end if end if C **** if(bVcon.and.NLVVL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last good IPS velocity value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3,A,F8.3)') & 'I =',NLVVL,' Year =',IYRSV(NLVVL),' DOY =',DOYSV8(NLVVL),' In-situ density =',VOBS(NLVVL) if(bVACE.and.NLVVVL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last good in-situ velocity value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3,A,F8.3)') & 'I =',NLVVVL,' Year =',IYRSV(NLVVVL),' DOY =',DOYSV8(NLVVVL),' In-situ density =',VOBS(NLVVVL) end if end if if(bTcon) then print *, ' ' call Say('smeiipstd0n_inp','I','Info','# good B observations ' & //cStr(:Int2Str(NLGG,cStr))) else call Say('smeiipstd0n_inp','I','Info','# good G observations ' & //cStr(:Int2Str(NLGG,cStr))) end if call Say('smeiipstd0n_inp','I','Info','# good V observations ' & //cStr(:Int2Str(NLVV,cStr))) if(NLVV.le.NLVmin) then bVcon = .FALSE. write(*,'(A,I5,A)') 'NLVV =',NLVV,', too few IPS values for a velocity deconvolution.' end if call Say('smeiipstd0n_inp','I','Info','# good in-situ D observations ' & //cStr(:Int2Str(NLGGG,cStr))) call Say('smeiipstd0n_inp','I','Info','# good in-situ V observations ' & //cStr(:Int2Str(NLVVV,cStr))) ffact = 15./48. ! about fifteen of 48 days total are being reconstructed redundant = float(NLG)/((nLng/Alng)*(nLat)*(nTG*ffact)) print *, ' ' write(*,'(A,F10.3,A)'), 'The D obs. are',redundant,' times more redundant than the digital res. requires.' print *, ' ' 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 DEN1AUO = DEN1AU DEN1AUN = DEN1AU BASEOBSO = 0.0 OBSBASEO = 0.0 DB = DEN1AU*((R1AU/RRD)**FALLOFFD) 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) C print *, 'DMAPV min max N',amin,amax,N call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,0.0) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,0.0) C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N',amin,amax,N C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,nTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP min max N',amin,amax,N end do end if print*, 'Before preliminary V MkShiftdn0n, bVcon =',bVcon if(bVcon) then ! BVJ added 4/30/2024 print*, 'bVcon =', bVcon, 'Go into preliminary V MkShiftdn0n' call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) end if 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,FALLOFFD,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_in(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM, & N_INV,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in) C call MkPostd(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, C & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8, C & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) if(NitV.eq.NiterT.and.bWritelos) call writelosmapvg_in(0,NLOSVP1,NLV,N_INV,nTV,NBSV,NCOFF,ITIM, & ITIM_in,XLON,XLON_in,XLAT,XLAT_in) C 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 Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DVfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,VFAC_in) call MkVModeltd_in(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX, & WTSV,DFAC,VFAC,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSV_in,DFAC_in,VFAC_in,PWV, & NLV,N_INV,NLOSV,NLOSVP1,VMAP,VLIM, & DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij,VWTij_in) C call MkVModeltd(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX, C & WTSV,DFAC,VFAC,PWV,NLV,NLOSV,NLOSVP1,VMAP,VLIM, C & DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij) do I=1,NLV+N_INV if(WTSV(1,I).eq.0.and.I.le.NLV) then NBSV(I) = 0 VMOD(I) = Speed end if if(WTSV_in(I-N_INV).eq.0.and.I.gt.NLV) then NBSV(I) = 0 VMOD(I) = Speed end if end do call FixModeltdn(1,NitV,NLV+N_INV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) C call FixModeltdn(1,NLV+N_INV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) C call FixModeltdn(1,NLV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) C if(bWriteMS.and.NitV.eq.NiterT) then call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV+N_INV,PWV,VRATIO, & NLmaxV,NSRV,SRCV,DOYSV8,XCEV,VOBS,VMOD,XLON,XLproj,XCtpr,NLOSVP1,FIX,VSSIG,NBSV) C call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV,PWV,VRATIO, C & 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_in(LON,LAT,ITIM,RPX,LON_in,LAT_in,ITIM_in,RPX_in,XEV,VWTij,VWTij_in, & FIX,NLV,N_INV,NLOSV,NLOSVP1, & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VSIGB,NVS,NBSV,VMAPSIG, & ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,WTPV_in, & Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp) C print *, 'Out of MkVmaptdn0n_in' C call MkVMaptdn0n(LON,LAT,ITIM,RPX,XEV,VWTij,FIX,NLV,NLOSV,NLOSVP1,XLON,XLAT, C & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VMOD,VSIGB,NVS,NBSV,VMAPSIG, C & ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV, C & 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 mkshiftdn0n. C if (.not.bDproxy.or.Nit.ge.1) then if (bGcon .or. bTcon .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 MkShiftdn0n DMAP min max N',amin,amax,N AEmax = AEmaxD/(RRS*RRS) AEmin = AEminD/(RRS*RRS) call FillMapOLS(2,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,Speed,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD,3,0.0,0.0) C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before 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 Mkshiftdn0n' call stopwatch(' ',' ') end if if (bGcon.or.bTcon) then C C Condition the next velocity map set for density (VMAPD) to use in mkshiftdn0n. 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 MkShiftdn0n VMAPD 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,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 MkShiftdn0n VMAPD min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTG,VMAPD,0,1.0*CONDT/aNdayG,0.,Dtmp) do N=1,NTG call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) end do end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After arranging VMAPD and before D Mkshiftdn0n' call stopwatch(' ',' ') end if if(bGcon.or.bTcon) then C 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) C else call MkShiftdn0n(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG,Speed, & nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After D Mkshiftdn0n' end if if (bGcon.or.bTcon) 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,FALLOFFD,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,FALLOFFD,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_in(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM, & N_ING,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in) C call MkPostd(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, C & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8, C & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) if(NitG.eq.NiterT.and.bWritelos) call writelosmapvg_in(1,NLOSGP1,NLG,N_ING,nTG,NBSG,NCOFF, & ITIM,ITIM_in,XLON,XLON_in,XLAT,XLAT_in) C if(NitG.eq.NiterT.and.bWritelos) call writelosmapvg(1,NLOSGP1,NLG,nTG,NBSG,NCOFF,ITIM,XLON,XLAT) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC) C write(*,'(A,5I6,I9)') 'Get4dval3 nLng, nLat, nMap, nT, NLOSP1, NL',nLng, nLat, nMap, nTG, NLOSGP1, NLG C call Get4Dval3(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, C & DDfact,NLG,NLOSGP1,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,N_ING,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in) C call Get4Dval3(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, C & DDfact,N_ING,1,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in) if(.not.bTcon) then ! G-level model-maker from density call MkGModeltd_in(XCbegG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,PWG,DEN1AU,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSD_in,DFAC_in, & NLG,N_ING,NLOSG,NLOSGP1,DMAP, & nLng,nLat,nTG,RRS,GMOD2,GWTij,GWTij_in,WTPD,GMWTi) C call MkGModeltdn(XCbegG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, C & WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP, C & nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi) NGval = 2 ! G-level fix else ! Brightness model-maker from density C C If in-situ density is available define a new DEN1AU called DEN1AUN if bFACE is yes C DEN1 = 0.0 ADEN1 = 0.0 if(bFACE) then do I=NLG+1,NLG+N_ING if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then ! Observed density DEN1 = DEN1 + GOBS(I) ADEN1 = ADEN1 + 1.0 end if end do DEN1AUO = DEN1AUN if(ADEN1 .ne.0.0) DEN1AUN = DEN1/ADEN1 ! The mean in-situ density value DEN1AUMN = DEN1AUN if(abs(DEN1AUN-DEN1AUO).gt.0.000001) then print *, ' ' write(*,'(A,F9.5,A,F9.5,A,I5,A)') & 'The previous mean in-situ density',DEN1AUO,' is now',DEN1AUN,' using',nint(ADEN1),' values' DEN1 = 0.0 ADEN1 = 0.0 do I=NLG+1,NLG+N_ING if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then ! Observed density DEN1 = DEN1 + (GOBS(I) - DEN1AUN)**2 ADEN1 = ADEN1 + 1.0 end if end do if(ADEN1 .ne.0.0) DENSQ = SQRT(DEN1/ADEN1) ! The mean in-situ density value write(*,'(A,F9.5,A,I5,A,F9.5,A)') & 'The RMS density is',DENSQ,' using',nint(ADEN1),' values, or',DEN1AUN-DENSQ,' from zero' end if C bBFIX = .FALSE. if(bBFIX) then ! If the base density needs to be modified to match observations if(Nit.eq.0.or.Nit.eq.(LimitoG+1)) then print *, 'base density set',Nit,LimitoG+1 BASEOBS = 0.0 ABASEOBS = 0.0 do I=1,NLG if(NBSG(I).ne.0) then C BASEOBS = BASEOBS + GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) C if(I.ge.2000.and.I.le.2020) then C print *, 'I, XEG, GOBS, factor', I, XEG(I),GOBS(I),((abs(XEG(I))**2.5)/(90.0**2.5)) C end if BASEOBS = BASEOBS + GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) ABASEOBS = ABASEOBS + 1.0 end if end do if(ABASEOBS.ne.0.0) BASEOBS = BASEOBS/ABASEOBS write(*,'(A,F8.5,A,I8,A)') & 'This indicates a brightness change to',BASEOBS,' from ',nint(ABASEOBS),' values' NEEDRMSF = 0 ! Throw out the below if(NEEDRMSF.eq.1) then if(abs(BASEOBS-BASEOBSO).gt.0.000001) then BASEOBSS = 0.0 ABASEOBSS = 0.0 BLOWESTT = 0.0 BLOWEST = -999999.0 ABLOWEST = 0.0 do I=1,NLG if(NBSG(I).ne.0) then BASEOBSS = BASEOBSS + (GOBS(I)*((abs(XEG(I))**2.5)/(90.0**2.5)) - BASEOBS)**2 ABASEOBSS = ABASEOBSS + 1.0 end if if(ABLOWEST.lt.10.0) then OBSLOWESTL = 999999.0 do J=1,NLG if(NBSG(J).ne.0) then OBSLOWEST = GOBS(J)*((abs(XEG(J))**2.5)/(90.0**2.5)) if(OBSLOWEST .gt. BLOWEST .and. OBSLOWEST .lt. OBSLOWESTL) then OBSLOWESTL = OBSLOWEST end if end if end do BLOWEST = OBSLOWESTL if(BLOWEST.gt.-999999.0) then ABLOWEST = ABLOWEST + 1.0 C print *, 'Lowest value', ABLOWEST ,BLOWEST BLOWESTT = BLOWESTT + BLOWEST end if end if end do if(ABLOWEST.ne.0) BLOWESTT = BLOWESTT/ABLOWEST if(ABASEOBSS .ne.0.0) BASEOBSS = SQRT(BASEOBSS/ABASEOBSS) ! The mean in-situ density value write (*,'(A,F10.5,A,I8,A)') & 'The RMS brightness is',BASEOBSS,' using',nint(ABASEOBSS),' values' write (*,'(A,F9.5,A,I4,A)') & 'The average lowest brightness is',BLOWESTT,' from',nint(ABLOWEST),' values' DENLOWEST = (DEN1AUN - abs(BLOWESTT/BASEOBSS)*DENSQ) write (*,'(A,F9.5,A)') & 'This implies a minimum brightness',abs(BLOWESTT/BASEOBSS),' times lower than the RMS value' write (*,'(A,F9.5,A)') ' or a lowest density of', DENLOWEST,'. Since no brightness should be 2 sigma -zero,' DEN1AU = (DEN1AUN - DENLOWEST)*2.0 write (*,'(A,F9.5)') ' this brightness variation implies a 1 AU mean density of at least', DEN1AU end if else DEN1AU = DEN1AUN end if write (*,'(A,F9.5)') ' The 1 AU mean density has now been set to:', DEN1AU call Mk_Base_in(bFACE,VMAPD,DMAP,nLng,nLat,nTG,nTminTS,nTmaxTS,RRV,DEN1AU,FALLOFFD,VEL1AU, & BadR4(),BadR4(),DBASE) end if end if end if call MkDHModeltd_in(XCbegG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,FALLOFFD,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSD_in,DFAC_in, & NBSG,NLG,N_ING,NLOSG,DISTG,NLOSGP1,DMAP,DBASE,BBASE, & nLng,nLat,nTG,RRS,GMOD2,GWTij,GWTij_in) C call MkDHModeltd(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, C & WTSG,DFAC,FALLOFFD,NLG,NLOSG,DISTG,NLOSGP1,DMAP,DBASE,BBASE, C & nLng,nLat,nTG,RRS,GMOD2,GWTij) NGVal = 3 ! Thomson scattering fix end if IBASEOBSO = 1 IMODBRIGHT = 0 IMODBRIGHT2 = 0 IMODBRIGHTB = 0 IMODBRIGHTB2 = 0 ipcheck = 0 if(ipcheck.eq.1) then print *,' ' PRINT *,' NLG =',NLG,' N_ING =',N_ING end if do I=1,NLG+N_ING if(.not.bTcon) then ! Non-SMEI part if(I.le.NLG) then GOBSS(I) = GOBS(I) else GOBSS(I) = (GOBS(I)/(DEN1AU*((1.0/RPX_in(I-NLG))**2)))**(2.0*PWG) ! Relative to ACE densities C print *, GOBSS(I),GOBS(I),DEN1AU,RPX_in(I-NLG),PWG end if end if ! SMEI part if(IBASEOBSO.eq.1) then ! Old version if(bTcon) then if(I.le.NLG) then if(WTSG(1,I).ne.0.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBS(I) + BBASE(I) ! Purpose of the B base is to up each LOS by mean end if else if(WTSD_in(I-NLG).ne.0.0.and.NBSG(I).ne.0) then GOBSS(I) = WTSD_in(I-NLG)*GOBS(I) + BBASE(I) ! Relative to weighted modeled ACE densities if(GOBSS(I).lt.0.0) GOBSS(I) = 0.001 ! No negative (or zero) in-situ densities allowed end if C if(WTSD_in(I-NLG).ne.0.0) then C GMODD_in(I-NLG) = (GMOD2(I) + BBASE(I))/WTSD_in(I-NLG) C end if end if end if end if C if(IBASEOBSO.eq.0) then ! If IBASEOBSO is 1, then don't do this C if(bTcon) then ! Thomson scattering 3D reconstruction is performed C if(I.le.NLG) then C if(WTSG(1,I).ne.0.0.and.NBSG(I).ne.0) then C GOBSS(I) = GOBSS(I) - BASEOBS ! Make sure the observed brightness has the interval mean C end if C else C if(WTSD_in(I-NLG).ne.0.0.and.NBSG(I).ne.0) then C GOBSS(I) = WTSD_in(I-NLG)*GOBS(I) + BBASE(I) ! to compare with the weighted modeled in-situ densities C end if C end if C end if C end if if(I.le.NLG) then ! This is simply a counter and setting of bad values if(WTSG(1,I).eq.0.0) then ! via NBSG(I) if WTSG is zero. NBSG(I) = 0 GMOD2(I) = 1.0 IMODBRIGHTB = IMODBRIGHTB + 1 else ! Put in 5/16/11 BVJ This is the total model bright GMOD2(I) = GMOD2(I) if(GMOD2(I).lt.0.0) then GMOD2(I) = 0.001 ! No negative (or zero) model brightness allowed IMODBRIGHT = IMODBRIGHT + 1 end if end if else if(WTSD_in(I-NLG).eq.0.0) then NBSG(I) = 0 GMOD2(I) = 1.0 IMODBRIGHTB2 = IMODBRIGHTB2 + 1 end if if(NBSG(I).ne.0) then GMOD2(I) = GMOD2(I) - BBASE(I) if(GMOD2(I).lt.0.0) then ! No negative (or zero) in-situ model densities GMOD2(I) = 0.001 IMODBRIGHT2 = IMODBRIGHT2 + 1 end if end if end if end do if(IMODBRIGHTB.ne.0) print *, IMODBRIGHTB,' model brightness weights are zero after MkDHModeltd_in.' if(IMODBRIGHTB2.ne.0) print *,IMODBRIGHTB2,' model in-situ weights are zero after MkDHModeltd_in.' if(IMODBRIGHT.ne.0) print *, IMODBRIGHT,' model brightness values are below zero after MkDHModeltd_in.' if(IMODBRIGHT2.ne.0) print *, IMODBRIGHT2,' model in-situ values are below zero after MkDHModeltd_in.' IOBSBRIGHT = 0 IOBSBRIGHT2 = 0 do I=1,NLG+N_ING if(IBASEOBSO.eq.1) then ! Old version if(NBSG(I).ne.0) then if(I.le.NLG) then C GGOBS(I) = GOBSS(I) + BBASE(I) ! Each iteration, use modified observed brightnesses C GGOBS(I) = GOBSS(I) - BASEOBS ! Put in 5/16/11 BVJ Makes sure the mean brightness during the interval is the mean GGOBS(I) = GOBSS(I) ! Remove BASEOBS change BVJ 4/29/2020 if(GGOBS(I).lt.0.0) then GGOBS(I) = 0.001 ! No negative (or zero) brightness values are allowed end if else GGOBS(I) = GOBSS(I) + DEN1AU*WTSD_in(I-NLG) ! The observed "brightness" at in-situ location of total density C GGOBS(I) = GOBSS(I) end if ! to compare with modeled brightness (of total density) end if end if if(NBSG(I).ne.0) then if(I.le.NLG) then GGOBS(I) = GOBSS(I) + BBASE(I) ! Adds the modeled mean density brightness to each LOS to increase the observed value if(NLmaxG.gt.1000000.and.ipcheck.eq.1) then if(I.eq.(100000)) print *,' Brightness # = I BBASE(I) GOBSS(I) GGOBS(I) GMOD2(I)' if(I.eq.(2500000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(2000000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(1500000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(1000000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-25000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-20000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-15000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-10000)) write(*,'(A,I12,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) end if if(NLmaxG.lt.1000000.and.ipcheck.eq.1) then if(I.eq.(100000)) print *,' Brightness # = I BBASE(I) GOBSS(I) GGOBS(I) GMOD2(I)' if(I.eq.(250000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(200000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(150000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(100000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-2500)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-2000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-1500)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-1000)) write(*,'(A,I8,F11.5,F11.5,2F11.5)') & ' Brightness value', I, BBASE(I), GOBSS(I), GGOBS(I), GMOD2(I) end if if(GGOBS(I).lt.0.0) then GGOBS(I) = 0.001 ! No negative (or zero) brightness values are allowed IOBSBRIGHT = IOBSBRIGHT + 1 end if else C GGOBS(I) = GOBSS(I) + BBASE(I) ! Adds the modeled mean in-situ density to each in-situ observed value C GGOBS(I) = GOBSS(I) + BBASE(I) - DEN1AU*WTSD_in(I-NLG) GGOBS(I) = GOBSS(I) if(I.eq.NLG+200.and.ipcheck.eq.1) & write(*,'(A,I8,4F13.3)') ' In-situ value ',I,BBASE(I),DEN1AU*WTSD_in(I-NLG),GGOBS(I),GMOD2(I) if(I.eq.NLG+200.and.ipcheck.eq.1) & write(*,'(A,I8,I11,3F13.3)')' In-situ original',I,NBSG(I),WTSD_in(I-NLG), GOBS(I),GOBSS(I) if(I.eq.NLG+400.and.ipcheck.eq.1) & write(*,'(A,I8,4F13.3)') ' In-situ value ',I,BBASE(I),DEN1AU*WTSD_in(I-NLG),GGOBS(I),GMOD2(I) if(I.eq.NLG+400.and.ipcheck.eq.1) & write(*,'(A,I8,I11,3F13.3)')' In-situ original',I,NBSG(I),WTSD_in(I-NLG), GOBS(I),GOBSS(I) if(GGOBS(I).lt.0.0) then GGOBS(I) = 0.001 ! No negative (or zero) in-situ densities are allowed IOBSBRIGHT2 = IOBSBRIGHT2 + 1 end if end if end if end do print *,' ' C if(IOBSBRIGHT.ne.0) write(*,'(I8,A)') IOBSBRIGHT, ' observed brightness values are below zero before FixModeltdn' C if(IOBSBRIGHT2.ne.0) write(*,'(I8,A)') IOBSBRIGHT2,' observed in-situ values are below zero before FixModeltdn' call FixModeltdn(NGVal,NitG,NLG+N_ING,GGOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) C call FixModeltdn(NGVal,NLG+N_ING,GGOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) C call FixModeltdn(NGVal,NLG,GGOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) C if(ipcheck.eq.1) then print *,' ' PRINT *,' NLG =',NLG end if do I=1,NLG+N_ING if(NBSG(I).ne.0) then if(I.le.NLG) then GGOBS(I) = GOBSS(I) + BBASE(I) ! Adds the modeled mean density brightness to each LOS to increase the observed value if(NLmaxG.gt.1000000.and.ipcheck.eq.1) then if(I.eq.(100000)) print *,'Brightness Fix I NBSG(I) FIX(I) GGOBS(I) GMOD2(I)' if(I.eq.(2500000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(2000000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(1500000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(1000000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-25000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-20000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-15000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-10000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) end if if(NLmaxG.lt.1000000.and.ipcheck.eq.1) then if(I.eq.(100000)) print *,'Brightness Fix I NBSG(I) FIX(I) GGOBS(I) GMOD2(I)' if(I.eq.(250000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(200000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(150000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(100000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-2500)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-2000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-1500)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) if(I.eq.(NLG-1000)) write(*,'(A,I8,I5,3F13.5)') & ' Brightness Fix', I, NBSG(I), FIX(I), GGOBS(I), GMOD2(I) end if else if(I.eq.NLG+200.and.ipcheck.eq.1) write(*,'(A,I8,I5,3F13.5)') ' In-situ Fix',I,NBSG(I),FIX(I),GGOBS(I),GMOD2(I) if(I.eq.NLG+400.and.ipcheck.eq.1) write(*,'(A,I8,I5,3F13.5)') ' In-situ Fix',I,NBSG(I),FIX(I),GGOBS(I),GMOD2(I) end if end if end do if(ipcheck.eq.1) then write(*,'(A,F10.5)') 'BASEOBS =', BASEOBS print *,' ' end if C do I=1,NLG+N_ING C GMOD2(I) = GMOD2(I) - BBASE(I) ! Model brightness minus base compares with GOBS in write C GMOD2(I) = GMOD2(I) - BBASE(I) + BASEOBS ! Put in 5/16/11 BVJ GMOD2(I) = GMOD2(I) - BBASE(I) ! Remove BASEOBS change BVJ 4/29/2020 end do if(bWriteMS.and.NitG.eq.NiterT) then call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG+N_ING,PWG,GRATIO, & NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GGOBS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG) C call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO, C & 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 if(.not.bTcon) then ! density map maker from g-level call MkDMaptdn0n_in(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in, & PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) C call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, C & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, C & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, C & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) else ! density map maker from brightness call MkDMaptdn0n_in(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in, & PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) C call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, C & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, C & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, C & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) end if if(bFACE.and.Nit.eq.NiterT+1) & write (*,'(A,F10.5,A)') 'The 1 AU density was set to',DEN1AU,' through the iterations' C 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 and this day has come in the C form of the above (see Mk_Base_in 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 C if(bTcon) then C call Mk_Base_in(bFACE,VMAPD,DMAP,nLng,nLat,nTG,nTminTS,nTmaxTS,RRV,DEN1AU,FALLOFF,VEL1AU, C & BadR4(),BadR4(),DBASE) C end if C write(*,'(A,F12.6)') 'The DMAP lower density 1 AU limit is',DB2 do N=1,NTG if (NitG .eq. NiterT) then if(bBFIX) then if(N.eq.1) then print *, 'Since the data are adjusted to the insitu, change the adjusted to the mean by' write (*,'(F10.5,A,F10.5)') (-DEN1AUMN+DEN1AU),' The mean density has been reset to',DEN1AUMN end if call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N min max ',N,amin,amax do I=1,nLng do J=1,nLat if(DMAP(I,J,N).ne.BADR4()) then DMAP(I,J,N) = DMAP(I,J,N) + (-DEN1AUMN-DEN1AU)/(RRG*RRG) end if end do end do C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'After base removed DMAP N,amin,amax',N,amin,amax end if call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) end if C call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG)) call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG)) end do if(bBFIX .and. NitG .eq. NiterT) DEN1AU = DEN1AUMN 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 mkshiftdn0n. 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 Mkshiftdn0n' call stopwatch(' ',' ') end if C C Condition the next density map set for velocity (DMAPV) to use in mkshiftdn0n 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 MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) end if end if end if if(bstopw) then call stopwatch(' ',' ') print *, 'After V Mkshiftdn0n' end if end do ! End of iteration loop C******************************************************************************************************** C******************************************************************************************************** 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 if(bVCon) 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 if(VOBS(I).gt.-9999.999.and.VOBS(I).lt.99999.0) then write (*,'(3A,I4,A,F8.1,A,F9.3)') ' Bad V',SRCV(I),' L.O.S. removed, # =', & I,', V =',VOBS(I),', elong. =',XEV(I) end if end if end do write (*,'(A,I4,A)') ' There were',J,' bad V L.O.S. removed.' J = 0 ! List bad IPS sources do I=NLV+1,NLV+N_INV if (NBSV(I) .eq. 0) then J = J+1 if(VOBS(I).gt.-99999.9.and.VOBS(I).lt.999999.9) then write (*,'(3A,I4,A,F8.1)') ' Bad in-situ velocity',SRCV(I),' removed, # =', & I,', V =',VOBS(I) end if end if end do write (*,'(A,I4,A)') ' There were',J,' bad in-situ velocity values removed.' J = 0 do I=1,NLG if (NBSG(I).eq.0) then J = J+1 if(GOBS(I).gt.-99.999.and.GOBS(I).lt.999.999) then 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 if end do write (*,'(A,I8,A)') ' There were',J,' bad brightness L.O.S. removed.' J = 0 ! List bad IPS sources do I=NLG+1,NLG+N_ING if (NBSG(I) .eq. 0) then J = J+1 if(GOBS(I).gt.-999.99.and.GOBS(I).lt.9999.99) then write (*,'(3A,I4,A,F7.2)') ' Bad in-situ density',SRCV(I),' removed, # =', & I,', D =',GOBS(I) end if end if end do write (*,'(A,I4,A)') ' There were',J,' bad in-situ density values removed.' C Here VMAPHOLE contains the original point-P map (bVcon=.FALSE.) (or is equal to VMAP) do N=1,NTV call arrR4getminmax(nLngLat,VMAPHOLE(1,1,N),amin,amax) write(*,'(A,I4,E12.4,F11.4)') 'Iter End VMAPHOLE N, amin and amax = ', N, amin, amax end do 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) 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 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 *,' ' do N=1,NTG call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) write(*,'(A,I4,E12.4,F11.4)'), 'After iter VMAPD N, amin and amax = ', N, amin, amax end do print *,' ' print *, ' first Extractdvd (EA file)',XCbegG(1,1),XCtbegG,XCstrt print *,' ' write(*,'(2(A,I5))') 'Before the first extractdvd (EA file), nLng =',nLng,' nLat =',nLat do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) write(*,'(A,I4,E12.4,F11.4)') , 'DMAP N, amin and amax = ', N, amin, amax end do do N=1,NTG call arrR4getminmax(-nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do NinterD = NinterDD C call Extractdn_intel(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat,nMap,nTG, C & NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) call Extractdvd(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, ! New BVJ 2020 & nMap,nTG,NinterD,FALLOFFD,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. C RRSCON = (RRS/R1AU)**FALLOFFD !not in ipstd BVJ 2020 so removed C do N=1,nTG C call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) C end do nTF = nTG Nit = Nit + 1 nMapn = nMap cPre = 'nv3d' C The first write3D_infotd3DM makes maps that have holes to be used to show Carrington maps print *, ' ' print *, 'Before first write3Dinfo nv3d' C call write3D_infotd3DM(0,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, C & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) call write3D_infotd3D_3(NOneOr1,bForecast,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT) C These (mode 2) write3D_infotd calls make 3D arrays with time error maps to show reconstruction reliability do I=1,4*nTmaxG ! Will not work for Extractd3D if XCbegG's are different New BVJ 2020 XCbe(1,I) = XCbegG(1,1) XCbe(2,I) = XCbegG(2,1) end do C Begin conditioning for error writes and high resolution checks of errors C ANMAPD and AWTMAPD produced in MkDMaptdn0n_in C Begin conditioning of Density checks 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() VDEVD(I,J,N) = ANMAPD(I,J,N) ! Won't work with VDEV if NTG not = to NTV BVJ 2020 if(AWTMAPD(I,J,N).eq.0.0) AWTMAPD(I,J,N) = badR4() GDEV(I,J,N) = AWTMAPD(I,J,N) end do end do end do call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, ! New BVJ 2020, just to be safe nTG,nTG & nTG,nTG,ConsTG,VDEVD,ANMAPD,Vtmpt,DDT,DDD) call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, ! New BVJ 2020, just to be safe nTG,nTG & nTG,nTG,ConsTG,GDEV,AWTMAPD,Vtmpt,DDT,DDD) PIRSQD = (CONDT/aNdayG)*3.14159*((CONRD/(20.0/float(NGRESS)))**2) ! normalize weights and LOS to temporal and spatial resolution PIRSQD = PIRSQD*(3.0/float(NGRESS)) ! normalize according to the density LOS increment PIRSQD10 = PIRSQD/10.0 do N=1,nTG if(N.eq.1) then print *,' ' write(*,'(A,I5)'), ' Check to see if ANMAMPD and AWTMAPD nominal nTG =', nTG end if 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),PIRSQD10,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 call arrR4getminmax(-nLngLat,ANMAPD(1,1,N),amin,amax) ! Check to see that density error1 limits OK write(*,'(A,I4,F12.8,F15.4)'), 'ANMAPD N, amin and amax = ', N, amin, amax call arrR4getminmax(-nLngLat,AWTMAPD(1,1,N),amin,amax) ! Check to see that density error2 wts OK write(*,'(A,I4,F12.8,F15.4)') , 'AWTMAPD N, amin and amax = ', N, amin, amax end do C Begin conditioning of V checks C ANMAPV and AWTMAPV produced in MkVMaptdn0n_in 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 C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, ! This is bad in ipstd_20n BVJ 2020 C & nTV,nTV,ConsTV,VDEV,ANMAPV,Vtmpt,VVT,VVV) ! Provides zeros in check for densities C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, C & nTV,nTV,ConsTV,GDEV,AWTMAPV,Vtmpt,VVT,VVV) call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,ANMAPV,ANDMAPV,Dtmpt,VVT,VVV) call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,AWTMAPV,AWTDMAPV,Dtmpt,VVT,VVV) C The digital temporal resolution in the error files has the same temporal spacing as in the density files, C but have been made from the velocity 3d reconstruction and so to show use velocity resolution criteria. PIRSQV = (CONVT/aNdayV)*3.14159*((CONRV/(20.0/float(NVRESS)))**2) ! normalize weights and LOS to temporal and spatial res. PIRSQV = PIRSQV*(3.0/float(NGRESS)) ! normalize according to the velocity LOS increment do N=1,nTG ! Changed to nTG 2/8/2011 BVJ if(N.eq.1) then print *,' ' write(*,'(2(A,I5))'), ' Check to see if ANMAMPV and AWTMAPV nominal nTV =', nTV, ' nTG =',nTG end if call arrR4getminmax(nLngLat,ANDMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,ANDMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small ANDMAPV holes call ArrR4TimesConstant(-nLngLat,ANDMAPV(1,1,N),PIRSQV,ANDMAPV(1,1,N)) end if call arrR4getminmax(nLngLat,AWTDMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,AWTDMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small AWTMAPV holes call ArrR4TimesConstant(-nLngLat,AWTDMAPV(1,1,N),PIRSQV,AWTMAPV(1,1,N)) end if call arrR4getminmax(-nLngLat,ANDMAPV(1,1,N),amin,amax) ! Check to see that velocity error1 limits OK write(*,'(A,I4,F15.5,F15.5)'), 'ANDMAPV N, amin and amax = ', N, amin, amax call arrR4getminmax(-nLngLat,AWTDMAPV(1,1,N),amin,amax) ! Check to see that velocity error2 wts OK write(*,'(A,I4,F15.9,F15.5)'), 'AWTDMAPV N, amin and amax = ', N, amin, amax end do if(bWrerr) then C NinterC = NinterCC cPre = 'er1_' print *, ' ' print *, 'Before the error1 (L.O.S. crossings) write' C call write3D_infotd3DM(2,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, C & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & 90.0,Scale,ANDMAPV,ANMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,ANDMAPV,ANMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) print *, ' ' print *, 'Before the error2 (projected weights) write' C NinterC = NinterCC cPre = 'er2_' C call write3D_infotd3DM(2,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, C & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & 90.0,Scale,AWTMAPV,AWTMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,AWTMAPV,AWTMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) end if print *, ' ' print *, 'Before first filling' call Timesmooth(nLng,nLat,nTV,VMAP,2,CONVT/aNdayV,0.75*CONVT/aNdayV,VZtmp) ! Fill in times one time away New BVJ 2020 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 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 do print *, ' ' print *, 'Before second filling' call Timesmooth(nLng,nLat,nTG,DMAP,2,CONDT/aNdayG,0.75*CONDT/aNdayG,Dtmp) ! Fill in times one time away New BVJ 2020 do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy new BVJ 2020 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(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD/2.0,0,0.0,0.0) ! Change 5/9/08 to smooth polar holes end do print *, ' ' print *, 'After second filling' if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1).and.XCtbegG.eq.XCtbegV) then call ArrR4Copy(nLngLat*nTG,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 The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views print *, ' ' print *, 'Before nv3f* write' do N=1,NTG call arrR4getminmax(-nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do print *, ' ' print *, 'Before second write3D_infotd nv3f file' cPre = 'nv3f' C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. C call write3D_infotd3DM(0,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, C & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) call write3D_infotd3D_3(NOneOr2,bForecast,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT) if(bWrite3D_HR) then C Must change to low value in current version ! New RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do 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. print *, 'Before nv3h* write' do N=1,NTG call arrR4getminmax(-nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do cPre = 'nv3h' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = NinterDD ! 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 Mo = 0 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. print *, ' ' write(*,'(A,2(F7.4,A),I4,A,I4)') & 'Before 2nd High Resolution write. Rs =',RRRS,' Rt =',RRRSM,' nMapb =',nMapb,' nMapHR =',nMapHR call write3D_infotd3DMBE_HR_3(Mo,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, & iYrBG,RRRS,dRR,Bbegin3,EEnd3,nLng,nLat,nMapb,nMapHR,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdener,bVdener,bVverer,bDverer,ERDLOS,ERVLOS,ANMAPD,ANDMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) C Must change back to regular value in current version RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do end if C if(bforecast) then if(bforecast.and..not.bmagex) then ! New BVJ 2020 C RRSCON = (R1AU/RRS)**FALLOFFD C do N=1,nTG C call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) C end do print *, 'Before 2nd forecast and not magnetic extractdvd' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do print *, ' Before 2nd Extractdvd',XCbegG(1,1),XCtbegG,XCstrt NinterD = NinterDD call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C RRSCON = (RRS/R1AU)**FALLOFFD C do N=1,nTG C call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) C end do end if ************************************************************************************************************************************* C Produce a set of nv3o or nv3d files that are completely filled as well as we can make them. C The third write3D_infotd makes maps that have no holes to be used to provide the best possible boundary for a 3D MHD model if(bWrite3Do.or.bWrite3D_HRb.or..not.bForecast) then C ixxxxxxx = 0 C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop print *, 'Before third third filling' print *, 'bWrite3Do.or.bWrite3D_HRb.or..not.bForecast. NTG=',NTG do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do print *, 'bWrite3Do.or.bWrite3D_HRb.or..not.bForecast. NTV=',NTV 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 third velocity 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 print *, 'After third velocity filling. There should be no zeros' 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 third density 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 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 do call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) ! Fill all the density times print *, 'After third density filling. There should be no bad values' do N=1,NTG call arrR4getminmax(-nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do print *, 'After third 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(nLngLat*nTG,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 print *, 'After third V & D filling. There should be no bad values in VMAPD' 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 Change the density back since the falloff at distance is produced in Extractdvd C End change if(bWrite3Do) then C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) ! positive should not processes bad values end do C End change print *, 'Before third write3D_infotd nv3o files',nLng,nLat print *, '(There should be no bad values)' do N=1,NTG call arrR4getminmax(-nLngLat,DMAP(1,1,N),amin,amax) ! negative process bad values (should be none) print *, 'DMAP N, amin and amax = ', N, amin, amax end do cPre = 'nv3o' NinterD = NinterDD C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_3(0,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) C Change the density back since the falloff at distance is produced in Extractdvd RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change end if if(bWrite3D_HRb) then C Must change to low value in current version Removed BVJ 4/30/2020 RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do print *,' ' print *, 'Before base write and after DMAP decrease',nLng,nLat do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) write(*,'(A,I5,2F12.5)') 'DMAP N, amin and amax = ', N, amin, amax end do C End change cPre = 'nv3b' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = NinterDD ! In original map coordinates (so nMapHR = 1.5 AU). if(nMapHRb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRb = nMap ! Make the setting of nMapHR idiot proof if(nMapbb.gt.nMapHRb) nMapbb = nMapHRb Mo = 0 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. print *, ' ' write(*,'(A,2(F7.4,A),I4,A,I4)') & 'Before 3rd High Res. Base write. Rs =',RRRRS,' Rt =',RRRRSM,' nMapbb =',nMapbb,' nMapHR =',nMapHRb call write3D_infotd3DMBE_HR_3(Mo,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, & iYrBG,RRRRS,dRR,Bbeginb3,EEndb3,nLng,nLat,nMapbb,nMapHRb,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdener,bVdener,bVverer,bDverer,ERDLOS,ERVLOS,ANMAPD,ANDMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) C Change the density back since the falloff at distance is produced elsewhere RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do print *,' ' print *, 'after base write and after DMAP increase',nLng,nLat do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do C End change end if end if C************************************************************************************************************************************* if ( bMagnetic ) then C inexc = 0 C if(inexc.eq.0) go to 4567 dTT = (XCtendG-XCtbegG)/(sngl((nTF-1)*(NinterD+1))) ! Time increment call T3D_set (T3D__TT , 0, sngl(XCtbegG)) ! Needed for header write in Write3D_bbtm_3. call T3D_set (T3D__DTT , 0, dTT ) call T3D_iset(T3D__NCOFF , 0, NCoff ) ! Needed to pass these values, and for the header write in Write3D_bbtm_3 call T3D_set (T3D__RR , 0, RRS ) call T3D_set (T3D__DRR , 0, dRR ) call T3D_iset(T3D__NLNG , 0, nLng ) call T3D_iset(T3D__NLAT , 0, nLat ) call T3D_iset(T3D__LAT , 0, -90 ) call T3D_iset(T3D__LAT+1 , 0, +90 ) call T3D_iset(T3D__NRAD , 0, nMap ) call T3D_iset(T3D__NTIM , 0, nTim ) C 4567 continue call T3D_iset(T3D__NCOFF , 0, NCoff ) ! Must pass these values to Bfield_Get_3 in Get_bbtm_3 call T3D_set (T3D__RR , 0, RRS ) call T3D_iset(T3D__NLNG , 0, nLng ) call T3D_iset(T3D__NLAT , 0, nLat ) call T3D_iset(T3D__NRAD , 0, nMap ) call T3D_iset(T3D__NTIM , 0, nTim ) call T3D_iset(T3D__LAT , 0, -90 ) call T3D_iset(T3D__LAT+1 , 0, +90 ) write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC call get_bbtm_3(XCbeGG,nLng,nLat,nTG,RRS,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN, & XCINTG,BRR2DT,BRC2DT,BTC2DT,BNC2DT,iBflag) if(bmagexx) then print *, ' ' print *, 'Before write3D_bbtm_HR_3 after HR base' NintLng = NintLn NintLat = NintLa NintHt = NintH Ninter = 3 ClipLng = 90.0 nMapbb = nMapHRb ! In original map coordinates (so nMapHRMb = 0.25 AU). if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb IradialBRR = 0 IradialBRC = 0 IradialBTC = 0 IradialBNC = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 to smooth polar holes IradialBRR = IradialBRR + 1 ! Case 1 end if call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 to smooth polar holes IradialBRC = IradialBRC + 1 ! Case 2 end if call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 to smooth polar holes IradialBTC = IradialBTC + 1 ! Case 3 end if call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 to smooth polar holes IradialBNC = IradialBNC + 1 ! Case 4 end if end do if(bforecast) then if(IradialBRR.ne.0.and.IradialBRR.ne.(nTF+1)) then ! Case 1 write(*,'(A,I4,A)') ' There are no radial magnetic fields for',(nTF+1)-IradialBRR,' BRR source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRR2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRR2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BRR source surface',N,' was completely filled by gridsphere' else if(amin.eq.BadR4().and.amax.eq.BadR4()) then write(*,'(A,I4,A)') 'BRR source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRR source surfaces could not be filled with gridsphere' end if if(IradialBRC.ne.0.and.IradialBRC.ne.(nTF+1)) then ! Case 2 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBRC,' BRC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BRC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BRC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRC source surfaces were not filled with gridsphere' end if if(IradialBTC.ne.0.and.IradialBTC.ne.(nTF+1)) then ! Case 3 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBTC,' BTC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BTC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BTC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BTC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BTC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BTC source surfaces were not filled with gridsphere' end if if(IradialBNC.ne.0.and.IradialBNC.ne.(nTF+1)) then ! Case 4 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBNC,' BNC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BNC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BNC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BNC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BNC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BNC source surfaces were not filled with gridsphere' end if end if write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BEEndm3FACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC C BFACTORRR = 1.0 C BFACTORRC = 1.0 C BFACTORTC = 1.0 C BFACTORNC = 0.0 C do N=1,nTF+1 C call ArrR4Constant(nLngLat,BFACTORRR,BRR2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORRC,BRC2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORTC,BTC2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORNC,BNC2DT(1,1,N)) C end do C write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC do N=1,nTF if(bFillmag) then call ArrR4Constant(nLngLat,VMAPD(1,1,N),VMAPF(1,1,N)) else call ArrR4Constant(nLngLat,VMAP(1,1,N),VMAPF(1,1,N)) end if end do C call write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG,iYrBG, C & RRS,dRR,nLng,nLat,nMapbb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & ClipLng,VMAPF,XCshift,DVfact, C & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR) call write3D_bbtmBE_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG,iYrBG, & RRS,dRR,nLng,nLat,Bbeginm3,EEndm3,nMapbb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & ClipLng,VMAPF,XCshift,DVfact, & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR) end if if(bmagex) then print *,' ' print *, 'Before final Extractdvdm_3 magnetic field write' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do write(*,'(A,F10.6,F12.7,F10.6)') ' Before Final Extractdvdm_3',XCbegG(1,1),XCtbegG,XCstrt C ifieldw2 = 1 ifieldw2 = 0 if(ifieldw2.ne.0) then print *, ' ' write (*,'(A,100(I4,E12.5))') 'BRR2DT', (L,BRR2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BRC2DT', (L,BRC2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BTC2DT', (L,BTC2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BNC2DT', (L,BNC2DT(10,4,L),L=1,ntG) end if print *, ' ' NinterD = 3 call Extractdvdm_3(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift, & DVfact,DDfact,TT3D,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar) end if end if C************************************************************************************************************************************* if(.not.bMagnetic) then NinterD = NinterDD print *,' ' print *, 'Before final Extractvd non magnetic field write' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) print *, 'DMAP N, amin and amax = ', N, amin, amax end do call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if call Say('smeiipstd','S','Stop','program successfully completed') stop end