C+ C NAME: C ipstd_20n4r_inp_mag3_v18_mhd_intel C PURPOSE: C Construct deconvolved synoptic maps for the IPS observations using C a "time-dependent" technique that projects all observations to a C single lower reference surface. Note - uneven time cadences are C allowed for velocity and density data sets. This version incorporates C 3-component velocities, magnetic fields, and density and temperature C into the mkshiftdn routines. C CATEGORY: C Data processing C CALLING SEQUENCE: C run ipstd_20n_inp_mag3_v18_mhd_intel 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_20n4r_inp_mag3_v18_mhd 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 = 121, ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) (61) (121) & 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.025) ! Resolution (AU) in radial direction (0.1) (0.05) (0.025) 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 & OFIXG (NLmaxG), ! Old Fix to G model to make OK & OFIXV (NLmaxV), ! Old Fix to V 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 NiterT & 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), & ND(4) 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) & BMAPR (nLng,nLat,nTmaxG), & BMAPT (nLng,nLat,nTmaxG), & BMAPN (nLng,nLat,nTmaxG), & 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), & XCsh (3), & PP (4), & Scale (4) C Addition for mkshiftdnma.f, the subroutine that provides an input of volumes and traceback from an external source. real VV3 (nLng,nLat,NMap,nTmaxG,3), ! Input local volumetric velocities & BB3 (nLng,nLat,NMap,nTmaxG,3), ! Input local volumetric magnetic fields & DD1 (nLng,nLat,NMap,nTmaxG), ! Input local volumetric densities & TT1 (nLng,nLat,NMap,nTmaxG), ! Input local volumetric temperatures & SSS1 (nLng,nLat), ! Scratch array 1 for volumetric magnetic fields & SSS2 (nLng,nLat), ! Scratch array 2 for volumetric magnetic fields & SSS3 (nLng,nLat), ! Scratch array 3 for volumetric magnetic fields & VVV3 (nLng,nLat,NMap,nTmaxG,3), ! Input external volumetric velocities & BBB3 (nLng,nLat,NMap,nTmaxG,3), ! Input external volumetric magnetic fields & DDD1 (nLng,nLat,NMap,nTmaxG), ! Input external volumetric densities & TTT1 (nLng,nLat,NMap,nTmaxG), ! Input external volumetric temperatures & RADS (NMap+1), ! Input external heights (should be same as tomographic heights) & XCshift3(nLng,nLat,nMap,nTmaxG,3), ! Map of accumulated shifts at heights (from 1 RR) (in long, lat, time) & XCshiftM(nLng,nLat,nTmaxG,3), ! Map of one shift at heights (from 1 RRMS to RR) (in long, lat, time) & Vratio3 (nLng,nLat,nMap,nTmaxG,3), ! Map of accumulated velocity ratios at heights (from 1 RR) & Bratio3 (nLng,nLat,nMap,nTmaxG,3), ! Map of accumulated magnetic field ratios at heights (from 1 RR) & Dratio1 (nLng,nLat,nMap,nTmaxG), ! Map of accumulated density ratios at heights (from 1 RR) & Tratio1 (nLng,nLat,nMap,nTmaxG), ! Map of accumulated temperature ratios at heights (from 1 RR) & XLLTtmp (nLng,nLat,nTmaxG,8) ! Scratch array 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, & bENLIL, bMSFLUKSS, bH3DMHD, bSmooth, & 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 & bMagnetic_ext /.TRUE./, ! Magnetic map available for iteration for an external source iteration & bMagnet /.FALSE./, ! Magnetic map check & bmagexx /.FALSE./, ! Produce 3D magnetic field data files & bmagex /.FALSE./, ! Extract magnetic field data as well as velocity and density data & bmagexx_ext /.FALSE./, ! Produce 3D magnetic field data files from and external source & bmagex_ext /.FALSE./, ! Extract magnetic field data as well as velocity and density data from external & 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 /.FALSE./, & bVdenerb /.FALSE./, C & bVdenerb /.TRUE./, C & bVvererb /.FALSE./, & bVvererb /.FALSE./, & bDvererb /.FALSE./, & bWrite3D_HR /.TRUE./, & bWrite3D_HRb /.FALSE./, & bWrite3D_ext /.FALSE./, & bWrite3Do /.TRUE./, & bWriteMS /.FALSE./, & bWritelos /.FALSE./, & bOneOr1 /.TRUE./, & bOneOr2 /.TRUE./, & bForecast, & bFcst /.FALSE./, & b3DMHD /.TRUE./, & bCont /.TRUE./, & bcheck /.FALSE./, & bcheck2 /.FALSE./, & bcheck3 /.FALSE./, & bzeroth /.TRUE./, & 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 logical BField_Choose ! Logical function 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'/ character cPrefixfnim(NTP)*2 /'M1','Vm','Em','M4','Jm','Sm','U7','Nm','Pm','Um','Ha','Hb','S1','S2','Mm', & '6m'/ character cPrefixfnic(NTP)*2 /'M1','Vc','Ec','M4','Jc','Sc','U7','Nc','Pc','Uc','Ha','Hb','S1','S2','Mm', & '6m'/ character cPrefixx(NTP)*3 /'m1_','v2_','e3_','m4_','j5_','s6_','u7_','n8_','p9_','uy_','ha_','hb_','s1_','s2_','mx_', & '67_'/ character cPrefixm(NTP)*3 /'m1m','v2m','e3m','m4m','j5m','s6m','u7m','n8m','p9m','uym','ham','hbm','s1m','s2m','mxm', & '67m'/ character CNAMEF*15 ! pLACED HERE TEMPORARILY 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 /48/, ! Fixed especially for the MHD runs & NTG /48/ ! Fixed especially for the MHD runs 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/, & NiterMHD /3/, & ireadmhd /10/, ! # number of mhd input files needed to be read for an attempt 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 & NBSVF /0/, ! # of bad velocity model values & NBSGF /0/ ! # of bad g-level model values 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/, & 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 / 30.0 /, ! Velocity lower limit for base maps & DUL / 200.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, & XCtt, & XCtimm 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, & cWild3DMHD*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' 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 in write3dinfo_td3d & XC3DT (nLng,nLat,nMap,4*nTmaxG,3), & XC3DTN (nLng,nLat,nMap,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 call ArrR4Constant(NLmaxV,1.0,OFIXV) ! Initialize OFIXV array call ArrR4Constant(NLmaxG,1.0,OFIXG) ! Initialize OFIXG 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 print*, ' ' call AskYN('Do you want to Iterate on a 3D-MHD file?',b3DMHD) if (b3DMHD) then ! A 3-D MHD external input is used at the end of the tomography program NETF = nTmaxG ! Total number of day times (to be set below if necessary) call AskI4('Minimum input files needed for an attempt?',ireadmhd) iYes = 2 ! iYes is 1 if a test write is made, 2 is a read. call AskI4('Read/Write test files: 0=No, 1=write, 2=read, 3=both?',iYes) Nit2 = 0 ! This says there has been no external input read yet. Nit2 = Nit2 + 1 with each external read. NiterT2 = 0 ! If zero, no external read input read yet. NiterT2 > 0 limits external reads to this number. call AskWhat('Use an external input from: ENLIL, MS-FLUKSS, H3D-MHD?$1$1',I) bENLIL = I .eq. 1 bMSFLUKSS = I .eq. 2 bH3DMHD = I .eq. 3 C C ENLIL is used as a tomography input C if(bENLIL) then C RAD1 = 0.1 + (1.0/320.0) ! ENLIL beginning height C RAD1 = 0.1 ! ENLIL beginning height - compatable with MS-FLUKSS print *, ' RAD1 - ENLIL lower Solar distance = ', RAD1 RADMS = 15.0*SUN__RAU ! This must be the same magnetic surface as the kinematic model MHDs = 1 ! ENLIL file read RRS = RAD1 ! Should be this print *, ' Start tomography at RSS = ', RRS C NETF = 47 ! Number of ENLIL time files available to read in call ForeignFile(-iVar-1,cVar,'ENL',cWild3DMHD) if (cWild3DMHD .eq. ' ') then cWild3DMHD = './mhd_files/ENLIL_XXXXXXXXXX_XXXXXXXX' end if end if C C End of ENLIL used as a tomography input C C******* C C MS-FLUKSS is used as a tomography input C if(bMSFLUKSS) then RAD1 = 0.1 ! MS-FLUKSS beginning height print *, ' RAD1 - MS-FLUKSS lower Solar distance = ', RAD1 RADMS = 15.0*SUN__RAU ! This must be the same magnetic surface as the kinematic model MHDs = 2 ! MS-FLUKSS file read RRS = RAD1 ! Should be this print *, ' Start tomography at RSS = ', RRS call ForeignFile(-iVar-1,cVar,'ENL',cWild3DMHD) if (cWild3DMHD .eq. ' ') then cWild3DMHD = './mhd_files/MS-FLUKSS_XXXXXXXXXX_XXXXXXXX' end if end if C C End of MS-FLUKSS used as a tomography input C bzeroth = .FALSE. call AskYN('Is this the zeroth Iteration for this MHD analysis?',bzeroth) C******* C New parameters used in all 3-D MHD inputs C dRRMS = RAD1 - RADMS RRMS = RADMS KPF = 1 ! The source surface should always begin at lowest MHD input level. FALLOFFT = 2.0 C & FALLOFFD /2.0/, C & FALLOFFV /0.0/, C & FALLOFFBR /2.0/, C & FALLOFFBT /1.0/, C & FALLOFFBN /1.34/, C FALLOFFN = 2.0 ! different name C FALLOFFBR = 2.0 ! same C FALLOFFBT = 1.0 ! same C FALLOFFBN = 1.0 ! different value VLLLL = 50.0 VLLUL = -50.0 end if print*, ' ' NiterT = 18 call AskI4('Number of kinematic model iterations?',NiterT) ! # iterations if (b3DMHD) call AskI4('Number of 3D-MHD model iterations?',NiterMHD) ! # 3D-MHD 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) XChi = XChi + 1.0 print *, 'XChi', XChi 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) 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(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 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 C if(bENLIL) XCbeg = 2056. if(bENLIL) XCbeg = 2114. if(bMSFLUKSS) XCbeg = 2114. if(bGCamb) XCbeg = 1884 call AskR4('Carrington rotation (0 = STOP)?'//cStr,XCbeg) if (XCbeg .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay') NintF = 1 nTF = nTG if(bForecast) then ! This never happens since here bforecast is not 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?$no',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.17 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?$no',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 BFACTORRC = 0.0 if(bWildMag(2).or.bWildMag(6).or.bWildMag(10).or.bWildMag(14)) BFACTORRC = 1.0 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 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.EQ.1.0.and.BFACTORRC.EQ.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?$yes',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======= do I=1,NLG+N_ING if(bGCamb) then XCG (NLG-I+1) = XCGG(I) YLG (NLG-I+1) = YLGG(I) DISTG(NLG-I+1) = DISTGG(I) GOBS (NLG-I+1) = GGOBS(I) DGOBS (NLG-I+1) = DGGOBS(I) GFREQ (NLG-I+1) = GFFREQ(I) GSIZE (NLG-I+1) = GSSIZE(I) XEG (NLG-I+1) = XEGG(I) XLSG (NLG-I+1) = XLSGG(I) XCEG (NLG-I+1) = XCEGG(I) XLLG (NLG-I+1) = XLLGG(I) XDLG (NLG-I+1) = XDLGG(I) DOYSG8(NLG-I+1) = DOYSGG8(I) IYRSG(NLG-I+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-I+1) = SRCGS end if if(bGNago.or.bGOoty.or.bGGen) then XCG(I) = XCGG(I) YLG(I) = YLGG(I) DISTG(I) = DISTGG(I) GOBS(I) = GGOBS(I) DGOBS(I) = DGGOBS(I) GFREQ(I) = GFFREQ(I) GSIZE(I) = GSSIZE(I) XEG(I) = XEGG(I) XLSG(I) = XLSGG(I) XCEG(I) = XCEGG(I) XLLG(I) = XLLGG(I) XDLG(I) = XDLGG(I) DOYSG8(I) = DOYSGG8(I) IYRSG(I) = IYRSGG(I) end if end do 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 if(.not.bForecast) then ! For non-forecast mode. Special for going 4 days into the future beyond the CR C NTV = NTV + 4 C NTG = NTG + 4 C ConsTMV = ConsTV*nTV C ConsTMD = ConsTG*nTG C NTV = NTV + 0 C NTG = NTG + 0 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 if(nTG.gt.NETF) then NTG = NETF NTV = NETF 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.or.WTSG(1,NLOSG/2).lt.0.00001) 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.or.WTSV(1,NLOSV/2).lt.0.00001) 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 if(.not.bzeroth) go to 9991 Nit = -1 NitV = -1 NitG = -1 9999 continue do while (Nit .lt. NiterT) Nit = Nit+1 print *, ' ' if(niterT2 .lt.1) then call Say('IPSD2020','I','Info','Kinematic Iteration '//cStr(:Int2Str(Nit,cStr))) else write(*,'(A)') '++++++++++++(This is an MHD iteration.)++++++++++++' call Say('IPSD2020','I','Info','MHD Internal Iteration '//cStr(:Int2Str(Nit,cStr))) call Say('IPSD2020','I','Info',' MHD New Iteration '//cStr(:Int2Str(Nit2,cStr))) end if C write(*,'(A,I3,4F7.3)') 'nit, CONRV, CONVT, CONSTV CONSTMV',nit, CONRV, CONVT, CONSTV, CONSTMV NVS = 0 if (Nit .eq. LimitoV .and. niterT2 .eq. 0) NVS = 1 if (Nit .eq. LimitoV .and. niterT2 .gt. 0) NVS = -1 NGS = 0 if (Nit .eq. LimitoG .and. niterT2 .eq. 0) NGS = 1 if (Nit .eq. LimitoG .and. niterT2 .gt. 0) 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) 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 MkVModeltd4r_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 C NBSV(I) = 0 ! Old way if(niterT2.eq.0) NBSV(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSV(I) = -2 ! MHD G setting VMOD(I) = Speed end if if(WTSV_in(I-NLV).eq.0.and.I.gt.NLV) then C NBSV(I) = 0 ! Old way if(niterT2.eq.0) NBSV(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSV(I) = -2 ! MHD G setting 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 ! Old way C if(niterT2.eq.0) NBSV(I) = -1 ! Kinematic G setting C if(niterT2.gt.0) NBSV(I) = -2 ! MHD G setting 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 ! Old way C if(niterT2.eq.0) NBSV(I) = -1 ! Kinematic G setting C if(niterT2.gt.0) NBSV(I) = -2 ! MHD G setting 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 ! Old way C if(niterT2.eq.0) NBSV(I) = -1 ! Kinematic G setting C if(niterT2.gt.0) NBSV(I) = -2 ! MHD G setting C end if C end do call FixModeltdn_mhdn(1,NLV+N_INV,VOBS,VMOD,PWV,NBSV,NBSVF,VSIG,FIX,OFIXV,VSSIG,VRATIO) C call FixModeltdn_mhd(1,NLV+N_INV,VOBS,VMOD,PWV,NBSV,NBSVF,VSIG,FIX,VSSIG,VRATIO) C 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), C & 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 if(NitV.gt.1) then 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 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(bCheck2.and..not.bVcon) then ! Do this as above for bCheck2 to run magnetic field analysis. write(*,'(2(A,F11.4))') 'Here 3 NEW DMAP(22,5,22) ',DMAP(22,5,22), ' NEW DMAPV(22,5,22)',DMAPV(22,5,22) do N=1,NTV 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 if(NiterT2.gt.0) then ! This says there is "valid" input volumetric data C iPrint = 1 ! if 1 prints out a lot iPrint = 0 iDo = 1 ! if iDo = 1 this shifts over over the xcshifts matrix used kinematic C iDo = 2 ! if iDo = 2 this shifts over the xcshifts and vratio and dratio values C iDo = 10 ! if 10, also smooths the input base value noi_pre = 0 ! Don't go through the shift value C noi_pre = 1 ! Go through the shift value if(noi_pre.eq.1) then call MkShiftdnma_pre(iPrint,iDo,bSmooth,XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG, & VMAPD,DMAP,VV3,BB3,DD1,TT1,XCshift,XCshift3,XCshiftM,DVfact,Vratio3,Bratio3,DDfact,Dratio1,Tratio1, & RRS,RRMS,dRR,dRRMS,KPF,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CONRV,CONRD,CONDT, & VLL,VUL,VLLLL,VLLUL,DLL/(RRG*RRG),DUL/(RRG*RRG),Vtmp,Dtmp,XLLTtmp) C end if end if if(NiterT2.eq.0) then ! This says to use the kinematic tomographic analysis 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 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), C & 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() C NBSG(I) = 0 ! Old setting type if(niterT2.eq.0) NBSG(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSG(I) = -2 ! MHD G setting 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).gt.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).gt.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).gt.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).gt.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() C NBSG(I) = 0 ! Old way if(niterT2.eq.0) NBSG(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSG(I) = -2 ! MHD G setting 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).gt.0) then GOBSS(I) = GOBSS(I)*DGOBSS end if if(GOBSS(I).le.GLIMP.and.GOBSS(I).gt.1.0.and.NBSG(I).gt.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).gt.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).gt.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 C NBSG(I) = 0 ! Old way if(niterT2.eq.0) NBSG(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSG(I) = -2 ! MHD G setting GMOD2(I) = 1.0 end if if(WTSD_in(I-NLG).eq.0.and.I.gt.NLG) then C NBSG(I) = 0 ! Old way if(niterT2.eq.0) NBSG(I) = -1 ! Kinematic G setting if(niterT2.gt.0) NBSG(I) = -2 ! MHD G setting GMOD2(I) = 1.0 end if end do call FixModeltdn_mhdn(2,NLG+N_ING,GOBSS,GMOD2,PWG,NBSG,NBSGF,GSIG,FIX,OFIXG,GSSIG,GRATIO) C call FixModeltdn_mhd(2,NLG+N_ING,GOBSS,GMOD2,PWG,NBSG,NBSGF,GSIG,FIX,GSSIG,GRATIO) C call FixModeltdn(2,NLG+N_ING,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), C & 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 if(nitG.gt.1) then 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 end if 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(NiterT2.gt.0) go to 9988 ! If > 0 this says there is "valid" input MHD volumetric data, so update even on last, C and this jumps around this caveat until the end of the iteration loop. if(Nit.ne.NiterT) then ! Update velocity shifts and vmap except on last 9988 continue C C If in-situ density is available define DEN1AU except on the last if bFACE is yes C if(bCheck2) go to 9997 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).gt.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).gt.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).gt.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 9997 continue 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) 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 if(NiterT2.gt.0) then ! If > 0 this says there is "valid" input MHD volumetric data C iPrint = 1 ! if 1 prints out a lot iPrint = 0 iDo = 1 ! if iDo = 1 this shifts over over the xcshifts matrix to be used as kinematic C iDo = 2 ! if iDo = 2 this shifts over the xcshifts and vratio and dratio values C iDo = 10 ! if 10, also smooths the input base value noi_pre = 0 ! Don't go through the shift value C noi_pre = 1 ! Go through the shift value if(noi_pre.eq.1) then call MkShiftdnma_pre(iPrint,iDo,bSmooth,XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG, & VMAP,DMAPV,VV3,BB3,DD1,TT1,XCshift,XCshift3,XCshiftM,DVfact,Vratio3,Bratio3,DDfact,Dratio1,Tratio1, & RRS,RRMS,dRR,dRRMS,KPF,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CONRV,CONRD,CONDT, & VLL,VUL,VLLLL,VLLUL,DLL/(RRG*RRG),DUL/(RRG*RRG),Vtmp,Dtmp,XLLTtmp) end if end if if(NiterT2.eq.0) then ! If = 0, this says to use the kinematic tomographic analysis call MkShiftdn0n(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV,Speed, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFFD,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) end if end if end if end if C print *, ' ' C print *, 'VMAP I,1-10,2 in iteration loop' C write(*,'(8F10.4)') (VMAP(I,1,2),I=1,15) C write(*,'(8F10.4)') (VMAP(I,6,2),I=1,15) C write(*,'(8F10.4)') (VMAP(I,10,2),I=1,15) C print *, 'VMAP I,1-10,30' C write(*,'(8F10.4)') (VMAP(I,1,30),I=1,15) C write(*,'(8F10.4)') (VMAP(I,6,30),I=1,15) C write(*,'(8F10.4)') (VMAP(I,10,30),I=1,15) C print *, 'VMAP I,1-10,40' C write(*,'(8F10.4)') (VMAP(I,1,40),I=1,15) C write(*,'(8F10.4)') (VMAP(I,6,40),I=1,15) C write(*,'(8F10.4)') (VMAP(I,10,40),I=1,15) C print *, ' ' C print *, 'DMAPV I,1-10,2' C write(*,'(8F10.4)') (DMAPV(I,1,2),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,6,2),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,10,2),I=1,15) C print *, 'DMAPV I,1-10,30' C write(*,'(8F10.4)') (DMAPV(I,1,30),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,6,30),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,10,30),I=1,15) C print *, 'DMAPV I,1-10,40' C write(*,'(8F10.4)') (DMAPV(I,1,40),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,6,40),I=1,15) C write(*,'(8F10.4)') (DMAPV(I,10,40),I=1,15) C print *, ' ' end do ! End of iteration loop C stop ! stop 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(bCheck2) go to 9998 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 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 print *, ' ' write (*,'(A,I4,A,I4,A)') 'MHD modeling removed',NBSVF,' velocity sources and',NBSGF,' g-level sources' print *, ' ' J = 0 ! List bad IPS sources NGOODSDVAL = 0 do I=1,NLV+N_INV if (NBSV(I) .le. 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).le.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 Now that the good and bad sources are written out, for the MHD analysis skip directly to the base writes with VMAP and DMAPV if(niterT2.gt.0) go to 9998 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 LOS 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 print *, ' ' write(*,'(A,13F10.4))') 'Here 12', (VMAP(I,1,2),I=1,13) print *, ' ' 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 write(*,'(A,3I10,2F15.8)') 'Before copy to VMAPD', nTG, nTV,nLngLatnTG, XCbegG(1,1),XCbeV(1,1) 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 following 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,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if C Must change back to regular value in current version prior 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***************************************************************************************************************************** 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 C do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP N, amin and amax = ', N, amin, amax C end do print *, 'Before forth filling of DMAP' print *, ' ' print *, 'VMAP I,1-10,2' write(*,'(8F10.4)') (VMAP(I,1,2),I=1,15) write(*,'(8F10.4)') (VMAP(I,6,2),I=1,15) write(*,'(8F10.4)') (VMAP(I,10,2),I=1,15) print *, 'VMAP I,1-10,30' write(*,'(8F10.4)') (VMAP(I,1,30),I=1,15) write(*,'(8F10.4)') (VMAP(I,6,30),I=1,15) write(*,'(8F10.4)') (VMAP(I,10,30),I=1,15) print *, 'VMAP I,1-10,40' write(*,'(8F10.4)') (VMAP(I,1,40),I=1,15) write(*,'(8F10.4)') (VMAP(I,6,40),I=1,15) write(*,'(8F10.4)') (VMAP(I,10,40),I=1,15) print *, ' ' 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 *, ' ' print *, 'After forth filling of DMAP' 9998 continue Print *, 'This is the value of nitert after 9998 continue', nitert if(niterT2.gt.0) then if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,DMAPV,DMAP) end if end if 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 print *, 'Copy array from VMAP to VMAPD through CopyVtoVDN' call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV) end if print *, ' ' print *, 'Before writes to filled nv3o or base volumetric data.' print *, ' ' print *, 'Check that base maps that are filled' print *, 'DMAP I,1-10,2' write(*,'(8F10.4)') (DMAP(I,1,2),I=1,55) write(*,'(8F10.4)') (DMAP(I,6,2),I=1,55) write(*,'(8F10.4)') (DMAP(I,10,2),I=1,55) print *, 'DMAP I,1-10,30' write(*,'(8F10.4)') (DMAP(I,1,30),I=1,55) write(*,'(8F10.4)') (DMAP(I,6,30),I=1,55) write(*,'(8F10.4)') (DMAP(I,10,30),I=1,55) print *, 'DMAP I,1-10,40' write(*,'(8F10.4)') (DMAP(I,1,40),I=1,55) write(*,'(8F10.4)') (DMAP(I,6,40),I=1,55) write(*,'(8F10.4)') (DMAP(I,10,40),I=1,55) print *, ' ' print *, 'VMAPD I,1-10,2' write(*,'(8F10.4)') (VMAPD(I,1,2),I=1,55) write(*,'(8F10.4)') (VMAPD(I,6,2),I=1,55) write(*,'(8F10.4)') (VMAPD(I,10,2),I=1,55) print *, 'VMAPD I,1-10,30' write(*,'(8F10.4)') (VMAPD(I,1,30),I=1,55) write(*,'(8F10.4)') (VMAPD(I,6,30),I=1,55) write(*,'(8F10.4)') (VMAPD(I,10,30),I=1,55) print *, 'VMAPD I,1-10,40' write(*,'(8F10.4)') (VMAPD(I,1,40),I=1,55) write(*,'(8F10.4)') (VMAPD(I,6,40),I=1,55) write(*,'(8F10.4)') (VMAPD(I,10,40),I=1,55) print *, ' ' print *, 'Before third write3D_infotd' Print *, 'This is the value of nitert after setting before bWrite3Do', nitert 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.or.bWRite3D_ext) then C Must change to low value in current version following 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 *, 'This is the value of nitert after setting before bWrite3D_HRb', nitert if(bWrite3D_HRb .and. niterT2 .eq. 0) then ! A kinematic D, V Base output 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 coord. (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. Mode 2 begins the 3D information one rotation before rotation center. if(.not.bForecast) then nar = 2 bFcst = .TRUE. else nar = 0 bFcst = .TRUE. end if call write3D_infotd3DM_HR_3(nar,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, & iYrBG,RRS,dRR,nLng,nLat,nMapb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,bFcst,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if if(bWrite3D_ext .and. niterT2 .gt. 0) then ! An MHD D, V Base output print *, ' ' print *, 'Before write3D_infotd of High Resolution MHD output' cPre = 'nv3m' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = 3 C nMapb = nMapHRb ! In original coordinates (so nMapHRMb = 0.25 AU). nMapm = nMapHRm if(nMapHRMm*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMm = nMapm ! Make the setting of nMapHR idiot proof if(nMapm.gt.nMapHRMm) nMapm = nMapHRMm C if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb ! Make the setting of nMapHR idiot proof C if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb write(*,'(A,2F10.7,2I6)') 'RRS, dRR, nMapm, nMapHRMm',RRS, dRR, nMapm, nMapHRMm C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. Mode 2 begins the 3D one rotation before the rotation center. if(.not.bForecast) then nar = 2 bFcst = .TRUE. else nar = 0 bFcst = .TRUE. end if C niterT = 17 Print *, 'This is the value of nitert before bWrite3D_infotd3DM_HR_3', nitert call write3D_infotd3DM_HR_3(nar,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, & iYrBG,RRS,dRR,nLng,nLat,nMapm,nMapHRMm,nMap,nTF,nTmaxG,nCar,JDCar,bFCst,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) stop C call write3D_infotd3DM_HR_3(nar,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, C & iYrBG,RRS,dRR,nLng,nLat,nMapb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,bFcst,XCbeg,xInc,NCC, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, C & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) C call write3D_infotd3DM_HR_3DMHD(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegg,XCtbegG,XCtendG, C & iYrBG,RRS,dRR,nLng,nLat,nMapm,nMapHRMm,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xIncd,NCC, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, C & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift3,DVfact,DDfact,V3DHR,D3DHR, C & DDD1,VVV3) ! BVJ uses nMapm,nMapHRMm,XCshift3, Writes DDD1, VVV3(i,j,k,n,1). VVV3(i,j,k,n,2) C & VVV3(i,j,k,n,3) C call write3D_infotd3DM_HR_3DMHD(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG, C & iYrBG,RRS,dRR,nLng,nLat,nMapm,nMapHRMm,nMap,nTF,nTmaxG,nCar,JDCar,bForecast,XCbeg,xInc,NCC, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, C & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift3,DVfact,DDfact,V3DHR,D3DHR) ! BVJ uses nMapm,nMapHRMm,XCshift3 end if 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 .or. bMagnetic_ext) then 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_set (T3D__DRR , 0, dRR ) 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, RRMS ) 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 ) call get_bbtm_3(XCbeGG,nLng,nLat,nTG,RRMS,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN, & XCINTG,BRR2DT,BRC2DT,BTC2DT,BNC2DT,iBflag) if( bmagexx .or. bmagexx_ext) then ! Magnetic field files output either from internal or external files if( niterT2 .eq. 0) then ! A magnetic field Base output write of kinematic files print *, ' ' print *, 'Before write3D_bbtm_HR_3 after HR base; niterT2 equals 0' NintLng = NintLn NintLat = NintLa NintHt = NintH Ninter = 3 C NintLng = 0 C NintLat = 0 C NintHt = 0 C Ninter = 0 ClipLng = 90.0 nMapbb = nMapHRb ! In original map coordinates (so nMapHRMb = 0.25 AU). if(nMapHRMb*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMb = nMapb ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHRMb) nMapb = nMapHRMb IradialBRR = 0 IradialBRC = 0 IradialBTC = 0 IradialBNC = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 smootsh polar holes IradialBRR = IradialBRR + 1 ! Case 1 end if call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 smooths polar holes IradialBRC = IradialBRC + 1 ! Case 2 end if call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 smooths polar holes IradialBTC = IradialBTC + 1 ! Case 3 end if call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD/2.0,2,0.0,0.0) ! Change 6/27/16 smooths polar holes IradialBNC = IradialBNC + 1 ! Case 4 end if end do if(bforecast) then if(IradialBRR.ne.0.and.IradialBRR.ne.(nTF+1)) then ! Case 1 write(*,'(A,I4,A)') ' There are no radial magnetic fields for',(nTF+1)-IradialBRR,' BRR source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRR2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRR2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRR2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRR2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BRR source surface',N,' was completely filled by gridsphere' else if(amin.eq.BadR4().and.amax.eq.BadR4()) then write(*,'(A,I4,A)') 'BRR source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRR source surfaces could not be filled with gridsphere' end if if(IradialBRC.ne.0.and.IradialBRC.ne.(nTF+1)) then ! Case 2 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBRC,' BRC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BRC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BRC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BRC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BRC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BRC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BRC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BRC source surfaces were not filled with gridsphere' end if if(IradialBTC.ne.0.and.IradialBTC.ne.(nTF+1)) then ! Case 3 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBTC,' BTC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BTC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BTC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BTC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BTC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BTC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BTC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BTC source surfaces were not filled with gridsphere' end if if(IradialBNC.ne.0.and.IradialBNC.ne.(nTF+1)) then ! Case 4 write(*,'(A,I4,A)') '.There are no radial magnetic fields for',(nTF+1)-IradialBNC,' BNC source surfaces.' print *, 'These are filled in the following code.' print *, ' ' call FillWholeT(2,Nit,NiterT,XCbeGG,nLng,nLat,nTF+1,ConsTMD,BNC2DT,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTF+1,ConsTG,BNC2DT,DDD) Iradial = 0 do N=1,nTF+1 call arrR4getminmax(nLngLat,BNC2DT(1,1,N),amin,amax) if(amax.ne.BadR4().and.amin.eq.BADR4()) then call GridSphere2D(ALng,nLng,nLat,1,BNC2DT(1,1,N),CONRD,4,0.0,0.0) write(*,'(A,I4,A)') 'BNC source surface',N,' was completely filled by gridsphere' else if(amax.ne.BadR4()) then write(*,'(A,I4,A)') 'BNC source surface',N,' was completely empty' Iradial = Iradial + 1 end if end if end do if(Iradial.ne.0) write(*,'(A,I4,A)') 'Error!', Iradial,' BNC source surfaces were not filled with gridsphere' end if end if call xc3dtshift_rrms(0,niterT2,nLng,nLat,nMap,nTF,RRS,dRR,RADMS,XCshift,XCshift3,XCshiftM) C call write3D_bbtm_HR_3(bForecast,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG, C & iYrBG,RRS,dRR,nLng,nLat,nMapbb,nMapHRMb,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC, C & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, C & ClipLng,VMAP,XCshift,DVfact, C & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR) if(.not.bForecast) then nar = 2 bFcst = .TRUE. else nar = 0 bFcst = .TRUE. end if call write3D_bbtm_HR_3DMHD(nar,bFcst,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG, & iYrBG,RRS,dRR,RRMS,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,XCshift3,XCshiftM,DVfact, ! Uses RRMS nMapHRm, nMapbm, XCshift3, XCshiftM & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR, & BBB3) ! Uses RRMS nMapHRm, nMapbm, XCshift3, XCshiftM, Writes BBB3 end if if( niterT2 .ne. 0) then ! A magnetic field Base output write of MHD files print *, ' ' print *, 'Before write3D_bbtm_HR_3DMHD after HR base; niterT2 not equal to 0' NintLng = NintLn NintLat = NintLa NintHt = NintH Ninter = 3 ClipLng = 90.0 nMapbm = nMapHRm ! In original map coordinates (so nMapHRMb = 0.25 AU). if(nMapHRMm*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHRMm = nMapbm ! Make the setting of nMapHR idiot proof if(nMapbm.gt.nMapHRMm) nMapbm = nMapHRMm C C The subroutine below (uncommented) provides an approximate value for XCshiftM if needed. C call xc3dtshift_rrms(0,niterT2,nLng,nLat,nMap,nTF,RRS,dRR,RADMS,XCshift,XCshift3,XCshiftM) if(.not.bForecast) then nar = 2 bFcst = .TRUE. else nar = 0 bFcst = .TRUE. end if write(*,'(A,3F10.7,2I6)') 'RRS, dRR, RRMS, nMapbm, nMapHRMm',RRS, dRR, RRMS, nMapbm, nMapHRMm call write3D_bbtm_HR_3DMHD(nar,bFcst,Nit,NiterT,NintLng,NintLat,NintHt,Ninter,NCoff,XCintDG,XCbeGG,XCtbegG,XCtendG, & iYrBG,RRS,dRR,RRMS,nLng,nLat,nMapbm,nMapHRMm,nMap,nTF,nTmaxG,nCar,JDCar,XCbeg,xInc,NCC, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & ClipLng,VMAP,XCshift3,XCshiftM,DVfact, ! Uses RRMS nMapHRm, nMapbm, XCshift3, XCshiftM & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,FALLOFFBN,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR, & BBB3) ! Uses RRMS nMapHRm, nMapbm, XCshift3, XCshiftM, Writes BBB3 end if end if if( bmagex .or. bmagex_ext) then ! A magnetic field extraction of in-situ data if( niterT2 .eq. 0) then ! A kinematic planet extraction of in-situ data write(*,'(A,F10.6,F12.7,F10.6)') ' Before Final Extractdvdm_3DMHD with fields',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 xc3dtshift_rrms(0,niterT2,nLng,nLat,nMap,nTG,RRS,dRR,RADMS,XCshift,XCshift3,XCshiftM) call Extractdvdm_3DMHD(bForecast,bExtractf,cPrefixx,NTP,RRS,dRR,RRMS,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift3,XCshiftM, & DVfact,DDfact,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar) ! Uses cPrefixm, RRMS, XCshift3, XCshiftM C call Extractdvdm_3(bForecast,bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, C & nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift, C & DVfact,DDfact,TT3D,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar) end if if( niterT2 .gt. 0) then ! An MHD planet extraction of in-situ data. NinterD = 3 C C The subroutine below (uncommented) provides an approximate value for XCshiftM if needed C call xc3dtshift_rrms(0,niterT2,nLng,nLat,nMap,nTG,RRS,dRR,RADMS,XCshift,XCshift3,XCshiftM) call Extractdvdm_3DMHD(bForecast,bExtractf,cPrefixm,NTP,RRS,dRR,RRMS,XCbeGG,XCtbegG,XCtendG,IYRBG,nLng,nLat, & nMap,nTG,NinterD,FALLOFFD,FALLOFFV,FALLOFFBR,FALLOFFBT,FALLOFFBN,VMAPD,DMAP,XCshift3,XCshiftM, & DVfact,DDfact,BRR2DT,BRC2DT,BTC2DT,BNC2DT,NCoff,XCstrt,nCar,JDCar) ! Uses cPrefixm, RRMS, XCshift3, XCshiftM end if end if end if if(.not. bMagnetic .and. niterT2 .eq. 0) then ! A kinematic, non-magnetic in-situ extraction output print *, ' Before Final Extractdvd, no fields',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 if(.not. bMagnetic .and. niterT2 .gt. 0) then ! An MHD, non-magnetic in-situ extraction output print *, ' Before Final Extractdvd, no fields',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdvd(bForecast,bExtractf,cPrefixfnic,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 C****************************************************************** C C This is the normal end of the v18 tomography program C C****************************************************************** C if( b3DMHD ) then 9991 continue ! The input read didn't work. After a pause try to input another input file externally C C Set up the tomography program to read and provide data from an external source (at the spatial, temporal cadence, C and height resolution of the original tomography program. C write(*,'(A,I4)') 'Set up an external set of files to read, NiterT2 =',NiterT2 9986 continue C call AskYN('First time: Do you want to continue?$no',bCont) if(.not.bCont) go to 9986 call AskYN('Final time: Do you really want to continue?$yes',bCont) if(.not.bCont) go to 9985 nTF = nTG C C Write out a test input file to be read externally or read in an external test file. C print *, 'Before Test External Write, iYes =', iYes C C Originally LLBegg came from the write of the 4d files, now this is not defined. Thus, put in some number. C C If niterT2 = 0, the kinematic iteration, then XCshift3 is filled with XCshift values and XCshiftM is filled C with an approximation by xc3dtshift_rrms.f. If niterT2 > 0 then XCshift3 and XCshiftM are filled by C MkShiftdnma_pre.f from the input MHD files. C LLfst = 1 LLEndd = nTG print *, NTG,NTV,NTP,LLfst,LLEndd if( iYes .eq. 0 ) print *, 'Dont write out any files from the modeling.' if(niterT2 .eq. 0 .and. iYes .eq. 1) print *, 'Write out test external files from the kinematic model.' if(niterT2 .gt. 0 .and. iYes .eq. 1) print *, 'Write out external files from the externally input model.' if( iYes .eq. 2 ) print *, 'Read in files from an external model.' write(*,'(3(A,F11.4))') & 'Before Ext VMAPD(22,5,22)',VMAPD(22,5,22),' VMAP(22,5,22)',VMAP(22,5,22),' DMAP(22,5,22)',DMAP(22,5,22) C call ExternalRWMHD(MHDs,bForecast,iYes,cWild3DMHD,nLng,nLat,nMap,nTF,nTmaxG,XCbegg,XCbeg,xInc,NCC, & LLfst,LLEndd,JDCar,nCar,RRS,dRR,RRMS,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,R1AU, & BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,iBFlag,BRR2DT,BRC2DT,BTC2DT,BNC2DT,BR3DHR,BT3DHR,BN3DHR, & 'null',Nit,NiterT,NCoff,XCintDG,XCtbegG,XCtendG,iYrBG, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & bDdenerb,bVdenerb,bVvererb,bDvererb,ERDLOSB,ERVLOSB,ANMAPD,ANMAPV, & Scale,VMAPD,VMAP,DMAP,VMHR,DMHR,XCshift3,XCshiftM,DVfact,DDfact,V3DHR,D3DHR, & ConsTMD,ConsTMV,ConstL,aNdayG,aNdayV, & DDD1,TTT1,VVV3,BBB3,DD1,TT1,VV3,BB3,iNum) do L=1,0 C do L=1,nTF do I=1,nLng do J=1,nLat do N=1,nMap if(N.eq.10.and.J.eq.5.and.I.eq.22) write(*,'(A,4I3,F6.2,F7.1,5F7.3,F7.2)') & 'After ExRead',I,J,N,L, DDD1(I,J,N,L),VVV3(I,J,N,L,1),VVV3(I,J,N,L,2),VVV3(I,J,N,L,3), & BBB3(I,J,N,L,1),BBB3(I,J,N,L,2),BBB3(I,J,N,L,3),TTT1(I,J,N,L) end do end do end do end do C stop write(*,'(3(A,F11.4))') 'After Ext VMAPD(22,5,22)',VMAPD(22,5,22),' VMAP(22,5,22)',VMAP(22,5,22),' DMAP(22,5,22)',DMAP(22,5,22) C C End of the write out of an input test file. C C End of read of an external input file C C If the external read was successful and the new 3D input has been set up, begin a new external read iteration. C iPrint = 1 do I=1, nLng do J=1, nLat do N=1, nMap do L=1, nTG if(iPrint.eq.1.and.N.eq.1.and.I.eq.22.and.J.eq.5.and.bzeroth) then write(*,'(A,I3,3(A,F10.4))') & 'L ',L,' Old VMAPD(22,5,L)',VMAPD(22,5,L),' VMAP(22,5,L)',VMAP(22,5,L), ' DMAP(22,5,L)',DMAP(22,5,L) end if do NX=1,3 if(iPrint.eq.1.and.bzeroth) then if(NX.eq.1.and.J.eq.5.and.I.eq.30.and.N.eq.2) & write(*,'(A,4I4,3F12.8)'), 'Here 1',I,J,L,NX,XCshift(I,J,2,L,NX),XCshift3(I,J,2,L,NX),XCshiftM(I,J,L,NX) if(NX.eq.3.and.J.eq.5.and.I.eq.30.and.N.eq.2) & write(*,'(A,4I4,3F12.8)'), 'Here 1',I,J,L,NX,XCshift(I,J,2,L,NX),XCshift3(I,J,2,L,NX),XCshiftM(I,J,L,NX) end if end do end do end do end do end do if(iNum.gt.ireadmhd) then write(*,'(A,I3)') ' There has been a successful 3-D MHD external read.' else write(*,'(A,I3)') ' The 3-D MHD external read was unsuccessful.' go to 9991 end if C C ************************************************************************************************* C Continue MHD iterations with a read of the input files C print *, ' ' call AskYN('Do you want to write out the MHD analysis input?$no',bcheck3) print *, ' ' if(bcheck3) then call write3D_input_3D_MHD_DV(cPre,Nit,NiterT,NinterC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR,R1AU, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,FALLOFFV,FALLOFFD,VVV3,DDD1,V3D,D3D,VM,DM) end if call AskYN('Do you want to write out the traced MHD input and stop?$no',bcheck2) if(bCheck2) bCont = .TRUE. print *, ' ' if(.not.bcheck2) call AskYN('Do you want to continue the MHD tomography iterations?$no',bCont) if(.not.bCont) go to 9987 C Nit = -1 ! Reset the internat MHD iterations NitV = -15 NitG = -1 Nit2 = Nit2 + 1 ! Begin a new MHD at iteration Nit2 = Nit2 + 1 NiterT = 18 bzeroth = .FALSE. NBSVF = 0 NBSGF = 0 if(.not.bcheck2) call AskI4('Number of MHD iterations per this input?',NiterT) ! # iterations if(bcheck2) NiterT = 1 if(NiterT.gt.0) then NiterT2 = 3 else NiterT2 = 0 end if PWRGF = PWRG PWRVF = PWRV if(.not.bcheck2) call AskR4('Radial G-level falloff to fit G density?',PWRG) PWRV = PWRG if(.not.bcheck2) call AskR4('Radial G-level falloff to fit V density?',PWRV) NNBSG = 0 do I=1,NLMAXG ! reset the g-level source values changed bad to good if(NBSG(I).eq.-1) then NBSG(I) = 1 NNBSG = NNBSG + 1 end if end do if(.not.bcheck2) write(*,'(A,I4,A)') 'You just reset',NNBSG,' bad density kinematic source values to good' if(PWRG.ne.PWRGF) then print *,' ' write (*,'(A,A,I10,F8.1,3F8.4,3F8.3,F8.1)') '1029, NBSG, GFREQ, GSIZE, XEG, DISTG WTSG(1), WTSG(10), WTSG(10), PWRG', & ' ',NBSG(1029),GFREQ(1029),GSIZE(1029),XEG(1029),DISTG(1029),WTSG(1,1029),WTSG(10,1029),WTSG(20,1029),PWRG do I=1,NLG C call MkLOSWeightsx(NLOSWG,dLOSG,NLOSG,NTGD,GFREQ(NTGB),GSIZE(NTGB),XEG(NTGB),DISTG(NTGB),WTSG(1,NTGB),PWRG) if(NBSG(I).ne.0) call MkLOSWeightsx(NLOSWG,dLOSG,NLOSG,1,GFREQ(I),GSIZE(I),XEG(I),DISTG(I),WTSG(1,I),PWRG) end do write (*,'(A,A,I10,F8.1,3F8.4,3F8.3,F8.1)') '1029, NBSG, GFREQ, GSIZE, XEG, DISTG WTSG(1), WTSG(10), WTSG(10), PWRG', & ' ',NBSG(1029),GFREQ(1029),GSIZE(1029),XEG(1029),DISTG(1029),WTSG(1,1029),WTSG(10,1029),WTSG(20,1029),PWRG print *,' ' end if NNBSV = 0 do I=1,NLMAXV ! reset the velocity source values changed bad to good if(NBSV(I).eq.-1) then NBSV(I) = 1 NNBSV = NNBSV + 1 end if end do if(.not.bcheck2) write(*,'(A,I4,A)') 'You just reset',NNBSV,' bad velocity kinematic source values to good' if(PWRV.ne.PWRVF) then print *,' ' write (*,'(A,A,I10,F8.1,3F8.4,3F8.2,F8.1)') '699, NBSV, VFREQ, VSIZE, XEV, DISTV WTSV(1), WTSV(10), WTSV(10), PWRV', & ' ',NBSV(699),VFREQ(699),VSIZE(699),XEV(699),DISTV(699),WTSV(1,699),WTSV(10,699),WTSV(20,699),PWRV do I=1,NLV C call MkLOSWeightsx(NLOSWV,dLOSV,NLOSV,NTVD,VFREQ(NTVB),VSIZE(NTVB),XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) if(NBSV(I).ne.0) call MkLOSWeightsx(NLOSWV,dLOSV,NLOSV,1,VFREQ(I),VSIZE(I),XEV(I),DISTV(I),WTSV(1,I),PWRV) end do write (*,'(A,A,I10,F8.1,3F8.4,3F8.2,F8.1)') '699, NBSV, VFREQ, VSIZE, XEV, DISTV WTSV(1), WTSV(10), WTSV(10), PWRV', & ' ',NBSV(699),VFREQ(699),VSIZE(699),XEV(699),DISTV(699),WTSV(1,699),WTSV(10,699),WTSV(20,699),PWRV print *,' ' end if if(NiterT.gt.4) LimitoV = NiterT/2 if(.not.bcheck2) call AskI4('Iteration at which to limit V-source deviations, -1=none?',LimitoV) if(NiterT.gt.4) LimitoG = NiterT/2 if(.not.bcheck2) call AskI4('Iteration at which to limit G-source deviations, -1=none?',LimitoG) C C Now ask if you want to write out nv3m*, magnetic, and e3 files at the end of the MHD run C if(bcheck2) bwrite3D_ext = .True. if(.not.bcheck2) call AskYN('Output 3D high resolution MHD base files?$yes',bWrite3D_ext)! Filled 3D MHD "m" file write if(bWrite3D_ext) then C RRRRS = RRS if(.not.bcheck2) call AskR4('Minimum high-resolution base height in AU?',RRRRS) if(RRRRS.le.RRS) then nMapHRm = 1 else nMapHRm = nint((RRRRS-RRS)/dRR) + 1 end if if(nMapHRm.gt.nMap) nMapHRm=nMap ! Limit maximum height to maximum available C RRRRSME = 0.17 if(bcheck2) RRRRSME = 1.5 call AskR4('Maximum high-resolution base height in AU?',RRRRSME) nMapHRMm = nint((RRRRSME-RRS)/dRR) + 1 if(nMapHRMm.gt.nMap) nMapHRMm = nMap ! Limit maximum height to maximum available if(nMapHRMm.lt.nMapHRm) nMapHRMm = nMapHRm ! Don't be stupid! bDdenerb = .FALSE. ! If overwritten this will place the base maps the same as done for kinematic bVvererb = .FALSE. ! If overwritten this will place the base maps the same as done for kinematic C end if C if (bMagnetic) call AskYN('Produce magnetic 3DMHD field files?$yes',bmagexx_ext) if (bMagnetic) call AskYN('Extract magnetic 3DMHD fields?$yes',bmagex_ext) C C Condition the inputs for the tomography program C C Provide values for DMAPV from the input MHD data file for density. C Assume the density and velocity files have the same resolution C print *, ' ' write(*,'(2(A,F11.4))') 'OLD VMAPD(22,5,22)',VMAPD(22,5,22),' OLD VMAP(22,5,22) ',VMAP(22,5,22) write(*,'(2(A,F11.4))') 'OLD DMAP(22,5,22) ',DMAP(22,5,22), ' OLD DMAPV(22,5,22)',DMAPV(22,5,22) write(*,'(A)') 'Someday other than radial velocities may need to be reconstructed.' do I=1,nLng do J=1,nLat if (bVcon) then do L=1,nTV VMAPD(I,J,L) = VV3(I,J,1,L,1) VMAP(I,J,L) = VV3(I,J,1,L,1) end do end if if (bGcon) then do L=1,nTV DMAP(I,J,L) = DD1(I,J,1,L) DMAPV(I,J,L) = DD1(I,J,1,L) if(iPrint.eq.1.and.N.eq.1.and.I.eq.22.and.J.eq.5) then write(*,'(A,I3,3(A,F10.4))') & 'L ',L,' Ext VMAPD(22,5,L)',VMAPD(22,5,L),' VMAP(22,5,L)',VMAP(22,5,L), ' DMAP(22,5,L)',DMAP(22,5,L) end if end do end if end do end do write(*,'(2(A,F11.4))') 'NEW VMAPD(22,5,22)',VMAPD(22,5,22), ' NEW VMAP (22,5,22)',VMAP(22,5,22) write(*,'(2(A,F11.4))') 'NEW DMAP (55,5,22)',DMAP(55,5,22), ' NEW DMAPV(22,5,22)',DMAPV(22,5,22) write(*,'(2(A,F11.4))') ' DD1(55,5,1,22) ',DD1(55,5,1,22), ' VV3(22,5,1,22,1) ',VV3(22,5,1,22,1) write(*,'(2(A,F11.4))') ' VV3(22,5,10,22,2)',VV3(22,5,10,22,2),' VV3(22,5,10,22,3) ',VV3(22,5,10,22,3) write(*,'(2(A,F11.4))') 'BAD DMAP (55,5,11)',DMAP(55,5,11), ' BAD DMAPV(55,5,22)',DMAPV(55,5,22) write(*,'(A,3F11.7)') 'XCshift (22,5,10,22,1),2),3)',XCshift(22,5,10,22,1),XCshift(22,5,10,22,2),XCshift(22,5,10,22,3) write(*,'(A,3F11.7)') 'XCshift3(22,5,10,22,1),2),3)',XCshift3(22,5,10,22,1),XCshift3(22,5,10,22,2),XCshift3(22,5,10,22,3) write(*,'(A,3F11.7)') 'XCshiftM(22,5,10,22,1),2),3)',XCshiftM(22,5,22,1),XCshiftM(22,5,22,2),XCshiftM(22,5,22,3) C C Check below (bcheck) to see that the input base files are input correctly and give good volumetric data with the C kinematic model upward traced XCshift values Note, this only works if the kinematic maodel is run and iterated first. C if(bWrite3D_ext.and..not.bcheck2) call AskYN('Produce files to check input densities and velocities?$no',bcheck) if(bWrite3D_ext.and.bcheck) then ! A kinematic D, V Base output C C Must change to low value in current version C RRSCON = (RRS/R1AU)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do print *, ' ' print *, 'Before test write3D_infotd of High Resolution base check' cPre = 'nv3c' 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. if(.not.bForecast) then nar = 2 bFcst = .TRUE. else nar = 0 bFcst = .TRUE. end if 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,ANMAPV, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) C print *, ' ' write(*,'(2(A,F11.4))') 'BAD DMAP(55,5,11) ',DMAP(55,5,11), ' NEW DMAPV(55,5,22)',DMAPV(55,5,22) C C Change the density back since the falloff at distance is produced elsewhere C RRSCON = (R1AU/RRS)**FALLOFFD do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do write(*,'(2(A,F11.4))') 'Here 0 NEW DMAP(22,5,22) ',DMAP(22,5,22), ' NEW DMAPV(22,5,22)',DMAPV(22,5,22) NinterD = 3 call Extractdvd(bForecast,bExtractf,cPrefixfnim,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 C End check (b check) to see that the input 3-D MHD files are translated correctly with kinematic model values of XCshift C C C If noi_pre = 1 below provide the traceback from the input 3-D MHD model C C print *, ' ' print *, ' Before MkShiftdnma_pre subroutine' print *, ' ' C C iPrin = 1 ! if 1 prints out a lot iPrin = 0 C iDo = 1 ! if iDo = 1 this shifts over over the xcshifts matrix to be used as kinematic iDo = 2 ! if iDo = 2 this shifts over the xcshifts and vratio and dratio values C iDo = 10 ! if 10, also smooths the input base value C noi_pre = 0 ! Don't go through the shift value noi_pre = 1 ! Go through the shift value if(noi_pre.eq.1) then bSMOOTH=.FALSE. if(bMSFLUKSS) bSMOOTH=.TRUE. call MkShiftdnma_pre(iPrin,iDo,bSmooth,XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG,VMAP,DMAPV, & VV3,BB3,DD1,TT1,XCshift,XCshift3,XCshiftM,DVfact,Vratio3,Bratio3,DDfact,Dratio1,Tratio1,RRS,RRMS, & dRR,dRRMS,KPF,FALLOFFD,FALLOFFT,FALLOFFBR,FALLOFFBT,FALLOFFBN,CONRV,CONRD,CONDT,VLL,VUL,VLLLL,VLLUL, & DLL/(RRG*RRG),DUL/(RRG*RRG),Vtmp,Dtmp,XLLTtmp) bSMOOTH=.FALSE. end if print *, ' ' print *, 'After MkShiftdnma_pre subroutine' if(ido.eq.0) print *, 'iDo = 0 , XCSHIFT, etc, not reset from original kinematic value.' if(ido.eq.1) print *, 'iDo = 1 , XCSHIFT, etc, reset to mhd value., Input base is VMAP, DMAP.' if(ido.eq.10) print *,'iDo = 10, XCSHIFT, etc, reset to mhd value., Input base is smoothed.' C write(*,'(A,3F11.7)')'After pre Vratio3(22,5,10,22,1),2),3)',Vratio3(22,5,10,22,1),Vratio3(22,5,10,22,2), C & Vratio3(22,5,10,22,3) C write(*,'(A,3F11.7)')'After pre Bratio3(22,5,10,22,1),2),3)',Bratio3(22,5,10,22,1),Bratio3(22,5,10,22,2), C & Bratio3(22,5,10,22,3) write(*,'(2(A,F11.4))') 'NEW VMAPD(22,5,22)',VMAPD(22,5,22),' NEW VMAP(22,5,22) ',VMAP(22,5,22) write(*,'(2(A,F11.4))') 'NEW DMAP(22,5,22) ',DMAP(22,5,22), ' NEW DMAPV(22,5,22)',DMAPV(22,5,22) write(*,'(A,3F11.7)') 'XCshift (22,5,10,22,1),2),3)',XCshift(22,5,10,22,1),XCshift(22,5,10,22,2),XCshift(22,5,10,22,3) write(*,'(A,3F11.7)') 'XCshift3(22,5,10,22,1),2),3)',XCshift3(22,5,10,22,1),XCshift3(22,5,10,22,2),XCshift3(22,5,10,22,3) write(*,'(A,3F11.7)') 'XCshiftM(22,5,22,1),2),3) ',XCshiftM(22,5,22,1), XCshiftM(22,5,22,2), XCshiftM(22,5,22,3) print *, ' ' write(*,'(A,F11.7)') 'DVfact (22,5,10,22) ', DVfact(22,5,10,22) write(*,'(A,3F11.7)') 'VRatio3(22,5,10,22,1),2),3) ', Vratio3(22,5,10,22,1),Vratio3(22,5,10,22,2),Vratio3(22,5,10,22,3) write(*,'(A,F11.7)') 'DDfact (22,5,10,22) ', DDfact(22,5,10,22) write(*,'(A,F11.7)') 'Dratio1(22,5,10,22) ', Dratio1(22,5,10,22) write(*,'(A,3F11.7)') 'Bratio3(22,5,10,22,1),2),3) ', Bratio3(22,5,10,22,1),Bratio3(22,5,10,22,2),Bratio3(22,5,10,22,3) write(*,'(A,F11.7)') 'Tratio1(22,5,10,22) ', Tratio1(22,5,10,22) C C End of external input precondition C print *, ' ' print *, 'DMAPV sample I,1-10,2' write(*,'(6F11.4)') (DMAPV(I,1,2),I=1,55) write(*,'(6F11.4)') (DMAPV(I,6,2),I=1,55) write(*,'(6F11.4)') (DMAPV(I,10,2),I=1,55) print *, 'DMAPV sample I,1-10,30' write(*,'(6F11.4)') (DMAPV(I,1,30),I=1,55) write(*,'(6F11.4)') (DMAPV(I,6,30),I=1,55) write(*,'(6F11.4)') (DMAPV(I,10,30),I=1,55) print *, 'DMAPV sample I,1-10,40' write(*,'(6F11.4)') (DMAPV(I,1,40),I=1,55) write(*,'(6F11.4)') (DMAPV(I,6,40),I=1,55) write(*,'(6F11.4)') (DMAPV(I,10,40),I=1,55) print *, ' ' print *, 'VMAP sample I,1-10,2' write(*,'(6F11.4)') (VMAP(I,1,2),I=1,55) write(*,'(6F11.4)') (VMAP(I,6,2),I=1,55) write(*,'(6F11.4)') (VMAP(I,10,2),I=1,55) print *, 'VMAP sample I,1-10,30' write(*,'(6F11.4)') (VMAP(I,1,30),I=1,55) write(*,'(6F11.4)') (VMAP(I,6,30),I=1,55) write(*,'(6F11.4)') (VMAP(I,10,30),I=1,55) print *, 'VMAP sample I,1-10,40' write(*,'(6F11.4)') (VMAP(I,1,40),I=1,55) write(*,'(6F11.4)') (VMAP(I,6,40),I=1,55) write(*,'(6F11.4)') (VMAP(I,10,40),I=1,55) print *, ' ' print *, 'XCshift longitude shift (in CR) sample at 1 AU I,1-10,10,2,1' write(*,'(6F11.4)') (XCshift(I,1,10,2,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,2,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,2,1),I=1,55) print *, 'XCshift longitude shift (in CR) sample at 1 AU I,1-10,10,30,1' write(*,'(6F11.4)') (XCshift(I,1,10,30,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,30,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,30,1),I=1,55) print *, 'XCshift longitude shift (in CR) sample at 1 AU I,1-10,10,40,1' write(*,'(6F11.4)') (XCshift(I,1,10,40,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,40,1),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,40,1),I=1,55) print *, ' ' print *, 'XCshift latitude shift (in deg) sample at 1 AU I,1-10,10,2,2' write(*,'(6F11.4)') (XCshift(I,1,10,2,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,2,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,2,2),I=1,55) print *, 'XCshift latitude shift (in deg) sample at 1 AU I,1-10,10,30,2' write(*,'(6F11.4)') (XCshift(I,1,10,30,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,30,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,30,2),I=1,55) print *, 'XCshift latitude shift (in deg) sample at 1 AU I,1-10,10,40,2' write(*,'(6F11.4)') (XCshift(I,1,10,40,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,40,2),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,40,2),I=1,55) print *, ' ' print *, 'XCshift time shift (in days) sample at 1 AU I,1-10,10,2,3' write(*,'(6F11.4)') (XCshift(I,1,10,2,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,2,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,2,3),I=1,55) print *, 'XCshift time shift (in days) sample at 1 AU I,1-10,10,30,3' write(*,'(6F11.4)') (XCshift(I,1,10,30,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,30,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,30,3),I=1,55) print *, 'XCshift time shift (in days) sample at 1 AU I,1-10,10,40,3' write(*,'(6F11.4)') (XCshift(I,1,10,40,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,6,10,40,3),I=1,55) write(*,'(6F11.4)') (XCshift(I,10,10,40,3),I=1,55) print *, ' ' if(bcheck2) then ! Write out input 3-D MHD inputs with the input traced matrix and stop Nit = -1 NitV = -1 NitG = -1 bGcon = .FALSE. bVcon = .FALSE. nMapm = nMap nMapbm = 1 nMapHRm = 1 nMapHRMm = nMap write(*,'(A,3F10.7,3I5)') 'RRS,dRR,RRMS,nMapm,nMapbm,nMapHRMm',RRS, dRR, RRMS, nMapm, nMapbm, nMapHRMm write(*,'(2(A,F11.4))') 'Here 1 NEW DMAP(22,5,22) ',DMAP(22,5,22), ' NEW DMAPV(22,5,22)',DMAPV(22,5,22) end if go to 9999 ! check the new initiation set from an externally read input file and stop. end if call Say('IPSD2020','S','Stop','program successfully completed') stop 9985 continue call Say('IPSD2020_MHD','S','Stop','program terminated before external read') stop 9987 continue call Say('IPSD2020_MHD','S','Stop','program terminated after external read') stop end