C+ C NAME: C ipstd_20n_inp_mag3_v16a 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 CATEGORY: C Data processing C CALLING SEQUENCE: C run ipstd_20n_inp_mag3_v16a C INPUTS: C Data file $CAM:G*. C OUTPUTS: C Rectangular grid of data points C FUNCTIONS/SUBROUTINES: C Local2UT,MAP_TZERO,READ_GO, C SetGipsy C INCLUDE: C include 'sun.h' C COMMON BLOCKS: C SIDE EFFECTS: C RESTRICTIONS: C PROCEDURE: C Latest 9/15/00 work - tried a copydtodv here and this filled C the density matrix a lot better when it was placed just before write3dinfo - C There was no apparent change in extractdv from the ea_file but maybe something C was different with ulysses. At this stage, this program does not give exactly C the same thing as the current version on the vax, but it is close. Prior I tried C using. Found an error 10/15/00 (A gridsphere2d was placed outside a do loop after C the first MKshiftdn0n. This version should now have no seam in each temporal value. C C > The modified Julian day MJD used for the IPS data is related to the C regular Julian day (as used in subroutine Julian) by C MJD = JD-2400000.5 C (i.e. the modified Julian day starts at midnight). C > Once the logical $CAM is set (during the first run of IPSD2020), the C only way to force IPSD2020 to look in another directory is to delete the C logical $CAM. C MODIFICATION HISTORY: C JUN-1992, Paul Hick (UCSD) C NOV-1995, Bernie Jackson (STEL) C March-1999 Bernie Jackson (UCSD) C June -1999 Bernie Jackson (UCSD) C August -1999 Bernie Jackson (UCSD) C Jan/Feb 2001 Bernie Jackson (UCSD) C Mar 2001 Bernie Jackson (UCSD) C- program ipstd_20n_inp_mag3_v16a parameter (NLmaxG = 30000, ! Max # G data points (13000) (400,000) & NLmaxV = 30000, ! Max # V data points (2000) & NGREST = 1, ! G temporal resolution factor better than 1 day & NVREST = 1, ! V temporal resolution factor better than 1 day & NGRESS = 1, ! Spatial G resolution factor better than 20 degrees & NVRESS = 1, ! Spatial V resolution factor better than 20 degrees & NF = 2, ! LOS resolution factor (Should be 1 if Spatial resolution is >= 2) & nTmaxG = 55*NGREST, ! Max # G times (50) (100) (150) & nTmaxV = 55*NVREST) ! Max # V times (50) (100) (150) parameter (NLOSG = 20*NF*NGRESS, ! # resolution elements along each G line of sight (40) (80) (125) & NLOSV = 20*NF*NVRESS, ! # resolution elements along each V line of sight (40) (80) (125) & NLOSGP1 = NLOSG+1, & NLOSVP1 = NLOSV+1, & dLOSG = 0.1/(NGRESS*NF), ! LOS G Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & dLOSV = 0.1/(NVRESS*NF), ! LOS V Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & xInc = 0.5) ! Carrington map "extra" increment parameter (nLng = 54*NGRESS+1, ! # Longitudes, resolution is 360/(nLng-1) (55) (109) (217) (325) & iLng = 18*NGRESS+1, ! # Longitudes to use in extracted map (19) (37) (73) (109) & nLat = 9*NGRESS+1, ! # Latitudes , resolution is 180/(nLat-1) (9)(10) (19) (37) (55) & nMap = 31, ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) & nLngLat = nLng*nLat, & nLngLatnTG =nLng*nLat*nTMaxG, & nLngLatnTV =nLng*nLat*nTMaxV) parameter (NintLn = 3, ! High resolution output parameters & NintLa = 3, & NintH = 3) parameter (dRR = 0.1) ! Resolution (AU) in radial direction (0.1) parameter (NTP = 16) ! number of body locations possible to extract. integer iXPG (NLmaxG), ! Array indices from daily G maps & iYPG (NLmaxG), & iMJDG (NLmaxG), ! Original MJD of daily G maps & IYRFG (NLmaxG), & IRECG (NLmaxG), & IYRSG (NLmaxG), ! Year of source observation & IYRSGG (NLmaxG), ! Year of source observation & NBSG (NLmaxG), & LON (NLOSGP1,NLmaxG), & LAT (NLOSGP1,NLmaxG), & ITIM (NLOSGP1,NLmaxG), & LON_in (NLmaxG), & LAT_in (NLmaxG), & ITIM_in (NLmaxG) real DISTG (NLmaxG), ! G distance from Sun to Earth & DISTGG (NLmaxG), ! G distance from Sun to Earth & XLSG (NLmaxG), ! Ecliptic longitude of Sun & XLSGG (NLmaxG), ! Ecliptic longitude of Sun & XLLG (NLmaxG), ! Longitude diff. (Earth-Sun) & XLLGG (NLmaxG), ! Longitude diff. (Earth-Sun) & XDLG (NLmaxG), ! Declination of Source & XDLGG (NLmaxG), ! Declination of Source & XCEG (NLmaxG), ! Carrington variable of Earth & XCEGG (NLmaxG), ! Carrington variable of Earth & XEG (NLmaxG), ! Source Elongation & XEGG (NLmaxG), ! Source Elongation & XCG (NLmaxG), ! Carrington variable & XCGG (NLmaxG), ! Carrington variable & YLG (NLmaxG), ! Heliographic latitude (degree) & YLGG (NLmaxG), ! Heliographic latitude (degree) & GGGOBS (NLmaxG), ! Dummy Solar wind observed G values & GOBS (NLmaxG), ! Solar wind observed G values & DGOBS (NLmaxG), ! Solar wind observed G value errors & GGOBS (NLmaxG), ! Dummy Solar wind observed G values & DGGOBS (NLmaxG), ! Dummy Solar wind observed G value errors & GOBSS (NLmaxG), ! Solar wind observed G values & GSSIZE (NLmaxG), ! Dummy Source g-level size & GSIZE (NLmaxG), ! Source g-level size & GFFREQ (NLmaxG), ! Dummy source g-level frequency & GFREQ (NLmaxG), ! Source g-level frequency & VVOBS (NLmaxV), ! Dummy Solar wind observed velocities & VOBS (NLmaxV), ! True Solar wind observed velocities & DVOBS (NLmaxV), ! True Solar wind observed velocity errors & VSIZE (NLmaxV), ! Source velocity size & VFREQ (NLmaxV), ! Source velocity frequency & GMOD2 (NLmaxG), ! Model G values & FIX (NLmaxG), ! Fix to model to make OK & GSSIG (NLmaxG), ! G-level source sigma & WTSG (NLOSG ,NLmaxG), & GWTij (NLOSG ,NLmaxG), & XLON (NLOSGP1,NLmaxG), & XLONP (NLOSGP1,NLmaxG), & XLAT (NLOSGP1,NLmaxG), & RPX (NLOSGP1,NLmaxG), & XLproj (NLOSGP1,NLMAXG), & VFAC (NLOSGP1,NLmaxG), ! New array (3/99) & DFAC (NLOSGP1,NLmaxG), & XLON_in (NLmaxG), & XLONP_in(NLmaxG), & XLAT_in (NLmaxG), & RPX_in (NLmaxG), & XLproj_in(NLmaxG), & VFAC_in (1,NLmaxG), & DFAC_in (1,NLmaxG) character SRCV (NLmaxV)*9, ! Velocity source name & SRCVLAST*9, & SRCVG (NLmaxV)*115, ! Full source info and g-value and velocity file & SRCGV (NLmaxG)*115, ! Full source info and g-value and velocity file & SRCVGsav(NLmaxV)*115, & SRCGVsav(NLmaxG)*115, & SRCVsav (NLmaxV)*9, & SRCGsav (NLmaxG)*9, & SRCGG (NLmaxG)*9, ! G-level source name & SRCGLAST*9, & SRCALAST*9, & cvgfiles(NLmaxG)*11, & cggfiles(NLmaxG)*11, & SRCGS*9 integer IYRFV (NLmaxV), & IRECV (NLmaxV), & NSRV (NLmaxV), & NSRG (NLmaxG), & IYRSV (NLmaxV), ! Year of source observation & NBSV (NLmaxV) real DISTV (NLmaxV), ! V distance from Sun to Earth & XLSV (NLmaxV), ! Ecliptic longitude of Sun & XLLV (NLmaxV), ! Longitude diff. (Earth-Sun) & XDLV (NLmaxV), ! Declination of Source & XCEV (NLmaxV), ! Carrington variable of Earth & XEV (NLmaxV), ! Source Elongation & XCV (NLmaxV), ! PLH & YLV (NLmaxV), ! PLH & VMOD (NLmaxV), ! Model velocities & VSSIG (NLmaxV), ! Velocity source sigma & VWTij (NLOSV,NLmaxV), & WTSV (NLOSV,NLmaxV), & VWTij_in(NLmaxV), & GWTij_in(NLmaxG), & WTSV_in (NLmaxV), & WTSD_in (NLmaxG), & WTPV_in (NLmaxV), & WTPD_in (NLmaxG) real XCinTV (nTmaxV), & XCinTG (nTmaxG), & DOYV (nTmaxV), & DOYG (nTmaxG), & XCbeV (2,nTmaxV), & XCbeGG (2,nTmaxG), & XCbe (2,4*nTmaxG) integer IYRV (nTmaxV), & IYRG (nTmaxG), & NMAPD (nLng,nLat,nTmaxG), ! Scratch integer arrays & NMAPV (nLng,nLat,nTmaxV) real VDEV (nLng,nLat,nTmaxV), & VPNT (nLng,nLat,nTmaxV), & PPGMAP (nLng,nLat,nTmaxG), ! Point-P G map (with holes) & VMAPHOLE(nLng,nLat,nTmaxV), & DM (nLng,nLat), & VM (nLng,nLat), & DMHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1), & VMHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1), & DMAPHOLE(nLng,nLat,nTmaxG), & DMAP (nLng,nLat,nTmaxG), ! Density map (no holes) & VMAP (nLng,nLat,nTmaxV), & VMAPD (nLng,nLat,nTmaxG), & DMAPV (nLng,nLat,nTmaxV), & GDEV (nLng,nLat,nTmaxG), & GPNT (nLng,nLat,nTmaxG), & DDfact (nLng,nLat,nMap,nTmaxG), & DVfact (nLng,nLat,nMap,nTmaxG), ! New array (3/99) & XCshift (nLng,nLat,nMap,nTmaxG,3), & V3D (nLng,nLat,nMap), & D3D (nLng,nLat,nMap), & V3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & D3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BR3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BT3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & BN3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & V3DT (nLng,nLat,nMap,4*nTmaxG), & D3DT (nLng,nLat,nMap,4*nTmaxG), & Vtmp (nLng,nLat,nTmaxG,3), ! Scratch real arrays & VtmpV (nLng,nLat,nTmaxV,3), & VZtmp (nLng,nLat,nTmaxV), & Dtmp (nLng,nLat,nTmaxG), & XLTtmp (nLng,nTmaxG,5), & XLtmp (nLng,nTmaxG), & XLtmpt (nLng,nTmaxG), & Dtmpt (nLng,nTmaxG), & Vtmpt (nLng,nTmaxG), & Vtp (nLng,nLat), & Dtp (nLng,nLat), & ANMAPD (nLng,nLat,nTmaxG), & AWTMAPD (nLng,nLat,nTmaxG), & DMAD (nLng,nLat,nTmaxG), & FIXMD (nLng,nLat,nTmaxG), & WTSMD (nLng,nLat,nTmaxG), & DMEANFD (nLng,nLat,nTmaxG), & DFIXM2D (nLng,nLat,nTmaxG), & WTPD (NLOSG,NLmaxG), & GMWTi (NLmaxG), & ANMAPV (nLng,nLat,nTmaxV), & AWTMAPV (nLng,nLat,nTmaxV), & VMAV (nLng,nLat,nTmaxV), & FIXMV (nLng,nLat,nTmaxV), & WTSMV (nLng,nLat,nTmaxV), & VMEANFV (nLng,nLat,nTmaxV), & VFIXM2V (nLng,nLat,nTmaxV), & WTPV (NLOSV,NLmaxV), & VVV (nLng), & DDD (nLng), & VVT (nTmaxV), & DDT (nTmaxG), & VRATION (200), & GRATION (200), C & Scale (2) & Scale (4) byte NSIDE(9*nLng*nLat) parameter (nCar = 2000) real*8 JDCar(nCar), & JD, JEpoch, & dLngSun,dLatSun,dDisSun real*8 DOYSG8 (NLmaxG), ! DOY (and fraction) of source obs. & DOYSGG8 (NLmaxV), ! DOY (and fraction) of source obs. & DOYSV8 (NLmaxV), ! DOY (and fraction) of source obs. & RADUMG8 (NLMAXG), ! Dummy array for glevel sort & RADUMV8 (NLMAXV), ! Dummy array for velocity sort & DOYSsav8, & dorbit real*8 XCtpr (NLOSGP1,NLMAXG), & XCtim (NLOSGP1,NLMAXG), & XCtpr_in(NLMAXG), & XCtim_in(NLMAXG), & XCintDV (nTmaxV), & XCintDG (nTmaxG) C------- C String written to output files. First three chars are also used in file name construction character cStrPPV(1)*80 /'PPV: Point P velocity map with holes'/, & cStrPPG(1)*80 /'PPG: Point P G-level map with holes'/, & cStrPPD(1)*80 /'PPD: 1 AU Point P density map with holes'/, & cStrDV (1)*80 /'DVV: Deconvolved velocity map'/, & cStrDD (1)*80 /'DDD: Deconvolved density map'/, & cStrDVs(1)*80 /'DVS: Deconvolved, smoothed velocity map'/, & cStrDDs(1)*80 /'DDS: Deconvolved, smoothed density map'/, & cStrVname*14 /'vts1941_00.txt'/, & cStrGname*14 /'gts1941_00.txt'/, & cSym*17 /'LIB__SIGNAL_QUIET'/ logical bGcon, bVcon, & bVNago, bVOoty, bVUCSD, bVEISC, bVGen, & bGUCSD, bGCamb, bGOoty, bGNago, bGGen, & bVACE, bDACE, bFACE, bGACE, & bVFACE,bVACE2, bDACE2, & bVWind, bDWind, & bVCELIAS,bDCELIAS, & bVmv2 /.FALSE./, ! Use mv^2 = constant & bGmv2 /.FALSE./, ! Use mv^2 = constant & bDaydV, bDaydG, & bOnly1, ! If .TRUE., use only one shift matrix (on density matrix) & bVproxy /.FALSE./, ! Velocity proxy map used for deconvolution & bDproxy /.FALSE./, ! Density proxy map used for deconvolution & bMagnetic /.TRUE./, ! Magnetic map available & bMagnet /.FALSE./, ! Magnetic map check & bmagexx /.TRUE./, ! Produce 3D magnetic field data files & bmagex /.TRUE./, ! Extract magnetic field data as well as velocity and density data & bAdjust /.TRUE./, ! Adjust the magnetic field BF factors for solar cycle and instrument types & bGWRITE /.FALSE./, ! Write out a G data file & bVWRITE /.FALSE./, ! Write out a V data file & bMirror /.FALSE./, & bjMirror /.FALSE./, & bWrite /.FALSE./, & bWrite3D1 /.FALSE./, & bWrite3D2 /.FALSE./, & bWrerr /.FALSE./, & bDdener /.TRUE./, & bVdener /.FALSE./, C & bVdener /.TRUE./, C & bVverer /.FALSE./, & bVverer /.TRUE./, & bDverer /.FALSE./, & bDdenerb /.TRUE./, & bVdenerb /.FALSE./, C & bVdenerb /.TRUE./, C & bVvererb /.FALSE./, & bVvererb /.TRUE./, & bDvererb /.FALSE./, & bWrite3D_HR /.TRUE./, & bWrite3D_HRb /.FALSE./, & bWrite3Do /.TRUE./, & bWriteMS /.FALSE./, & bWritelos /.FALSE./, & bOneOr1 /.TRUE./, & bOneOr2 /.TRUE./, & bForecast, & bSource /.TRUE./, ! Write out good sources & bExtract(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,4*.FALSE./, & bEWRITE, ! Write out an Earth file & bMWRITE, ! Write out a Mercury file & bVEWRITE, ! Write out a Venus file & bMAWRITE, ! Write out a Mars file & bUWRITE, ! Write out a Ulysses file & bSAWRITE, ! Write out a STEREO A file & bSBWRITE, ! Write out a STEREO B file & bMSWRITE, ! Write out a MESSENGER file & b67WRITEF, ! Write out a comet 67P/CG file & bExtractf(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,4*.FALSE./, & bEWRITEF, ! Write out an Earth filled file & bMWRITEF, ! Write out a Mercury filled file & bV2WRITE, ! Write out a Venus filled file & bMAWRITEF, ! Write out a Mars filled file & bUWRITEF, ! Write out a Ulysses filled file & bSAWRITEF, ! Write out a STEREO A filled file & bSBWRITEF, ! Write out a STEREO B filled file & bMXWRITEF, ! Write out a MESSENGER filled file & bC67WRITEF ! Write out a comet 67P/CG filled file equivalence (bExtract( 3), bEWRITE), & (bExtract( 1), bMWRITE), & (bExtract( 2), bVEWRITE), & (bExtract( 4), bMAWRITE), & (bExtract(10), bUWRITE), & (bExtract(13), bSAWRITE), & (bExtract(14), bSBWRITE), & (bExtract(15), bMSWRITE), & (bExtract(16), b67WRITE), & (bExtractf( 3), bEWRITEF), & (bExtractf( 1), bMWRITEF), & (bExtractf( 2), bV2WRITEF), & (bExtractf( 4), bMAWRITEF), & (bExtractf(10), bUWRITEF), & (bExtractf(13), bSAWRITEF), & (bExtractf(14), bSBWRITEF), & (bExtractf(15), bMXWRITEF), & (bExtractf(16), bC67WRITEF) character cPrefix (NTP)*2 /'Me','Ve','Ea','Ma','Ju','Sa','Ur','Ne','Pl','Ul','H1','H2','Sa','Sb','MS', & 'CG'/ character cPrefixf(NTP)*2 /'m1','v2','e3','m4','j5','s6','u7','n8','p9','uy','ha','hb','s1','s2','mx', & '67'/ integer iEorWG /0/, & iEorWV /0/, c & MJD, c & MJDfrst /48180/, ! MJD for earliest data file c & MJDlast /49622/, ! MJD for latest data file & MJDfrst / 46100/, & MJDlast / 60000/, & LimitoV /-1/, & LimitoG /-1/, & NLOSWV /1/, & NLOSWG /1/, & NTV /44/, & NTG /44/ parameter iVarp = 27 ! Number of potential data files accessed via a parameter list (6 IPS + 5 in-situ + 16 mag) parameter iSyst = 7 ! total number of IPS systems accessed at one time integer iSys (iSyst) ! System number 1-UCSD, 2-Eiscat, 3-Ooty, 4-Nagoya, 5-MEXART, 6-General integer NiterT /18/, C & iOffUT / 0/, & iOffUT / 7/, ! Changed 5/2/11 BVJ & NAV / 0/, & NC / 1/, ! # rotations averaged & nDUMP / 0/, ! Threshold for #pnts/bin & nFILL / 3/, ! Min. # neighbours for fill & NmidHRV /0/, & NmidHRG /0/, & NLVmin /50/, ! Default minimum number of velociy lines of sight & NLGmin /50/ ! Default minimum number of g-level lines of sight real YLbeg /-90./, & YLend / 90./, & XCadj(3) /3*0./, & XCadjB(3) /3*0./, & XCmargin /.5/, & XElowG /11.5/, & XElowV /11.5/, & XElowNG /11.5/, & XElowNV /11.5/, & XElowOG /11.5/, & XElowOV /11.5/, & XElowKG /11.5/, & XElowKV /11.5/, & XElowMG /21.0/, & XElowMV /21.0/, & XElowPG /26.1/, & XElowPV /26.1/, & XEhighG /80./, & XEhighV /180./, & XRlimG /90./, & XRlimV /90./, & GPOWER /0./, & R1AU /1.0/, & DEN1AU /5.0/, & ANLIMITV /1.2/, & ANLIMITG /1.2/, & ERDLOS /1.2/, & ERVLOS /1.2/, & ERDLOSB /1.2/, & ERVLOSB /1.2/, & aNdayV /1.0/, & aNdayG /1.0/ real VMAPSIG /1.0/, & DMAPSIG /1.0/, & GSIG /1.0/, & GSIGB /3.0/, & VSIG /1.0/, & VSIGB /3.0/, & FCONSV /1.0/, & FCONSG /1.0/, & CONRV /14./, & CONRD /14./, & CONRO /8.0/, & CONVT /0.65/, & CONDT /0.65/, & CONOT /0.45/, & ConsTV /1.0/, & ConsTG /1.0/, & VLIMM /1200./, ! Velocity upper limit for sources & VLIML /201./, ! Velocity lower limit for sources & VLIM /1200./, ! Velocity upper limit for tomography & VLIL /100./, ! Velocity lower limit for tomography & GLIMP /4.0/, ! G-level upper limit for sources & GLIMOF /0.5/, ! Ooty g-level input excursion multiplier & GLIMNF /1.0/, ! Nagoya g-level input excursion multiplier & GLIMKF /1.0/, ! KSWC g-level input excursion multiplier & GLIMMF /0.2/, ! MEXART g-level excursion input multiplier & GLIMPF /1.0/, ! Pushchino g-level input excursion multiplier & GLIMF /1.0/, ! General g-level input excursion multiplier & GLIMM /0.2/, ! G-level lower limit for sources & DMAX /200./, ! Density upper limit for in-situ inputs & DMIN /0.1/, ! Density lower limit for in-situ inputs & VVMAX /1500./, ! Velocity upper limit for in-situ inputs & VVMIN /100./ ! Velocity lower limit for in-situ inputs real RRS / -15. /, ! Reference distance for deconvolution & RRX / 1.0 /, ! Reference distance for Map extraction & VUL / 1500.0 /, ! Velocity upper limit for maps & VLL / 100.0 /, ! Velocity lower limit for maps & DUL / 500.0 /, ! Density upper limit for maps & DLL / 0.0 / ! Density lower limit for maps real DAV /1./, & Speed /400./, & FALLOFFD /2.0/, & FALLOFFV /0.0/, & FALLOFFBR /2.0/, & FALLOFFBT /1.0/, & FALLOFFBN /1.34/, & PWG /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.20/, ! Radial power of density & PWRV /0.20/, ! Radial power of velocity & AEmaxV /800.0/, ! Anti-Earth maximum velocity & AEmaxD /10.0/, ! Anti-Earth maximum density at 1 AU & AEminV /200.0/, ! Anti-Earth minimum velocity & AEminD /2.0/, ! Anti-Earth minimum density at 1 AU & AEangV /135./, ! Angular distance from Earth to begin V limit & AEangD /135./ ! Angular distance from Earth to begin D limit real*8 MJDref, & MJDcntr, & XCtbegV, & XCtbegG, & XCtendV, & XCtendG character cCam*4, & cMon*3, & cStrID*11, & cVar(27)*256, & ccVar*256, & cArg*100, & cStr*80, & cPre*4 character cWildNagoya*100, & cWildOoty*100, & cWildGen*80, & cWildGens(6)*100, ! No of IPS systems general files & cWildACE*80, & cWildACE2*80, & cWildWind*80, & cWildVCELIAS*80, & cWildDCELIAS*80, & cWild*80 external iReadNagoyan8, iProcessNagoyan8 external iReadOotyn8, iProcessOotyn8 external iReadGen8_n, iProcessGen8_n external EARTH include 'dirspec.h' include 'sun.h' include 't3d_array_3.h' include 'b3d_param_3.h' C include 't3d_array.h' character cWildMAG(B3D__NN)*256 logical bWildMAG(B3D__NN) integer Flt2Str c**TJD**added for magnetic time dependant version real Btmp(nLng,nLat,3), !scratch c & XCbegMAT(nTmaxG), !3-Car time frame c & XCendMAT(nTmaxG), & XCbegMAT(nTmaxG), !3-Car time frame & XCendMAT(nTmaxG), & XCbegROI(nTmaxG), !ROI time frame & XCendROI(nTmaxG), & TrueVel(nLng,nLat,nMap,nTmaxG), !real velocity array & BR3D(nLng,nLat,nMap), !scratch & BT3D(nLng,nLat,nMap), !scratch & BN3D(nLng,nLat,nMap), !scratch & BRR2DT(nLng,nLat,nTmaxG), !scratch & BRC2DT(nLng,nLat,nTmaxG), !scratch & BTC2DT(nLng,nLat,nTmaxG), !scratch & BNC2DT(nLng,nLat,nTmaxG), !scratch & TT3D (4*nTmaxG), !holds XC corresponding to 6 hour int. time steps output by write3dinfo_td3d & XC3DT (nLng,nLat,nMap,4*nTmaxG,3), & BRcheck (nLng,nLat), & Vcheck5 (nLng,nLat), & Vcheck7 (nLng,nLat) call T3D_set_mode(T3D__MODE_TIME, 1) iVar = iVarp ! Up to 6 IPS pointers, 5 in-situ pointers, 16 magnetic field input = 27 total call ForeignArg(' ',iVar,cVar,cArg) call ArrI4Constant(iSyst,0,iSys) ! Initialize iSys array to 0 (the total possible number of IPS systems, iSyst = 7) do I=1,iVar write(*,'(A,I3,1X,A64)') 'I, cVar(I)', I, cVar(I) ccVar = cVar(I) if(ccVar(1:9) .eq.'UCSD=UCSD') iSys(1) = 1 if(ccVar(1:8) .eq.'gen=UCSD') then call ForeignFile(1,ccVar,'gen',cwildGens(1)) iSys(1) = 1 end if if(ccVar(1:13).eq.'EISCAT=EISCAT') iSys(2) = 2 if(ccVar(1:10).eq.'gen=EISCAT') then call ForeignFile(1,ccVar,'gen',cwildGens(2)) iSys(2) = 2 end if if(ccVar(1:9) .eq.'ooty=ooty') iSys(3) = 3 if(ccVar(1:8) .eq.'gen=ooty') then call ForeignFile(1,ccVar,'gen',cwildGens(3)) iSys(3) = 3 end if if(ccVar(1:13).eq.'nagoya=nagoya') iSys(4) = 4 if(ccVar(1:10).eq.'gen=nagoya') then call ForeignFile(1,ccVar,'gen',cwildGens(4)) iSys(4) = 4 end if if(ccVar(1:8) .eq.'gen=KSWC') then call ForeignFile(1,ccVar,'gen',cwildGens(5)) iSys(5) = 5 end if if(ccVar(1:10).eq.'gen=MEXART') then call ForeignFile(1,ccVar,'gen',cwildGens(6)) iSys(6) = 6 end if if(ccVar(1:8).eq.'gen=push') then call ForeignFile(1,ccVar,'gen',cwildGens(7)) iSys(7) = 7 end if end do C------- C Calculate start times for Carrington rotations between 1985 and now. C Rotation NCoff+I starts at time JDCar(I). aNdayV = aNdayV/NVREST aNdayG = aNdayG/NGREST ConsTMV = ConsTV*nTV ConsTMD = ConsTG*nTG ConstL = nLng/3 DGOBSS = 1.0 Scale(1) = 0.0 Scale(2) = 1.0 Scale(3) = 0.0 Scale(4) = 1.0 I = iDeleteSymbol(cSym,2) C iYr = 1997 C iYr = 1985 iYr = 1985 Doy = 1. call MAP_TZERO(EARTH,iYr,Doy,.01,nCar,JDCar) NCoff = N_CARRINGTON(iYr,Doy) call ArrR4Constant(NLmaxV,327.0,VFREQ) ! Initialize VFREQ array ! Nagoya or Ooty default call ArrR4Constant(NLmaxG,327.0,GFREQ) ! Initialize GFREQ array ! Nagoya or Ooty defaut call ArrR4Constant(NLmaxV,0.1,VSIZE) ! Initialize VSIZE array ! Nagoya or Ooty default call ArrR4Constant(NLmaxG,0.1,GSIZE) ! Initialize GSIZE array ! Nagoya or Ooty default call ArrR4Constant(NLmaxV,5000.,XCEV) ! Initialize XCEV array call ArrR4Constant(NLmaxG,5000.,XCEG) ! Initialize XCEGG array call ArrR4Constant(NTmaxV,5000.,XCintV) ! Initialize XCintV array call ArrR4Constant(NTmaxG,5000.,XCintG) ! Initialize XCintG array C------ C iOffUT is an offset time (in hours) between the local computer time and UT. This is used only if the program C is operated in forecast mode. Note that the local time is assumed to be the same as the C computer time. The default offset is 7 hours (the difference between PDT and GST). call AskYN('Forecast mode?$no',bForecast) ! Regular or forecast mode if (bForecast) then call AskI4('Offset between local computer time and UT (hours)?',iOffUT) NLVmin = NLVmin/2 NLGmin = NLGmin/2 end if C niterT = 2 call AskI4('Number of iterations?',NiterT) ! # iterations call AskR4('Deconvolution distance (AU; 0=solar surface)?',RRS) if (RRS .lt. 0) RRS = -RRS*SUN__RAU ! Neg. RRS are in units of solar radii RRS = max(RRS,sngl(SUN__RAU)) ! Keep RRS > 1 solar radius ! Deconvolve velocity, G-level or both? call AskWhat('Deconvolve: velocity, G-level, both?$1$3',I) C------- C A proxy map for g-values and/or velocity can be used to replace the point-P C map determined from the observations. C The map should be in files DMAP.DAT and VMAP.DAT. They must be ascii files of dimension C nLng, nLat. C C If bVmv2=TRUE then bVcon=FALSE C If bGmv2=TRUE then bGcon=FALSE bVcon = I .eq. 1 .or. I .eq. 3 bGcon = I .eq. 2 .or. I .eq. 3 if (bVcon) then ! V deconvolution: select IPS station call AskWhat('Use velocity data from: UCSD, EISCAT, Ooty, Nagoya, Gen,?$1$4',I) bVUCSD = I .eq. 1 bVEISC = I .eq. 2 bVOoty = I .eq. 3 bVNago = I .eq. 4 bVGen = I .eq. 5 noSys = 0 if(bVUCSD) XElowV = 35.0 if(bVEISC) XElowV = 3.0 if(bVOoty) XElowV = 11.5 if(bVNago) XElowV = 11.5 do ii=1, iSyst if(iSys(ii).ne.0) noSys = noSys + 1 end do if(noSys.eq.0) then Print*, 'You asked for a General input velocity IPS file but did not specify its path. Please specify the path.' stop end if if(bVUCSD) then Print *, 'Sorry: The UCSD Velocity reconstructions are not yet supported. Try another.' stop end if bVmv2 = .FALSE. if(bVEISC) NLOSWV = 4 if(bVNago) NLOSWV = 1 if(bVOoty) NLOSWV = 1 if(bVGen) NLOSWV = 1 call AskI4('Minimum # of velocity lines of sight to use in the reconstruction?',NLVmin) call AskYN('Use in-situ velocity in the 3D reconstructions?$yes',bVFACE) bVACE = .FALSE. bVACE2 = .FALSE. bVWind = .FALSE. bVCELIAS = .FALSE. if(bVFACE) then call AskWhat('Use in-situ velocity from: ACE_L0, ACE_L2, Wind, CELIAS$1$3',J) bVACE = J .eq. 1 bVACE2 = J .eq. 2 bVWind = J .eq. 3 bVCELIAS = J .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 bVEISC = .FALSE. bVNago = .FALSE. bVOoty = .FALSE. call AskYN('Do you want to make a velocity map using mv^2 = con?$no',bVmv2) end if if(.not.bVFACE) bVACE = .FALSE. if (bGcon) then ! G deconvolution: select IPS station if(bVOoty) call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$3',I) if(bVNago) call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$4',I) if(bVGen) call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$5',I) if(.not.bVUCSD.and..not.bVNago.and..not.bVOoty.and..not.bVGen) & call AskWhat('Use G-level data from: UCSD, Cambridge, Ooty, Nagoya, Gen,?$1$5',I) bGUCSD = I .eq. 1 bGCamb = I .eq. 2 bGOoty = I .eq. 3 bGNago = I .eq. 4 bGGen = I .eq. 5 noSys = 0 do ii=1, iSyst if(iSys(ii).ne.0) noSys = noSys + 1 end do if(noSys.eq.0) then Print*, 'You asked for a General input density IPS file but did not specify its path. Please specify the path.' stop end if if(bGUCSD) then Print *, 'Sorry: The UCSD density reconstructions are not supported. Try another.' stop end if bGmv2 = .FALSE. cCam = cEnvi//'CAM' if(bGGen ) NLOSWG = 1 if(bGNago) NLOSWG = 1 if(bGOoty) NLOSWG = 1 if(bGcamb) NLOSWG = 2 if(bGUCSD) NLOSWG = 1 call AskI4('Minimum # of density lines of sight to use in the reconstruction?',NLGmin) else ! No g deconvolution bGUCSD = .FALSE. bGCamb = .FALSE. bGOoty = .FALSE. bGNago = .FALSE. bGgen = .FALSE. call AskYN('Do you want to make a density map using mv^2 = con?',bGmv2) end if if(bGcon) call AskYN('Use in-situ density in the 3D reconstructions?$yes',bFACE) bDACE = .FALSE. bDACE2 = .FALSE. bDWind = .FALSE. bDCELIAS = .FALSE. if(bGcon.and.bFACE) then I = 0 if(bVACE) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) if(bVACE2) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$2',I) if(bVWind) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$3',I) if(bVCELIAS) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$4',I) C if(I.eq.0) call AskWhat('Use in-situ density 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.bGcon.and.bFACE) then if(bVACE) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$1',I) if(bVACE2) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$2',I) if(bVWind) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$3',I) if(bVCELIAS) call AskWhat('Use in-situ density from: ACE_L0, ACE_L2, Wind, CELIAS$1$4',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(bFACE) call AskR4('Minimum density for fit in-situ density values?',DMIN) if(bFACE) call AskR4('Maximum density for fit in-situ density values?',DMAX) DLIM = DMAX/(RRS*RRS) if(bDACE) call AskYN('Fit mean g-level to ACE_L0?$yes',bGACE) if(bDACE2) call AskYN('Fit mean g-level to ACE_L2?$yes',bGACE) if(bDWind) call AskYN('Fit mean g-level to Wind?$yes',bGACE) if(bDCELIAS) call AskYN('Fit mean g-level to CELIAS?$yes',bGACE) if (bGCamb) call SetGipsy(cCam,0,iEdt,MJDfrst,MJDlast) if (bVNago .or. bGNago .or. bVEISC) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya) if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty ) C if (bVGen .or. bGGen ) call ForeignFile(iVar,cVar,'gen' ,cWildGen ) if (bVACE .or .bDACE) then call ForeignFile(-iVar-1,cVar,'ace',cWildACE) if (cWildACE .eq. ' ') then cWildACE = '/zshare/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 = '/zshare/dat/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 = '/zshare/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 = '/zshare/dat/insitu/celias_????.hravg' end if if (cWildVCELIAS .eq. ' '.and.bForecast) then cWildVCELIAS = '/zshare/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 = '/zshare/dat/insitu/celias_????.hravg' end if if (cWildDCELIAS .eq. ' '.and.bForecast) then cWildDCELIAS = '/zshare/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(bVEISC) NmidHRV = -1 if(bVNago) NmidHRV = -9 if(bVGen) NmidHRV = -9 if(bVOoty) NmidHRV = -6 call AskI4('Number of hours from V UT midnight?',NmidHRV) end if call AskI4('Number of velocity time values?',NTV) ! # NTV time itervals if(NTV.gt.NTmaxV-2) then NTV = NTmaxV-2 call AskI4('Number of velocity time values?',NTV) ! Ask again end if end if if(bGcon) then bDaydG = .TRUE. call AskYN('Deconvolve G from one day to the next$yes',bDaydG) if(bDaydG) then call AskR4('Number of days to combine for g-level?',aNdayG) if(bGNago) NmidHRG = -9 if(bGGen) NmidHRG = -9 if(bGOoty) NmidHRG = -6 call AskI4('Number of hours from g UT midnight?',NmidHRG) end if call AskI4('Number of g-level time values?',NTG) ! # NTG time itervals if(NTG.gt.NTmaxG-2) then NTG = NTmaxG-2 call AskI4('Number of g-level time values?',NTG) ! Ask again end if end if C I = iFilePath('/home/bjackson/dat/nagoya/',1,'yearly','nagoya.'//cWildChar(:4),cWildNagoya) C I = iFilePath(cEnvi//'DAT',1,'Ooty' ,'Ooty.' //cWildChar(:4),cWildOoty ) C------- C Calculate range of Carrington variables to work with: [XClo,XChi] C XClo is based on the first data file available (MJDfrst). C XChi is based on the current time. call Julian(11,iYr,Doy,dble(MJDfrst),JEpoch) XClo = NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XClo = XClo-.5 C call Local2UT(iOffUT,iYr,Doy) call Local2UT(NmidHRV,iYr,Doy) ! Changed 5/2/11 BVJ XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo) RRV = RRS RRG = RRS if (bVcon) then call AskR4('# of V-points per bin that constitute a deconvolution?',ANLIMITV) if(NiterT.gt.8) LimitoV = NiterT/2 call AskI4('Iteration at which to limit V-source deviations, -1=none?',LimitoV) C------- C Fix the reference distance used for the V map formed. C Get E or W or EW; Set the elongation angles call AskWhat('All velocity data, East data, West data ?$0',iEorWV) if (iEorWV .eq. 2) iEorWV = -1 if(bVGen) then if(iSys(3).eq.3) call AskR4('Minimum elongation for Ooty V-data?$0$11.5$',XElowOV) if(iSys(4).eq.4) call AskR4('Minimum elongation for Nagoya V-data?$11.5$',XElowNV) if(iSys(5).eq.5) call AskR4('Minimum elongation for KSWC V-data?$11.5$',XElowKV) if(iSys(6).eq.6) call AskR4('Minimum elongation for MEXART V-data?$21.0$',XElowMV) if(iSys(7).eq.7) call AskR4('Minimum elongation for Pushchino V-data?$21.0$',XElowPV) else call AskR4('Minimum elongation for V-data??$0$180$',XElowV) end if call AskR4('Maximum elongation for V-data?$'//cStr(:Flt2Str(XElowV,1,cStr))//'$180$',XEhighV) call AskR4('Anti-Earth maximum velocity at 1 AU?$'//cStr(:Flt2Str(600.0,1,cStr))//'$1000.0$',AEmaxV) call AskR4('Anti-Earth minimum velocity at 1 AU?$'//cStr(:Flt2Str(0.,1,cStr))//'$400.0$',AEminV) call AskR4('Angular distance from Earth to begin V limit?$'//cStr(:Flt2Str(90.,1,cStr))//'$180.0$',AEangV) call AskR4('Maximum RA relative to Sun for V-data?$0$180$',XRlimV) RRV = RRS call AskR4('V-map reference distance in AU (0=solar surface)?',RRV) if (RRV .lt. 0) RRV = -RRV*SUN__RAU RRV = max(RRV,sngl(SUN__RAU)) call AskYN('Do you want to use a proxy map for velocity?',bVproxy) call AskR4('Maximum upper limit for tomography velocity?',VLIM) call AskR4('Maximum lower limit for tomography velocity?',VLIL) call AskR4('Maximum upper limit for source velocity?',VLIMM) if(bVNago) VLIML = VLIL if(bVGen) VLIML = VLIL call AskR4('Maximum lower limit for source velocity?',VLIML) end if if (bGcon) then call AskR4('# of G-points per bin that constitute a deconvolution?',ANLIMITG) if(NiterT.gt.8) LimitoG = NiterT/2 call AskI4('Iteration at which to limit G-source deviations, -1=none',LimitoG) C Get E or W or EW; Set the elongation angles call AskWhat('All G-level data, East data, West data ?$0',iEorWG) if (iEorWG .eq. 2) iEorWG = -1 if(bGNago) XElowG = XElowV if(bGGen) XElowG = XElowV if(bGOoty) XElowG = XElowV if(bGCAMB) XElowG = 30.0 if(bGUCSD) XElowG = XElowG if(bGGen) then if(iSys(3).eq.3) call AskR4('Minimum elongation for Ooty G-data?$0$11.5$',XElowOG) if(iSys(4).eq.4) call AskR4('Minimum elongation for Nagoya G-data?$11.5$',XElowNG) if(iSys(5).eq.5) call AskR4('Minimum elongation for KSWC G-data?$11.5$',XElowKG) if(iSys(6).eq.6) call AskR4('Minimum elongation for MEXART G-data?$21.0$',XElowMG) if(iSys(7).eq.7) call AskR4('Minimum elongation for Pushchino G-data?$21.0$',XElowPG) else call AskR4('Minimum elongation for G-data?$0$180$',XElowG) end if if(bGNago) XEhighG = 180. if(bGGen) XEhighG = 180. if(bGOoty) XEhighG = 180. if(bGCAMB) XEhighG = 80.0 if(bGUCSD) XEhighG = 80.0 call AskR4('Maximum elongation for G-data?$'//cStr(:Flt2Str(XElowG,1,cStr))//'$180$',XEhighG) call AskR4('Anti-Earth maximum density at 1 AU?$'//cStr(:Flt2Str(5.,1,cStr))//'$30.$',AEmaxD) call AskR4('Anti-Earth minimum density at 1 AU?$'//cStr(:Flt2Str(0.,1,cStr))//'$4.0$',AEminD) call AskR4('Angular distance from Earth to begin D limit?$'//cStr(:Flt2Str(90.,1,cStr))//'$180.0$',AEangD) call AskR4('Maximum RA relative to Sun for G-data?$0$180$',XRlimG) RRG = RRS call AskR4('G-map reference distance in AU (0=solar surface)?',RRG) if (RRG .lt. 0) RRG = -RRG*SUN__RAU RRG = max(RRG,sngl(SUN__RAU)) call AskYN('Do you want to use a proxy map for density?',bDproxy) call AskR4('Maximum limit for G deconvolution?',GLIMP) call AskR4('Minimum limit for G deconvolution?',GLIMM) NDUP = 0 if(bGGen) then if(iSys(3).eq.3) call AskR4('Factor by which to modify Ooty g-level excursions?',GLIMOF) if(iSys(4).eq.4) call AskR4('Factor by which to modify Nagoya g-level excursions?',GLIMNF) if(iSys(5).eq.5) call AskR4('Factor by which to modify KSWC g-level excursions?',GLIMKF) if(iSys(6).eq.6) call AskR4('Factor by which to modify MEXART g-level excursions?',GLIMMF) if(iSys(6).eq.6) call AskI4('Duplicate MEXART g-level excursions?',NDUP) if(iSys(7).eq.7) call AskR4('Factor by which to modify Pushchino g-level excursions?',GLIMPF) else call AskR4('Factor by which to modify g-level excursions?',GLIMF) end if end if call AskR4('G-level time hole filter?',ConsTG) call AskR4('Velocity time hole filter?',ConsTV) CONRD = CONRD/float(NGRESS) call AskR4(' G-level spatial filter?',CONRD) if(bVEISC) CONRV = CONRD if(bVNago) CONRV = CONRD if(bVGen) CONRV = CONRD if(bVOoty) CONRV = CONRD if(CONRD.ne.CONRV) CONRV = CONRD call AskR4('Velocity spatial filter?',CONRV) CONDT = CONDT/float(NGREST) call AskR4(' G-level temporal filter?',CONDT) if(bVEISC) CONVT = CONDT if(bVNago) CONVT = CONDT if(bVGen) CONVT = CONDT if(bVOoty) CONVT = CONDT if(CONVT.ne.CONDT) CONVT = CONDT call AskR4('Velocity temporal filter?',CONVT) if(bGNago) PWG = 0.35 ! Returned back to fit data 5/30/2011 BVJ C if(bGNago) PWG = 0.20 ! Changed 4/28/2011 to fit data BVJ if(bGGen) PWG = 0.35 if(bGOoty) PWG = 0.35 if(bGCamb) PWG = 0.15 if(bGUCSD) PWG = 0.15 call AskR4(' G-level power to fit density?',PWG) if(bGcon) call AskR4('1 AU average density to fit base?',DEN1AU) if(bVEISC) PWV = 0.35 if(bVNago) PWV = 0.35 if(bVGen) PWV = 0.35 if(bVOoty) PWV = 0.35 if(PWV.ne.PWG) PWV = PWG call AskR4('G-level power to fit velocity?',PWV) call AskR4('Radial G-level falloff to fit G density?',PWRG) if(bVEISC) PWRV = PWRG if(bVNago.and.bGNago) PWRV = PWRG if(bVGen.and.bGGen) PWRV = PWRG if(PWRV.ne.PWRG) PWRV = PWRG call AskR4('Radial G-level falloff to fit V density?',PWRV) call AskR4('Solar wind speed for G-level traceback?',Speed) RRD = RRG call AskR4('Map Extraction distance (AU; 0=solar surface)?',RRX) if (RRX .lt. 0) RRX = -RRX*SUN__RAU ! Neg. RRX are in units of solar radii RRX = max(RRX,sngl(SUN__RAU)) ! Keep RRX > 1 solar radius C------- Prompts for time if (bForecast) then ! Forecast mode call Local2UT(iOffUT,iYr,Doy) ! calls the local computer time clock. With the offset calculates the Doy = universal time call Julian(10,iYr,Doy,JD,JEpoch) ! Current time in MJD write (*,'(/,10X,A,I6,A,I6,A,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('IPSD2020','I','StopPlay','StopPlay') C------- Make sure iYr,Doy match MJDref C Get heliographic coordinates for center of plot (Earth position) C MJDref = MJDref + 7 call Julian(11,iYr,Doy,MJDref,JEpoch) Idoy = int(Doy) HR = (Doy-float(Idoy))*24.0 IH = int(HR) IMN = nint((HR-float(IH))*60.0) call DATE_DOY(1,iYr,cMon,iM,iDay,iDoy) write(*,'(A,I2.2,A,I2.2,A,I2.2,A,A,I5,A,F8.3,A)') & ' (This MJD is ',IH,':',IMN,' UT ',iDay,' ',cMon,iYr,' -- DOY',Doy,')' call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) C------- Get modified Carrington variable for center of plot C [Note that EarthLng=(1-frac(XCEarth))*360] XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XCbeg = XCEarth-.5 call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) call T3D_set(T3D__TNOW, 0, XCbeg+0.5) XCtest1 = XCbeg-XCmargin ! Search limits XCtest2 = XCbeg+NC+XCmargin XCend= XCbeg + NC XCstrt = XCbeg XCbeg= XCbeg - xInc XCend= XCend + xInc MJDcntr = MJDref NC = 0 NCC = NC else ! Normal mode write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$' XCbeg = max(XClo,NCoff+XCbeg) C XCbeg = 1964.6 ! Bastille-day event C XCbeg = 2055.85 ! Used for EISCAT data run C XCbeg = 2061. C XCbeg = 2147. C XCbeg = 2153. C XCbeg = 2156.7 XCbeg = 2056. if(bGCamb) XCbeg = 1884 call AskR4('Carrington rotation (0 = STOP)?'//cStr,XCbeg) if (XCbeg .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay') NintF = 1 if(bForecast) then nTF = nTG - 1 NintF = 0 end if call AskI4('# Intermediate rotations interpolated',NintF) call AskI4('# Rotations averaged',NC) NC = min(NC,1+int(XChi-XCbeg) ) ! Stay within available data range NCC = NC XCbeg = XCbeg-NCoff call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) C------- Info only: start time of Carrington rotation call Julian(1,iYr,Doy,JDCar(int(XCbeg)),JEpoch) iDoy = Doy call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+int(XCbeg), & ' starts on date ',iDay,'-',cMon,'-',iYr,' ( Doy ',iDoy,' )' C------- Get modified Carrington variable for center of plot C Get heliographic coordinates for center of plot (Earth position) XCEarth = XCbeg+.5 ! Center of plot I = XCEarth ! Rotation containing plot center MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I)) ! Julian day for plot center call Julian(1,iYr,Doy,MJDcntr,JEpoch) ! iYr,Doy for plot center MJDcntr = MJDcntr-SUN__MJDtoJD ! Modified Julian day for plot center call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) ! Heliographic coord. for plot center XCtest1 = XCbeg-XCmargin ! Search limits XCtest2 = XCbeg+NC+XCmargin XCend= XCbeg + NC XCstrt = XCbeg XCbeg= XCbeg - xInc XCend= XCend + xInc I = XCtest2 MJDref = JDCar(I)-SUN__MJDtoJD+(XCtest2-I)*(JDCar(I+1)-JDCar(I)) ! Start MJD for data search end if C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCbeg call AskYN('Do you want to output model-source files?$no',bWriteMS) call AskYN('Do you want to output l.o.s. source surface maps?$no',bWritelos) call AskYN('Do you want one shift matrix per iteration?$no',bOnly1) call AskYN('Output good source files?$yes',bSource) call AskYN('Output three D error files?$no',bWrerr) NOneOr1 = -2 NOneOr2 = -2 NinterC = 3 ! Default first extract matrix NinterD = 3 ! Default forecast extract matrix call AskYN('Output 3D files with holes?$no',bWrite3D1) ! 3D 1 file write if(bWrite3D1) then call AskI4('Number of intermediate rotation densities with holes?',NinterC) ! Affects first extract call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr1) NOneOr1 = 0 if(bOneOr1) NOneOr1 = 1 end if call AskYN('Output 3D mostly filled files?$no',bWrite3D2) ! 3D "f" file write if(bWrite3D2) then call AskI4('Number of intermediate filled rotation densities?',NinterD) ! Affects forecast extract call AskYN('Remove a further 1/r beyond 1 AU?$no',bOneOr2) NOneOr2 = 0 if(bOneOr2) NOneOr2 = 1 end if nMapb = 1 nMapHRb = 1 nMapHRMb = nMap call AskYN('Output 3D high resolution files?$no',bWrite3D_HR)! 3D "h" file write if(bWrite3D_HR) then RRRS = RRS call AskR4('Minimum high-resolution height in AU?',RRRS) if(RRRS.le.RRS) then nMapb = 1 else nMapb = nint((RRRS-RRS)/dRR) + 1 end if if(nMapb.gt.nMap) nMapb = nMap ! Limit minimum height to maximum available RRRSM = 3.0 call AskR4('Maximum high-resolution height in AU?',RRRSM) nMapHR = nint((RRRSM-RRS)/dRR) + 1 if(nMapHR.gt.nMap) nMapHR = nMap ! Limit maximum height to maximum available if(nMapHR.lt.nMapb) nMapHR = nMapb ! Don't be stupid! call AskYN('Do you want the density error to limit HR densities?$no',bDdener) if(bDdener) then call ASKR4('What density L.O.S. crossing limit?',ERDLOS) call AskYN('Do you want the density error to limit HR velocities?$no',bVdener) end if call AskYN('Do you want the velocity error to limit HR velocities?$no',bVverer) if(bVverer) then bVdener = .FALSE. call ASKR4('What velocity L.O.S. crossing limit?',ERVLOS) call AskYN('Do you want the velocity error to limit HR densities?$no',bDverer) if(bDverer) bDdener = .FALSE. end if end if call AskYN('Output completely filled 3D files?$no',bWrite3Do) ! Completely filled 3D "o" file write call AskYN('Output a 3D high resolution base?$no',bWrite3D_HRb)! Completely filled 3D base "b" file write if(bWrite3D_HRb) then RRRRS = RRS call AskR4('Minimum high-resolution base height in AU?',RRRRS) if(RRRRS.le.RRS) then nMapHRb = 1 else nMapHRb = nint((RRRRS-RRS)/dRR) + 1 end if if(nMapHRb.gt.nMap) nMapHRb=nMap ! Limit maximum height to maximum available RRRRSM = 0.25 call AskR4('Maximum high-resolution base height in AU?',RRRRSM) nMapHRMb = nint((RRRRSM-RRS)/dRR) + 1 if(nMapHRMb.gt.nMap) nMapHRMb = nMap ! Limit maximum height to maximum available if(nMapHRMb.lt.nMapHRb) nMapHRMb = nMapHRb ! Don't be stupid! call AskYN('Do you want the den error to limit HR base densities?$no',bDdenerb) if(bDdenerb) then call ASKR4('What base density L.O.S. crossing limit?',ERDLOSB) call AskYN('Do you want the den error to limit HR base velocities?$no',bVdenerb) end if call AskYN('Do you want the vel error to limit HR base velocities?$no',bVvererb) if(bVvererb) then bVdenerb = .FALSE. call ASKR4('What base velocity L.O.S. crossing limit?',ERVLOSB) call AskYN('Do you want the vel error to limit HR base densities?$no',bDvererb) if(bDvererb) bDdenerb = .FALSE. end if end if C C Do you want to include a magnetic map at the reference surface height? C If yes, then read the source surface map (at 2.5Rs) at C the resolution of the reference arrays. bMagnetic = Bfield_choose_3(iVarp,cVar) call AskYN('Do you want to make magnetic source surface maps?$yes',bMagnetic) if (bMagnetic) then bMagnet = .TRUE. bMagnetic = .FALSE. do J=1, B3D__NN bWildMag(J) = .FALSE. call ForeignFile(-iVarp-1,cVar,B3D__STRCMD(J),cWildMag(J)) if (itrim(cWildMag(J)) .gt. 0) then bMagnetic = .TRUE. bWildMag(J) = .TRUE. end if end do if (bMagnet.and..not.bMagnetic) then print *, "You requested magnetic maps but didn't specify one correct path to find them." stop end if if (bMagnetic) call AskYN('Do you want to produce magnetic field 3D files?$yes',bmagexx) if (bMagnetic) call AskYN('Do you want to extract magnetic fields?$yes',bmagex) BFACTORRR = 0.0 if(bWildMag(1).or.bWildMag(5).or.bWildMag(9).or.bWildMag(13)) BFACTORRR = 1.0 ! radial open component factor BFACTORRC = 0.0 if(bWildMag(2).or.bWildMag(6).or.bWildMag(10).or.bWildMag(14)) BFACTORRC = 1.0 ! radial closed component factor BFACTORTC = 0.0 if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.1.0) then BFACTORRC = 2.0 if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 ! tangential closed component factor end if if(BFACTORRR.EQ.1.0.and.BFACTORRC.EQ.0.0) then BFACTORRR = 2.0 if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 end if if((BFACTORRR+BFACTORRC).GE.1.0) then if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 1.0 end if if(BFACTORRR.EQ.0.0.and.BFACTORRC.EQ.0.0) then if(bWildMag(3).or.bWildMag(7).or.bWildMag(11).or.bWildMag(15)) BFACTORTC = 2.0 end if BFACTORNC = 0.0 if(bWildMag(4).or.bWildMag(8).or.bWildMag(12).or.bWildMag(16)) BFACTORNC = 2.0 if (bMagnetic.and.BFACTORRR.ne.0.0) call AskR4('Multiplicative factor for CSSS fields?',BFACTORRR) if (bMagnetic.and.BFACTORRC.ne.0.0) call AskR4('Multiplicative factor for radial closed fields?',BFACTORRC) if (bMagnetic.and.BFACTORTC.ne.0.0) call AskR4('Multiplicative factor for tangential closed fields?',BFACTORTC) if (bMagnetic.and.BFACTORNC.ne.0.0) call AskR4('Multiplicative factor for normal closed fields?',BFACTORNC) if (bMagnetic.and.(BFACTORRR.ne.0.0.or.BFACTORRC.ne.0.0.or.BFACTORTC.ne.0.0.or.BFACTORNC.ne.0.0)) & call AskYN('Adjust field factors for possible known effects?$no',bAdjust) if(bAdjust) call Adjust_BF(bWildMag,B3D__NN,iYr,Doy,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC) end if call AskYN('Output an Earth parameter file?$no',bEWRITE) ! Earth extraction call AskYN('Output a Mercury parameter file?$no',bMWRITE) ! Mercury extraction call AskYN('Output a Venus parameter file?$no',bVEWRITE) ! Venus extraction call AskYN('Output a Mars parameter file?$no',bMAWRITE) ! Mars extraction if(iYr.lt.2009) then call AskYN('Output a Ulysses parameter file?$no',bUWRITE) ! Ulysses extraction end if if(iYr.ge.2007) then call AskYN('Output a STEREO A parameter file?$no',bSAWRITE) ! STEREO A extraction call AskYN('Output a STEREO B parameter file?$no',bSBWRITE) ! STEREO B extraction end if if(iYr.ge.2007) then call AskYN('Output a MESSENGER parameter file?$no',bMSWRITE) ! MESSENGER extraction end if call AskYN('Output a Comet 67P parameter file?$no',b67WRITE) ! Comet 67P/GC extraction 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 Venus parameter filled file?$no',bV2WRITEF) ! Venus extraction call AskYN('Output a Mars parameter filled file?$no',bMAWRITEF) ! Mars filled extraction if(iYr.lt.2009) then call AskYN('Output a Ulysses parameter filled file?$no',bUWRITEF) ! Ulysses filled extraction end if if(iYr.ge.2007) then call AskYN('Output a STEREO A parameter filled file?$no',bSAWRITEF) ! STEREO A filled extraction call AskYN('Output a STEREO B parameter filled file?$no',bSBWRITEF) ! STEREO B filled extraction end if if(iYr.ge.2007) then call AskYN('Output a MESSENGER parameter filled file?$no',bMXWRITEF) ! MESSENGER filled extraction end if call AskYN('Output a Comet 67P filled parameter file?$no',bC67WRITEF) ! Comet 67P/GC filled extraction C------- C The range of the Carrington variable plotted is [XCbeg,XCend]. The target range C [XCtest1,XCtest2] includes an extra safety margin to make sure that all relevant C data are read from file. NLV = 0 NLG = 0 do I=1,NLmaxV IYRSV(I) = 0 end do do I=1,NLmaxG IYRSGG(I) = 0 end do if (bVNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVEISC) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEVsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVGen) then do ii=1, iSyst if(iSys(ii).ne.0) then if(ii.eq.3) XElowV = XElowOV if(ii.eq.4) XElowV = XElowNV if(ii.eq.5) XElowV = XElowKV if(ii.eq.6) XElowV = XElowMV if(ii.eq.7) XElowV = XElowPV NLVB = NLV + 1 call ReadVIPSn8_n(iSys(ii),iReadGen8_n,iProcessGen8_n,cWildGens(ii), & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,DVOBS,GGOBS,DGGOBS,VFREQ,VSIZE, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,VFsav,VSsav) NLVE = NLV C if(ii.eq.4) then ! This is STELAB velocity DATA C do III=NLVB+10,NLVB+20 C print *, 'III, VFREQ(III)', III, VFREQ(III) C end do C end if if(ii.eq.6) then ! This is MEXART Velocity DATA print *, ' ' print *, 'MEXART speed data are input' IVMEX = 0 do III=NLVB,NLVE write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)') & 'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III) IVMEX = IVMEX + 1 end do write(*,'(A,I4,A)') 'There were ',IVMEX,' MEXART IPS V sources found' end if if(ii.eq.7) then ! This is Pushchino Velocity DATA print *, ' ' print *, 'Pushchino speed data are input' IVPUSH = 0 do III=NLVB,NLVE write(*,'(A,F5.0,A,F5.2,A,F7.0,A,F6.0,A,F6.2)') & 'FREQ =',VFREQ(III),' Size =',VSIZE(III),' VOBS =',VOBS(III),' DVOBS =',DVOBS(III),' Elong = ',XEV(III) IVPUSH = IVPUSH + 1 end do write(*,'(A,I4,A)') 'There were ',IVPUSH,' Pushchino IPS V sources found' end if end if end do C Sort the IPS velocity array input parameters C do i=1, 200 C write (*,'(A,I5,A10,I5,F8.4)') 'unsorted velocity array ',i, SRCV(I), IYRSV(I),DOYSV8(I) C end do call sort_IPS_data(NLV,NLMAXV,RADUMV8,cvgfiles,SRCVG,SRCV,IYRFV,IRECV, & XCV,YLV,DISTV,VOBS,DVOBS,VFREQ,VSIZE,XEV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8) C do i=1,200 C write (*,'(A,I5,A10,I5,F8.4)') 'sorted velocity array ',i, SRCV(I), IYRSV(I),DOYSV8(I) C end do end if print *, ' ' if (bGCamb) call ReadGIPS8('$CAM',nCar,JDCar, & iEdt,.FALSE.,bForecast, & XCtest1,XCtest2,MJDref,MJDfrst,MJDlast, & RRG,Speed,GPOWER,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG, & iXPG,iYPG,iMJDG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,GGOBS, & 1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav8,XDSsav,XLSsav, & XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGGen) then do ii=1, iSyst if(iSys(ii).ne.0) then if(ii.eq.3) XElowG = XElowOG if(ii.eq.4) XElowG = XElowNG if(ii.eq.5) XElowG = XElowKG if(ii.eq.6) XElowG = XElowMG if(ii.eq.7) XElowG = XElowPG NLGB = NLG + 1 call ReadVIPSn8_n(-iSys(ii),iReadGen8_n,iProcessGen8_n,cWildGens(ii), & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG,cggfiles,SRCGV,SRCGG,SRCGVsav,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,DVOBS,GGOBS,DGGOBS,GFFREQ,GSSIZE, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,DVVsav,GGsav,DGGsav,GFsav,GSSsav) NLGE = NLG do III=NLGB,NLGE if(ii.eq.3) then ! Ooty GGOBS(III) = (GGOBS(III)-1.0)*GLIMOF + 1.0 if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMOF end if if(ii.eq.4) then ! Nagoya GGOBS(III) = (GGOBS(III)-1.0)*GLIMNF + 1.0 if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMNF end if if(ii.eq.5) then ! KSWC GGOBS(III) = (GGOBS(III)-1.0)*GLIMKF + 1.0 if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMKF end if if(ii.eq.6) then ! MEXART GGOBS(III) = (GGOBS(III)-1.0)*GLIMMF + 1.0 if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMMF end if if(ii.eq.7) then ! Pushchino GGOBS(III) = (GGOBS(III)-1.0)*GLIMPF + 1.0 if(DGGOBS(III).gt.0.0) DGGOBS(III) = DGGOBS(III)*GLIMPF end if end do if(ii.eq.6) then ! This is MEXART g-level DATA print *, ' ' print *, 'MEXART g-level data are input' IGMEX = 0 do III=NLGB,NLGE write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') & 'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III) IGMEX = IGMEX + 1 end do write(*,'(A,I4,A)') 'There were ',IGMEX,' MEXART IPS g-level sources found' end if if(ii.eq.7) then ! This is Pushchino g-level DATA print *, ' ' print *, 'Pushchino g-level data are input' IGPUSH = 0 do III=NLGB,NLGE write(*,'(A,F5.0,A,F5.2,A,F7.3,A,F6.3,A,F6.2)') & 'FREQ =',GFFREQ(III),' Size =',GSSIZE(III),' GGOBS =',GGOBS(III),' DGGOBS =',DGGOBS(III),' Elong = ',XEGG(III) IGPUSH = IGPUSH + 1 end do write(*,'(A,I4,A)') 'There were ',IGPUSH,' Pushchino IPS g-level sources found' end if end if end do C Sort the IPS g-level array input parameters call sort_IPS_data(NLG,NLMAXG,RADUMG8,cggfiles,SRCGV,SRCGG,IYRFG,IRECG, & XCGG,YLGG,DISTGG,GGOBS,DGGOBS,GFFREQ,GSSIZE,XEGG,XLSGG,XCEGG,XLLGG,XDLGG,IYRSGG,DOYSGG8) C do i=1, 200 C write (*,'(A,I5,A10,I5,F8.4)') 'sorted g-level array ',i,SRCGV(I),IYRSGG(I),DOYSGG8(I) C end do end if print *, ' ' 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 print *, ' ' N_INV = 0 if (bVACE) call ReadACE8_v(cWildACE,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVACE2) call ReadACE28_v(cWildACE2,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVWind) call ReadWind8_v(cWildWind,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVCELIAS.and..not.bForecast) call ReadCELIAS8_v(cWildVCELIAS,nCar,JDCar,XCtest1,XCtest2,VVMAX,VVMIN,NCoff,RRV, & NLmaxV,NLV,N_INV,cvgfiles,SRCVG,SRCV, & IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS) if (bVCELIAS.and.bForecast) 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,SRCGV,SRCGG, & IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS) if (bDACE2) call ReadACE28_d(cWildACE2,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG, & IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS) if (bDWind) call ReadWind8_d(cWildWind,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG, & IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS) if (bDCELIAS.and..not.bForecast) call ReadCELIAS8_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG, & IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS) if (bDCELIAS.and.bForecast) call ReadCELIAS8R_d(cWildDCELIAS,nCar,JDCar,XCtest1,XCtest2,DMAX,DMIN,NCoff,RRG, & NLmaxG,NLmaxV,NLG,N_ING,cggfiles,SRCGV,SRCGG, & IYRSGG,DOYSGG8,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,GGOBS) call ArrR4Zero(NLmaxV,VSSIG) call ArrI4Constant(NLmaxV,1 ,NBSV ) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLmaxG,1 ,NBSG ) if(bforecast) write(*,'(A,F10.4,A,I5)') 'The Forecast cut-off time is DOY',Doy,' of',iYr write(*,'(A,4I10)') 'NLG, NLV, N_ING, N_INV',NLG,NLV,N_ING,N_INV LASTVVALUE = 0 ALASTVDOY = 0.0 LASTViYr = 0 do I=1,NLV C if (I.lt.100) print *, I, SRCV(I), SRCGG(I) if (VOBS(I) .gt. VLIMM) then NBSV(I) = 0 VOBS(I) = BadR4() end if if (VOBS(I) .lt. VLIML) then NBSV(I) = 0 VOBS(I) = BadR4() end if if(bforecast) then ! Provides a sharp cut-off in the velocity forecast lines of sight. if(IYRSV(I).ge.iYr.and.DOYSV8(I).gt.dble(Doy).or.IYRSV(I).eq.0) then if(LASTVVALUE.eq.0) then write(*,'(A,A,A,F10.4,A,I5)') 'The last valid IPS velocity LOS was ',SRCVLAST,' on',ALASTVDOY,' of',LASTViYr LASTVVALUE = 1 end if NBSV(I) = 0 VOBS(I) = BadR4() end if if(I.le.NLV.and.NBSV(I).ne.0) then ALASTVDOY = sngl(DOYSV8(I)) LASTViYr = IYRSV(I) SRCVLAST = SRCV(I) end if end if end do LASTGVALUE = 0 ALASTGDOY = 0.0 LASTGiYr = 0 do I=1,NLG if(GLIMF.ne.1.0) then ! g-level site specific multiplicative factor not specified GGOBS(I) = (GGOBS(I)-1.0)*GLIMF + 1.0 if(DGGOBS(I).gt.0.0) DGGOBS(I) = DGGOBS(I)*GLIMF end if if (GGOBS(I) .gt. GLIMP) then NBSG(I) = 0 GGOBS(I) = BadR4() end if if (GGOBS(I) .lt. GLIMM) then NBSG(I) = 0 GGOBS(I) = BadR4() end if if(bforecast) then ! Provides a sharp cut-off in the g-level forecast lines of sight. C print *, I, NLmaxG, IYRSGG(I), iYr, DOYSGG8(I), dble(Doy), LASTGVALUE if(IYRSGG(I).ge.iYr.and.DOYSGG8(I).gt.dble(Doy).or.IYRSGG(I).eq.0) then if(LASTGVALUE.eq.0) then write(*,'(A,A,A,F10.4,A,I5)') 'The last valid IPS g-level LOS was ',SRCGLAST,' on',ALASTGDOY,' of',LASTGiYr LASTGVALUE = 1 end if NBSG(I) = 0 GGOBS(I) = BadR4() end if if(I.le.NLG.and.NBSG(I).ne.0) then ALASTGDOY = sngl(DOYSGG8(I)) LASTGiYr = IYRSGG(I) SRCGLAST = SRCGG(I) end if end if end do LASTAVALUE = 0 ALASTADOY = 0.0 LASTAiYr = 0 do I=NLV+1,NLmaxV if(bforecast) then ! Provides a sharp cut-off in the in-situ V forecast measurements. if(IYRSV(I).ge.iYr.and.DOYSV8(I).gt.dble(Doy)) then if(LASTAVALUE.eq.0) then write(*,'(A,A,A,F10.4,A,I5)'), 'The last valid in-situ velocity was ',SRCALAST,' on',ALASTADOY,' of',LASTAiYr LASTAVALUE = 1 end if NBSV(I) = 0 VOBS(I) = BadR4() end if if(I.le.NLV+N_INV.and.NBSV(I).ne.0) then ALASTADOY = sngl(DOYSV8(I)) LASTAiYr = IYRSV(I) SRCALAST = SRCV(I) end if end if end do LASTAVALUE = 0 ALASTADOY = 0.0 LASTAiYr = 0 NINSITUBEND = 0 do I=NLG+1,NLmaxG if(bforecast) then ! Provides a sharp cut-off in the in-situ D forecast measurements. if(IYRSGG(I).ge.iYr.and.DOYSGG8(I).gt.dble(Doy)) then if(LASTAVALUE.eq.0) then write(*,'(A,A,A,F10.4,A,I5)') 'The last valid in-situ density was ',SRCALAST,' on',ALASTADOY,' of',LASTAiYr LASTAVALUE = 1 end if NBSG(I) = 0 NINSITUBEND = NINSITUBEND + 1 C write(*,*) 'GGOBS(I)', GGOBS(I), ' IYRSGG(I)',IYRSGG(I), ' DOYSGG8(I)', DOYSGG8(I) GGOBS(I) = BadR4() end if if(I.le.NLG+N_ING.and.NBSG(I).ne.0) then ALASTADOY = sngl(DOYSGG8(I)) LASTAiYr = IYRSGG(I) SRCALAST = SRCV(I) end if end if if (GGOBS(I) .gt. DMAX) then NBSG(I) = 0 GGOBS(I) = BadR4() end if if (GGOBS(I) .lt. DMIN) then NBSG(I) = 0 GGOBS(I) = BadR4() end if end do if(NINSITUBEND.gt.0) write(*,'(A,I5)') 'The number of in-situ density points beyond the cut-off =',NINSITUBEND if(bFACE) write(*,'(A,F8.3)') 'The 1 AU upper limit of the in-situ density is:',DMAX if(bFACE) write(*,'(A,F8.3)') 'The 1 AU lower limit of the in-situ density is:',DMIN C======= II = 0 III = 0 do I=1,NLG+N_ING II = II + 1 if(bGCamb) then XCG (NLG-II+1) = XCGG(I) YLG (NLG-II+1) = YLGG(I) DISTG(NLG-II+1) = DISTGG(I) GOBS (NLG-II+1) = GGOBS(I) DGOBS (NLG-II+1) = DGGOBS(I) GFREQ (NLG-II+1) = GFFREQ(I) GSIZE (NLG-II+1) = GSSIZE(I) XEG (NLG-II+1) = XEGG(I) XLSG (NLG-II+1) = XLSGG(I) XCEG (NLG-II+1) = XCEGG(I) XLLG (NLG-II+1) = XLLGG(I) XDLG (NLG-II+1) = XDLGG(I) DOYSG8(NLG-II+1) = DOYSGG8(I) IYRSG(NLG-II+1) = IYRSGG(I) XLLGD5 = XLLG(NLG-I+1)/5.0 XDLGD5 = XDLG(NLG-I+1)/5.0 if(XLLGD5.gt.0.) then XLLGD5 = XLLGD5 + 0.5 else XLLGD5 = XLLGD5 - 0.5 end if if(XLLGD5.gt.0.) then XDLGD5 = XDLGD5 + 0.5 else XDLGD5 = XDLGD5 - 0.5 end if IXLLGD5 = XLLGD5 IXDLGD5 = XDLGD5 write(SRCGS(1:5),'(I5)') IXLLGD5*5 write(SRCGS(6:9),'(I4)') IXDLGD5*5 SRCGG(NLG-II+1) = SRCGS end if if(bGNago.or.bGOoty.or.bGGen) then XCG(II) = XCGG(I) YLG(II) = YLGG(I) DISTG(II) = DISTGG(I) GOBS(II) = GGOBS(I) DGOBS(II) = DGGOBS(I) GFREQ(II) = GFFREQ(I) GSIZE(II) = GSSIZE(I) XEG(II) = XEGG(I) XLSG(II) = XLSGG(I) XCEG(II) = XCEGG(I) XLLG(II) = XLLGG(I) XDLG(II) = XDLGG(I) DOYSG8(II) = DOYSGG8(I) IYRSG(II) = IYRSGG(I) if(bGGen.and.GFFREQ(I).eq.140.0.and.NDUP.ne.0) then ! MEXART and other data do IDUP=1,NDUP-1 III = III + 1 II = II + 1 XCG(II) = XCGG(I) YLG(II) = YLGG(I) DISTG(II) = DISTGG(I) GOBS(II) = GGOBS(I) DGOBS(II) = DGGOBS(I) GFREQ(II) = GFFREQ(I) GSIZE(II) = GSSIZE(I) XEG(II) = XEGG(I) XLSG(II) = XLSGG(I) XCEG(II) = XCEGG(I) XLLG(II) = XLLGG(I) XDLG(II) = XDLGG(I) DOYSG8(II) = DOYSGG8(I) IYRSG(II) = IYRSGG(I) end do end if end if end do NLG = NLG + III ALng = 1.*(nLng-1)/(iLng-1) C print *, 'Spot check of IPS data - 30', IYRSG(30),IYRSV(30),DOYSG8(30),DOYSV8(30),XEG(30),XEV(30),GOBS(30),VOBS(30) if(bForecast) then ! For forecast mode NTV = nint((float(NTV) + 3.0)/2.0) + 7 NTG = nint((float(NTG) + 3.0)/2.0) + 7 end if print *, ' ' write(*,'(A,2I5)') ' Before MkTimes NTV, NTG ',NTV, NTG write(*,'(2(A,I9))') ' NLV =',NLV,' NLG =',NLG print *, ' ' NLVVV = 0 do I=1,NLV if (VOBS(I) .ge. VLIML .and. VOBS(I) .le. VLIMM) NLVVV = NLVVV + 1 end do do I=NLV+1,NLV+N_INV if (VOBS(I) .ge. VLIML .and. VOBS(I) .le. VLIMM) NLVVV = NLVVV + 1 end do if (NLVVV .eq. 0) then print *, 'There are no valid IPS or in-situ velocity measurements. Set bVcon to false.' bVcon = .FALSE. end if NLGGG = 0 do I=1,NLG if (GOBS(I) .ge. GLIMM .and. GOBS(I) .le. GLIMP) NLGGG = NLGGG + 1 end do do I=NLG+1,NLG+N_ING if (GOBS(I) .ge. DMIN .and. GOBS(I) .le. DMAX) NLGGG = NLGGG + 1 end do if (NLGGG .eq. 0) then print *, 'There are no valid g-level or in-situ density measurements. Set bGcon to false.' bGcon = .FALSE. end if if(.not.bVcon.and..not.bGcon) then print *, ' STOP!!!!!!!' print *, ' There are no sources for either a velocity' print *, ' or a g-level deconvolution during this interval' stop end if print *, 'Before MkTimes MJDcntr =', MJDcntr print *, ' ' ItoffG = 0 ItoffV = 0 IdaybegV = 365 IdaybegG = 365 IYRBV = IYRSV(1) IYRBG = IYRSG(1) if(NLV.eq.0) bVcon = .FALSE. if(bVcon) & call MkTimes(bForecast,MJDcntr,NmidHRV,aNdayV,NTmaxV,ALng,nLng,nCar,JDCar, & 0,NTV,IYRBV,IdaybegV,XCintDV,XCintV,XCbeV,XCtbegV,XCtendV,IYRV,DOYV) if(NLG.eq.0) bGcon = .FALSE. if(bGcon) & call MkTimes(bForecast,MJDcntr,NmidHRG,aNdayG,NTmaxG,ALng,nLng,nCar,JDCar, & 0,NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG) C print *, 'SRCVG(30)',SRCVG(30) C print *, 'Spot check of IPS data - 30', IYRSG(30),IYRSV(30),DOYSG8(30),DOYSV8(30),XEG(30),XEV(30) C print *, 'SRCVG(409)',SRCVG(409) C print *, 'Spot check of IPS data - 409', IYRSG(409),IYRSV(409),DOYSG8(409),DOYSV8(409),XEG(409),XEV(409) write (*,'(A,4I5)') ' After MkTimes IYRBV, IYRBG, IdaybegV, IdaybegG',IYRBV,IYRBG,IdaybegV,IdaybegG if(.not.bGcon) then nTG = nTV IYRBG = IYRBV ItoffG = ItoffV XCtbegG = XCtbegV XCtendG = XCtendV do I=1,nTmaxG XCintG(I) = XCintV(I) XCintDG(I) = XCintDV(I) XCbegG(1,I) = XCbeV(1,I) XCbegG(2,I) = XCbeV(2,I) end do end if if(.not.bVcon) then nTV = nTG IYRBV = IYRBG ItoffV = ItoffG XCtbegV = XCtbegG XCtendV = XCtendG do I=1,nTmaxV XCintV(I) = XCintG(I) XCintDV(I) = XCintDG(I) XCbeV(1,I) = XCbegG(1,I) XCbeV(2,I) = XCbegG(2,I) end do end if if(bOnly1) then nTV = nTG ItoffV = ItoffG XCtbegV = XCtbegG XCtendV = XCtendG do I=1,nTmaxV ! The temporal spaces for XCintV, etc. better be the same as for XCintG XCintDV(I) = XCintDG(I) XCintV(I) = XCintG(I) XCbeV(1,I) = XCbegG(1,I) XCbeV(2,I) = XCbegG(2,I) end do end if write(*,'(A,2F12.5,I12)') ' XCtbegV, XCtendV, NTV', XCtbegV, XCtendV, NTV write(*,'(A,2F12.5,I12)') ' XCtbegG, XCtendG, NTG', XCtbegG, XCtendG, NTG write(*,'(/A)') ' XCbeV =' write(*,'(8F9.6)') XCbeV write(*,'(/A)') ' XCbeGG =' write(*,'(8F9.6)') XCbeGG write(*,'(/A)') ' XCintV =' write(*,'(8F9.6)') XCintV write(*,'(/A)') ' XCintG =' write(*,'(8F9.6)') XCintG write(*,'(/A)') ' XCintDV =' write(*,'(8F9.4)') XCintDV write(*,'(/A)') ' XCintDG =' write(*,'(8F9.4)') XCintDG XCtbegF = XCtbegG ! For forecast mode XCtendF = XCtendG ! For forecast mode if(NLVVV.lt.NLVmin.and.N_INV.lt.NLVmin) then bVcon = .FALSE. NLV = 0 write(*,'(A,I3,A)') 'There are too few (',NLVVV+N_INV,') sources for a velocity deconvolution. NLV set to zero.' print *, ' ' end if if(NLGGG.lt.NLGmin.and.N_ING.lt.NLGmin) then bGcon = .FALSE. NLG = 0 write(*,'(A,I3,A)') 'There are too few (',NLGGG+N_ING,') sources for a g-level deconvolution. NLG set to zero.' print *, ' ' end if if(.not.bVcon.and..not.bGcon) then print *, ' STOP!!!!!!!' write(*,'(2(A,I9))') 'NLV =',NLV,' NLG =', NLG print *, ' There are too few sources for either a velocity' print *, ' or a g-level deconvolution' print *, ' You will need to set the lower limit of either the velocity' print *, ' or g-level LOS (or both) to a smaller number.' stop end if if(bGcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP) if(bVcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),VMAP) NLmax = max(NLmaxG,NLmaxV) NLVBEG = 0 NLGBEG = 0 NLVEND = NLV NLGEND = NLG print *, 'NLVEND =', NLVEND, 'NLGEND =',NLGEND NTVV = ItoffV NTGG = ItoffG NTVE = 0 NTGE = 0 K = 0 Nbadvel = 0 Nbadden = 0 do KKK=1,NLmax+nTG ! Loop before deconvolution to set up VMAP, DMAP K = K + 1 IVSW = 0 ! Switch turn on for NTVV value of K if(K.le.NLVEND+1) then ! If K is more than the source # + 1 don't turn on switch C print *, ' Before V switch loop K, NTVV, NLmax', K,NTVV,NLmax,XCEV(K),IYRSV(NTVV+2),DOYV(NTVV+2),XCintV(NTVV+2) if(bVcon.and.(XCEV(K).gt.XCintV(NTVV+2).and.XEV(K).le.XEhighV).or.(XEV(K).le.XEhighV.and.K.ge.NLV)) then NTVV = NTVV + 1 NTVB = NTVE + 1 NTVE = K - 1 if(NTVB.ne.NTVE+1) IVSW = 1 if(K.le.NLVEND) K = K - 1 NTVD = NTVE-NTVB+1 if(NTTV.eq.1) NLVBEG = NTVD if(NTVD.ne.0) & write (*,'(A,I4,A,2I6,2F10.3)') 'V loop ',NTVV,' - S#, K, XCEV, XCintV', NTVD,K,NCoff+XCEV(K),NCoff+XCintV(NTVV+2) if(NTVV.eq.NTV) NLVEND = NTVE end if end if if(IVSW.eq.1) then C print *, ' V switch ',IVSW,' loop',NTVV if (NLVEND .eq. 0) call Say('IPSD2020','E','NoV','no velocity IPS data') if (bVproxy) then ! Read proxy point-P V map from file I = iReadProxyMapN(1,NCoff+XCintV(NTVV),nLng,nLat,VMAP(1,1,NTVV)) call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP',NTVV,' min and max after read is', amin, amax end if call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),aminV,amaxV) if (amaxV.eq.BadR4()) then ! Try to calculate a proxy point-P V map call MapGrid(XCbeV(1,NTVV),XCbeV(2,NTVV),YLbeg,YLend,NC, & -NTVD,XCV(NTVB),YLV(NTVB),VOBS(NTVB),NAV,DAV,nLng,nLat, & VMAP(1,1,NTVV),VDEV(1,1,NTVV),VPNT(1,1,NTVV), & Vmin,Vmax,VDmin,VDmax) C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintV(NTVV) end if ! Calculate weights, WTSV, along all lines of sight C call MkLOSWeightsm(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) call MkLOSWeightsx(NLOSWV,dLOSV,NLOSV,NTVD,VFREQ(NTVB),VSIZE(NTVB),XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAPHOLE(1,1,NTVV)) if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,NTVV),CONRV,3,0.0,0.0) end if end if call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,NTVV,' min and max is', amin, amax if(amax.ne.BadR4()) Nbadvel = Nbadvel + 1 if(NTVV.eq.ItoffV+1) then ! Fill earlier arrays do N=1,ItoffV amaxV = BadR4() if(bVproxy) then I = iReadProxyMapN(1,NCoff+XCintV(N),nLng,nLat,VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),aminV,amaxV) print *, 'The VMAP' ,N,' min and max after read is', aminV, amaxV end if if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,N,' min and max is', amin, amax if(amax.ne.BadR4()) then call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAP(1,1,N)) call Say('IPSD2020','I','Info','Copy first source array into data file ' & //cStr(:Int2Str(N,cStr))) end if end if end do end if end if IGSW = 0 ! Switch turn on for NTGG value of K if(K.le.NLGEND+1) then ! If K is more than the source # +1 don't turn on switch if(K.gt.0) then C print *, ' Before G switch loop K, NTGG, NLmax', K,NTGG,NLmax,XCEG(K),IYRSG(NTGG+2),DOYG(NTGG+2),XCintG(NTGG+2) if(bGcon.and.(XCEG(K).gt.XCintG(NTGG+2).and.XEG(K).le.XEhighG).or.(XEG(K).le.XEhighG.and.K.ge.NLG)) then NTGG = NTGG + 1 NTGB = NTGE + 1 NTGE = K - 1 if(NTGB.ne.NTGE+1) IGSW = 1 if(K.le.NLGEND) K = K - 1 NTGD = NTGE-NTGB+1 if(NTGG.eq.1) NLVBEG = NTGD if(NTGD.ne.0) & write (*,'(A,I4,A,2I6,2F10.3)') 'G loop ',NTGG,' - S#, K, XCEG, XCintG', NTGD,K,NCoff+XCEG(K),NCoff+XCintG(NTGG+2) if(NTGG.eq.NTG) NLGEND = NTGE end if end if end if if(IGSW.eq.1) then if(NTGB.ne.0) then if (NLGEND .eq. 0) call Say('IPSD2020','E','NoG','no G-level IPS data') if (bDproxy) then ! Read proxy density map from file I = iReadProxyMapN(2,NCoff+XCinTG(NTGG),nLng,nLat,DMAP(1,1,NTGG)) end if call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),aminDMAP,amaxDMAP) if (amaxDMAP.eq.badR4()) then ! Try to calculate a point-P synoptic G map call MapGrid(XCbeGG(1,NTGG),XCbeGG(2,NTGG),YLbeg,YLend,NC, & -NTGD,XCG(NTGB),YLG(NTGB),GOBS(NTGB),NAV,DAV,nLng,nLat, & DMAP(1,1,NTGG),GDEV(1,1,NTGG),GPNT(1,1,NTGG), & Gmin,Gmax,GDmin,GDmax) C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCintG(NTGG) end if ! Calculate weights, WTSG, along all lines of sight C call MkLOSWeightsm(NLOSWG,dLOSG,NLOSG,NTGD,XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) call MkLOSWeightsx(NLOSWG,dLOSG,NLOSG,NTGD,GFREQ(NTGB),GSIZE(NTGB),XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),PPGMAP(1,1,NTGG)) ! Convert G to density if(amaxDMAP.eq.badR4()) then call ConvertG2D(0,nLngLat,PPGMAP(1,1,NTGG),DMAP(1,1,NTGG),RRG,DEN1AU,PWG) call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,NTGG),CONRD,3,0.0,0.0) end if end if call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) print *, 'The DMAP',NTGG,' min and max is', amin, amax if(amax.ne.BadR4()) Nbadden = Nbadden + 1 C======== call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAPHOLE(1,1,NTGG)) ! Dummy density map for print out call ArrR4TimesConstant(-nLngLat,DMAPHOLE(1,1,NTGG),RRS*RRS,DMAPHOLE(1,1,NTGG)) ! n*R^2 with holes possible if(NTGG.eq.ItoffG+1) then ! Fill earlier arrays do N=1,ItoffG amaxD = BadR4() if(bDproxy) then I = iReadProxyMapN(2,NCoff+XCintG(N),nLng,nLat,DMAP(1,1,N)) call arrR4getminmax(nLngLat,DMAP(1,1,N),aminD,amaxD) print *, 'The DMAP',N,' min and max after read is', aminD, amaxD end if if(amaxD.eq.BadR4()) then call arrR4getminmax(nLngLat,DMAP(1,1,NTGG),amin,amax) print *, 'The DMAP',N,' min and max is', amin, amax if(amax.ne.BadR4()) then call ArrR4Copy(nLngLat,DMAP(1,1,NTGG),DMAP(1,1,N)) call Say('IPSD2020','I','Info','Copy first source array into data file ' & //cStr(:Int2Str(N,cStr))) end if end if end do end if end if end if end do if(Nbadvel .eq. 0 .and. .not. bVproxy) then ! If all else fails, fill the velocity arrays with a valid values do N=1,nTV do I=1,Nlng do J=1,Nlat VMAP(I,J,N) = 400.0 end do end do end do print *, ' ' print *, 'The velocity arrays all needed to be filled with valid values.' print *, ' ' end if if(Nbadden .eq. 0 .and. .not. bDproxy) then ! If all else fails, fill the density arrays with a valid values do N=1,nTG do I=1,Nlng do J=1,Nlat DMAP(I,J,N) = DEN1AU*R1AU/RRG end do end do end do print *, ' ' print *, 'The density arrays all needed to be filled with valid values.' print *, ' ' end if if(bVFACE) then do I=1,N_INV WTSV_in(I) = 10.0 end do end if if(bFACE) then do I=1,N_ING WTSD_in(I) = 10.0 end do end if if(.not.bVcon) then do N=1,NTV call ArrR4Constant(nLngLat,Speed,VMAP(1,1,N)) call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N)) call ArrR4Copy(nLngLat,VMAP(1,1,N),DMAPV(1,1,N)) end do call arrR4getminmax(nLngLat,VMAP(1,1,NTV),amin,amax) print *, 'All',NTV,' VMAPs beginning with .not. bVcon are', amax end if if(.not.bGcon) then do N=1,NTG call ArrR4Constant(nLngLat,1.,DMAP(1,1,N)) call ArrR4Copy(nLngLat,DMAP(1,1,N),PPGMAP(1,1,N)) call ConvertG2D(0,nLngLat,PPGMAP(1,1,N),DMAP(1,1,N),RRG,DEN1AU,PWG) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) call ArrR4Copy(nLngLat,DMAP(1,1,N),VMAPD(1,1,N)) end do call arrR4getminmax(nLngLat,DMAP(1,1,NTG),amin,amax) print *, 'All',NTG,' DMAPs beginning with .not. bGcon are', amax end if NLGG = 0 NLGGL = 0 do L=1,NLG if(L.lt.NLGBEG.or.L.gt.NLGEND) then GOBS(L) = BadR4() NBSG(L) = 0 end if if(GOBS(L) .ne. BadR4() .and. GOBS(L) .ne.0.0) 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) .ne. BadR4()) then NLVV = NLVV + 1 NLVVL = L end if end do if(NLGG.lt.NLGmin) then print *, ' ' write(*,'(A,I3,A)') 'There are too few g-level (',NLGG,') remote sources for a good deconvolution.' write(*,'(A)') 'All remote sources set bad.' do L=1,NLG NBSG(L) = 0 end do NLGG = 0 end if NLVV = 0 NLVVL = 0 do L=1,NLV if(L.lt.NLVBEG.or.L.gt.NLVEND) then VOBS(L) = BadR4() NBSV(L) = 0 end if if(VOBS(L) .ne. BadR4() .and. VOBS(L) .ne. 0.0) then NLVV = NLVV + 1 NLVVL = L end if end do if(NLVV.lt.NLVmin) then print *, ' ' write(*,'(A,I3,A)') 'There are too few velocity (',NLVV,') remote sources for a good deconvolution.' write(*,'(A)') 'All remote sources set bad.' do L=1,NLV NBSV(L) = 0 end do NLVV = 0 end if NLGGG = 0 NLGGGL = 0 do L=NLG+1,NLG+N_ING if(GOBS(L) .ne. BadR4() .and. GOBS(L) .ne. 0.0) then NLGGG = NLGGG + 1 NLGGGL = L end if end do NLVVV = 0 NLVVVL = 0 do L=NLV+1,NLV+N_INV if(VOBS(L) .ne. BadR4() .and. VOBS(L) .ne. 0.0) then NLVVV = NLVVV + 1 NLVVVL = L end if end do if(bGcon.and.NLGGL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last IPS g-value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') & 'I =',NLGGL,' Year =',IYRSG(NLGGL),' DOY =',DOYSG8(NLGGL),' IPS g-level =',GOBS(NLGGL) end if if(bFACE.and.NLGGGL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last good in-situ density value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') & 'I =',NLGGGL,' Year =',IYRSG(NLGGGL),' DOY =',DOYSG8(NLGGGL),' In-situ density =',GOBS(NLGGGL) end if if(bVcon.and.NLVVL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last good IPS velocity value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') & 'I =',NLVVL,' Year =',IYRSV(NLVVL),' DOY =',DOYSV8(NLVVL),' IPS velocity =',VOBS(NLVVL) end if if(bVFACE.and.NLVVVL.ne.0) then print *, ' ' write (*,'(1X,A)') 'The last good in-situ velocity value was:' write (*,'(1X,A,I10,A,I5,A,F10.5,A,F8.3)') & 'I =',NLVVVL,' Year =',IYRSV(NLVVVL),' DOY =',DOYSV8(NLVVVL),' In-situ velocity =',VOBS(NLVVVL) end if do N=1,NTV do I=1,Nlng do J=1,Nlat VMAP(I,J,N) = 400.0 end do end do end do print *, ' ' call Say('ipsdt_20n_inp','I','Info',' # good g-level observations ' & //cStr(:Int2Str(NLGG,cStr))) call Say('ipsdt_20n_inp','I','Info',' # good V observations ' & //cStr(:Int2Str(NLVV,cStr))) call Say('ipsdt_20n_inp','I','Info',' # good in-situ D observations ' & //cStr(:Int2Str(NLGGG,cStr))) call Say('ipsdt_20n_inp','I','Info',' # good in-situ V observations ' & //cStr(:Int2Str(NLVVV,cStr))) 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) if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,DMAP,DMAPV) else call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD) end if if(bVproxy) call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV) if(bDproxy) then ! Do this for use after velocity iteration if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,nTV,nTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV) end if call FillWholeT(1,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD) end if if(.not.bVproxy) then call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV) do N=1,NTV C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N',amin,amax,N C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,0.0) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,0.0) C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N',amin,amax,N C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) end do end if call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N', amin,amax,N call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,1,0.0,0.0) end do C C Comment out call MkDensity(nLng,nLat,nMap,XCshift,DDfact,VMAP,RRV,dRR) C Nit = -1 ! Why three iterations counters?? NitV = -1 NitG = -1 do while (Nit .lt. NiterT) Nit = Nit+1 print *, ' ' call Say('IPSD2020','I','Info','Iteration '//cStr(:Int2Str(Nit,cStr))) C write(*,'(A,I3,4F7.3)') 'nit, CONRV, CONVT, CONSTV CONSTMV',nit, CONRV, CONVT, CONSTV, CONSTMV NVS = 0 if (Nit .eq. LimitoV) NVS = 1 NGS = 0 if (Nit .eq. LimitoG) NGS = 1 if (bVcon) then ! V deconvolution NitV = NitV+1 ! # V iterations ! bGmv2=.TRUE., implies bGcon=.FALSE. if (bGmv2) then RRD = RRV call Veleq3(nLng,nLat,NTV,VMAP,SPEEDEQ) print *, 'Solar equatorial mv^2 = constant value SPEEDEQ =',SPEEDEQ do N=1,NTV call Mk_V2DN(nLngLat,VMAP(1,1,N),DMAPV(1,1,N),SPEEDEQ,DEN1AU,RRV) end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,DMAPV,DMAP) else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,NTV,NTG,ConsTV,DMAPV,VMAPD,Dtmpt,VVT,VVV) call ArrR4Copy(nLngLatnTG,VMAPD,DMAP) end if if(NitV .eq. NiterT) then do N=1,NTG call ArrR4Copy(nLngLat,VMAPD(1,1,N),DMAPHOLE(1,1,N)) end do end if end if call MkPostd_in(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM, & N_INV,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in) C print *, 'Spot check after Mkpostd_V - 1,30', ITIM(1,1),ITIM(1,30),IYRSV(1),IYRSV(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30) if(NitV.eq.NiterT.and.bWritelos) & call writelosmapvg_in(0,NLOSVP1,NLV,N_INV,nTV,NBSV,NCOFF,ITIM,ITIM_in,XLON,XLON_in,XLAT,XLAT_in) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DVfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,VFAC) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DVfact,N_INV,XCtim_in,XLONP_in,XLAT_in,RPX_in,VFAC_in) call MkVModeltd_in(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX, & WTSV,DFAC,VFAC,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSV_in,DFAC_in,VFAC_in,PWV, & NLV,N_INV,NLOSV,NLOSVP1,VMAP,VLIM, & DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij,VWTij_in) do I=1,NLV+N_INV C if(I.le.NLV) print *, 'Before NLV',I,VMOD(I),VOBS(I),WTSV(20,I),NBSV(I),VWTij(20,I) C if(I.gt.NLV) print *, 'Through NLV',I,VMOD(I),VOBS(I),WTSV_in(I-NLV),NBSV(I),VWTij_in(I-NLV) if(WTSV(1,I).eq.0.and.I.le.NLV) then NBSV(I) = 0 VMOD(I) = Speed end if if(WTSV_in(I-NLV).eq.0.and.I.gt.NLV) then NBSV(I) = 0 VMOD(I) = Speed end if end do C if(nit.eq.18) stop C if(VMOD(I).eq.BadR4().or.VOBS(I).eq.BadR4().and.I.gt.NLV) then C NBSV(I) = 0 C else C if(I.le.NLV) print *, 'Before NLV',I,VMOD(I),VOBS(I),WTSV(20,I),NBSV(I),VWTij(20,I) C if(VMOD(I).gt.2000..or.VMOD(I).lt.100.) print *, I,VMOD(I),VOBS(I) C if(I.gt.NLV) print *, 'Through NLV',I,VMOD(I),VOBS(I),WTSV_in(I-NLV),NBSV(I),VWTij_in(I-NLV) C NBSV(I) = 1 C end if C end do C if(I.gt.367.and.I.le.NLV) then C if(WTSV(1,I).eq.0.and.I.le.NLV) then C print *, 'Before the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I) C NBSV(I) = 0 C end if C if(WTSV_in(I-NLV).eq.0.and.I.gt.NLV) then C print *, 'Before the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I) C NBSV(I) = 0 C end if C end do call FixModeltdn(1,NLV+N_INV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) C do I=1,NLV+N_INV C if(I.gt.367.and.I.le.NLV) then C print *, 'After the fix is:', I,VOBS(I),VMOD(I),NBSV(I),VSIG,FIX(I),VSSIG(I) C end if C end do if(NitV .eq. 0 .or. NitV .eq. NiterT) then do I=1,NLV+N_INV,24 if(I.le.NLV) then C print *, 'Observed and model V:', I,' ',SRCV(I),VOBS(I),VMOD(I),FIX(I),NBSV(I) else C print *, 'Observed and model V:', I,' ',SRCV(I),VOBS(I),VMOD(I),FIX(I),NBSV(I),XLON_in(I-NLV),XLAT_in(I-NLV),XCtim_in(I-NLV) end if end do end if if(bWriteMS.and.NitV.eq.NiterT) then call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV+N_INV,PWV,VRATIO, & NLmaxV,NSRV,SRCV,DOYSV8,XCEV,VOBS,VMOD,XLON,XLproj,XCtpr,NLOSVP1,FIX,VSSIG,NBSV) end if VRATION(NitV+1) = VRATIO ! Put in 7/9/00 if(NitV.ge.2.and.VRATION(NitV+1).gt.(VRATION(NitV)+(VRATION(NitV+1)*0.05)).and. & VRATION(NitV+1).gt.(VRATION(NitV-1)+(VRATION(NitV+1)*0.10))) then print *, 'THE SUMMED VELOCITY RATIO IS NO LONGER CONVERGING' C NitV = NiterT end if if(NVREST.eq.1) then call MkVMaptdn0n_in(LON,LAT,ITIM,RPX,LON_in,LAT_in,ITIM_in,RPX_in,XEV,VWTij,VWTij_in, & FIX,NLV,N_INV,NLOSV,NLOSVP1, & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VSIGB,NVS,NBSV,VMAPSIG, & ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,WTPV_in, & Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp) end if if(NVREST.gt.1) then call MkVMaptdn0n_inm(LON,LAT,ITIM,RPX,LON_in,LAT_in,ITIM_in,RPX_in,XEV,VWTij,VWTij_in, & FIX,NLV,N_INV,NLOSV,NLOSVP1, & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VSIGB,NVS,NBSV,VMAPSIG, & ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV,WTPV_in, & Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp) end if do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C if(amax.eq.BadR4()) then C print *, 'VMAP is bad before map determination N =', N,' Iter =', Nit C end if C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'In V determination VMAP min max',amin,amax if (NitV .eq. NiterT) call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N)) call ArrR4SetMinMax(-nLngLat,VMAP(1,1,N),VLL,VUL) end do end if if (.not.bDproxy.or.Nit.ge.1) then if (bGcon .or. bVmv2) then ! Update density and density shifts even on the last call FillWholeT(0,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTG,ConsTG,DMAP,DDD) do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before Grids MkShiftdn0n DMAP min max N',amin,amax,N AEmax = AEmaxD/(RRS*RRS) AEmin = AEminD/(RRS*RRS) call FillMapOLS(2,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,Speed,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD,3,0.0,0.0) C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkShiftdn0n DMAP min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTG,DMAP, 0,1.0*CONDT/aNdayG,0.,Dtmp) do N=1,NTG call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) end do end if if (bGcon) then C if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then C call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) C else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV, & nLng,nLat,NTV,NTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV) C end if do N=1,NTG C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before FillWholeT VMAPD min max',amin,amax end do call FillWholeT(1,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,VMAPD,DDD) do N=1,NTG C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before Grids MkShiftdn0n VMAPD min max N',amin,amax,N call FillMapL(2,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) call FillMapL(0,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAPD(1,1,N),CONRD,3,0.0,0.0) C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before Timesmooth and MkShiftdn0n VMAPD min max N',amin,amax,N end do call Timesmooth(nLng,nLat,nTG,VMAPD,0,1.0*CONDT/aNdayG,0.,Dtmp) do N=1,NTG call FillMapL(3,XCbeGG(1,N),XCbeGG(2,N),nLng,nLat,ConstL,VMAPD(1,1,N)) C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before MkShiftdn0n VMAPD min max N',amin,amax,N end do end if end if if(bGcon) then call MkShiftdn0n(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG,Speed, & nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C print *, VMAPD (30,5,20),VMAPD (31,5,20),VMAPD (30,6,20),VMAPD (31,6,20),VMAPD (30,5,21),VMAPD (31,5,21), C & DMAP (30,5,20),DMAP (31,5,20),DMAP (30,6,20),DMAP (31,6,20),DMAP (30,5,21),DMAP (31,5,21), C & XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11), C & DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11), C & DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11) end if if (bGcon) then ! G deconvolution NitG = NitG+1 ! # G iterations ! bVmv2=.TRUE., implies bVcon=.FALSE. if (bVmv2) then RRV = RRD call Deneq3(nLng,nLat,NTG,DMAP,RRD,DENEQ) print *, 'Solar equatorial mv^2 = constant value DENEQ =',DENEQ do N=1,NTG call Mk_D2VN(nLngLat,DMAP(1,1,N),VMAPD(1,1,N),DENEQ,Speed,RRD) end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAPD,VMAP) else call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTG,VMAPD,DMAPV,Vtmpt,DDT,DDD) call ArrR4Copy(nLngLatnTV,DMAPV,VMAP) end if if(NitG .eq. NiterT) then do N=1,NTV call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAPHOLE(1,1,N)) end do end if end if call MkPostd_in(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM, & N_ING,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in) C call MkPostd_in(XCbeV,XCtbegV,XCtendV,IdaybegV,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, C & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV8, C & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM, C & N_INV,XLON_in,XLONP_in,XLproj_in,XLAT_in,RPX_in,XCtpr_in,XCtim_in,LON_in,LAT_in,ITIM_in) C call MkPostd(XCbeGG,XCtbegG,XCtendG,IdaybegG,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, C & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG8, C & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) C print *, 'Spot check after Mkpostd_G - 1,30', ITIM(1,1),ITIM(1,30),IYRSG(1),IYRSG(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30) if(NitG.eq.NiterT.and.bWritelos) & call writelosmapvg_in(1,NLOSGP1,NLG,N_ING,nTG,NBSG,NCOFF,ITIM,ITIM_in,XLON,XLON_in,XLAT,XLAT_in) C if(NitG.eq.NiterT.and.bWritelos) call writelosmapvg(1,NLOSGP1,NLG,nTG,NBSG,NCOFF,ITIM,XLON,XLAT) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,N_ING,XCtim_in,XLONP_in,XLAT_in,RPX_in,DFAC_in) C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkGModeltdn DMAP min max N',amin,amax,N C end do call MkGModeltd_in(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,PWG,DEN1AU,XCtpr_in,XLON_in,XLproj_in,RPX_in,WTSD_in,DFAC_in, & NLG,N_ING,NLOSG,NLOSGP1,DMAP, & nLng,nLat,nTG,RRS,GMOD2,GWTij,GWTij_in,WTPD,GMWTi) C call MkGModeltdn(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, C & WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP, C & nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi) do I=1,NLG+N_ING if(I.le.NLG) then if(GOBS(I).eq.badR4()) then GOBSS(I) = badR4() NBSG(I) = 0 else GOBSS(I) = GOBS(I) if(bFACE) then if(DGOBSS.le.1.0) then if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)-(1.0-DGOBSS) end if else if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)-(1.0-DGOBSS) end if if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if end if end if end if else if(GOBS(I).eq.badR4()) then GOBSS(I) = badR4() NBSG(I) = 0 else GOBSS(I) = (GOBS(I)/(DEN1AU*((1.0/RPX_in(I-NLG))**2)))**(2.0*PWG) if(bFACE) then if(DGOBSS.le.1.0) then if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)-(1.0-DGOBSS) end if else if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)-(1.0-DGOBSS) end if if(GOBSS(I).ge.GLIMM.and.GOBSS(I).le.1.0.and.NBSG(I).ne.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if end if end if end if C print *, GOBSS(I),GOBS(I),DEN1AU,RPX_in(I-NLG),PWG end if if(WTSG(1,I).eq.0.and.I.le.NLG) then NBSG(I) = 0 GMOD2(I) = 1.0 end if if(WTSD_in(I-NLG).eq.0.and.I.gt.NLG) then NBSG(I) = 0 GMOD2(I) = 1.0 end if end do call FixModeltdn(2,NLG+N_ING,GOBSS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) C call FixModeltdn(2,NLG,GOBSS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) if(NitG .eq. 0 .or. NitG .eq. NiterT) then do I=1,NLG+N_ING,24 if(I.le.NLG) then C print *, 'Observed and model g:', I,' ',SRCGG(I),GOBSS(I),GMOD2(I),FIX(I),NBSG(I) else C print *, 'Observed and model g:', I,' ',SRCGG(I),GOBSS(I),GMOD2(I),FIX(I),NBSG(I),XLON_in(I-NLG),XLAT_in(I-NLG),XCtim_in(I-NLG) end if end do end if if(bWriteMS.and.NitG.eq.NiterT) then call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG+N_ING,PWG,GRATIO, & NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GOBSS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG) end if C if(bWriteMS.and.NitG.eq.NiterT) then C call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO, C & NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GOBSS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG) C end if GRATION(NitG+1) = GRATIO ! Put in 7/9/00 if(NitG.ge.2.and.GRATION(NitG+1).gt.(GRATION(NitG)+(GRATION(NitG+1)*0.05)).and. & GRATION(NitG+1).gt.(GRATION(NitG-1)+(GRATION(NitG+1)*0.10))) then print *, 'THE SUMMED DENSITY RATIO IS NO LONGER CONVERGING' C Nit = NiterT end if C print *, 'Into MkDmaptdn0n_in',nLng,nLat,nTG,NLOSG,NLG if(NGREST.eq.1) then call MkDMaptdn0n_in(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in, & PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) end if if(NGREST.gt.1) then call MkDMaptdn0n_inm(LON,LAT,ITIM,GWTij,LON_in,LAT_in,ITIM_in,GWTij_in, & PWG,FIX,NLG,N_ING,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,DLIM,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD,WTPD_in, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) end if C print *, 'Out of MkDmaptdn0n_in' C call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, C & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, C & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, C & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'After MkDmaptd0 DMAP min max',amin,amax if (NitG .eq. NiterT) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG)) end do end if if(Nit.ne.NiterT) then ! Update velocity shifts and vmap except on last C C If in-situ density is available define DEN1AU except on the last if bFACE is yes C DEN1 = 0.0 ADEN1 = 0.0 if(bFACE) then do I=NLG+1,NLG+N_ING if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then DEN1 = DEN1 + GOBS(I) ADEN1 = ADEN1 + 1.0 end if end do DEN1AUO = DEN1AU DEN1AU = DEN1/ADEN1 if(abs(DEN1AU-DEN1AUO).gt.0.000001) print *, 'The old 1 AU density',DEN1AUO,' has been changed to',DEN1AU end if if(bFACE.and.Nit.eq.NiterT) print *, 'The final 1 AU density has been set to',DEN1AU if(bGACE) then DGOBSSO = DGOBSS DGOBSS = 0.0 ADGOBSS = 0.0 C do I=1,NLG C if(GOBS(I).le.GLIMP.and.GOBS(I).ge.GLIMM.and.NBSG(I).ne.0) then C DGOBSS = DGOBSS + GOBS(I) C ADGOBSS = ADGOBSS + 1.0 C end if C end do do I=NLG+1,NLG+N_ING if(GOBS(I).le.DMAX.and.GOBS(I).ge.DMIN.and.NBSG(I).ne.0) then DGOBSS = DGOBSS + (GOBS(I)/(DEN1AU*((1.0/RPX_in(I-NLG))**2)))**(2.0*PWG) ADGOBSS = ADGOBSS + 1.0 end if end do DGOBSS = DGOBSS/ADGOBSS if(abs(DGOBSS-DGOBSSO).gt.0.000001) print *, 'The old G-level mean',DGOBSSO,' has been changed to',DGOBSS end if if(bGACE.and.Nit.eq.NiterT) print *, 'The final G-level mean of 1.0 has been set to',DGOBSS if (bVcon.or.bGmv2) then call FillWholeT(1,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) call FillMaptN(0,XCbeV,nLng,nLat,nTV,ConsTV,VMAP,VVV) C if(niTV.eq.2) then C print *, 'Line 2828 iteration 2 stop' C stop C end if do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'After MkDmaptdn0 VMAP min max N',amin,amax,Ncd IPS call FillMapOLS(2,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,Speed,ConstL,AEmaxV,AEminV,AEangV,VMAP) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,3,0.0,0.0) end do call Timesmooth(nLng,nLat,nTV,VMAP, 0,1.0*CONVT/aNdayV,0.,VZtmp) ! prior to 2016 BVJ C call Timesmooth(nLng,nLat,nTV,VMAP, 2,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) end do end if if (.not.bOnly1) then if (bVcon) then C if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then C call ArrR4Copy(nLngLatnTG,DMAP,DMAPV) C else call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD) C end if call FillWholeT(0,NitV,NiterT,XCbeV,nLng,nLat,nTV,ConsTMV,DMAPV,VVV) do N=1,NTV C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'After MkDmaptdn0 DMAPV min max N',amin,amax,N call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call FillMapL(0,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) call GridSphere2D(ALng,nLng,nLat,1,DMAPV(1,1,N),CONRV,3,0.0,0.0) end do call Timesmooth(nLng,nLat,nTV,DMAPV,0,1.0*CONVT/aNdayV,0.,VZtmp) ! prior to 2016 BVJ C call Timesmooth(nLng,nLat,nTV,DMAPV,2,1.0*CONVT/aNdayV,0.,VZtmp) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) end do call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C print *, VMAP (30,5,20),VMAP (31,5,20),VMAP (30,6,20),VMAP (31,6,20),VMAP (30,5,21),VMAP (31,5,21), C & DMAPV (30,5,20),DMAPV (31,5,20),DMAPV (30,6,20),DMAPV (31,6,20),DMAPV (30,5,21),DMAPV (31,5,21), C & XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11), C & DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11), C & DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11) end if end if end if end do ! End of iteration loop C ********************************************************************* C The below is needed at least for header writes in the NV3 files C C call T3D_set (T3D__CLIPLNG , 0, 90.0 ) C call T3D_iset(T3D__ITERATIONS , 0, NiterT ) C call T3D_iset(T3D__NCOFF , 0, NCoff ) C call T3D_set (T3D__RR , 0, RRS ) C call T3D_set (T3D__DRR , 0, dRR ) C call T3D_set (T3D__TT , 0, TT ) C call T3D_set (T3D__DTT , 0, dTT ) C call T3D_iset(T3D__NLNG , 0, nLng ) C call T3D_iset(T3D__NLAT , 0, nLat ) C call T3D_iset(T3D__NRAD , 0, nMap ) C call T3D_iset(T3D__NTIM , 0, nTim ) C call T3D_set (T3D__PWN_V , 0, PWV ) C call T3D_set (T3D__PWN_G , 0, PWG ) C call T3D_set (T3D__PWR_V , 0, PWRV ) C call T3D_set (T3D__PWR_G , 0, PWRG ) C call T3D_set (T3D__D1AU , 0, DEN1AU ) C call T3D_set (T3D__SMOOTH_V , 0, CONRV ) C call T3D_set (T3D__SMOOTH_D , 0, CONRD ) C call T3D_set (T3D__FILL_V , 0, CONSTV ) C call T3D_set (T3D__FILL_D , 0, CONSTG ) C call T3D_set (T3D__SMOOTH_TIME_V , 0, CONVT ) C call T3D_set (T3D__SMOOTH_TIME_D , 0, CONDT ) call T3D_iset(T3D__NL_V , 0, NLV ) call T3D_iset(T3D__NL_G , 0, NLG ) call T3D_iset(T3D__NLOS_V , 0, NLOSV ) call T3D_iset(T3D__NLOS_G , 0, NLOSG ) call T3D_set (T3D__DLOS_V , 0, dLOSV ) call T3D_set (T3D__DLOS_G , 0, dLOSG ) call T3D_set (T3D__BINX_V , 0, ANLIMITV ) call T3D_set (T3D__BINX_D , 0, ANLIMITG ) C ********************************************************************* if(bSource) then ! Write out good sources used if(NLG.eq.NLV) then print *, 'NLG equals NLV for writegoodsource. Both =',NLG istation = 1 if(bGNago) istation = 1 if(bGOoty) istation = 2 C call writegoodsourcegv(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NBSG,NLV,N_INV) if(bGgen) call writegoodsources(NLMAXV,SRCVG,NBSV,NBSG,NLV,N_INV) if(.not.bGgen) call writegoodsourcegv(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NBSG,NLV,N_INV) else if(NLG.ne.NLV.and.bGgen) then Print *, 'Sorry, this sourcewrite option for V is not yet supported' else print *, 'NLG does not equal NLV for writegoodsource. NLG =',NLG,' NLV =',NLV istation = 1 if(bGNago) istation = 1 if(bGOoty) istation = 2 call writegoodsourceg(istation,NLMAXG,cggfiles,SRCGV,NBSG,NLG,N_ING) end if if(NLG.ne.NLV.and.bVgen) then Print *, 'Sorry, this sourcewrite option for g is not yet supported' else print *, 'Out of sourcewrite for g' if(bVEISC) istation = 1 if(bVNago) istation = 1 if(bVOoty) istation = 2 call writegoodsourcev(istation,NLMAXV,cvgfiles,SRCVG,NBSV,NLV,N_INV) end if end if end if J = 0 ! List bad IPS sources NGOODSDVAL = 0 do I=1,NLV+N_INV if (NBSV(I) .eq. 0) then J = J+1 if(VOBS(I).eq.BadR4()) VOBS(I) = -9999.9 write (*,'(3A,I5,A,F7.1)') ' Bad V ', SRCV(I),' source removed, # =', & I,', elongation =',XEV(I) end if if(I.gt.NLV.and.NBSV(I).eq.1) then if(NGOODSDVAL.eq.0) FIRSTDOY = sngl(DOYSV8(I)) NGOODSDVAL = NGOODSDVAL + 1 FINALDOY = sngl(DOYSV8(I)) end if end do write (*,'(A,I8,A)') ' There were',J,' bad V sources removed.' if(N_INV.ne.0) then write (*,'(A,I9,2(A,F6.1))') & ' There are',NGOODSDVAL," good in-situ velocities. 1'st DOY =",FIRSTDOY ,' Final DOY =',FINALDOY end if J = 0 NGOODSDVAL = 0 do I=1,NLG+N_ING if (NBSG(I).eq.0) then J = J+1 if(GOBSS(I).eq.BadR4()) GOBSS(I) = -99.999 write (*,'(3A,I5,A,F7.3)') ' Bad G ',SRCGG(I),' source removed, # =', & I,', elongation =',XEG(I) end if if(I.gt.NLG.and.NBSG(I).eq.1) then if(NGOODSDVAL.eq.0) FIRSTDOY = sngl(DOYSG8(I)) NGOODSDVAL = NGOODSDVAL + 1 FINALDOY = sngl(DOYSG8(I)) end if end do write (*,'(A,I8,A)') ' There were',J,' bad G sources removed.' if(N_ING.ne.0) then write (*,'(A,I9,2(A,F6.1))') & ' There are',NGOODSDVAL," good in-situ densities. 1'st DOY =",FIRSTDOY ,' Final DOY =',FINALDOY end if C Here VMAPHOLE contains the original point-P map (bVcon=.FALSE.) ( or is equal to VMAP do N=1,NTV call ArrR4Copy(nLngLat,VMAPHOLE(1,1,N),DMAPV(1,1,N)) ! Place into a dummy array end do C Changed 9/21/01 BVJ C if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then ! Do this if smoother V maps needed C call ArrR4Copy(nLngLatnTV,DMAPV,VMAP) ! In any case, replace dummy array C else C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, C & nTV,nTV,ConsTV,DMAPV,VMAP,Vtmpt,VVT,VVV) C end if call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,DMAPV,VMAP,Vtmpt,VVT,VVV) do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam C XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region (Change 1/18 BVJ) C XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/27.275) ! Mid region of interest on map call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,Speed,ConstL,AEmaxV,AEminV,AEangV,VMAP) C call FillMapOLS(XCbeV(1,N),XCbeV(2,N),XCmidV,nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 180 deg bad spots call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy end if end do do N=1,NTG call ArrR4Copy(nLngLat,DMAPHOLE(1,1,N),VMAPD(1,1,N)) ! Place into a dummy array end do C Changed below 9/21/01 BVJ C if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then ! Do this if smoother D maps needed C call ArrR4Copy(nLngLatnTG,VMAPD,DMAP) ! Replace dummy array C else C call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, C & nTG,nTG,ConsTG,VMAPD,DMAP,Vtmpt,DDT,DDD) C end if call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,VMAPD,DMAP,Vtmpt,DDT,DDD) do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam C XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region (Changed 1/18/02) C XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/27.275) ! Mid region of interest on Map AEmax = AEmaxD/(RRS*RRS) AEmin = AEminD/(RRS*RRS) call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,Speed,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy end if end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV) end if C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,nLng,nLat, C & NTG,NTV,ConsTG,DMAP,DMAPV,Vtmpt,DDT,DDD) print *, ' Before the first Extractdvd',XCbegG(1,1),XCtbegG,XCstrt C***** Take data out at different spacecraft (and earth) locations (bExtract tells where) NinterD = 3 call Extractdvd(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C call Extractdnn(bForecast,bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat,nMap,nTG,NinterD,FALLOFFD, C & VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C***** Write to disk of results interpolated to the even latitudes, longitudes C and times over the heights (specified as nMapn) of the XCshift model. nTF = nTG nMapn = nMap Nit = Nit + 1 cPre = 'nv3d' C The first write3D_infotd makes maps that have holes to be used to show Carrington maps call write3D_infotd3D_3(NOneOr1,cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT) do I=1,4*nTmaxG ! Will not work for Extractd3D if XCbegG's are different XCbe(1,I) = XCbegG(1,1) XCbe(2,I) = XCbegG(2,1) end do print *, 'Out of write3D_infotd3D_3' C if(bWrerr) then ! Changed BVJ 18/5/11 C These (mode 2) write3D_infotd calls make 3D arrays with time error maps to show reconstruction reliability do N=1,nTG do I=1,nLng do J=1,nLat if(ANMAPD(I,J,N).eq.0.0) ANMAPD(I,J,N) = badR4() VDEV(I,J,N) = ANMAPD(I,J,N) if(AWTMAPD(I,J,N).eq.0.0) AWTMAPD(I,J,N) = badR4() GDEV(I,J,N) = AWTMAPD(I,J,N) end do end do end do call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,VDEV,ANMAPD,Vtmpt,DDT,DDD) call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,GDEV,AWTMAPD,Vtmpt,DDT,DDD) PIRSQD = (CONDT/aNdayG)*3.14159*((CONRD/(20.0/float(NGRESS)))**2) ! normalize weights and lines of sight to spatial resolution do N=1,nTG call arrR4getminmax(nLngLat,ANMAPD(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,ANMAPD(1,1,N),NSIDE,Zmin,Zmax) ! Remove small ANMAPD holes call ArrR4TimesConstant(-nLngLat,ANMAPD(1,1,N),PIRSQD,ANMAPD(1,1,N)) end if call arrR4getminmax(nLngLat,AWTMAPD(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,AWTMAPD(1,1,N),NSIDE,Zmin,Zmax) ! Remove small AWTMAPD holes) call ArrR4TimesConstant(-nLngLat,AWTMAPD(1,1,N),PIRSQD,AWTMAPD(1,1,N)) end if end do do N=1,nTV do I=1,nLng do J=1,nLat if(ANMAPV(I,J,N).eq.0.0) ANMAPV(I,J,N) = badR4() VDEV(I,J,N) = ANMAPV(I,J,N) if(AWTMAPV(I,J,N).eq.0.0) AWTMAPV(I,J,N) = badR4() GDEV(I,J,N) = AWTMAPV(I,J,N) end do end do end do call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,VDEV,ANMAPV,Vtmpt,VVT,VVV) call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,GDEV,AWTMAPV,Vtmpt,VVT,VVV) PIRSQV = (CONVT/aNdayV)*3.14159*((CONRV/(20.0/float(NVRESS)))**2) ! normalize weights and lines of sight to spatial resolution do N=1,nTV call arrR4getminmax(nLngLat,ANMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,ANMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small ANMAPV holes call ArrR4TimesConstant(-nLngLat,ANMAPV(1,1,N),PIRSQV,ANMAPV(1,1,N)) end if call arrR4getminmax(nLngLat,AWTMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,AWTMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small AWTMAPV holes call ArrR4TimesConstant(-nLngLat,AWTMAPV(1,1,N),PIRSQV,AWTMAPV(1,1,N)) end if end do if(bWrerr) then ! Changed BVJ 18/5/11 NinterCC = 0 cPre = 'er1_' print *, 'Before the error1 (L.O.S. crossings) write' call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,ANMAPV,ANMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) print *, 'Before the error2 (projected weights) write' NinterCC = 0 cPre = 'er2_' call write3D_infotd3DM_3(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,AWTMAPV,AWTMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) end if if(.not. bForecast) then print *, ' Before Extractd3D',XCbegG(1,1),XCtbegG,XCstrt call Extractd3D(bExtract,RRS,dRR,XCbe,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterC+1,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar) end if call Timesmooth(nLng,nLat,nTV,VMAP,2,CONVT/aNdayV,0.75*CONVT/aNdayV,VZtmp) ! Fill in times one time away do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy SpeedV = 0.0 call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam end if call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV/2.0,0,0.0,0.0) ! Change 5/9/08 to smooth polar holes end do call Timesmooth(nLng,nLat,nTG,DMAP,2,CONDT/aNdayG,0.75*CONDT/aNdayG,Dtmp) ! Fill in times one time away do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy AEmax = AEmaxD/(RRS*RRS) AEmin = AEmin/(RRS*RRS) SpeedD = 0.0 call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD/2.0,0,0.0,0.0) ! Change 5/9/08 to smooth polar holes end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV) end if C The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views C The first write3D_infotd makes maps that have holes to be used to show Carrington maps cPre = 'nv3f' call write3D_infotd3D_3(NOneOr2,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT) C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change if(bWrite3D_HR) then print *, 'Before second write3D_infotd at High Resolution' cPre = 'nv3h' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = 3 ! In original map coordinates (so nMapHR = 1.5 AU). if(nMapHR*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHR = nMap ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHR) nMapb = nMapHR print *, RRRS, RRRSM, nMapb, nMapHR C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR_3(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapb,nMapHR,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdener,bVdener,bVverer,bDverer,ERDLOS,ERVLOS,ANMAPD,ANDMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if C Must change back to regular value in current version RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change if(bForecast.and..not.bmagex) then C print *, 'DATA EXTRACTED FROM FINAL MATRIX FOR FORECAST' print *, ' Before 2nd Extractdvd',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if ************************************************************************************************************************************* C Produce a set of nv3o files that are completely filled as well as we can make them. C The third write3D_infotd makes maps that have no holes to be used to provide the best possible boundary for a 3D MHD model if(bWrite3Do.or.bWrite3D_HRb.or..not.bForecast) then C ixxxxxxx = 0 C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop print *, 'Before third filling' do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP N, amin and amax = ', N, amin, amax if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP N, amin and amax = ', N, amin, amax if(amax.ne.BadR4()) then call GridFill(0,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Fill ALL velocity holes in any map with a hole C call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy C C SpeedV = 0.0 C call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP) C call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam end if end if end do call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) ! Fill all the velocity times do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP N, amin and amax = ', N, amin, amax end do print *, 'Before forth filling' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(0,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Fill ALL density holes in any map with a hole C AEmax = AEmaxD/(RRS*RRS) C AEmin = AEmin/(RRS*RRS) (needs fixing) C AEmin = AEminD/(RRS*RRS) (fixed 2/29/2008 BVJ) C SpeedD = 0.0 C call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP) C call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy C call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if end if end do call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) ! Fill all the density times print *, 'After forth filling' if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1).and.XCtbegG.eq.XCtbegV) then print *, 'Copy array from VMAP to VMAPD directly' call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else C print *, 'VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78) call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV) C print *, 'VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78) end if C end if !ixxxxxxx loop print *, 'Before third write3D_infotd' C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'DMAP N, amin and amax = ', N, amin, amax C end do C do N=1,NTG C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'VMAPD N, amin and amax = ', N, amin, amax C end do if(bWrite3Do) then C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change cPre = 'nv3o' NinterD = 3 C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_3(0,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) C end if !ixxxxxxx loop C Change the density back since the falloff at distance is produced in Extractdvd RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change end if if(bWrite3D_HRb) then C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change print *, ' ' print *, 'Before third write3D_infotd of High Resolution base' cPre = 'nv3b' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = 3 nMapb = nMapHRb ! In original map coordinates (so nMapHRMb = 0.25 AU). if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb print *, RRRRS, RRRRSM, nMapb, nMapHRMb C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR_3(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANDMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) C Change the density back since the falloff at distance is produced elsewhere RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change end if end if C************************************************************************************************************************************* if ( bMagnetic ) then C inexc = 0 C if(inexc.eq.0) go to 4567 dTT = (XCtendG-XCtbegG)/(sngl((nTF-1)*(NinterD+1))) ! Time increment call T3D_set (T3D__TT , 0, sngl(XCtbegG)) ! Needed for header write in Write3D_bbtm_3. call T3D_set (T3D__DTT , 0, dTT ) call T3D_iset(T3D__NCOFF , 0, NCoff ) ! Needed to pass these values, and for the header write in Write3D_bbtm_3 call T3D_set (T3D__RR , 0, RRS ) call T3D_set (T3D__DRR , 0, dRR ) call T3D_iset(T3D__NLNG , 0, nLng ) call T3D_iset(T3D__NLAT , 0, nLat ) call T3D_iset(T3D__LAT , 0, -90 ) call T3D_iset(T3D__LAT+1 , 0, +90 ) call T3D_iset(T3D__NRAD , 0, nMap ) call T3D_iset(T3D__NTIM , 0, nTim ) C 4567 continue call T3D_iset(T3D__NCOFF , 0, NCoff ) ! Must pass these values to Bfield_Get_3 in Get_bbtm_3 call T3D_set (T3D__RR , 0, RRS ) call T3D_iset(T3D__NLNG , 0, nLng ) call T3D_iset(T3D__NLAT , 0, nLat ) call T3D_iset(T3D__NRAD , 0, nMap ) call T3D_iset(T3D__NTIM , 0, nTim ) call T3D_iset(T3D__LAT , 0, -90 ) call T3D_iset(T3D__LAT+1 , 0, +90 ) write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC call get_bbtm_3(XCbeGG,nLng,nLat,nTG,RRS,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN, & XCINTG,BRR2DT,BRC2DT,BTC2DT,BNC2DT,iBflag) if(bmagexx) then print *, ' ' print *, 'Before write3D_bbtm_HR_3 after HR base' NintLng = NintLn NintLat = NintLa NintHt = NintH Ninter = 3 ClipLng = 90.0 nMapbb = nMapHRb ! In original map coordinates (so nMapHRMb = 0.25 AU). if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb 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 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 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 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 end if end do write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC C BFACTORRR = 1.0 C BFACTORRC = 1.0 C BFACTORTC = 1.0 C BFACTORNC = 0.0 C do N=1,nTF+1 C call ArrR4Constant(nLngLat,BFACTORRR,BRR2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORRC,BRC2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORTC,BTC2DT(1,1,N)) C call ArrR4Constant(nLngLat,BFACTORNC,BNC2DT(1,1,N)) C end do C write(*,'(A,4F6.3)') 'BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC',BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC call write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG,iYrBG, & RRS,dRR,nLng,nLat,nMapbb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & ClipLng,VMAP,XCshift,DVfact, & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR) end if if(bmagex) then write(*,'(A,F10.6,F12.7,F10.6)') ' Before Final Extractdvdm_3',XCbegG(1,1),XCtbegG,XCstrt C ifieldw2 = 1 ifieldw2 = 0 if(ifieldw2.ne.0) then print *, ' ' write (*,'(A,100(I4,E12.5))') 'BRR2DT', (L,BRR2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BRC2DT', (L,BRC2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BTC2DT', (L,BTC2DT(10,4,L),L=1,ntG) print *, ' ' print *, ' ' write (*,'(A,100(I4,E12.5))') 'BNC2DT', (L,BNC2DT(10,4,L),L=1,ntG) end if print *, ' ' NinterD = 3 call Extractdvdm_3(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift, & DVfact,DDfact,TT3D,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar) end if end if if(.not. bMagnetic) then print *, ' Before Final Extractdvd',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdvd(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if call Say('IPSD2020','S','Stop','program successfully completed') stop end