C+ C NAME: C ipsdt 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: program ipsdt C INPUTS: C Data file $CAM:G*. C OUTPUTS: C Rectangular grid of data points C CALLS: 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 MKshiftd. 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- parameter (NLmaxG = 15000) ! Max # G data points parameter (NLmaxV = 15000) ! Max # V data points parameter (NGREST = 1) ! G temporal resolution factor better than 1 day parameter (NVREST = 1) ! V temporal resolution factor better than 1 day parameter (NGRESS = 1) ! Spatial G resolution factor better than 20 degrees parameter (NVRESS = 1) ! Spatial V resolution factor better than 20 degrees parameter (NF = 2) ! LOS resolution factor (Should be 1 if Spatial resolution is >= 2) parameter (nTmaxG = 50*NGREST) ! Max # G times parameter (nTmaxV = 50*NVREST) ! Max # V times 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) ! Resolution (AU) in radial direction (0.1) parameter (dRR = 0.1/(NGREST*NF/2)) 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 DOYSG (NLmaxG), ! DOY (and fraction) of source obs. & DOYSGG (NLmaxG), ! DOY (and fraction) of source obs. & 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 & SRCVS*9,SRCVsav(NLmaxV)*9, & SRCG (NLmaxG)*9, ! G-level source name & SRCGS*9,SRCGsav(NLmaxG)*9 integer IYRFV (NLmaxV), & IRECV (NLmaxV), & NSRV (NLmaxV), & NSRG (NLmaxG), & IYRSV (NLmaxV), ! Year of source observation & NBSV (NLmaxV) real DOYSV (NLmaxV), ! DOY (and fraction) of source obs. & 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), & 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), & 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), & Vtp (nLng,nLat), & Dtp (nLng,nLat), & ANMAPD (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), & 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 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, & 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./, & bWriteMS /.FALSE./, & bWritelos /.FALSE./, & bOneOr1 /.TRUE./, & bOneOr2 /.TRUE./, & bForecast, & bExtract(12) /2*.FALSE.,.TRUE.,9*.FALSE./, & bEWRITE, ! Write out an Earth file & bMWRITE, ! Write out a Mercury file & bUWRITE ! Write out a Ulysses file logical BField_Choose ! Logical function equivalence (bExtract( 3), bEWRITE), & (bExtract(1), bMWRITE), & (bExtract(10), bUWRITE) 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 /30./, & XElowV /11.5/, & XEhighG /80./, & XEhighV /180./, & XRlimG /90./, & XRlimV /90./, & GPOWER /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./, & CONVT /0.65/, & CONDT /0.65/, & ConsTV /1.0/, & ConsTG /1.0/, & VLIM /1200./, ! Velocity upper limit for sources & VLIL /100./, ! Velocity lower limit for sources & 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 iReadNagoyan, iProcessNagoyan external iReadOotyn, iProcessOotyn 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: Nagoya, Ooty?$1$1',I) bVNago = I .eq. 1 bVOoty = I .eq. 2 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 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 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 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 limit for V deconvolution?',VLIM) end if if (bGcon) then if(bGCamb) PWG = 0.15 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 call AskR4('Minimum elongation for G-data$0$180$',XElowG) if(bGNago) XEhighG = 180. 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) call AskR4(' G-level spatial filter?',CONRD) if(bVNago.and.bGNago) CONRV = CONRD call AskR4('Velocity spatial filter?',CONRV) call AskR4(' G-level temporal filter?',CONDT) if(bVNago.and.bGNago) CONVT = CONDT call AskR4('Velocity temporal filter?',CONVT) call AskR4(' G-level power to fit density?',PWG) if(bVNago.and.bGNago) 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 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) NOneOr1 = -2 NOneOr2 = -2 NinterC = 3 ! Default first extract matrix NinterD = 3 ! Default forecast extract matrix call AskYN('Output Three D 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 Three D filled files?$yes',bWrite3D2) ! 3D 2 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 an Earth parameter file$yes',bEWRITE) ! Earth extraction call AskYN('Output a Mercury parameter file$no',bMWRITE) ! Mercury extraction call AskYN('Output a Ulysses parameter file$no',bUWRITE) ! Ulysses extraction 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 ReadVIPSn(iReadNagoyan,iProcessNagoyan,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,SRCV,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVOoty) call ReadVIPSn(iReadOotyn, iProcessOotyn,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV,SRCV,SRCVsav, & IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav) if (bGCamb) call ReadGIPS('$CAM',nCar,JDCar, & iEdt,.FALSE.,bForecast, & XCtest1,XCtest2,MJDref,MJDfrst,MJDlast, & RRG,Speed,GPOWER,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG, C & NLmaxG,NLG,SRCG,SRCGsav, & iXPG,iYPG,iMJDG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VVOBS,GGOBS, & 1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav,XDSsav,XLSsav, & XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGOoty) call ReadVIPSn(iReadOotyn, iProcessOotyn,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, c & RRV,XElowG,XEhighG,iEorWG,XRlimG, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,SRCG,SRCGsav, c & NLmaxV,NLV,SRCG,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, c & IYRFG,IRECG,IYRSGG,DOYSGG,DISTV,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGNago) call ReadVIPSn(iReadNagoyan,iProcessNagoyan,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG,SRCG,SRCGsav, & IYRFG,IRECG,IYRSGG,DOYSGG,DISTGG,XLSGG,XLLGG,XDLGG,XCEGG,XEGG,XCGG,YLGG,VOBS,GGOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,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), SRCG(I) if (VOBS(I) .gt. VLIM) then NBSV(I) = 0 VOBS(I) = BadR4() end if if (VOBS(I) .lt. VLIL) 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) DOYSG(NLG-I+1) = DOYSGG(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 SRCG(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) DOYSG(I) = DOYSGG(I) IYRSG(I) = IYRSGG(I) end if end do ALng = 1.*(nLng-1)/(iLng-1) if(bForecast) then ! For forecast mode NTV = (NTV + 1.5)/2.0 + 6 NTG = (NTG + 1.5)/2.0 + 6 end if 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, & 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, & NTG,IYRBG,IdaybegG,XCintDG,XCintG,XCbeGG,XCtbegG,XCtendG,IYRG,DOYG) 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,90.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,90.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,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,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,90.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,90.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 MkShiftdn(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,90.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,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,Idaybeg,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) 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_Scomp(1,NCoff,XCintV,NitV,cStrVname,nTmaxV,NTV,NLV,PWV,VRATIO, & NLmaxV,NSRV,SRCV,DOYSV,XCEV,VOBS,VMOD,XLON,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 call MkVMaptdn(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,VMAV,FIXMV,WTSMV,VMEANFV,VFIXM2V,WTPV, & Alng,ConstL,CONRV,CONVT,aNdayV,VZtmp) 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 MkShiftd DMAP min max N',amin,amax,N AEmax = AEmaxD/(1.d0*RRS*RRS) AEmin = AEminD/(1.d0*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,90.0) C call arrR4getminmax(nLngLat,DMAP(1,1,N),amin,amax) C print *, 'Before MkShiftd 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,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 MkShiftd 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,90.0) C call arrR4getminmax(nLngLat,VMAPD(1,1,N),amin,amax) C print *, 'Before Timesmooth and MkShiftd 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 MkShiftd VMAPD min max N',amin,amax,N end do end if end if if(bGcon) then call MkShiftdn(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,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,Idaybeg,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG, & XLON,XLONP,XLproj,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) if(NitG.eq.NiterT.and.bWritelos) call writelosmap(NLOSGP1,NLG,nTG,NBSG,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_Scomp(2,NCoff,XCintG,NitG,cStrGname,nTmaxG,NTG,NLG,PWG,GRATIO, & NLmaxG,NSRG,SRCG,DOYSG,XCEG,GOBS,GMOD2,XLON,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 MkDMaptdn(LON,LAT,ITIM,GWTij,PWG,FIX,NLG,NLOSG,NLOSGP1, & nLng,nLat,nTG,DMAP,ANLIMITG,GSSIG,GSIGB,NGS,NBSG,DMAPSIG, & ANMAPD,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 MkDmaptd 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 MkDmaptdn 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,90.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,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 MkDmaptdn 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,90.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 MkShiftdn(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 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 (*,'(A,I4,A,F7.1,A,F7.1)') ' Bad V source removed, # =', & I,', V =',VOBS(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 (*,'(A,I5,A,F7.3,A,F7.3)') ' Bad G source removed, # =', & I,', G = ',GOBS(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,VVT,VVV) C end if call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,DMAPV,VMAP,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,DDT,DDD) C end if call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,VMAPD,DMAP,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/(1.d0*RRS*RRS) AEmin = AEminD/(1.d0*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,VVT,VVV) end if C call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegG,XCtendG,nLng,nLat, C & NTG,NTV,ConsTG,DMAP,DMAPV,DDT,DDD) print *, ' Before Extractd',XCbegG(1,1),XCtbegG,XCstrt C***** Take data out at different spacecraft (and earth) locations (bExtract tells where) call Extractd(bExtract,RRS,dRR,XCbeGG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,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,(NinterC+1)*nTG ! 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(.not. bForecast) then C 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 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/(1.d0*RRS*RRS) AEmin = AEminD/(1.d0*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 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,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) if(bForecast) then C print *, 'DATA EXTRACTED FROM FINAL MATRIX FOR FORECAST' call Extractd3D(bExtract,RRS,dRR,XCbe,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterD+1,V3DT,D3DT,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_bbtt(XCbeGG,LLbeg,(NinterD+1)*nTF,TT3D,XC3DT,V3DT,BR2DT,BR3D,BT3D) end if call Say('IPSD2020','S','Stop','program successfully completed') end