C+ C NAME: C ipsg2s C PURPOSE: C Deconvolve IPS observations; (quasi-) time-independent. C CATEGORY: program ipsg2s C CALLING SEQUENCE: C run $exe:ipsg2s C INPUTS: C IPS data files C OUTPUTS: C t3d* files in current directory C CALLS: C Say, itrim, LocFirst, ForeignArg, ForeignFile, ReadVIPSValidData, ReadVIPS C ReadGIPS, AskYN, AskR4 C SetGrid, IPSProcessing, IPSConstraints, T3D_set, T3D_iset, T3D_iget, T3D_show, C T3D_get_grid, LOSPosition, LOSProjection, LOSReach, LOSClean, LOSWeights C LOSIntegralG, LOSIntegralV, LOSTweak, Str2Str, Flt2Str, ArrR4TimesArrR4 C ArrR4Constant, ArrR4Copy, ArrR4Function, ArrR4Mask, MapGrid, SanityCheck C CvG2D, CvD2G, MkD2V, MkV2D, iArrI4Total C SW_Model_Kinematic, BuildSourceSurface, GridSphere2D, GridFill, Write3D_nv, Write3D_bb C Write3D_nv_UT, Write3D_bb_UT, WriteLOSY C WriteLOSD, StopWatch, T3D_marker C iReadNagoya, iProcessNagoya, iReadOoty, iProcessOoty C INCLUDE: include 'dirspec.h' include 't3d_param.h' include 't3d_array.h' include 't3d_index.h' C EXTERNAL: external iReadNagoya external iProcessNagoya ! Passed to ReadVIPS external iReadOoty external iProcessOoty C PROCEDURE: C > !!! Normalized densities are used throughout. C These are defined as n*r^2 with n the density in cm^-3 and r the heliocentric C distance in AU. Typical values would be in the range 1-10. C > !!! Velocities are in km/s 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 IPSD), the C only way to force IPSD to look in another directory is to delete it. C C II = iand(I__MODE,TOM__MOD): C C LOSReach: C uses I__MODE to stop execution when LOS segments fall outside C the range [XCbegMAT,XCendMAT] (only if II = 0). C C SW_Model_Kinematic: C uses I__MODE three times to deal with mod(XCendMAT-XCbegMAT) issues (if II != 0) C and twice to avoid array index errors (if II = 0). C The 2nd of these latter is the only place in the NOMOD scheme where a C rather ugly kludge is needed. C C LOSIntegralG, LOSIntegralV C passes I__MODE to Get3DValue C C LOSProjection C passes I__MODE to Get3DValue C C BuildSourceSurface C uses I__MODE to map XC values into the range [XCbegMAT,XCendMAT] (II != 0 only). C C POTENTIAL MODIFICATIONS/IMPROVEMENTS: C C The first and last grid point in longitude and latitude now are located C on the edges of the interval covered. It makes more sense to put the grid C at bin centers instead. E.g. at a 10 degree resolution the latitude grid C now is located at -90,-80,...,80,90. It should be located at -85,-75,75,85. C This is probably fairly easy to fix in the Fortran tomography program. C It may be more difficult to adapt the IDL visualization programs. C C Switch all Carrington variables (time and longitude) to double precision C MODIFICATION HISTORY: C JUN-1992, Paul Hick (UCSD) C point-P synoptic map program C NOV-1995, Bernie Jackson (STEL) C conversion to tomography program C APR-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C subtantial rewrite. Program now uses g^2 rather than density. C- !------- ! Line of sight resolution parameter (NLG0 = 12000) ! Max # G data points parameter (NLV0 = 3000) ! Max # V data points parameter (NLOS0 = 40) ! # resolution elements along each line of sight parameter (dLOS0 = 0.05) ! Resolution along los in AU parameter (NLmax0 = NLG0) ! Use the larger of NLG0 and NLV0 to define NLmax0 parameter (NLTot0 = NLOS0*NLmax0) ! Max # los elements !------- ! We set up the grid size separately here. They are used to declare ! array sizes, and are also the default grid sizes used. The grid size ! can be modified by specifying the appropriate command line argument. ! As long as the total # of grid points is not bigger than nSST0 (at ! the source surface) and n3DT0 (in the 3D volume) everything should ! work. parameter (nLng0 = 37) ! # Longitudes, resolution is 360/(nLng-1) parameter (nLat0 = 19) ! # Latitudes , resolution is 180/(nLat-1) parameter (nRad0 = 31) ! # Radial heights. Range covered is [0,(nRad-1)*dRR] parameter (nTim0 = 1) ! Max # long/lat elements in source surface map parameter (nSS0 = (3*(nLng0-1)+1)*nLat0) parameter (n3D0 = nSS0 *nRad0) ! Max # long/lat/rad elements parameter (nSST0 = nSS0 *nTim0) ! Max # long/lat/tim elements in source surface map parameter (n3DT0 = nSST0*nRad0) ! Max # long/lat/rad/tim elements in whole volume real XCbegMAT(nTim0) real XCendMAT(nTim0) real XCbegROI(nTim0) real XCendROI(nTim0) !------- ! Arrays for IPS velocity observations (max. of NLmaxV los observations) ! The first two array are needed to retrieve the original data from the ! input file and write them to a separate file with the model IPS data ! attached. integer IYRFV (NLV0) ! Year of file where los observation is located integer IRECV (NLV0) ! Record number in file integer IYRSV (NLV0) ! Year of los observation integer NBSV (NLV0) ! Initialized by LOSReach ! IPS V data arrays real DOYSV (NLV0) ! Day of year (incl. fraction for time of day) of observation real DISTV (NLV0) ! Sun-Earth distance at time of obs. real XLSV (NLV0) ! Geocentric ecliptic longitude of Sun real XLLV (NLV0) ! Geocentric ecliptic longitude diff.lng(los)-lng(Sun) real XDLV (NLV0) ! Geocentric ecliptic latitude los real XCEV (NLV0) ! Carrington variable of sub-Earth point real XEV (NLV0) ! Los elongation (deg) (>0: East of Sun; <0: West of Sun) real XCV (NLV0) ! Carrington variable of point-P traced back to source surface at speed VOBS real YLV (NLV0) ! Heliographic latitude of point-P (deg) real VOBS (NLV0) ! Observed IPS velocities real VMDL (NLV0) ! Model IPS velocities at end real VMDL0 (NLV0) ! Model IPS velocities after iteration 0 real PPlosV (LOS__N,NLOS0,NLV0) real VPlosV (NLOS0,NLV0) ! sin(Chi) = Vperp/Vrad real WWlosV (NLOS0,NLV0) !------- ! Arrays for IPS G-level observations (max. of NLmaxG los observations) ! The first two array are needed to retrieve the original data from the ! input file and write them to a separate file with the model IPS data ! attached. integer iXPG (NLG0) ! Array indices from daily G maps integer iYPG (NLG0) integer iMJDG (NLG0) ! Original MJD of daily G maps integer IYRFG (NLG0) ! Year of yearly IPS file integer IRECG (NLG0) ! Record # in yearly file equivalence (iMJDG, IYRFG) equivalence (iXPG , IRECG) ! IPS G data arrays integer IYRSG (NLG0) ! Year of source observation integer NBSG (NLG0) ! Initialized by LOSReach real DOYSG (NLG0) ! Day of year (incl. fraction for day of year) of obs. real DISTG (NLG0) ! Sun-Earth distance at time of obs. real XLSG (NLG0) ! Geocentric ecliptic longitude of Sun real XLLG (NLG0) ! Geocentric ecliptic longitude diff.lng(los)-lng(Sun) real XDLG (NLG0) ! Geocentric ecliptic latitude los real XCEG (NLG0) ! Carrington variable of sub-Earth point real XEG (NLG0) ! Los elongation (deg) (>0: East of Sun; <0: West of Sun) real XCG (NLG0) ! Carrington variable of point-P traced back to source surface at speed V0 real YLG (NLG0) ! Heliographic latitude of point-P (deg) real GOBS (NLG0) ! Observed IPS g-levels real G2OBS (NLG0) ! Observed IPS g^2-levels real G2MDL (NLG0) ! Model IPS G-levels at end real G2MDL0 (NLG0) ! Model IPS G-levels after iteration 0 real PPlosG (LOS__N,NLOS0,NLG0) real WWlosG (NLOS0,NLG0) real M2Mean (NLG0) !------- ! Scratch space for line of sight manipulations. ! Note that some of the arrays are equivalenced, so be careful!! integer iLOSTmp (NLTot0,3) real rLOSTmp (NLTot0,3) equivalence (rLOSTmp, iLOSTmp) real GMDL (NLmax0) real VSUM (NLmax0) equivalence (GMDL,VSUM) ! Arrays shared by IPS V and G real PPprj (PRJ__N,NLTot0) ! Lng, lat, time of projected point on los real Weight (NLTot0) real FIX (NLmax0) ! Fix to model to make OK !------- ! Source surface arrays (includes scratch arrays) real VMap (nSST0) ! Source surface velocity map real VFin (nSST0) real GMap (nSST0) ! Source surface density map real GFin (nSST0) real BINXV (nSST0) ! # los crossings in V (scratch array) real BINXG (nSST0) ! # los crossings in G (scratch array) real TmpSST (nSST0*3) ! Scratch array for SW_Model_Kinematic real TmpSST1 (nSST0) real TmpSST2 (nSST0) equivalence (TmpSST1, TmpSST( 1)) equivalence (TmpSST2, TmpSST(nSST0+1)) real Gtmp (nSS0,2) ! Scratch arrays for SW_Model_Kinematic, BuildSourceSurface real Vtmp (nSS0,2) byte Btmp (nSS0,9) ! Scratch array equivalence (Btmp,Vtmp) !------- ! Volume array (includes scratch arrays) real XC3D (n3DT0*3) ! 3D 'traceback' array real V3D (n3DT0) ! Velocities real G3D (n3DT0) ! g^2 real Tmp3D1 (n3D0 ) ! Scratch array real Tmp3D2 (n3D0 ) ! Scratch array logical bGcon logical bVcon logical bVNago logical bVOoty logical bGNago logical bGOoty logical bGCamb logical bMV2 ! Use mv^2 = constant logical bBoundV /.TRUE./ ! Limit V-map boundaries? logical bBoundG /.TRUE./ ! Limit D-map boundaries? logical bForecast logical bRemove logical bPointP logical bMod logical bDiff logical bTrackG integer iEorWG / 0/ integer iEorWV / 0/ integer NAV / 0/ integer NitTotal /18/ integer NitBadV /-1/ integer NitBadG /-1/ integer nFill / 4/ ! Min. # neighbours for fill integer NT / 0/ real BINXVmin / 4.0/ ! Minimum required # los crossings at source surface real BINXGmin / 10.0/ ! (passed to BuildSourceSurface) real YLbeg /-90.0/ ! Heliographic latitude limits (MapGrid) real YLend / 90.0/ real XElowG / 30.0/ ! Limits on los elongations (ReadVIPS,ReadGIPS) real XElowV / 5.0/ real XEhighG / 80.0/ real XEhighV /180.0/ real XRlimG / 90.0/ ! Limits on hour angle difference with Sun (ReadVIPS,ReadGIPS) real XRlimV / 90.0/ real GPOWER / 0.0/ ! Determines radial dependence of G (ReadGIPS) real ClipLng / 0.0/ ! Needs to be a non-zero for NOMOD version (90?) real V0 /400.0/ ! Passed to ReadGIPS; constant V map if no IPS V data real PWN / 0.3/ ! Power for normalized density dependence real PWR / 1.4/ ! Power for radial dependence (deviation from r^-2) real SMOOTH_D / 6.0/ ! Gaussian width (deg) for angular smoothing (GridSphere2D) real SMOOTH_V / 7.5/ real SMOOTH_TIME / 0.0/ real D1AU / 1.0/ ! Density at 1 AU real GSIGB / 3.0/ real VSIGB / 3.0/ real DAV / 1.0/ ! Passed to MapGrid parameter (nCar = 1000) real*8 JDCar(nCar) real*8 MJDref character cCam*4 character cStr*80 character cSay*6 /'ipsg2s'/ character cVar(2)*100 character cArg*100 character cWildNagoya*50 character cWildOoty *50 intrinsic sqrt integer Flt2Str integer Str2Str real Scale(4) /0.0,1.0,0.0,1.0/ iVar = 2 call ForeignArg(' ',iVar,cVar,cArg) Bad = BadR4() I__MODE = TOM__NORATIO !------- ! Check command line for switches (the grid spacings will be checked later) if (itrim(cArg) .gt. 0) then call Say(cSay,'I','CmdLine',cArg) if (LocFirst(cSwitch//'pointp',cArg) .ne. 0) I__MODE = I__MODE+TOM__POINTP ! Point-P map as initial condition if (LocFirst(cSwitch//'mod' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__MOD ! Mod vs nomod version if (LocFirst(cSwitch//'diff' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__DIFF ! Differences instead of ratios if (LocFirst(cSwitch//'mass' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__MASS ! Goes to SW_Model_Kinematic if (LocFirst(cSwitch//'inclg' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__INCLG ! Goes to LOSIntegralG if (LocFirst(cSwitch//'inclv' ,cArg) .ne. 0) I__MODE = I__MODE+TOM__INCLV ! Goes to LOSIntegralV if (LocFirst(cSwitch//'trackg',cArg) .ne. 0) I__MODE = I__MODE+TOM__TRACKG ! Tracks G-convergence in addition to G^2 end if bPointP = iand(I__MODE,TOM__POINTP) .ne. 0 bMod = iand(I__MODE,TOM__MOD ) .ne. 0 bDiff = iand(I__MODE,TOM__DIFF ) .ne. 0 bTrackG = iand(I__MODE,TOM__TRACKG) .ne. 0 I = 0 if (.not. bMod ) I = I+Str2Str('in quasi-time-dependent mode',cStr(I+1:)) if ( bMod ) I = I+Str2Str('in time-independent mode' ,cStr(I+1:)) if (.not. bDiff) I = I+Str2Str('#using ratios of model and observed los values#.',cStr(I+1:)) if ( bDiff) I = I+Str2Str('#using differences of model and observed los values#.',cStr(I+1:)) call Say(cSay,'I','Running',cStr) call IPSProcessing(bForecast,NitTotal,bVcon,bVNago,bVOoty,bGcon,bGNago,bGCamb,bGOoty,bMV2) if (bVNago .or. bGNago) call ForeignFile(iVar,cVar,'nagoya',cWildNagoya) if (bVOoty .or. bGOoty) call ForeignFile(iVar,cVar,'ooty' ,cWildOoty ) ! If bMV2=.TRUE. then either bGcon or bVcon is .FALSE. and the other is .TRUE. if (bMV2) then ! Either bGcon or bVcon is .FALSE. if (bVcon) SMOOTH_D = SMOOTH_V if (bGcon) SMOOTH_V = SMOOTH_D else ! Both bGcon and bVcon are .TRUE. if (bVNago .and. bGNago) SMOOTH_D = SMOOTH_V ! G and V from same source end if call AskR4('Power of g-dependence on normalized density' ,PWN) call AskR4('Power of g-dependence on heliocentric distance',PWR) if (bVcon) call IPSConstraints(bForecast,TOM__V,NitTotal,bBoundV,BINXVmin,NitBadV,iEorWV,XElowV,XEhighV,XRlimV) if (bGcon) then if (bVNago .and. bGNago) then bBoundG = bBoundV BINXGmin = BINXVmin NitBadG = NitBadV iEorWG = iEorWV XElowG = XElowV XEhighG = XEhighV XRlimG = XRlimV else call IPSConstraints(bForecast,TOM__G,NitTotal,bBoundG,BINXGmin,NitBadG,iEorWG,XElowG,XEhighG,XRlimG) end if call AskR4('Solar wind speed for G-level traceback',V0) end if !------- ! Establish defaults for numerical grid (updated in SetGrid) nLng = nLng0 nLat = nLat0 nRad = nRad0 nTim = nTim0 !------- ! Set up the t3d array. The array is saved internally. From here on access to ! its elements should be done only using T3D_iget, T3D_iset, T3D_get, T3D_set. call T3D_iset(T3D__MODE ,0,I__MODE) !------- ! Store grid dimensions. ! Then check whether arrays are big enough to handle them. call T3D_iset(T3D__NLNG ,0,nLng ) call T3D_iset(T3D__NLAT ,0,nLat ) call T3D_iset(T3D__NRAD ,0,nRad ) call T3D_iset(T3D__NTIM ,0,nTim ) call T3D_set (T3D__PWN_V ,0,PWN ) call T3D_set (T3D__PWN_G ,0,PWN ) call T3D_set (T3D__PWR_V ,0,PWR ) call T3D_set (T3D__PWR_G ,0,PWR ) call T3D_set (T3D__D1AU ,0,D1AU ) call T3D_set (T3D__SMOOTH_V ,0,SMOOTH_V) call T3D_set (T3D__SMOOTH_D ,0,SMOOTH_D) call T3D_set (T3D__FILL_V ,0,SMOOTH_V) call T3D_set (T3D__FILL_D ,0,SMOOTH_D) call T3D_set (T3D__SMOOTH_TIME_V,0,SMOOTH_TIME) call T3D_set (T3D__SMOOTH_TIME_D,0,SMOOTH_TIME) call T3D_iset(T3D__ITERATION ,0,0 ) call T3D_iset(T3D__ITERATIONS ,0,NitTotal) call T3D_set (T3D__BINX_V ,0,BINXVmin) call T3D_set (T3D__BINX_D ,0,BINXGmin) call T3D_set (T3D__CLIPLNG ,0,ClipLng) call T3D_set (T3D__SCALE ,4,Scale) ! Set time limits on data selection call SetGrid(bForecast, nCar, JDCar, XCbegMAT, XCendMAT, XCbegROI, XCendROI, & XCtst1, XCtst2, MJDref, bGCamb, MJDfrst, MJDlast, cCam, iEdt, NT) XCrange = XCendMAT(1)-XCbegMAT(1) !------- ! Pick up constants set by SetGrid call T3D_iget(T3D__MODE ,0,I__MODE) call T3D_iget(T3D__NCOFF,0,NCoff ) call T3D_get_grid(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1) call T3D_show() nSS1 = nLng*nLat ! # points at source surface at given time nSST = nSS1*nTim ! # points at source surface at all times combined n3DT = nSST*nRad ! # point in reconstruction space (spatial + temporal) if (nSS1 .gt. nSS0 ) call Say(cSay,'F','Too','many bins at source surface; nLng*nLat too large') if (nSST .gt. nSST0) call Say(cSay,'F','Too','many bins at source surface; nLng*nLat*nTim too large') if (n3DT .gt. n3DT0) call Say(cSay,'F','Too','many voxels in volume; nLng*nLat*nRad*nTim too large') !------- ! The range of the Carrington variable needed is [XCbegMAT,XCendMAT]. The target range ! [XCtst1,XCtst2] includes an extra safety margin to make sure that all relevant ! data are read from file. NLmaxG = NLG0 NLmaxV = NLV0 NLOS = NLOS0 dLOS = dLOS0 call T3D_set (T3D__DLOS,0,dLOS ) call T3D_iset(T3D__NLOS,0,NLOS ) ! Only one of (bVNago, bVOoty) is .TRUE. call ReadVIPSValidData(LOS__CHECKV,Bad,Bad,Bad,Bad) ! Check for valid V data if (bVNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya, & nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref, & RR,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.,XCtst1,XCtst2,NCoff,MJDref, & RR,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) ! Only one of (bGCamb, bGOoty, bGNago) is .TRUE. ! NOTE: the dummy array rLOSTmp is needed to avoid overwriting the array VOBS call ReadVIPSValidData(LOS__CHECKG,Bad,Bad,Bad,Bad) ! Check for valid G data if (bGCamb) call ReadGIPS(cCam,nCar,JDCar, & iEdt,.FALSE.,bForecast, & XCtst1,XCtst2,MJDref,MJDfrst,MJDlast, & RR,V0,GPOWER,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG, & iXPG,iYPG,iMJDG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS, & 1,iXPsav,iYPsav,iMJDsav,IYRSsav,DOYSsav,XDSsav,XLSsav, & XLLsav,XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGNago) call ReadVIPS(iReadNagoya,iProcessNagoya,cWildNagoya, & nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref, & RR,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG, & IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) if (bGOoty) call ReadVIPS(iReadOoty ,iProcessOoty ,cWildOoty, & nCar,JDCar,.FALSE.,XCtst1,XCtst2,NCoff,MJDref, & RR,XElowG,XEhighG,iEorWG,XRlimG,NLmaxG,NLG, & IYRFG,IRECG,IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,rLOSTmp,GOBS, & 1,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) call T3D_iset(T3D__NL_V,0,NLV ) call T3D_iset(T3D__NL_G,0,NLG ) call StopWatch(' ',' ') ! Start stopwatch !------- ! The nRad entry in the t3d array may be modified by LOSReach. ! Save the value here so we can restore it later. nRadFull = nRad !------- ! Process the line-of-sight data: ! - LOSPosition : Calculate positions of all segments ! - LOSReach : Calculate range of positions of all segments; lines of sight ! that fall outside the numerical grid are flagged as bad in NBS ! - LOSClean : iscard lines-of-sight flagged by LOSReach ! - LOSWeights : Calculate line-of-sight weights for all segments if (bVcon) then call LOSPosition(TOM__V,IYRSV,DOYSV,DISTV,XLSV,XCEV,XLLV,XDLV,PPlosV,VPlosV) call LOSReach(TOM__V,XCbegMAT,XCendMAT,PPlosV,NBSV) call LOSClean(TOM__V,.FALSE.,NBSV,iXPsav,iYPsav,iMJDsav, & IYRFV,IRECV, IYRSV,DOYSV,DISTV,XLSV,XLLV,XDLV,XCEV,XEV,XCV,YLV,VOBS,PPlosV) IPSV = 1 if (bVNago) IPSV = 1 ! Calculate weights, WWlosV, along all lines of sight call LOSWeights(TOM__V,IPSV,XEV,DISTV,PPlosV,WWlosV,rLOSTmp) call Say(cSay,'I',' ','.#.') end if if(bGcon) then call LOSPosition(TOM__G,IYRSG,DOYSG,DISTG,XLSG,XCEG,XLLG,XDLG,PPlosG,rLOSTmp) call LOSReach(TOM__G,XCbegMAT,XCendMAT,PPlosG,NBSG) call LOSClean(TOM__G,bGCamb,NBSG,iXPG,iYPG,iMJDG, & IYRFG,IRECG, IYRSG,DOYSG,DISTG,XLSG,XLLG,XDLG,XCEG,XEG,XCG,YLG,GOBS,PPlosG) IPSG = 2 ! Default: Cambridge if (bGNago) IPSG = 1 ! Calculate weights, WWlosG, along all lines of sight call LOSWeights(TOM__G,IPSG,XEG,DISTG,PPlosG,WWlosG,M2Mean) call Say(cSay,'I',' ','.#.') end if call T3D_iget(T3D__NL_V,0,NLV ) ! Updated by LOSClean call T3D_iget(T3D__NL_G,0,NLG ) ! Updated by LOSClean call T3D_iget(T3D__NRAD,0,nRad) ! Updated by LOSReach if (bGcon) call ArrR4TimesArrR4(NLG,GOBS,GOBS,G2OBS) ! Switch to g^2 !------- ! At this point all lines of sight to be used in the reconstruction have ! been collected. Before starting the reconstruction set up a standard ! for the 'goodness of fit', based on a 'flat' source surface velocity ! and density map (these are saved internally by LOSTweak) call ArrR4Constant(nSST,V0,VFin) call ArrR4Constant(nSST,1.,GFin) call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VFin,GFin,XC3D,V3D,G3D,TmpSST) if (NLV .gt. 0) then call LOSIntegralV(TOM__V,XCbegMAT,XCendMAT,V3D,G3D,PPlosV,WWlosV,VPlosV,VMDL,Weight,VSUM,rLOSTmp) call LOSTweak(TOM__V,VOBS,VMDL,NBSV,FIX,FIXMEAN,FIXSTDV) end if if (NLG .gt. 0) then call LOSIntegralG(TOM__G,XCbegMAT,XCendMAT,G3D,M2Mean,PPlosG,WWlosG,G2MDL,Weight,rLOSTmp) if (bTrackG) then call ArrR4Function(NLG,sqrt,G2MDL,GMDL) call LOSTweak(TOM__G,GOBS,GMDL,NBSG,FIX,FIXMEAN,FIXSTDV) end if call LOSTweak(TOM__G,G2OBS,G2MDL,NBSG,FIX,FIXMEAN,FIXSTDV) end if !------- ! Set up the initial maps at the source surface. ! The default is constant V, constant normalized density call ArrR4Constant(nSST,1.,GMap) call ArrR4Constant(nSST,V0,VMap) !------- ! Replace the default by a 'point-P' map. The point-P map may have holes in it ! so call GridSphere2D to fill the holes (it also smooths the map) if (bPointP) then NC = 1 if (bVcon) then call MapGrid(XCbegMAT(1),XCendMAT(1),YLbeg,YLend,NC, NLV,XCV,YLV,VOBS, & NAV,DAV,nLng,nLat,VMap,Gtmp,Vtmp,Zmin,Zmax,ZDmin,ZDmax) call GridSphere2D(XCrange,nLng,nLat,nTim,VMap,SMOOTH_V,3,0.0,ClipLng) N = 1 do I=2,nTim ! Use same V source surface map for all times N = N+nSS1 call ArrR4Copy(nSS1,VMap,VMap(N)) end do end if if (bGcon) then call MapGrid(XCbegMAT(1),XCendMAT(1),YLbeg,YLend,NC,NLG,XCG,YLG,G2OBS, & NAV,DAV,nLng,nLat,GMap,Gtmp,Vtmp,Zmin,Zmax,ZDmin,ZDmax) call GridSphere2D(XCrange,nLng,nLat,nTim,GMap,SMOOTH_D,3,0.0,ClipLng) N = 1 do I=2,nTim ! Use same G source surface map for all times N = N+nSS1 call ArrR4Copy(nSS1,GMap,VMap(N)) end do end if end if call LOSSanityCheck(bVcon,bGcon) !------- ! Iteration loop do Nit=1,NitTotal call Say(cSay,'I',' ','.#Iteration '//cStr(:Int2Str(Nit,cStr))) call T3D_iset(T3D__ITERATION,0,Nit) if (bVcon) then ! V deconvolution ! bMV2=.TRUE., implies bGcon=.FALSE. if (bMV2) then call MkV2D(nSST,VMap,GMap,V0) ! Calculate nr^2 call CvD2G(PWN,-nSST,GMap,GMap) ! nr^2 -> g^2 end if call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST) call LOSIntegralV(TOM__V,XCbegMAT,XCendMAT,V3D,G3D,PPlosV,WWlosV,VPlosV,VMDL,Weight,VSUM,rLOSTmp) if (Nit .eq. 1) call ArrR4Copy(NLV,VMDL,VMDL0) call LOSProjection(TOM__V,XCbegMAT,XCendMAT,XC3D,PPlosV,PPprj) call LOSTweak(TOM__V,VOBS,VMDL,NBSV,FIX,FIXMEAN,FIXSTDV) bRemove = Nit .eq. NitBadV c bRemove = Nit/50*50 .eq. Nit call BuildSourceSurface(TOM__V,bRemove,bBoundV,BINXVmin,NBSV,FIX,FIXMEAN,FIXSTDV, & PPprj,Weight,XCbegMAT,XCendMAT,VMap,VSIGB,Vtmp,Gtmp,BINXV,iLOSTmp) if (Nit .eq. NitTotal) call ArrR4Copy(nSST,VMap,VFin) call GridSphere2D(XCrange,nLng,nLat,nTim,VMap,SMOOTH_V,3,0.0,ClipLng) call Say(cSay,'I',' ','.') end if if (bGcon) then ! G deconvolution ! bMV2=.TRUE., implies bVcon=.FALSE. if (bMV2) then call CvG2D(PWN,-nSST,GMap,VMap) ! In: GMap = g^2 -> Out: VMap = nr^2 call MkD2V(nSST,VMap,VMap,V0) ! nr^2 --> Velocity end if call SW_Model_Kinematic(TOM__G,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST) call LOSIntegralG(TOM__G,XCbegMAT,XCendMAT,G3D,M2Mean,PPlosG,WWlosG,G2MDL,Weight,rLOSTmp) if (Nit .eq. 1) call ArrR4Copy(NLG,G2MDL,G2MDL0) call LOSProjection(TOM__G,XCbegMAT,XCendMAT,XC3D,PPlosG,PPprj) if (bTrackG) then call ArrR4Function(NLG,sqrt,G2MDL,GMDL) call LOSTweak(TOM__G,GOBS,GMDL,NBSG,FIX,FIXMEAN,FIXSTDV) end if call LOSTweak(TOM__G,G2OBS,G2MDL,NBSG,FIX,FIXMEAN,FIXSTDV) bRemove = Nit .eq. NitBadG c bRemove = Nit/50*50 .eq. Nit call BuildSourceSurface(TOM__G,bRemove,bBoundG,BINXGmin,NBSG,FIX,FIXMEAN,FIXSTDV, & PPprj,Weight,XCbegMAT,XCendMAT,GMap,GSIGB,Vtmp,Gtmp,BINXG,iLOSTmp) if (Nit .eq. NitTotal) call ArrR4Copy(nSST,GMap,GFin) call GridSphere2D(XCrange,nLng,nLat,nTim,GMap,SMOOTH_D,3,0.0,ClipLng)! Fill gaps and smooth call Say(cSay,'I',' ','.') end if end do call Say(cSay,'I',' ','.#.#Finished iterating#.#.') nRad = nRadFull call T3D_iset(T3D__NRAD,0,nRad) call T3D_iset(T3D__ITERATION,0,Nit) !------- ! The velocity map VFin and density map GFin (as produced by the last ! BuildSourceSurface calls) contain bad values indicating where the source surface ! maps are invalid. GMap and VMap have gone to one more GridSphere2D, ! filling in all bad values and smoothing the array. ! Use VMap and GMap to create the final shift array XC3D. Output from ! SW_Model_Kinematic (XC3D, V3D, G3D) does not contain any bad values (the 'bad ! value' information is stored in the source surface maps VMap, GMap). call SW_Model_Kinematic(TOM__V,XCbegMAT,XCendMAT,VMap,GMap,XC3D,V3D,G3D,TmpSST) ! In the original program there was an inconsistency here: XC3D, which ! is based on VMap and GMap, was combined with VFin and GFin to find ! the final velocity array. Here I'm using VMap and GMap after copying ! the 'holes' from VFin and GFin. To go back to the old treatment, ! comment the next 5 lines !======= call ArrR4Mask(nSST,VFin,Bad,Bad,Bad,0.,1.,VMap)! Copy bad values of VFin into VMap call ArrR4Mask(nSST,GFin,Bad,Bad,Bad,0.,1.,GMap)! Copy bad values of GFin into GMap call ArrR4Copy(nSST,VMap,VFin) ! Copy VMap into VFin call ArrR4Copy(nSST,GMap,GFin) ! Copy GMap into GFin !======= N = -nSS1+1 do I=1,nTim ! Fill in small holes N = N+nSS1 call GridFill(nFill,nLng,nLat,VFin(N),Btmp,Zmin,Zmax) call GridFill(nFill,nLng,nLat,GFin(N),Btmp,Zmin,Zmax) end do !------- ! Write3D_nv_UT does not modify any of its input arguments call Write3D_nv_UT(NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,VFin,GFin,XC3D,V3D,G3D,Tmp3D1,Tmp3D2) !------- ! Write the t3d* file. ! Write3D_nv returns velocity in V3D and normalized density in GFin and G3D. ! Output arrays contain bad values based on the bad values in the ! input arrays VFin and GFin. call Write3D_nv(XCbegMAT,XCendMAT,XCbegROI,XCendROI,VFin,GFin,XC3D,V3D,G3D,BINXV,BINXG) !------- ! Write the b3d* file. Note that V3D is the full 3D velocity matrix call Write3D_bb ( XCbegMAT,XCendMAT,XCbegROI,XCendROI,XC3D,V3D,Tmp3D1,Tmp3D2,Gtmp) call Write3D_bb_UT(NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,XC3D,V3D,Tmp3D1,Tmp3D2,TmpSST1,TmpSST2,Btmp) if (bVcon) then I = Str2Str('los'//char(TOM__V),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,VMDL0,VMDL,XEV,NBSV) if (bVOoty) call WriteLOSY(cStr,cWildOoty ,NLV,IYRFV,IRECV,VOBS,VMDL0,VMDL,XEV,NBSV) end if if (bGcon) then I = Str2Str('los'//char(TOM__G),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,G2OBS,G2MDL0,G2MDL,XEG,NBSG) if (bGOoty) call WriteLOSY(cStr,cWildOoty ,NLG,IYRFG,IRECG,G2OBS,G2MDL0,G2MDL,XEG,NBSG) if (bGCamb) call WriteLOSD(cStr,cCam,iEdt,NLG,iMJDG,iXPG,iYPG,G2OBS,G2MDL0,G2MDL,XEG,NBSG) end if call Say(cSay,'I','Removed',cStr(:Int2Str(NLV-iArrI4Total(NLV,NBSV,I),cStr))//' bad V sources') call Say(cSay,'I','Removed',cStr(:Int2Str(NLG-iArrI4Total(NLG,NBSG,I),cStr))//' bad G sources') call StopWatch(' ',' ') call Say(cSay,'S','Stop','Program successfully completed') end