C+ C NAME: C ipstd_20n C PURPOSE: C Construct deconvolved synoptic maps for the IPS observations using C a "time-dependent" technique that projects all observations to a C single lower reference surface. Note - uneven time cadences are C allowed for velocity and density data sets. C CATEGORY: C Data processing C CALLING SEQUENCE: C run ipstd_20n 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 extractd from the ea_file but maybe something C was different with ulysses. At this stage, this program does not give exactly C the same thing as the current version on the vax, but it is close. Prior I tried C using. Found an error 10/15/00 (A gridsphere2d was placed outside a do loop after C the first MKshiftd0. This version should now have no seam in each temporal value. C C > The modified Julian day MJD used for the IPS data is related to the C regular Julian day (as used in subroutine Julian) by C MJD = JD-2400000.5 C (i.e. the modified Julian day starts at midnight). C > Once the logical $CAM is set (during the first run of IPSD2020), the C only way to force IPSD2020 to look in another directory is to delete the C logical $CAM. C MODIFICATION HISTORY: C JUN-1992, Paul Hick (UCSD) C NOV-1995, Bernie Jackson (STEL) C March-1999 Bernie Jackson (UCSD) C June -1999 Bernie Jackson (UCSD) C August -1999 Bernie Jackson (UCSD) C Jan/Feb 2001 Bernie Jackson (UCSD) C Mar 2001 Bernie Jackson (UCSD) C- program ipstd_20n parameter (NLmaxG = 30000, ! Max # G data points (13000) (400,000) & NLmaxV = 30000, ! Max # V data points (2000) & NGREST = 1, ! G temporal resolution factor better than 1 day & NVREST = 1, ! V temporal resolution factor better than 1 day & NGRESS = 1, ! Spatial G resolution factor better than 20 degrees & NVRESS = 1, ! Spatial V resolution factor better than 20 degrees & NF = 2, ! LOS resolution factor (Should be 1 if Spatial resolution is >= 2) & nTmaxG = 55*NGREST, ! Max # G times (50) (100) (150) & nTmaxV = 55*NVREST) ! Max # V times (50) (100) (150) parameter (NLOSG = 20*NF*NGRESS, ! # resolution elements along each G line of sight (40) (80) (125) & NLOSV = 20*NF*NVRESS, ! # resolution elements along each V line of sight (40) (80) (125) & NLOSGP1 = NLOSG+1, & NLOSVP1 = NLOSV+1, & dLOSG = 0.1/(NGRESS*NF), ! LOS G Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & dLOSV = 0.1/(NVRESS*NF), ! LOS V Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & xInc = 0.5) ! Carrington map "extra" increment parameter (nLng = 54*NGRESS+1, ! # Longitudes, resolution is 360/(nLng-1) (55) (109) (217) (325) & iLng = 18*NGRESS+1, ! # Longitudes to use in extracted map (19) (37) (73) (109) & nLat = 9*NGRESS+1, ! # Latitudes , resolution is 180/(nLat-1) (9)(10) (19) (37) (55) & nMap = 31, ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) & nLngLat = nLng*nLat, & nLngLatnTG =nLng*nLat*nTMaxG, & nLngLatnTV =nLng*nLat*nTMaxV) parameter (NintLn = 3, ! High resolution output parameters & NintLa = 3, & NintH = 3) parameter (dRR = 0.1) ! Resolution (AU) in radial direction (0.1) parameter (NTP = 14) ! number of body locations possible to extract. integer iXPG (NLmaxG), ! Array indices from daily G maps & iYPG (NLmaxG), & iMJDG (NLmaxG), ! Original MJD of daily G maps & IYRFG (NLmaxG), & IRECG (NLmaxG), & IYRSG (NLmaxG), ! Year of source observation & IYRSGG (NLmaxG), ! Year of source observation & NBSG (NLmaxG), & LON (NLOSGP1,NLmaxG), & LAT (NLOSGP1,NLmaxG), & ITIM (NLOSGP1,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) & GOBS (NLmaxG), ! Solar wind observed G values & GGOBS (NLmaxG), ! Solar wind observed G values c & VOBS (NLmaxV), ! Solar wind observed velocities & VOBS (NLmaxG), ! Solar wind observed velocities & GMOD2 (NLmaxG), ! Model G values & FIX (NLmaxG), ! Fix to model to make OK & GSSIG (NLmaxG), ! G-level source sigma & WTSG (NLOSG ,NLmaxG), & GWTij (NLOSG ,NLmaxG), & XLON (NLOSGP1,NLmaxG), & XLONP (NLOSGP1,NLmaxG), & XLAT (NLOSGP1,NLmaxG), & RPX (NLOSGP1,NLmaxG), & XLproj (NLOSGP1,NLMAXG), & VFAC (NLOSGP1,NLmaxG), ! New array (3/99) & DFAC (NLOSGP1,NLmaxG) character SRCV (NLmaxV)*9, ! Velocity source name & 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 & 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 & VVOBS (NLmaxG), & VMOD (NLmaxV), ! Model velocities & VSSIG (NLmaxV), ! Velocity source sigma & VWTij (NLOSV,NLmaxV), & WTSV (NLOSV,NLmaxV) real XCinTV (nTmaxV), & XCinTG (nTmaxG), & DOYV (nTmaxV), & DOYG (nTmaxG), & XCbeV (2,nTmaxV), & XCbeGG (2,nTmaxG), & XCbe (2,4*nTmaxG) integer IYRV (nTmaxV), & IYRG (nTmaxG), & NMAPD (nLng,nLat,nTmaxG), ! Scratch integer arrays & NMAPV (nLng,nLat,nTmaxV) real VDEV (nLng,nLat,nTmaxV), & VPNT (nLng,nLat,nTmaxV), & PPGMAP (nLng,nLat,nTmaxG), ! Point-P G map (with holes) & VMAPHOLE(nLng,nLat,nTmaxV), & DM (nLng,nLat), & VM (nLng,nLat), & DMHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1), & VMHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1), & DMAPHOLE(nLng,nLat,nTmaxG), & DMAP (nLng,nLat,nTmaxG), ! Density map (no holes) & VMAP (nLng,nLat,nTmaxV), & VMAPD (nLng,nLat,nTmaxG), & DMAPV (nLng,nLat,nTmaxV), & GDEV (nLng,nLat,nTmaxG), & GPNT (nLng,nLat,nTmaxG), & DDfact (nLng,nLat,nMap,nTmaxG), & DVfact (nLng,nLat,nMap,nTmaxG), ! New array (3/99) & XCshift (nLng,nLat,nMap,nTmaxG,3), & V3D (nLng,nLat,nMap), & D3D (nLng,nLat,nMap), & V3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & D3DHR (nLng*(NintLn+1)+1,nLat*(NintLa+1)+1,nMap*(NintH+1)+1), & 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), & Scale (2) 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. & DOYSsav8, & dorbit real*8 XCtpr (NLOSGP1,NLMAXG), & XCtim (NLOSGP1,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, & bGCamb, bGOoty, bGNago, & 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 & 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./, & bWrite3D_HR /.TRUE./, & bWrite3D_HRb /.FALSE./, & bWrite3Do /.TRUE./, & bWriteMS /.FALSE./, & bWritelos /.FALSE./, & bOneOr1 /.TRUE./, & bOneOr2 /.TRUE./, & bForecast, & bSource /.TRUE./, ! Write out good sources & bExtract(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,2*.FALSE./, & bEWRITE, ! Write out an Earth file & bMWRITE, ! Write out a Mercury file & bMAWRITE, ! Write out a Mars file & bUWRITE, ! Write out a Ulysses file & bSAWRITE, ! Write out a STEREO A file & bSBWRITE, ! Write out a STEREO B file & bExtractf(NTP) /2*.FALSE.,.TRUE.,9*.FALSE.,2*.FALSE./, & bEWRITEF, ! Write out an Earth filled file & bMWRITEF, ! Write out a Mercury filled file & bMAWRITEF, ! Write out a Mars filled file & bUWRITEF, ! Write out a Ulysses filled file & bSAWRITEF, ! Write out a STEREO A filled file & bSBWRITEF ! Write out a STEREO B filled file logical BField_Choose ! Logical function equivalence (bExtract( 3), bEWRITE), & (bExtract( 1), bMWRITE), & (bExtract( 4), bMAWRITE), & (bExtract(10), bUWRITE), & (bExtract(13), bSAWRITE), & (bExtract(14), bSBWRITE), & (bExtractf( 3), bEWRITEF), & (bExtractf( 1), bMWRITEF), & (bExtractf( 4), bMAWRITEF), & (bExtractf(10), bUWRITEF), & (bExtractf(13), bSAWRITEF), & (bExtractf(14), bSBWRITEF) character cPrefix (NTP)*2 /'Me','Ve','Ea','Ma','Ju','Sa','Ur','Ne','Pl','Ul','H1','H2','Sa','Sb'/ character cPrefixf(NTP)*2 /'m1','v2','e3','m4','j5','s6','u7','n8','p9','uy','ha','hb','s1','s2'/ integer iEorWG /0/, & iEorWV /0/, c & MJD, c & MJDfrst /48180/, ! MJD for earliest data file c & MJDlast /49622/, ! MJD for latest data file & MJDfrst / 46100/, & MJDlast / 51800/, & LimitoV /-1/, & LimitoG /-1/, & NLOSWV /1/, & NLOSWG /1/, & NTV /44/, & NTG /44/ integer NiterT /18/, & iOffUT / 0/, & NAV / 0/, & NC / 1/, ! # rotations averaged & nDUMP / 0/, ! Threshold for #pnts/bin & nFILL / 3/, ! Min. # neighbours for fill & NmidHRV /0/, & NmidHRG /0/ real YLbeg /-90./, & YLend / 90./, & XCadj(3) /3*0./, & XCadjB(3) /3*0./, & XCmargin /.5/, & XElowG /11.5/, & XElowV /11.5/, & XEhighG /80./, & XEhighV /180./, & XRlimG /90./, & XRlimV /90./, & GPOWER /0./, & R1AU /1.0/, & DEN1AU /5.0/, & ANLIMITV /1.2/, & ANLIMITG /1.2/, & aNdayV /1.0/, & aNdayG /1.0/ real VMAPSIG /1.0/, & DMAPSIG /1.0/, & GSIG /1.0/, & GSIGB /3.0/, & VSIG /1.0/, & VSIGB /3.0/, & FCONSV /1.0/, & FCONSG /1.0/, & CONRV /14./, & CONRD /14./, & CONRO /8.0/, & CONVT /0.65/, & CONDT /0.65/, & CONOT /0.45/, & ConsTV /1.0/, & ConsTG /1.0/, & VLIMM /1200./, ! Velocity upper limit for sources & VLIML /201./, ! Velocity lower limit for sources & VLIM /1200./, ! Velocity upper limit for tomography & VLIL /100./, ! Velocity lower limit for tomography & GLIMP /4.0/, ! G-level upper limit for sources & GLIMF /1.0/, ! G-level input multiplier & GLIMM /0.2/ ! G-level lower limit for sources & real RRS / -15.0 /, ! Reference distance for deconvolution & RRX / 1.0 /, ! Reference distance for Map extraction & VUL / 1500.0 /, ! Velocity upper limit for maps & VLL / 100.0 /, ! Velocity lower limit for maps & DUL / 500.0 /, ! Density upper limit for maps & DLL / 0.0 / ! Density lower limit for maps real DAV /1./, & Speed /400./, & FALLOFF /2.0/, & PWG /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.35/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.20/, ! Radial power of density & PWRV /0.20/, ! Radial power of velocity & AEmaxV /800.0/, ! Anti-Earth maximum velocity & AEmaxD /10.0/, ! Anti-Earth maximum density at 1 AU & AEminV /200.0/, ! Anti-Earth minimum velocity & AEminD /2.0/, ! Anti-Earth minimum density at 1 AU & AEangV /135./, ! Angular distance from Earth to begin V limit & AEangD /135./ ! Angular distance from Earth to begin D limit real*8 MJDref, & MJDcntr, & XCtbegV, & XCtbegG, & XCtendV, & XCtendG character cCam*4, & cMon*3, & cStrID*11, & cVar(3)*100, & cArg*100, & cStr*80, & cPre*4 character cWildNagoya*80, & cWildOoty*80 external iReadNagoyan8, iProcessNagoyan8 external iReadOotyn8, iProcessOotyn8 external EARTH include 'dirspec.h' include 'sun.h' include 't3d_array.h' 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 & BR2DT(nLng,nLat,nTmaxG), !scratch & TT3D (4*nTmaxG), !holds XC corresponding to 6 hour int. time steps output by write3dinfo_td3d & XC3DT (nLng,nLat,nMap,4*nTmaxG,3), & BRcheck (nLng,nLat), & Vcheck5 (nLng,nLat), & Vcheck7 (nLng,nLat) call T3D_set_mode(T3D__MODE_TIME, 1) iVar = 3 call ForeignArg(' ',iVar,cVar,cArg) C------- C Calculate start times for Carrington rotations between 1985 and now. C Rotation NCoff+I starts at time JDCar(I). aNdayV = aNdayV/NVREST aNdayG = aNdayG/NGREST ConsTMV = ConsTV*nTV ConsTMD = ConsTG*nTG ConstL = nLng/3 I = iDeleteSymbol(cSym,2) 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,5000.,XCEV) ! Initialize XCEV array call ArrR4Constant(NLmaxG,5000.,XCEG) ! Initialize XCEGG array call ArrR4Constant(NTmaxV,5000.,XCintV) ! Initialize XCintV array call ArrR4Constant(NTmaxG,5000.,XCintG) ! Initialize XCintG array C------ C iOffUT is an offset time (in hours) between local time and UT. This is used only if the program C is operated in forecast mode. Note that the local time is assumed to be the same as the C computer time. The default offset is 7 hours (the difference between PST and GST). call AskYN('Forecast mode?$no',bForecast) ! Regular or forecast mode if (bForecast) call AskI4('Offset between local time and UT (hours)?',iOffUT) call AskI4('Number of iterations?',NiterT) ! # iterations call AskR4('Deconvolution distance (AU; 0=solar surface)?',RRS) if (RRS .lt. 0) RRS = -RRS*SUN__RAU ! Neg. RRS are in units of solar radii RRS = max(RRS,sngl(SUN__RAU)) ! Keep RRS > 1 solar radius ! Deconvolve velocity, G-level or both? call AskWhat('Deconvolve: velocity, 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, Ooty, Nagoya$1$3',I) bVUCSD = I .eq. 1 bVOoty = I .eq. 2 bVNago = I .eq. 3 if(bVUCSD) then Print *, 'Sorry: The UCSD Velocity reconstructions are not yet supported. Try another.' stop end if bVmv2 = .FALSE. if(bVNago) NLOSWV = 1 if(bVOoty) NLOSWV = 1 else ! No V deconvolution bVNago = .FALSE. bVOoty = .FALSE. call AskYN('Do you want to make a velocity map using mv^2 = con$no',bVmv2) end if if (bGcon) then ! G deconvolution: select IPS station if(bVOoty) call AskWhat('Use G-level data from: Cambridge, Ooty, Nagoya$1$2',I) if(bVNago) call AskWhat('Use G-level data from: Cambridge, Ooty, Nagoya$1$3',I) if(.not.bVNago.and..not.bVOoty) call AskWhat('Use G-level data from: Cambridge, Ooty, Nagoya$1$3',I) bGCamb = I .eq. 1 bGOoty = I .eq. 2 bGNago = I .eq. 3 bGmv2 = .FALSE. cCam = cEnvi//'CAM' if(bGNago) NLOSWG = 1 if(bGOoty) NLOSWG = 1 if(bGcamb) NLOSWG = 2 else ! No g deconvolution bGCamb = .FALSE. bGOoty = .FALSE. bGNago = .FALSE. call AskYN('Do you want to make a density map using mv^2 = con',bGmv2) end if if (bGCamb) call SetGipsy(cCam,iOffUT,iEdt,MJDfrst,MJDlast) if (bVNago .or. bGNago ) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya) if (bVOoty .or. bGOoty ) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty ) if(bVcon) then bDaydV = .TRUE. call AskYN('Deconvolve V from one day to the next$yes',bDaydV) if(bDaydV) then call AskR4('Number of days to combine for velocity?',aNdayV) if(bVNago) NmidHRV = -9 if(bVOoty) NmidHRV = -6 call AskI4('Number of hours from V UT midnight?',NmidHRV) end if call AskI4('Number of velocity time values?',NTV) ! # NTV time itervals if(NTV.gt.NTmaxV-2) then NTV = NTmaxV-2 call AskI4('Number of velocity time values?',NTV) ! Ask again end if end if if(bGcon) then bDaydG = .TRUE. call AskYN('Deconvolve 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(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 call Local2UT(iOffUT,iYr,Doy) XChi = max(NCoff+XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar)-1.,XClo) RRV = RRS RRG = RRS if (bVcon) then call AskR4('# of V-points per bin that constitute a deconvolution',ANLIMITV) if(NiterT.gt.8) LimitoV = NiterT/2 call AskI4('Iteration at which to limit V-source deviations, -1=none',LimitoV) C------- C Fix the reference distance used for the V map formed. C Get E or W or EW; Set the elongation angles call AskWhat('All velocity data, East data, West data ?$0',iEorWV) if (iEorWV .eq. 2) iEorWV = -1 call AskR4('Minimum elongation for V-data$0$180$',XElowV) call AskR4('Maximum elongation for V-data$'//cStr(:Flt2Str(XElowV,1,cStr))//'$180$',XEhighV) call AskR4('Anti-Earth maximum velocity at 1 AU$'//cStr(:Flt2Str(600.0,1,cStr))//'$1000.0$',AEmaxV) call AskR4('Anti-Earth minimum velocity at 1 AU$'//cStr(:Flt2Str(0.,1,cStr))//'$400.0$',AEminV) call AskR4('Angular distance from Earth to begin V limit$'//cStr(:Flt2Str(90.,1,cStr))//'$180.0$',AEangV) call AskR4('Maximum RA relative to Sun for V-data$0$180$',XRlimV) RRV = RRS call AskR4('V-map reference distance in AU (0=solar surface)?',RRV) if (RRV .lt. 0) RRV = -RRV*SUN__RAU RRV = max(RRV,sngl(SUN__RAU)) call AskYN('Do you want to use a proxy map for velocity?',bVproxy) call AskR4('Maximum upper limit for tomography velocity?',VLIM) call AskR4('Maximum lower limit for tomography velocity?',VLIL) call AskR4('Maximum upper limit for source velocity?',VLIMM) if(bVNago) VLIML = VLIL call AskR4('Maximum lower limit for source velocity?',VLIML) end if if (bGcon) then 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 = 11.5 if(bGOoty) XElowG = 11.5 if(bGCAMB) XElowG = 30.0 call AskR4('Minimum elongation for G-data$0$180$',XElowG) if(bGNago) XEhighG = 180. if(bGOoty) XEhighG = 180. if(bGCAMB) 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) call AskR4('Factor by which to modify G-level?',GLIMF) end if call AskR4('G-level time hole filter?',ConsTG) call AskR4('Velocity time hole filter?',ConsTV) if(bVOoty) CONRD = CONRD call AskR4(' G-level spatial filter?',CONRD) if(bVNago) CONRV = CONRD if(bVOoty) CONRV = CONRO if(CONRD.ne.CONRV) CONRV = CONRD call AskR4('Velocity spatial filter?',CONRV) if(bVOoty) CONDT = CONDT call AskR4(' G-level temporal filter?',CONDT) if(bVNago) CONVT = CONDT if(bVOoty) CONVT = CONOT if(CONVT.ne.CONDT) CONVT = CONDT call AskR4('Velocity temporal filter?',CONVT) if(bGNago) PWG = 0.35 if(bGOoty) PWG = 0.35 if(bGCamb) PWG = 0.15 call AskR4(' G-level power to fit density?',PWG) if(bVNago) 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(bVNago.and.bGNago) 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) ! 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) call Julian(11,iYr,Doy,MJDref,JEpoch) call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) C------- Get modified Carrington variable for center of plot C [Note that EarthLng=(1-frac(XCEarth))*360] XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XCbeg = XCEarth-.5 call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) 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 else ! Normal mode write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$' XCbeg = max(XClo,NCoff+XCbeg) XCbeg = 1964.6 if(bGCamb) XCbeg = 1884 call AskR4('Carrington rotation (0 = STOP)'//cStr,XCbeg) if (XCbeg .le. 0) call Say('IPSD2020','I','StopPlay','StopPlay') NintF = 1 if(bForecast) then nTF = nTG - 1 NintF = 0 end if call AskI4('# Intermediate rotations interpolated',NintF) call AskI4('# Rotations averaged',NC) NC = min(NC,1+int(XChi-XCbeg) ) ! Stay within available data range XCbeg = XCbeg-NCoff call AdjustJDCar(XCbeg, nCar, JDCar, NCoff) C------- Info only: start time of Carrington rotation call Julian(1,iYr,Doy,JDCar(int(XCbeg)),JEpoch) iDoy = Doy call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy) write (*,'(/,10X,A,I4,A,I2,3A,I4,A,I3,A,/)') 'Carrington rotation ',NCoff+int(XCbeg), & ' starts on date ',iDay,'-',cMon,'-',iYr,' ( Doy ',iDoy,' )' C------- Get modified Carrington variable for center of plot C Get heliographic coordinates for center of plot (Earth position) XCEarth = XCbeg+.5 ! Center of plot I = XCEarth ! Rotation containing plot center MJDcntr = JDCar(I)+(XCEarth-I)*(JDCar(I+1)-JDCar(I)) ! Julian day for plot center call Julian(1,iYr,Doy,MJDcntr,JEpoch) ! iYr,Doy for plot center MJDcntr = MJDcntr-SUN__MJDtoJD ! Modified Julian day for plot center call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun) EarthLng = mod(sngl(dLngSun)+180.,360.) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) ! Heliographic coord. for plot center XCtest1 = XCbeg-XCmargin ! Search limits XCtest2 = XCbeg+NC+XCmargin XCend= XCbeg + NC XCstrt = XCbeg XCbeg= XCbeg - xInc XCend= XCend + xInc I = XCtest2 MJDref = JDCar(I)-SUN__MJDtoJD+(XCtest2-I)*(JDCar(I+1)-JDCar(I)) ! Start MJD for data search end if C write (cStrID,'(I2.2,A,F8.3)') nint(10*RRS),'_',NCoff+XCbeg call AskYN('Do you want to output model-source files$no',bWriteMS) call AskYN('Do you want to output l.o.s. source surface maps$no',bWritelos) call AskYN('Do you want one shift matrix per iteration$no',bOnly1) call AskYN('Output good source files$yes',bSource) call AskYN('Output three D error files$yes',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?$yes',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 call AskYN('Output a 3D high resolution file$yes',bWrite3D_HR)! 3D "h" file write if(bWrite3D_HR) then nMapb = 1 nMapHR = 31 call AskI4('Minimum high resolution height number - 1-31?',nMapb) call AskI4('Maximum high resolution height number - 1-31?',nMapHR) if(nMapHR.gt.nMap) then nMapHR=nMap nMapb=nMapHR end if end if call AskYN('Output completely filled 3D files$yes',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 nMapHRb = 5 call AskI4('Minimum height number - 1-31?',nMapHRb) if(nMapHRb.gt.nMap) then nMapHRb=nMaps nMapb=nMapHRb end if end if call AskYN('Output an Earth parameter file$yes',bEWRITE) ! Earth extraction call AskYN('Output a Mercury parameter file$no',bMWRITE) ! Mercury extraction call AskYN('Output a Mars parameter file$no',bMAWRITE) ! Mars extraction call AskYN('Output a Ulysses parameter file$no',bUWRITE) ! Ulysses extraction if(iYr.ge.2007) then call AskYN('Output a STEREO A parameter file$yes',bSAWRITE) ! STEREO A extraction call AskYN('Output a STEREO B parameter file$yes',bSBWRITE) ! STEREO B extraction end if call AskYN('Output an Earth parameter filled file$yes',bEWRITEF) ! Earth filled extraction call AskYN('Output a Mercury parameter filled file$no',bMWRITEF) ! Mercury filled extraction call AskYN('Output a Mars parameter filled file$no',bMAWRITEF) ! Mars filled extraction call AskYN('Output a Ulysses parameter filled file$no',bUWRITEF) ! Ulysses filled extraction if(iYr.ge.2007) then call AskYN('Output a STEREO A parameter filled file$yes',bSAWRITEF) ! STEREO A filled extraction call AskYN('Output a STEREO B parameter filled file$yes',bSBWRITEF) ! STEREO B filled extraction end if C Do you want to include a magnetic map at the reference surface height? C If yes, then read the source surface map (at 2.5Rs) at C the resolution of the reference arrays. call AskYN('Do you want to make a magnetic map?',bMagnetic) if (bMagnetic) bMagnetic = BField_Choose(iVar,cVar) 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 if (bVNago) call ReadVIPSn8(iReadNagoyan8,iProcessNagoyan8,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVOoty) call ReadVIPSn8(iReadOotyn8, iProcessOotyn8,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,cvgfiles,SRCVG,SRCV,SRCVGsav,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV8,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, & XDLsav,XCEVsav,XEsav,XCsav,YLsav,VVsav,GGsav) 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) call ArrR4Zero(NLmaxV,VSSIG) call ArrI4Constant(NLmaxV,1 ,NBSV ) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLmaxG,1 ,NBSG ) do I=1,NLmaxV C if (I.lt.100) print *, I, SRCV(I), 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 end do do I=1,NLmaxG GGOBS(I) = GGOBS(I)*GLIMF 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 end do C======= do I=1,NLG 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) 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) then XCG(I) = XCGG(I) YLG(I) = YLGG(I) DISTG(I) = DISTGG(I) GOBS(I) = GGOBS(I) XEG(I) = XEGG(I) XLSG(I) = XLSGG(I) XCEG(I) = XCEGG(I) XLLG(I) = XLLGG(I) XDLG(I) = XDLGG(I) DOYSG8(I) = DOYSGG8(I) IYRSG(I) = IYRSGG(I) end if end do 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 = (NTV + 1.5)/2.0 + 6 NTG = (NTG + 1.5)/2.0 + 6 end if print *, 'MJDcntr', MJDcntr ItoffG = 0 ItoffV = 0 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) print *, 'SRCVG(30)',SRCVG(30) print *, 'Spot check of IPS data - 30', IYRSG(30),IYRSV(30),DOYSG8(30),DOYSV8(30),XEG(30),XEV(30) print *, ' 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 print *, 'XCtbegV, XCtendV, New NTV', XCtbegV, XCtendV, NTV print *, 'XCtbegG, XCtendG, New NTG', XCtbegG, XCtendG, NTG print *, 'XCbeV =' , XCbeV print *, 'XCintV =', XCintV print *, 'XCbeGG =', XCbeGG print *, 'XCintG =', XCintG print *, 'XCintDV =', XCintDV print *, 'XCintDG =', XCintDG XCtbegF = XCtbegG ! For forecast mode XCtendF = XCtendG ! For forecast mode if(bGcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),DMAP) if(bVcon) call ArrR4Constant(nLng*nLat*nTG,BadR4(),VMAP) NLmax = max(NLmaxG,NLmaxV) NTVV = ItoffV NTGG = ItoffG NTVE = 0 NTGE = 0 K = 0 do KKK=1,NLmax+nTG ! Loop before deconvolution to set up VMAP, DMAP K = K + 1 IVSW = 0 ! Switch turn on for NTVV value of K if(K.le.NLV+1) then ! If K is more than the source # + 1 don't turn on switch 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) then NTVV = NTVV + 1 NTVB = NTVE + 1 NTVE = K - 1 if(NTVB.ne.NTVE+1) IVSW = 1 if(K.le.NLV) K = K - 1 NTVD = NTVE-NTVB+1 write (*,'(A,I4,A,2I6,2F10.3)') 'V loop ',NTVV,' - S#, K, XCEV, XCintV', NTGD,K,NCoff+XCEV(K),NCoff+XCintV(NTVV+2) if(NTVV.eq.NTV) NLV = NTVE end if end if if(IVSW.eq.1) then C print *, ' V switch ',IVSW,' loop',NTVV if (NLV .eq. 0) call Say('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 call MkLOSWeights(NLOSWV,dLOSV,NLOSV,NTVD,XEV(NTVB),DISTV(NTVB),WTSV(1,NTVB),PWRV) call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAPHOLE(1,1,NTVV)) if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) if(amax.ne.BadR4()) then call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,NTVV),CONRV,3,0.0,0.0) end if end if call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,NTVV,' min and max is', amin, amax if(NTVV.eq.ItoffV+1) then ! Fill earlier arrays do N=1,ItoffV amaxV = BadR4() if(bVproxy) then I = iReadProxyMapN(1,NCoff+XCintV(N),nLng,nLat,VMAP(1,1,N)) call arrR4getminmax(nLngLat,VMAP(1,1,N),aminV,amaxV) print *, 'The VMAP' ,N,' min and max after read is', aminV, amaxV end if if(amaxV.eq.BadR4()) then call arrR4getminmax(nLngLat,VMAP(1,1,NTVV),amin,amax) print *, 'The VMAP' ,N,' min and max is', amin, amax if(amax.ne.BadR4()) then call ArrR4Copy(nLngLat,VMAP(1,1,NTVV),VMAP(1,1,N)) call Say('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.NLG+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) then NTGG = NTGG + 1 NTGB = NTGE + 1 NTGE = K - 1 if(NTGB.ne.NTGE+1) IGSW = 1 if(K.le.NLG) K = K - 1 NTGD = NTGE-NTGB+1 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) NLG = NTGE end if end if end if if(IGSW.eq.1) then if(NTGB.ne.0) then if (NLG .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 call MkLOSWeights(NLOSWG,dLOSG,NLOSG,NTGD,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 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(.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 do L=1,NLG if(GOBS(L) .ne. BadR4()) NLGG = NLGG + 1 end do NLVV = 0 do L=1,NLV if(VOBS(L) .ne. BadR4()) NLVV = NLVV + 1 end do call Say('IPSD2020','I','Info','# good g-level observations ' & //cStr(:Int2Str(NLGG,cStr))) call Say('IPSD2020','I','Info','# good V observations ' & //cStr(:Int2Str(NLVV,cStr))) 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 MkShiftdn0(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP min max N',amin,amax,N C call arrR4getminmax(nLngLat,DMAPV(1,1,N),amin,amax) C print *, 'DMAPV min max N', amin,amax,N call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV,1,0.0,0.0) end do C C Comment out call MkDensity(nLng,nLat,nMap,XCshift,DDfact,VMAP,RRV,dRR) C Nit = -1 ! Why three iterations counters?? NitV = -1 NitG = -1 do while (Nit .lt. NiterT) Nit = Nit+1 print *, ' ' call Say('IPSD2020','I','Info','Iteration '//cStr(:Int2Str(Nit,cStr))) NVS = 0 if (Nit .eq. LimitoV) NVS = 1 NGS = 0 if (Nit .eq. LimitoG) NGS = 1 if (bVcon) then ! V deconvolution NitV = NitV+1 ! # V iterations ! bGmv2=.TRUE., implies bGcon=.FALSE. if (bGmv2) then 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(XCbeV,XCtbegV,XCtendV,IdaybegG,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) C print *, 'Spot check after Mkpostd_V - 1,30', ITIM(1,1),ITIM(1,30),IYRSV(1),IYRSV(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30) if(NitV.eq.NiterT.and.bWritelos) call writelosmapvg(0,NLOSVP1,NLV,nTV,NBSV,NCOFF,ITIM,XLON,XLAT) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(1,XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DVfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,VFAC) call MkVModeltd(XCbeV,XCtbegV,XCtendV,XEV,XCtpr,XLON,XLproj,RPX, & WTSV,DFAC,VFAC,PWV,NLV,NLOSV,NLOSVP1,VMAP,VLIM, & DMAPV,nLng,nLat,nTV,RRS,VMOD,VWTij) call FixModeltdn(1,NLV,VOBS,VMOD,PWV,NBSV,VSIG,FIX,VSSIG,VRATIO) if(bWriteMS.and.NitV.eq.NiterT) then call writeM_Scomp8(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV,PWV,VRATIO, & NLmaxV,NSRV,SRCV,DOYSV8,XCEV,VOBS,VMOD,XLON,XLproj,XCtpr,NLOSVP1,FIX,VSSIG,NBSV) end if VRATION(NitV+1) = VRATIO ! Put in 7/9/00 if(NitV.ge.2.and.VRATION(NitV+1).gt.(VRATION(NitV)+(VRATION(NitV+1)*0.05)).and. & VRATION(NitV+1).gt.(VRATION(NitV-1)+(VRATION(NitV+1)*0.10))) then print *, 'THE SUMMED VELOCITY RATIO IS NO LONGER CONVERGING' C NitV = NiterT end if C print *, 'Into MkVmap',nLng,nLat,nTV,NLOSV,NLV call MkVMaptdn0n(LON,LAT,ITIM,RPX,XEV,VWTij,FIX,NLV,NLOSV,NLOSVP1,XLON,XLAT, & nLng,nLat,nTV,VMAP,ANLIMITV,VLIM,VSSIG,VMOD,VSIGB,NVS,NBSV,VMAPSIG, & ANMAPV,AWTMAPV,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV, & Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp) C print *, 'Out of MkVmap' do N=1,NTV C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C if(amax.eq.BadR4()) then C print *, 'VMAP is bad before map determination N =', N,' Iter =', Nit C end if C call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'In V determination VMAP min max',amin,amax if (NitV .eq. NiterT) call ArrR4Copy(nLngLat,VMAP(1,1,N),VMAPHOLE(1,1,N)) call ArrR4SetMinMax(-nLngLat,VMAP(1,1,N),VLL,VUL) end do end if if (.not.bDproxy.or.Nit.ge.1) then if (bGcon .or. bVmv2) then ! Update density and density shifts even on the last call FillWholeT(0,Nit,NiterT,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) call FillMaptN(0,XCbeGG,nLng,nLat,nTG,ConsTG,DMAP,DDD) do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before Grids MkShiftd0 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 MkShiftd0 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 MkShiftd0 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 MkShiftd0 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 MkShiftd0 VMAPD min max N',amin,amax,N end do end if end if if(bGcon) then call MkShiftdn0(XCbeGG,XCintG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG,nTmaxG, & nCar,JDCar,NCoff,VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,Vtmp,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C print *, VMAPD (30,5,20),VMAPD (31,5,20),VMAPD (30,6,20),VMAPD (31,6,20),VMAPD (30,5,21),VMAPD (31,5,21), C & DMAP (30,5,20),DMAP (31,5,20),DMAP (30,6,20),DMAP (31,6,20),DMAP (30,5,21),DMAP (31,5,21), C & XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11), C & DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11), C & DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11) end if if (bGcon) then ! G deconvolution NitG = NitG+1 ! # G iterations ! bVmv2=.TRUE., implies bVcon=.FALSE. if (bVmv2) then RRV = RRD call Deneq3(nLng,nLat,NTG,DMAP,RRD,DENEQ) print *, 'Solar equatorial mv^2 = constant value DENEQ =',DENEQ do N=1,NTG call Mk_D2VN(nLngLat,DMAP(1,1,N),VMAPD(1,1,N),DENEQ,Speed,RRD) end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAPD,VMAP) else call CopyDtoDVN(1,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG, & nLng,nLat,NTG,NTV,ConsTG,VMAPD,DMAPV,Vtmpt,DDT,DDD) call ArrR4Copy(nLngLatnTV,DMAPV,VMAP) end if if(NitG .eq. NiterT) then do N=1,NTV call ArrR4Copy(nLngLat,DMAPV(1,1,N),VMAPHOLE(1,1,N)) end do end if end if call MkPostd(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) C print *, 'Spot check after Mkpostd_G - 1,30', ITIM(1,1),ITIM(1,30),IYRSG(1),IYRSG(30),XLON(1,1),XLON(1,30),XLAT(1,1),XLAT(1,30) if(NitG.eq.NiterT.and.bWritelos) call writelosmapvg(1,NLOSGP1,NLG,nTG,NBSG,NCOFF,ITIM,XLON,XLAT) call Get4Dval(1,XCbeGG,XCtbegG,XCtendG,RRG,dRR,nLng,nLat,nMap,nTG, & DDfact,NLOSGP1*NLG,XCtim,XLONP,XLAT,RPX,DFAC) C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkGModeltdn DMAP min max N',amin,amax,N C end do call MkGModeltdn(XCbeGG,XCtbegG,XCtendG,XCtpr,XLON,XLproj,RPX, & WTSG,DFAC,PWG,DEN1AU,NLG,NLOSG,NLOSGP1,DMAP, & nLng,nLat,nTG,RRS,GMOD2,GWTij,WTPD,GMWTi) call FixModeltdn(2,NLG,GOBS,GMOD2,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) if(bWriteMS.and.NitG.eq.NiterT) then call writeM_Scomp8(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO, & NLmaxG,NSRG,SRCGG,DOYSG8,XCEG,GOBS,GMOD2,XLON,XLproj,XCtpr,NLOSGP1,FIX,GSSIG,NBSG) end if GRATION(NitG+1) = GRATIO ! Put in 7/9/00 if(NitG.ge.2.and.GRATION(NitG+1).gt.(GRATION(NitG)+(GRATION(NitG+1)*0.05)).and. & GRATION(NitG+1).gt.(GRATION(NitG-1)+(GRATION(NitG+1)*0.10))) then print *, 'THE SUMMED DENSITY RATIO IS NO LONGER CONVERGING' C Nit = NiterT end if call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,AWTMAPD,DMAD,FIXMD,WTSMD,DMEANFD,DFIXM2D,WTPD, & Alng,ConstL,CONRD,CONDT,aNdayG,Dtmp) do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'After MkDmaptd0 DMAP min max',amin,amax if (NitG .eq. NiterT) call ArrR4Copy(nLngLat,DMAP(1,1,N),DMAPHOLE(1,1,N)) call ArrR4SetMinMax(-nLngLat,DMAP(1,1,N),DLL/(RRG*RRG),DUL/(RRG*RRG)) end do end if if(Nit.ne.NiterT) then ! Update velocity shifts and vmap except on last 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) 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) do N=1,NTV call FillMapL(3,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,DMAPV(1,1,N)) end do call MkShiftdn0(XCbeV,XCintV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV,nTmaxV, & nCar,JDCar,NCoff,VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,FALLOFF,CONRV,CONRD,VLL,VUL, & NSIDE,VtmpV,Dtmp,XLTtmp,XLtmp,XLtmpt,Dtmpt) C print *, VMAP (30,5,20),VMAP (31,5,20),VMAP (30,6,20),VMAP (31,6,20),VMAP (30,5,21),VMAP (31,5,21), C & DMAPV (30,5,20),DMAPV (31,5,20),DMAPV (30,6,20),DMAPV (31,6,20),DMAPV (30,5,21),DMAPV (31,5,21), C & XCshift(30,5,20,11),XCshift(31,5,20,11),XCshift(30,6,20,11),XCshift(31,6,20,11),XCshift(30,5,21,11),XCshift(31,5,21,11), C & DVfact (30,5,20,11),DVfact (31,5,20,11),DVfact (30,6,20,11),DVfact (31,6,20,11),DVfact (30,5,21,11),DVfact (31,5,21,11), C & DDfact (30,5,20,11),DDfact (31,5,20,11),DDfact (30,6,20,11),DDfact (31,6,20,11),DDfact (30,5,21,11),DDfact (30,5,21,11) end if end if end if end do ! End of iteration loop if(bSource) then ! Write out good sources used istation = 1 if(bGNago) istation = 1 if(bGOoty) istation = 2 call writegoodsourceg(istation,cggfiles,SRCGV,GOBS,NBSG,NLG) if(bVNago) istation = 1 if(bVOoty) istation = 2 call writegoodsourcev(istation,cvgfiles,SRCVG,VOBS,NBSV,NLV) end if J = 0 ! List bad IPS sources do I=1,NLV if (NBSV(I) .eq. 0) then J = J+1 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 end do write (*,'(A,I4,A)') ' There were',J,' bad V sources removed.' J = 0 do I=1,NLG if (NBSG(I).eq.0) then J = J+1 if(GOBS(I).eq.BadR4()) GOBS(I) = -99.999 write (*,'(3A,I5,A,F7.3)') ' Bad G ',SRCGG(I),' source removed, # =', & I,', elongation =',XEG(I) end if end do write (*,'(A,I4,A)') ' There were',J,' bad G sources removed.' C Here VMAPHOLE contains the original point-P map (bVcon=.FALSE.) ( or is equal to VMAP do N=1,NTV call ArrR4Copy(nLngLat,VMAPHOLE(1,1,N),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 Extractdn',XCbegG(1,1),XCtbegG,XCstrt C***** Take data out at different spacecraft (and earth) locations (bExtract tells where) NinterD = 3 call Extractdn(bExtract,cPrefix,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) C***** Write to disk of results interpolated to the even latitudes, longitudes C and times over the heights (specified as nMapn) of the XCshift model. 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(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' if(bWrerr) then C These (mode 2) write3D_infotd calls make 3D arrays with time error maps to show reconstruction reliability do N=1,nTG do I=1,nLng do J=1,nLat if(ANMAPD(I,J,N).eq.0.0) ANMAPD(I,J,N) = badR4() VDEV(I,J,N) = ANMAPD(I,J,N) if(AWTMAPD(I,J,N).eq.0.0) AWTMAPD(I,J,N) = badR4() GDEV(I,J,N) = AWTMAPD(I,J,N) end do end do end do call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,VDEV,ANMAPD,Vtmpt,DDT,DDD) call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,GDEV,AWTMAPD,Vtmpt,DDT,DDD) PIRSQD = (CONDT/aNdayG)*3.14159*((CONRD/(20.0/float(NGRESS)))**2) ! normalize weights and lines of sight to spatial resolution do N=1,nTG call arrR4getminmax(nLngLat,ANMAPD(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,ANMAPD(1,1,N),NSIDE,Zmin,Zmax) ! Remove small ANMAPD holes call ArrR4TimesConstant(-nLngLat,ANMAPD(1,1,N),PIRSQD,ANMAPD(1,1,N)) end if call arrR4getminmax(nLngLat,AWTMAPD(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,AWTMAPD(1,1,N),NSIDE,Zmin,Zmax) ! Remove small AWTMAPD holes) call ArrR4TimesConstant(-nLngLat,AWTMAPD(1,1,N),PIRSQD,AWTMAPD(1,1,N)) end if end do do N=1,nTV do I=1,nLng do J=1,nLat if(ANMAPV(I,J,N).eq.0.0) ANMAPV(I,J,N) = badR4() VDEV(I,J,N) = ANMAPV(I,J,N) if(AWTMAPV(I,J,N).eq.0.0) AWTMAPV(I,J,N) = badR4() GDEV(I,J,N) = AWTMAPV(I,J,N) end do end do end do call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,VDEV,ANMAPV,Vtmpt,VVT,VVV) call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,GDEV,AWTMAPV,Vtmpt,VVT,VVV) PIRSQV = (CONVT/aNdayV)*3.14159*((CONRV/(20.0/float(NVRESS)))**2) ! normalize weights and lines of sight to spatial resolution do N=1,nTV call arrR4getminmax(nLngLat,ANMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,ANMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small ANMAPV holes call ArrR4TimesConstant(-nLngLat,ANMAPV(1,1,N),PIRSQV,ANMAPV(1,1,N)) end if call arrR4getminmax(nLngLat,AWTMAPV(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,AWTMAPV(1,1,N),NSIDE,Zmin,Zmax) ! Remove small AWTMAPV holes call ArrR4TimesConstant(-nLngLat,AWTMAPV(1,1,N),PIRSQV,AWTMAPV(1,1,N)) end if end do NinterCC = 0 cPre = 'er1_' print *, 'Before the error1 (L.O.S. crossings) write' call write3D_infotd3DM(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(2,cPre,Nit,NiterT,NinterCC,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,AWTMAPV,AWTMAPD,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) end if if(.not. bForecast) then print *, ' Before Extractd3D',XCbegG(1,1),XCtbegG,XCstrt call Extractd3D(bExtract,RRS,dRR,XCbe,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterC+1,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar) end if call Timesmooth(nLng,nLat,nTV,VMAP,2,CONVT/aNdayV,0.75*CONVT/aNdayV,VZtmp) ! Fill in times one time away do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy SpeedV = 0.0 call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP) call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam end if call GridSphere2D(ALng,nLng,nLat,1,VMAP(1,1,N),CONRV/2.0,0,0.0,0.0) ! Change 5/9/08 to smooth polar holes end do call Timesmooth(nLng,nLat,nTG,DMAP,2,CONDT/aNdayG,0.75*CONDT/aNdayG,Dtmp) ! Fill in times one time away do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy AEmax = AEmaxD/(RRS*RRS) AEmin = AEmin/(RRS*RRS) SpeedD = 0.0 call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP) call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if call GridSphere2D(ALng,nLng,nLat,1,DMAP(1,1,N),CONRD/2.0,0,0.0,0.0) ! Change 5/9/08 to smooth polar holes end do if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTV,VMAP,VMAPD,Dtmpt,VVT,VVV) end if C The second write3D_infotd makes maps that have fewer holes to be used to show as remote observer views C The first write3D_infotd makes maps that have holes to be used to show Carrington maps cPre = 'nv3f' call write3D_infotd3D(NOneOr2,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,LLbeg,LLend,TT3D,XC3DT,V3DT,D3DT) C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFF 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 C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapb,nMapHR,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if C Must change back to regular value in current version RRSCON = (R1AU/RRS)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change if(bForecast) then C print *, 'DATA EXTRACTED FROM FINAL MATRIX FOR FORECAST' print *, ' Before 2nd Extractdn',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdn(bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if ************************************************************************************************************************************* C Produce a set of nv3o files that are completely filled as well as we can make them. C The third write3D_infotd makes maps that have no holes to be used to provide the best possible boundary for a 3D MHD model if(bWrite3Do.or.bWrite3D_HRb.or..not.bForecast) then C ixxxxxxx = 0 C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop print *, 'Before third filling' do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP N, amin and amax = ', N, amin, amax if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small velocity holes call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) C print *, 'VMAP N, amin and amax = ', N, amin, amax if(amax.ne.BadR4()) then call GridFill(0,nLng,nLat,VMAP(1,1,N),NSIDE,Zmin,Zmax) ! Fill ALL velocity holes in any map with a hole C call FillMapL(5,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove 90 deg discrepancy C C SpeedV = 0.0 C call FillMapOLS(1,XCbeV,XCintV,nLng,nLat,nTV,N,NTmaxV,RRS,SpeedV,ConstL,AEmaxV,AEminV,AEangV,VMAP) C call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam end if end if end do call FillWholeT(1,0,0,XCbeV,nLng,nLat,nTV,ConsTMV,VMAP,VVV) ! Fill all the velocity times do N=1,NTV call arrR4getminmax(nLngLat,VMAP(1,1,N),amin,amax) print *, 'VMAP N, amin and amax = ', N, amin, amax end do print *, 'Before forth filling' do N=1,NTG call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(2,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Remove more small density holes call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) if(amax.ne.BadR4()) then call GridFill(0,nLng,nLat,DMAP(1,1,N),NSIDE,Zmin,Zmax) ! Fill ALL density holes in any map with a hole C AEmax = AEmaxD/(RRS*RRS) C AEmin = AEmin/(RRS*RRS) (needs fixing) C AEmin = AEminD/(RRS*RRS) (fixed 2/29/2008 BVJ) C SpeedD = 0.0 C call FillMapOLS(1,XCbegG,XCintG,nLng,nLat,nTG,N,NTmaxG,RRS,SpeedD,ConstL,AEmax,AEmin,AEangD,DMAP) C call FillMapL(5,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 90 deg discrepancy C call FillMapL(2,XCbegG(1,N),XCbegG(2,N),nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove seam end if end if end do call FillWholeT(0,0,0,XCbeGG,nLng,nLat,nTG,ConsTMD,DMAP,DDD) ! Fill all the density times print *, 'After forth filling' if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1).and.XCtbegG.eq.XCtbegV) then print *, 'Copy array from VMAP to VMAPD directly' call ArrR4Copy(nLngLatnTG,VMAP,VMAPD) else C print *, 'VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,39),VMAPD(50,9,77),VMAPD(50,9,78) call CopyVtoVDN(1,XCbeGG,XCtbegG,XCtendG,XCtbegV,XCtendV,nLng,nLat, & NTV,NTG,ConsTG*aNdayG,VMAP,VMAPD,Dtmpt,VVT,VVV) C print *, 'VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78)',VMAP(50,9,40),VMAPD(50,9,77),VMAPD(50,9,78) end if C end if !ixxxxxxx loop print *, 'Before third write3D_infotd' C do N=1,NTG C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'DMAP N, amin and amax = ', N, amin, amax C end do C do N=1,NTG C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'VMAPD N, amin and amax = ', N, amin, amax C end do C Must change to low value in current version RRSCON = (RRS/R1AU)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change if(bWrite3Do) then cPre = 'nv3o' NinterD = 3 C if(ixxxxxxx.eq.1) then ! ixxxxxxx loop C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM(0,cPre,Nit,NiterT,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) C end if !ixxxxxxx loop end if if(bWrite3D_HRb) then print *, 'Before third write3D_infotd of High Resolution base' cPre = 'nv3b' NintLng = NintLn NintLat = NintLa NintHt = NintH NinterD = 3 nMapb = nMapHRb ! In original map coordinates (so nMapHR = 1.5 AU). if(nMapHR*(NintHt+1)+1.gt.nMap*(NintHt+1)+1) nMapHR = nMap ! Make the setting of nMapHR idiot proof if(nMapb.gt.nMapHR) nMapb = nMapHR C Mode 1 places a r-3 falloff in the 3D matrix following 1 AU. call write3D_infotd3DM_HR(0,cPre,Nit,NiterT,NintLng,NintLat,NintHt,NinterD,NCoff,XCintDG,XCbegG,XCtbegG,XCtendG,iYrBG,RRS,dRR, & nLng,nLat,nMapb,nMapHR,nMap,nTF,nTmaxG,nCar,JDCar, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VMHR,DMHR,XCshift,DVfact,DDfact,V3DHR,D3DHR) end if end if C************************************************************************************************************************************* C Change the density back since the falloff at distance is produced in Extractdn RRSCON = (R1AU/RRS)**FALLOFF do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRSCON,DMAP(1,1,N)) end do C End change if(.not.bForecast) then print *, ' Before 2nd Extractdn',XCbegG(1,1),XCtbegG,XCstrt NinterD = 3 call Extractdn(bExtractf,cPrefixf,NTP,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD,FALLOFF,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) end if if ( bMagnetic ) then !dTT = (XCtendG-XCtbegG)/(nTF-1) ! Time increment (write3D_infotd3D divides by nTF. Why??) (Old) !TT = XCbeGG(1,1)+XCtbegG ! TT is start time (Old) dTT = (XCintG(nTF) - XCintG(1))/(nTF-1) print *, XCintG, dTT TT = XCintG(1) + dTT/2.0 print *, TT call T3D_set (T3D__CLIPLNG , 0, 90.0 ) call T3D_iset(T3D__ITERATIONS , 0, NiterT ) call T3D_iset(T3D__NCOFF , 0, NCoff ) call T3D_set (T3D__RR , 0, RRS ) call T3D_set (T3D__DRR , 0, dRR ) call T3D_set (T3D__TT , 0, TT ) call T3D_set (T3D__DTT , 0, dTT ) 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_set (T3D__PWN_V , 0, PWV ) call T3D_set (T3D__PWN_G , 0, PWG ) call T3D_set (T3D__PWR_V , 0, PWRV ) call T3D_set (T3D__PWR_G , 0, PWRG ) call T3D_set (T3D__D1AU , 0, DEN1AU ) call T3D_set (T3D__SMOOTH_V , 0, CONRV ) call T3D_set (T3D__SMOOTH_D , 0, CONRD ) call T3D_set (T3D__FILL_V , 0, CONSTV ) call T3D_set (T3D__FILL_D , 0, CONSTG ) call T3D_set (T3D__SMOOTH_TIME_V , 0, CONVT ) 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 ) call Write3D_bbt(XCbeGG,LLbeg,(NinterD+1)*nTF,TT3D,XC3DT,V3DT,BR2DT,BR3D,BT3D) end if call Say('IPSD2020','S','Stop','program successfully completed') end