C+ C NAME: C ipsd C PURPOSE: C Tomographic reconstruction based on IPS velocity and g-level C observations. The reconstruction is time-independent, but does C not assume that a mod 360 degree heliosphere exists. C CATEGORY: C Data processing C CALLING SEQUENCE: program ipsd C INPUTS: C IPS data files C OUTPUTS: C ASCII files C CALLS: C ForeignArg, ForeignFile, Say, MAP_TZERO, N_CARRINGTON, AskYN, AskI4, AskR4 C AskR8, AskWhat, Julian, XMAP_SC_POS, SunNewcomb, ECLIPTIC_HELIOGRAPHIC C AdjustJDCar, Local2UT, SetGipsy, DATE_DOY, ReadVIPS, ReadGIPS, FLINT8 C iReadProxyMap, MapGrid, GridFill, WR2DARR, MkLOSWeights, iArrI4Total, ArrR4Bad C ArrR4Zero, ArrI4Zero, ArrR4Constant, ArrR4Copy, ArrR4SetMinMax, ArrR4TimesConstant C ConvertG2D, FillMap, ExtractMap, WriteLOSY, WriteLOSD, MkPos, MkVMap, MkDMap C GridSphere2D, MkShift, Mk_V2D, Mk_D2V, Get3Dval, FixModel, MkVModel, MkGModel C T3D_set_mode, T3D_set, T3D_iset, T3D_marker, Write3D_nv, Write3D_bb, ExtractInsitu C BField_Choose, Str2Str, Flt2Str, ReadVIPSValidData C INCLUDE: include 'dirspec.h' include 'sun.h' include 't3d_array.h' include 't3d_index.h' C EXTERNAL: external iReadNagoya external iProcessNagoya external iReadOoty external iProcessOoty external EARTH C PROCEDURE: C > The input IPS data are organized in either yearly files (Nagoya and Oooty C data) or daily (Cambridge g-level data). Additional input are magnetic C field data at the source surface (radial magnetic field only.) C > The Cambridge data are located using the logical $CAM. This is set C interactively during the first run of the program. The definition is stored C in the file LOGFIL.TXT in the user home directory (Linux) or the root of the C C-drive (Windows). Once define, the only way to force the program to look in C another directory is remove or modify the definition from LOGFIL.TXT. C C > All other data (Nagoya, Ooty, UCSD and magnetic field data) are accessed C by specifying their location as a command line argument: C C IPS data: C nagoya=nagoya,/home/bjackson/dat/nagoya/yearly C ooty=ooty,/home/bjackson/dat/nagoya/yearly C Magnetic data (only one of these should be selected: C wso=wso(4)_(3)_(1).dat,$DAT,map,xmdi C wso_noaa=wso(4)_(3)_(1).dat,$DAT,map,xmdi C C > Compile using g77 as: C C f77 -w -fdollar-ok -ffixed-line-length-none -I$for/h -I$for/h/linux C ipsd.f -L$for -lgen -o$exe/ipsd C C The main features of this compilation are: C - the environment variable $for and $exe need to be defined C - generic include files are expected in $for/h C - platform dependent include files are expected in $for/h/linux (for linux) C or $for/h/win (for Win NT). C - object libraries are expected in $for C - the principal object library is libgen.a (Linux) or forgen.lib (Win NT) C - the executable is created in $exe 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 MODIFICATION HISTORY: C JUN-1992, Paul Hick (UCSD/CASS) C NOV-1995, Bernie Jackson (STEL) C MAR-1999, Bernie Jackson (UCSD) C NOV-2001, Tamsen Dunn (UCSD) C Added magnetic field and prepped for win/linux C MAY-2002, Paul Hick (UCSD/CASS) C Modified tail of program writing the 3D data files C and the timeseries extraction. Main change involves reading C of source surface magnetic maps, and subsequent writing C of 3D magnetic field. C AUG-2002, Paul Hick (UCSD/CASS) C Added check for type of magnetic source maps. C APR-2003, Paul Hick (UCSD/CASS) C If no V (or G) deconvolution is requested and no V (or G) data are available C then the requested deconvolution is now switched off, and deconvolution C continues using only G (or V) data (instead of aborting the program). C SEP-2003, Paul Hick (UCSD/CASS) C Bug fix. For G-deconvolution only, with momentum conservation set. The C source surface velocity map passed to the magnetic field was a constant C velocity map, instead of the map derived from momentum conservation. C- parameter (NLmaxG = 15000) ! Max # G data points parameter (NLmaxV = 4000) ! Max # V data points parameter (NLOS = 40) ! # resolution elements along each line of sight parameter (NTEST = 30) ! # Number of points inside map parameter (NLOSP1 = NLOS+1) parameter (DLOS = 0.05) ! Resolution along los in units of observer-Sun distance parameter (XCinc = 1.00) ! Carrington map "extra" increment integer iXPG (NLmaxG) ! Array indices from daily G mapsC integer iYPG (NLmaxG) integer iMJDG (NLmaxG) ! Original MJD of daily G maps integer IYRFG (NLmaxG) integer IRECG (NLmaxG) integer IYRSG (NLmaxG) ! Year of source observation integer NBSG (NLmaxG) integer LON (NLOSP1,NLmaxG) integer LAT (NLOSP1,NLmaxG) real DOYSG (NLmaxG) ! DOY (and fraction) of source obs. real DISTG (NLmaxG) ! G distance from Sun to Earth real XLSG (NLmaxG) ! Ecliptic longitude of Sun real XLLG (NLmaxG) ! Longitude diff. (Earth-Sun) real XDLG (NLmaxG) ! Declination of Source real XCEG (NLmaxG) ! Carrington variable of Earth real XEG (NLmaxG) ! Source Elongation real XCG (NLmaxG) ! Carrington variable real YLG (NLmaxG) ! Heliographic latitude (degree) real GOBS (NLmaxG) ! Solar wind observed G values !real VOBS (NLmaxV) ! Solar wind observed velocities real VOBS (NLmaxG) ! Solar wind observed velocities real GMDL (NLmaxG) ! Model G values real FIX (NLmaxG) ! Fix to model to make OK real GSSIG (NLmaxG) ! G-level source sigma real WTSG (NLOS ,NLmaxG) real GWTij (NLOS ,NLmaxG) real XLON (NLOSP1,NLmaxG) real XLONP (NLOSP1,NLmaxG) real XLAT (NLOSP1,NLmaxG) real RPX (NLOSP1,NLmaxG) real VFAC (NLOSP1,NLmaxG) ! New array (3/99) real DFAC (NLOSP1,NLmaxG) integer IYRFV (NLmaxV) integer IRECV (NLmaxV) integer IYRSV (NLmaxV) ! Year of source observation integer NBSV (NLmaxV) real DOYSV (NLmaxV) ! DOY (and fraction) of source obs. real DISTV (NLmaxV) ! V distance from Sun to Earth real XLSV (NLmaxV) ! Ecliptic longitude of Sun real XLLV (NLmaxV) ! Longitude diff. (Earth-Sun) real XDLV (NLmaxV) ! Declination of Source real XCEV (NLmaxV) ! Carrington variable of Earth real XEV (NLmaxV) ! Source Elongation real XCV (NLmaxV) ! PLH real YLV (NLmaxV) ! PLH real VVOBS (NLmaxG) real VMDL (NLmaxV) ! Model velocities real VSSIG (NLmaxV) ! Velocity source sigma real VWTij (NLOS,NLmaxV) real WTSV (NLOS,NLmaxV) parameter (nLng = 109) ! # Longitudes, resolution is 360/(nLng-1) parameter (iLng = 37) ! # Longitudes to use in extracted map parameter (nLat = 19) ! # Latitudes , resolution is 180/(nLat-1) parameter (nRad = 31) ! # Radial heights. Range covered is [0,(nRad-1)*dRR] parameter (nLngLat = nLng*nLat) parameter (dRR = 0.1) ! Resolution (AU) in radial direction real VDEV (nLng,nLat) real VPNT (nLng,nLat) real VMAPSAVE(nLng,nLat) real DMAP (nLng,nLat) ! Density map (no holes) real VMAP (nLng,nLat) real VMAPX (iLng,nLat) real DMAPX (iLng,nLat) real DMAPSAVE(nLng,nLat) real GDEV (nLng,nLat) real GPNT (nLng,nLat) real DDfact (nLng,nLat,nRad) real DVfact (nLng,nLat,nRad) ! New array (3/99) real XCshift (nLng,nLat,nRad) real Vtmp (nLng,nLat,2) ! Scratch arrays real Dtmp (nLng,nLat) real XLtmp (nLng) real ANMAPM (nLng,nLat) real AMA (nLng,nLat) real FIXM (nLng,nLat) real WTSM (nLng,nLat) real AMEANF (nLng,nLat) real AFIXM2 (nLng,nLat) real VWTP (NLOS,NLmaxV) real DWTP (NLOS,NLmaxV) real GMWTij (NLOS,NLmaxG) real GMWTi (NLmaxG) !real VV3D (nLng,nLat,nRad) !real DD3D (nLng,nLat,nRad) real BR3D (nLng,nLat,nRad) real BT3D (nLng,nLat,nRad) integer*1 NSIDE (9*nLng*nLat) parameter (nCar = 1000) real*8 JDCar(nCar) real*8 JD real*8 JEpoch real*8 dLngSun real*8 dLatSun real*8 dDisSun real*8 FLINT8 !------- ! 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'/ character cStrPPG(1)*80 /'PPG: Point P G-level map with holes'/ character cStrDV (1)*80 /'DVV: Deconvolved velocity map'/ character cStrDD (1)*80 /'DDD: Deconvolved density map'/ character cStrDVs(1)*80 /'DVS: Deconvolved, smoothed velocity map'/ character cStrDDs(1)*80 /'DDS: Deconvolved, smoothed density map'/ logical bGcon logical bVcon logical bVNago logical bVOoty logical bGCamb logical bGOoty logical bGNago logical bVmv2 logical bGmv2 ! Use mv^2 = constant logical bVproxy /.FALSE./ ! Velocity map used for deconvolution logical bGproxy /.FALSE./ ! G-level map used for deconvolution logical bMagnetic ! Magnetic field extraction ? logical BField_Choose ! Logical function logical bVWRITE /.FALSE./ ! Write out a V data file logical bWrite /.FALSE./ logical bMirror /.FALSE./ logical bForecast logical bExtract(12) /.TRUE.,.FALSE.,.TRUE.,9*.FALSE./ integer iEorWG /0/ integer iEorWV /0/ !integer MJD !integer MJDfrst /48180/ ! MJD for earliest data file !integer MJDlast /49622/ ! MJD for latest data file integer MJDfrst / 46100/ integer MJDlast / 51800/ integer NLOSWV /1/ integer NLOSWG /1/ integer LimitoV /-1/ integer LimitoG /-1/ integer nIterT /18/ integer iOffUT / 8/ integer NAV / 0/ integer NC / 1/ ! # rotations averaged integer nDUMP / 0/ ! Threshold for #pnts/bin integer nFILL / 3/ ! Min. # neighbours for fill real YLbeg /-90./ real YLend / 90./ real XCmargin /.5/ real XElowG /30./ real XElowV /5./ real XEhighG /80./ real XEhighV /180./ real XRlimG /90./ real XRlimV /90./ real GPOWER /0./ real DEN1AU /5.0/ real VMAPSIG /1.0/ real DMAPSIG /1.0/ real GSIG /1.0/ real GSIGB /3.0/ real VSIG /1.0/ real VSIGB /3.0/ real FCONSV /1.0/ real FCONSG /1.0/ real CONRD /6.0/ real CONRV /7.5/ real VLOW /100.0/ real VHIGH /820.0/ real GLOW /0.2/ real GHIGH /4.0/ real ANLIMITV /4./ real ANLIMITG /10./ real RRS / -2.5/ ! Reference distance for deconvolution real RRX / 1.0/ ! Reference distance for Map extraction real VUL /1500.0/ real VLL / 100.0/ real DAV /1./ real Speed /400./ real PWG /0.3/ !/0.5//1.0//2.0/ ! Power of density real PWV /0.3/ !/0.5//1.0//2.0/ ! Power of density real PWR /0.5/ real*8 MJDref real*8 MJDcntr character cSay*4 /'ipsd'/ character cMon*3 character cStrID*11 character cStrI3D*9 character cStrIDX*12 character cArg*100 character cStr*256 character cCambridge*4 character cWildNagoya*256 character cWildOoty *256 parameter (nVar = 3) character cVar(nVar)*100 integer Flt2Str integer Str2Str !------- ! Identify program as corotating tomography call T3D_set_mode(T3D__MODE_TIME, 0) iVar = nVar call ForeignArg(' ',iVar,cVar,cArg) !------- ! Calculate start times for Carrington rotations between 1985 and now. ! Rotation NCoff+I starts at time JDCar(I). iYr = 1985 Doy = 1. call MAP_TZERO(EARTH,iYr,Doy,.01,nCar,JDCar) NCoff = N_CARRINGTON(iYr,Doy) !------ ! iOffUT is an offset time (in hours) between local time and UT. This is used only if the program ! is operated in forecast mode. Note that the local time is assumed to be the same as the ! computer time. The default offset is 8 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 ! Deconvolution distance 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) !------- ! A proxy map for g-values and/or velocity can be used to replace the point-P ! map determined from the observations. ! The map should be in files DMAP.DAT and VMAP.DAT. They must be ascii files of dimension ! nLng, nLat. ! If bVmv2=TRUE then bVcon=bVproxy=FALSE ! If bGmv2=TRUE then bGcon=bGproxy=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$yes',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. if(bGNago) NLOSWG = 1 if(bGOoty) NLOSWG = 1 if(bGCamb) NLOSWG = 2 cCambridge = cEnvi//'cam' 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(cCambridge,iOffUT,iEdt,MJDfrst,MJDlast) if (bVNago .or. bGNago) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya) if (bVOoty .or. bGOoty) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty ) !------- ! Calculate range of Carrington variables to work with: [XClo,XChi] ! XClo is based on the first data file available (MJDfrst). ! 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) !------- ! Fix the reference distance used for the V map formed. ! 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('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('Replace point-P V-values by proxy map',bVproxy) call AskR4('Maximum limit for V deconvolution',VHIGH) end if if (bGcon) then !if(bGNago) PWG = 0.3 !if(bGCamb) PWG = 0.3 !if(bGNago) PWR = 1.4 if(bGCamb) PWG = 0.15 if(bGNago) ANLIMITG = 4 . if(bGNago) CONRD = CONRV 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) ! 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 = 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('Maximum RA relative to Sun for G-data$0$180$',XRlimG) 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',bGproxy) call AskR4('Maximum limit for G deconvolution',GHIGH) call AskR4('Minimum limit for G deconvolution',GLOW) call AskR4('G-level spatial filter',CONRD) CONRV = CONRD call AskR4('Velocity spatial filter',CONRV) call AskR4('G-level power to fit density',PWG) PWV = PWG call AskR4('G-level power to fit velocity',PWV) call AskR4('Radial G-level power',PWR) call AskR4('Solar wind speed for G-level traceback',Speed) end if RRD = RRG if (bWrite) then 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 end if !------- ! 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.4)') 'First MJD =',MJDfrst,' Last MJD =',MJDlast,' Current MJD =',JD MJDref = JD ! MJD at center of plot (Earth position) write (cStr,'(A,I5,A,F10.4,A)') '$',MJDfrst,'$',JD,'$0$' call AskR8('Modified Julian Day (0 = STOP)'//cStr,MJDref) if (MJDref .le. 0) call Say(cSay,'I','StopPlay','StopPlay') !------- ! Make sure iYr,Doy match MJDref ! 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.0,360.0) EarthLat = -dLatSun call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,EarthLng,EarthLat) !------- ! Get modified Carrington variable for center of plot ! [Note that EarthLng=(1-frac(XCEarth))*360] XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) XCbeg = XCEarth-0.5 call AdjustJDCar(XCbeg, nCar,JDCar, NCoff) XCtest1 = XCbeg - XCmargin ! Search limits XCtest2 = XCbeg + NC+XCmargin XCend = XCbeg + NC XCstrt = XCbeg XCbeg = XCbeg - XCinc XCend = XCend + XCinc MJDcntr = MJDref else ! Normal mode write (cStr,'(A,F7.2,A,F7.2,A)') '$',XClo,'$',XChi,'$0$' XCbeg = max(XClo,NCoff+XCbeg) XCbeg = 1941 call AskR4('Carrington rotation (0 = STOP)'//cStr,XCbeg) if (XCbeg .le. 0) call Say(cSay,'I','StopPlay','StopPlay') 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) !------- ! 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,' )' !------- ! Get modified Carrington variable for center of plot ! Get heliographic coordinates for center of plot (Earth position) XCEarth = XCbeg+0.5 ! Center of plot MJDcntr = FLINT8(1,nCar,JDCar,dble(XCEarth),0.0d0) ! 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.0,360.0) 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 - XCinc XCend = XCend + XCinc ! Start MJD for data search MJDref = FLINT8(1,nCar,JDCar,dble(XCtest2),0.d0)-SUN__MJDtoJD end if write (cStrID ,'(I2.2,A ,F8.3)') nint(10*RRS), '_',NCoff+XCbeg+XCinc write (cStrI3D,'( A ,F8.3)') '_',NCoff+XCbeg+XCinc write (cStrIDX,'(I2.2,A2,F8.3)') nint(10*RRX), 'B_',NCoff+XCbeg+XCinc if (bVcon) call AskYN('Save line-of-sight information to file',bVWRITE) call AskYN('Output an Mercury parameter file$no',bExtract( 1)) ! Mercury extraction call AskYN('Output an Earth parameter file$yes' ,bExtract( 3)) ! Earth extraction call AskYN('Output a Ulysses parameter file$no' ,bExtract(10)) ! Ulysses extraction !------- ! Do you want to include a magnetic map at the reference surface height? ! If yes, then read the source surface map at the resolution of the reference arrays. bMagnetic = BField_Choose(iVar,cVar) if (bMagnetic) call AskYN('Do you want to make a 3D magnetic field',bMagnetic) !------- ! The range of the Carrington variable plotted is [XCbeg,XCend]. The target range ! [XCtest1,XCtest2] includes an extra safety margin to make sure that all relevant ! data are read from file. !call ReadVipsLOScheck(XCbeg,XCend,NLOS,NTEST,DLOS) Bad = BadR4() call ReadVIPSValidData(LOS__CHECKV,VLOW,VHIGH,GLOW,GHIGH) ! Reject V <= 0 NLV = 0 NLG = 0 if (bVNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV, & IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bVOoty) call ReadVIPS(iReadOoty, iProcessOoty,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRV,XElowV,XEhighV,iEorWV,XRlimV, & NLmaxV,NLV, & IYRFV,IRECV,IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEVsav,XEVsav,XCVsav,YLVsav,VVsav,GGsav) call ReadVIPSValidData(LOS__CHECKG,VLOW,VHIGH,GLOW,GHIGH) ! Reject G = 0.0 if (bGCamb) call ReadGIPS(cCambridge,nCar,JDCar, & iEdt,.FALSE.,bForecast, & XCtest1,XCtest2,MJDref,MJDfrst,MJDlast, & RRG,Speed,GPOWER,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxG,NLG, & iXPG,iYPG,iMJDG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,VVOBS,GOBS, & 1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav,XDSsav,XLSsav, & XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGOoty) call ReadVIPS(iReadOoty, iProcessOoty,cWildOoty, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG, & IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,VVOBS,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya, & nCar,JDCar,.FALSE.,XCtest1,XCtest2,NCoff,MJDref, & RRG,XElowG,XEhighG,iEorWG,XRlimG, & NLmaxV,NLG, & IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,VVOBS,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) ! Check whether data for deconvolution are available. Switch of deconvolution ! if no data present. if (bVcon .and. NLV .eq. 0) then call Say(cSay,'W','NoV','no velocity IPS data; V deconvolution switched off') bVcon = .FALSE. end if if (bGcon .and. NLG .eq. 0) then call Say(cSay,'W','NoG','no g-level IPS data; G deconvolution switched off') bGcon = .FALSE. end if if (.not. bVcon .and. .not. bGcon) call Say(cSay,'E','sorry','no IPS data') !------ ! The list of T3D_set and T3D_iset calls store parameters in an internal array. ! This information can be pulled out with the T3D_iget and T3D_get subroutines. ! (see Write3D_nv and Write3D_bb) ! Some entries in the internal array only apply to the time-dependent ! tomography. Here they are filled with 'sensible' values. TT = 0.5*(XCbeg+XCend) ! Time assigned to reconstruction dTT = 0.0 call T3D_iset(T3D__NLNG , 0, nLng) ! # longitudes call T3D_iset(T3D__NLAT , 0, nLat) ! # latitudes call T3D_iset(T3D__NRAD , 0, nRad) ! # radial distances call T3D_iset(T3D__NTIM , 0, 1) ! # times (=1) call T3D_iset(T3D__NCOFF , 0, NCoff) ! Carrington offset call T3D_iset(T3D__ITERATION , 0, 0) ! Current iteration call T3D_iset(T3D__ITERATIONS , 0, nIterT) ! Total # iterations call T3D_iset(T3D__ITIME , 0, 1 ) ! Current time index call T3D_set (T3D__TIME , 0, TT ) ! Current time call T3D_set (T3D__TNOW , 0, TT ) ! Forecast time call T3D_set (T3D__XCMAT , 0, XCbeg) call T3D_set (T3D__XCMAT+1 , 0, XCend) !------ ! Define 'region of interest'. This is a range of 1 rotation (360 degrees in ! heliographic latitude extracted from the 3-rotation range [XCbeg,XCend]. ! The region of interest is needed for visualization purposes. call T3D_set (T3D__XCROI , 0, XCbeg+XCinc) call T3D_set (T3D__XCROI+1 , 0, XCend-XCinc) call T3D_set (T3D__LAT , 0, -90.0) ! Range of latitudes covered call T3D_set (T3D__LAT+1 , 0, 90.0) call T3D_set (T3D__TT , 0, TT ) ! First time call T3D_set (T3D__DTT , 0, dTT) ! Time step call T3D_set (T3D__RR , 0, RRS) ! Source surface distancd in AU call T3D_set (T3D__DRR , 0, dRR) ! Radial resolution in AU call T3D_set (T3D__PWN_V , 0, PWV) call T3D_set (T3D__PWN_G , 0, PWG) call T3D_set (T3D__PWR_V , 0, PWR) call T3D_set (T3D__PWR_G , 0, PWR) call T3D_set (T3D__D1AU , 0, DEN1AU) ! Density at 1 AU (cm^-3) call T3D_set (T3D__SMOOTH_V , 0, CONRV) call T3D_set (T3D__SMOOTH_D , 0, CONRD) call T3D_set (T3D__FILL_V , 0, CONRV) call T3D_set (T3D__FILL_D , 0, CONRD) call T3D_set (T3D__SMOOTH_TIME_V , 0, 0.0) call T3D_set (T3D__SMOOTH_TIME_D , 0, 0.0) call T3D_set (T3D__CLIPLNG , 0, 90.0) call T3D_iset(T3D__NL_V , 0, NLV) ! # velocity lines of sight call T3D_iset(T3D__NL_G , 0, NLG) ! # g-level lines of sight call T3D_set (T3D__DLOS_V , 0, dLOS) call T3D_set (T3D__DLOS_G , 0, dLOS) call T3D_iset(T3D__NLOS_V , 0, NLOS) call T3D_iset(T3D__NLOS_G , 0, NLOS) call T3D_set (T3D__BINX_V , 0, ANLIMITV) call T3D_set (T3D__BINX_D , 0, ANLIMITG) if (bVcon) then if (bVproxy) then ! Read proxy point-P V map from file I = iReadProxyMap(1,NCoff+XCbeg,nLng,nLat,VMAP) else ! Calculate proxy point-P V map call MapGrid(XCbeg,XCend,YLbeg,YLend,NC, & NLV,XCV,YLV,VOBS,NAV,DAV,nLng,nLat,VMAP,VDEV,VPNT,Vmin,Vmax,VDmin,VDmax) if (bWrite) call WR2DARR(1,nLng,nLat,VMAP,cStrPPV(1)(:3)//cStrID,'(F7.2)',bMirror,1,cStrPPV) end if ! Calculate weights, WTSG, along all lines of sight call MkLOSWeights(NLOSWV,dLOS,NLOS,NLV,XEV,DISTV,WTSV,PWR) else ! No V deconvolution: use constant V point-P map call ArrR4Constant(nLngLat,Speed,VMAP) end if if(bGcon) then if (bGproxy) then ! Read proxy point-P synoptic G map from file I = iReadProxyMap(3,NCoff+XCbeg,nLng,nLat,DMAP) else ! Calculate point-P synoptic G map call MapGrid(XCbeg,XCend,YLbeg,YLend,NC, & NLG,XCG,YLG,GOBS,NAV,DAV,nLng,nLat,DMAP,GDEV,GPNT,Gmin,Gmax,GDmin,GDmax) if (bWrite) call WR2DARR(1,nLng,nLat,DMAP,cStrPPG(1)(:3)//cStrID,'(F9.3)',bMirror,1,cStrPPG) end if ! Calculate weights, WTSG, along all lines of sight call MkLOSWeights(NLOSWG,dLOS,NLOS,NLG,XEG,DISTG,WTSG,PWR) else ! No G deconvolution: use constant G point-P map call ArrR4Constant(nLngLat,1.,DMAP) end if call ArrR4Copy(nLngLat,VMAP,VMAPSAVE) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,VMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRV,3,0.0,90.0) ! Convert G to density call ConvertG2D(0,nLngLat,DMAP,DMAP,RRG,DEN1AU,PWG) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,DMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAP,CONRD,3,0.0,90.0) call ArrR4Copy(nLngLat,DMAP,DMAPSAVE) ! Density map for V deconvolution call MkShift(XCBeg,XCEnd,nLng,nLat,nRad,VMAP,DMAP,XCshift, & DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,NSIDE,Vtmp,Dtmp,XLtmp) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRV,1,0.0,90.0) call ArrR4Zero(NLmaxV,VSSIG) call ArrR4Zero(NLmaxG,GSSIG) call ArrI4Constant(NLV,1,NBSV) call ArrI4Zero(NLmaxV-NLV,NBSV(NLV+1)) ! Probably unnecessary since these are never accessed call ArrI4Constant(NLG,1,NBSG) call ArrI4Zero(NLmaxG-NLG,NBSG(NLG+1)) ! Probably unnecessary since these are never accessed Nit = -1 ! Why three iterations counters?? NitV = -1 NitG = -1 do while (Nit .lt. nIterT) Nit = Nit+1 call Say(cSay,'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=bGproxy=.FALSE. if (bGmv2 .and. (NitV .eq. 0 .or. NitV .eq. LimitoV+2)) then call Mk_V2D(nLngLat,VMAP,DMAPSAVE,DEN1AU,RRV) CONRD = CONRV RRD = RRV call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAPSAVE,CONRD,1,0.0,90.0) end if call MkPos(XCbeg,XCend,RRS,dRR,DLOS,nLng,nLat,nRad,XCshift, & NLV,NLOSP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV, & XLON,XLONP,XLAT,RPX,LON,LAT) call Get3Dval(XCbeg,XCend,RRV,dRR,nLng,nLat,nRad,DDfact,NLOSP1*NLV,XLONP,XLAT,RPX,DFAC) call Get3Dval(XCbeg,XCend,RRV,dRR,nLng,nLat,nRad,DVfact,NLOSP1*NLV,XLONP,XLAT,RPX,VFAC) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,DMAPSAVE) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAPSAVE,CONRD,3,0.0,90.0) call MkVModel(XCbeg,XCend,XEV,XLON,XLAT,RPX,WTSV,DFAC,VFAC, & PWV,NLV,NLOS,NLOSP1,VMAP,VHIGH,DMAPSAVE,nLng,nLat, & RRD,RRV,RRS,VMDL,VWTij) call FixModel(1,NLV,VOBS,VMDL,PWV,NBSV,VSIG,FIX,VSSIG,VRAT) call MkVMap(LON,LAT,RPX,XEV,VWTij,FIX,NLV,NLOS,NLOSP1,DMAPSAVE,RRD, & XLON,XLAT,DFAC,nLng,nLat,PWV,VMAP,ANLIMITV,VHIGH,VSSIG,VMDL,VSIGB,NVS, & NBSV,VMAPSIG,ANMAPM,AMA,FIXM,WTSM,AMEANF,AFIXM2,VWTP,CONRV) call ArrR4Copy(nLngLat,VMAP,VMAPSAVE) call ArrR4SetMinMax(-nLngLat,VMAP,VLL,VUL) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,VMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRV,3,0.0,90.0) ! Note that we suddenly switch back to DMAP instead of DMAPSAVE. Why????????? call MkShift(XCbeg,XCend,nLng,nLat,nRad,VMAP,DMAP,XCshift, & DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,NSIDE,Vtmp,Dtmp,XLtmp) end if if (bGcon) then ! G deconvolution NitG = NitG+1 ! # G iterations ! bVmv2=.TRUE., implies bVcon=bVproxy=.FALSE. if (bVmv2 .and. (NitG .eq. 0 .or. NitG .eq. LimitoG+2)) then call Mk_D2V(nLngLat,DMAP,VMAP,DEN1AU,RRD) RRV = RRD call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRD,1,0.0,90.0) call ArrR4Copy(nLngLat,VMAP,VMAPSAVE)! PPH: bugfix SEP-2003 end if call MkPos(XCbeg,XCend,RRS,dRR,DLOS,nLng,nLat,nRad,XCshift, & NLG,NLOSP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG, & XLON,XLONP,XLAT,RPX,LON,LAT) !------- ! ??? Seems to me this should be done before entering the iteration loop ! Put in the following statements if pseudo observed velocity values are to be ! deconvolved. Also make sure that NLmaxV = NLmaxG in the parameter list ! of the main program and other subroutines that use this parameter. ! !if (NitG .eq. 0) then ! call MkVObs(XCbeg,XCend,VMAP,nLng,nLat,XLON,XLAT,NLG,NLV,NLOSP1, !& XLSG,XCEG,XLLG,XDLG,XEG,IYRSG,DOYSG,VOBS,XLSV,XCEV,XLLV,XDLV,XEV,IYRSV,DOYSV) ! NitV = -1 ! bVcon = .TRUE. !end if call Get3Dval(XCbeg,XCend,RRG,dRR,nLng,nLat,nRad,DDfact,NLOSP1*NLG,XLONP,XLAT,RPX,DFAC) call MkGModel(XCbeg,XCend,XLON,XLAT,RPX,WTSG,DFAC,PWG,DEN1AU,NLG, & NLOS,NLOSP1,DMAP,nLng,nLat,RRD,GMDL,GWTij,GMWTij,GMWTi) call FixModel(2,NLG,GOBS,GMDL,PWG,NBSG,GSIG,FIX,GSSIG,GRATIO) call MkDMap(LON,LAT,GWTij,PWG,FIX,NLG,NLOS,NLOSP1, & nLng,nLat,DMAP,ANLIMITG,GSSIG,GSIGB,NGS, & NBSG,DMAPSIG,ANMAPM,AMA,FIXM,WTSM,AMEANF,AFIXM2,DWTP,CONRD) call ArrR4Copy(nLngLat,DMAP,DMAPSAVE) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,DMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAP,CONRD,3,0.0,90.0) call MkShift(XCbeg,XCend,nLng,nLat,nRad,VMAP,DMAP,XCshift, & DVfact,DDfact,RRS,dRR,CONRV,CONRD,VLL,VUL,NSIDE,Vtmp,Dtmp,XLtmp) end if end do ! End of iteration loop do I=1,NLV ! List bad IPS sources if (NBSV(I) .eq. 0) write (*,'(A,I5,2(A,F7.1))') ' Bad V los, # =',I,', V =',VOBS(I),', elong =',XEV(I) end do write (*,'(A,2(I4,A))') ' There were',NLV-iArrI4Total(NLV,NBSV,I),'/',NLV,' bad V sources removed.' do I=1,NLG if (NBSG(I).eq. 0) write (*,'(A,I5,2(A,F7.3))') ' Bad G los, # =',I,', G = ',GOBS(I),', elong =',XEG(I) end do write (*,'(A,2(I4,A))') ' There were',NLG-iArrI4Total(NLG,NBSG,I),'/',NLG,' bad G sources removed.' ! Here VMAPSAVE contains the original point-P map (bVcon=.FALSE.) ( or is equal to VMAP ! except for a pass through ArrR4MinMax ????) call ArrR4Copy(nLngLat,VMAPSAVE,VMAP) call ArrR4Copy(nLngLat,DMAPSAVE,DMAP) write (*,'(A,F7.2,A)') ' Map begins at',(1-XCbeg+int(XCbeg))*360,' degrees.' ! Extract from unsmoothed maps call GridFill(2,nLng,nLat,VMAP,NSIDE,Zmin,Zmax) ! Remove small velocity holes call GridFill(2,nLng,nLat,DMAP,NSIDE,Zmin,Zmax) ! Remove small density holes if (bWrite) then ! Write data to file call ArrR4TimesConstant(-nLngLat,DMAP,RRS*RRS,DMAP) ! n*R^2 unsmoothed call ExtractMap(RRX,RRS,dRR,XCbeg,XCend,nLng,nLat,nRad,VMAP,DMAP, & XCshift,DVfact,DDfact,XCstrt,iLng,VMAPX,DMAPX) if (bVcon) then ! V unsmoothed call WR2DARR(1,nLng,nLat,VMAP,cStrDV(1)(:3)//cStrID,'(F7.2)',bMirror,1,cStrDV) call WR2DARR(1,iLng,nLat,VMAPX,cStrDV(1)(:3)//cStrIDX,'(F7.2)',bMirror,1,cStrDV) call GridSphere2D(XCend-XCbeg,iLng,nLat,1,VMAPX,CONRV,0,0.0,0.0) ! V smoothed call WR2DARR(1,iLng,nLat,VMAPX,cStrDVs(1)(:3)//cStrIDX,'(F7.2)',bMirror,1,cStrDV) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRV,0,0.0,0.0) ! V smoothed call WR2DARR(1,nLng,nLat,VMAP,cStrDVs(1)(:3)//cStrID,'(F7.2)',bMirror,1,cStrDVs) end if if (bGcon) then call WR2DARR(1,nLng,nLat,DMAP,cStrDD(1)(:3)//cStrID,'(F10.4)',bMirror,1,cStrDD) call WR2DARR(1,iLng,nLat,DMAPX,cStrDD(1)(:3)//cStrIDX,'(F10.4)',bMirror,1,cStrDD) call GridSphere2D(XCend-XCbeg,iLng,nLat,1,DMAPX,CONRD,0,0.0,0.0) ! n*R^2 smoothed call WR2DARR(1,iLng,nLat,DMAPX,cStrDDs(1)(:3)//cStrIDX,'(F10.4)',bMirror,1,cStrDD) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAP,CONRD,0,0.0,0.0) ! n*R^2 smoothed call WR2DARR(1,nLng,nLat,DMAP,cStrDDs(1)(:3)//cStrID,'(F10.4)',bMirror,1,cStrDDs) end if end if ! End write data to file if (bVWRITE) then if (bVcon) then I = Str2Str('losv',cStr) I = Flt2Str(NCoff+TT,-4,cStr(I+1:)) if (dTT .ne. 0) I = T3D_marker(cStr,icount) if (bVNago) call WriteLOSY(cStr,cWildNagoya,NLV,IYRFV,IRECV,VOBS,VMDL,VMDL,XEV,NBSV) if (bVOoty) call WriteLOSY(cStr,cWildOoty ,NLV,IYRFV,IRECV,VOBS,VMDL,VMDL,XEV,NBSV) end if if (bGcon) then I = Str2Str('losg',cStr) I = Flt2Str(NCoff+TT,-4,cStr(I+1:)) if (dTT .ne. 0) I = T3D_marker(cStr, icount) if (bGNago) call WriteLOSY(cStr,cWildNagoya,NLG,IYRFG,IRECG,GOBS,GMDL,GMDL,XEG,NBSG) if (bGOoty) call WriteLOSY(cStr,cWildOoty ,NLG,IYRFG,IRECG,GOBS,GMDL,GMDL,XEG,NBSG) if (bGCamb) call WriteLOSD(cStr,cCam,iEdt,NLG,iMJDG,iXPG,iYPG,GOBS,GMDL,GMDL,XEG,NBSG) end if end if call ArrR4Copy(nLngLat,VMAPSAVE,VMAP) call ArrR4Copy(nLngLat,DMAPSAVE,DMAP) call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,VMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,VMAP,CONRV,3,0.0,90.0) call ArrR4TimesConstant(-nLngLat,DMAP,RRS*RRS,DMAP) ! n*R^2 unsmoothed call FillMap(0,XCbeg,XCend,nLng,nLat,ConstL,DMAP) call GridSphere2D(XCend-XCbeg,nLng,nLat,1,DMAP,CONRD,3,0.0,90.0) call GridFill(2,nLng,nLat,VMAPSAVE,NSIDE,Zmin,Zmax) ! Remove small velocity holes call GridFill(2,nLng,nLat,DMAPSAVE,NSIDE,Zmin,Zmax) ! Remove small density holes call ArrR4TimesConstant(-nLngLat,DMAPSAVE,RRS*RRS,DMAPSAVE)! n*R^2 unsmoothed !------- ! For the final solution the iteration should be set larger than ! the total # iterations (the only result is that the iteration ! number is NOT encoded in the output file name). call T3D_iset(T3D__ITERATION , 0, nIterT+1) ! Current iteration !------- ! DMAPSAVE is a normalized to density (to 1AU). DDfact are normalized ! density ratios. DVfact are velocity ratios. Write3Dinfo_nv returns ! normalized density in DDfact and velocity in DVfact (i.e. they are ! NOT RATIOS ANYMORE). call Write3D_nv(VMAPSAVE,DMAPSAVE,XCshift,DVfact,DDfact) if (bMagnetic) then ! Get magnetic field data !------- ! Write3D_nv will try to find source surface maps to read and ! calculate the 3D magnetic field. The magnetic field is written to ! file, and is returned in BR3D and BT3D for use by ExtractInsitu. call Write3D_bb(XCshift,DVfact,BR3D,BT3D) else ! Set magnetic field data 'bad' call ArrR4Bad(nLng*nLat*nRad,BR3D) call ArrR4Bad(nLng*nLat*nRad,BT3D) end if !------- ! Extract the insitu time series at objects set in bExtract call ExtractInsitu(bExtract,DVfact,DDfact,BR3D,BT3D,nCar,JDCar) call Say(cSay,'S','Stop','program successfully completed') end