C+ C NAME: C IPSD2020 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 PURPOSE: C Construct deconvolved synoptic maps for the IPS observations using C a "time-dependent" technique that projects all observations to a C single lower reference surface. Note - uneven time cadences are C allowed for velocity and density data sets. C CATEGORY: C Data processing C CALLING SEQUENCE: C run ips2020 C INPUTS: C Data file $CAM:G*. C OUTPUTS: C Rectangular grid of data points C FUNCTIONS/SUBROUTINES: C Local2UT,MAP_TZERO,READ_GO, C SetGipsy C INCLUDE: C include 'sun.h' C COMMON BLOCKS: C SIDE EFFECTS: C RESTRICTIONS: C PROCEDURE: C > The modified Julian day MJD used for the IPS data is related to the C regular Julian day (as used in subroutine Julian) by C MJD = JD-2400000.5 C (i.e. the modified Julian day starts at midnight). C > Once the logical $CAM is set (during the first run of IPSD2020), the C only way to force IPSD2020 to look in another directory is to delete the C logical $CAM. C MODIFICATION HISTORY: C JUN-1992, Paul Hick (UCSD) C NOV-1995, Bernie Jackson (STEL) C March-1999 Bernie Jackson (UCSD) C June -1999 Bernie Jackson (UCSD) C August -1999 Bernie Jackson (UCSD) C Jan/Feb 2001 Bernie Jackson (UCSD) C Mar 2001 Bernie Jackson (UCSD) C- program IPSD2020 parameter (NLmaxG = 12200, ! Max # G data points (13000) (400,000) & NLmaxV = 1000, ! Max # V data points (2000) & NGREST = 1, ! G temporal resolution factor better than 1 day & NVREST = 1, ! V temporal resolution factor better than 1 day & NGRESS = 1, ! Spatial G resolution factor better than 20 degrees & NVRESS = 1, ! Spatial V resolution factor better than 20 degrees & NF = 2, ! LOS resolution factor (Should be 1 if Spatial resolution is >= 2) & nTmaxG = 50*NGREST, ! Max # G times (50) (100) (150) & nTmaxV = 50*NVREST) ! Max # V times (50) (100) (150) parameter (NLOSG = 20*NF*NGRESS, ! # resolution elements along each G line of sight (40) (80) (125) & NLOSV = 20*NF*NVRESS, ! # resolution elements along each V line of sight (40) (80) (125) & NLOSGP1 = NLOSG+1, & NLOSVP1 = NLOSV+1, & dLOSG = 0.1/(NGRESS*NF), ! LOS G Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & dLOSV = 0.1/(NVRESS*NF), ! LOS V Resolution in units of observer-Sun distance (0.05) (0.025) (0.016) & xInc = 0.5) ! Carrington map "extra" increment parameter (nLng = 54*NGRESS+1, ! # Longitudes, resolution is 360/(nLng-1) (55) (109) (217) (325) & iLng = 18*NGRESS+1, ! # Longitudes to use in extracted map (19) (37) (73) (109) & nLat = 9*NGRESS+1, ! # Latitudes , resolution is 180/(nLat-1) (9)(10) (19) (37) (55) & nMap = 31, ! # Radial heights. Range covered is [0,(nMap-1)*dRR) (31) & nLngLat = nLng*nLat, & nLngLatnTG =nLng*nLat*nTMaxG, & nLngLatnTV =nLng*nLat*nTMaxV) parameter (dRR = 0.1) ! Resolution (AU) in radial direction (0.1) 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), & XCtpr (NLOSGP1,NLMAXG), & XCtim (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) 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,nTmaxG), & VM (nLng,nLat,nTmaxG), & 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), & V3D (nLng,nLat,nMap), & D3D (nLng,nLat,nMap), & V3DT (nLng,nLat,nMap,nTmaxG,4), & D3DT (nLng,nLat,nMap,nTmaxG,4), & Vtmp (nLng,nLat,nTmaxG,3), ! Scratch real arrays & VZtmp (nLng,nLat,nTmaxV), & Dtmp (nLng,nLat,nTmaxG), & XLTtmp (nLng,nTmaxG,4), & 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 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 /.FALSE./, ! 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./, & bWrite3D /.FALSE./, & bWriteMS /.FALSE./, & 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 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 / 7/, & 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 /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 /13./, & CONRD /13./, & CONVT /0.5/, & CONDT /0.5/, & 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 / 0.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./, & PWG /0.34/, !/0.5//1.0//2.0/ ! Power of density & PWV /0.34/, !/0.5//1.0//2.0/ ! Power of density & PWRG /0.34/ , ! Radial power of density & PWRV /0.34/ ! Radial power of velocity real*8 MJDref, & MJDcntr character cCam*4, & cMon*3, & cStrID*11, & cVar(2)*100, & cArg*100, & cStr*80, & cPre*4 character cWildNagoya*50, & cWildOoty*50 external iReadNagoyan, iProcessNagoya n external iReadOotyn, iProcessOotyn external EARTH include 'dirspec.h' include 'sun.h' integer Flt2Str iVar = 2 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 iYr = 1985 C iYr = 1990 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('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 = 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) 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) 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 = 1965 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 one shift matrix per iteration$no',bOnly1) call AskYN('Output a Three D file$yes',bWrite3D) ! 3D file write if(bWrite3D) then NinterC = 3 call AskI4('Number of intermediate rotation densities with holes?',NinterC) NinterD = 0 call AskI4(' Number of intermediate filled rotation densities?',NinterD) 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) 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. 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) NBSV(I) = 0 if (VOBS(I) .lt. VLIL) NBSV(I) = 0 end do do I=1,NLmaxG GGOBS(I) = GGOBS(I)*GLIMF if (GGOBS(I) .gt. GLIMP) NBSG(I) = 0 if (GGOBS(I) .lt. GLIMM) NBSG(I) = 0 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 if(bVcon) & call MkTimesvd(bDaydV,bForecast,MJDcntr,NmidHRV,aNdayV,XCtest1,XCtest2,xInc,NTV,NTmaxV, & nCar,JDCar,ALng,XCintV,XCbeV,XCtbegV,XCtendV,XCtioffV,ItoffV,IYRV,DOYV) if(bGcon) & call MkTimesvd(bDaydG,bForecast,MJDcntr,NmidHRG,aNdayG,XCtest1,XCtest2,xInc,NTG,NTmaxG, & nCar,JDCar,ALng,XCintG,XCbeGG,XCtbegG,XCtendG,XCtioffG,ItoffG,IYRG,DOYG) if(.not.bGcon) then nTG = nTV ItoffG = ItoffV XCtbegG = XCtbegV XCtendG = XCtendV do I=1,nTmaxG XCintG(I) = XCintV(I) XCbegG(1,I) = XCbeV(1,I) XCbegG(2,I) = XCbeV(2,I) end do end if if(.not.bVcon) then nTV = nTG ItoffV = ItoffG XCtbegV = XCtbegG XCtendV = XCtendG do I=1,nTmaxV XCintV(I) = XCintG(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 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, XCintV, ItoffV print *, XCbeGG, XCintG, ItoffG 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 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 call Say('IPSD2020','I','Info','# g-level observations ' & //cStr(:Int2Str(NLG,cStr))) call Say('IPSD2020','I','Info','# V observations ' & //cStr(:Int2Str(NLV,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 MkShiftd(XCbeV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV, & VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,CONRV,CONRD,CONVT,aNdayV,VLL,VUL, & NSIDE,Vtmp,Dtmp,VZtmp,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,RRS,dRR,dLOSV,nLng,nLat,nMap,nTV, & XCshift,NLV,NLOSVP1,DISTV,XLSV,XCEV,XLLV,XDLV,IYRSV,DOYSV, & XLON,XLONP,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) call Get4Dval(XCbeV,XCtbegV,XCtendV,RRV,dRR,nLng,nLat,nMap,nTV, & DDfact,NLOSVP1*NLV,XCtim,XLONP,XLAT,RPX,DFAC) call Get4Dval(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,XLAT,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, DMAPV,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 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 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 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 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)) end do end if end if if(bGcon) then call MkShiftd(XCbeGG,XCtbegG,XCtendG,ALng,nLng,nLat,nMap,nTG, & VMAPD,DMAP,XCshift,DVfact,DDfact,RRS,dRR,CONRV,CONRD,CONVT,aNdayV,VLL,VUL, & NSIDE,Vtmp,Dtmp,VZtmp,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,RRS,dRR,dLOSG,nLng,nLat,nMap,nTG, & XCshift,NLG,NLOSGP1,DISTG,XLSG,XCEG,XLLG,XDLG,IYRSG,DOYSG, & XLON,XLONP,XLAT,RPX,XCtpr,XCtim,LON,LAT,ITIM) call Get4Dval(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,XLAT,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,N cd IPS 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 MkShiftd(XCbeV,XCtbegV,XCtendV,ALng,nLng,nLat,nMap,nTV, & VMAP,DMAPV,XCshift,DVfact,DDfact,RRS,dRR,CONRV,CONRD,CONVT,aNdayV,VLL,VUL, & NSIDE,Vtmp,Dtmp,VZtmp,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 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 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 if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then ! Do this if smoother V maps needed call CopyDtoDVN(2,XCbeV,XCtbegV,XCtendV,XCtbegV,XCtendV,nLng,nLat, & nTV,nTV,ConsTV,DMAPV,VMAP,VVT,VVV) else call ArrR4Copy(nLngLatnTV,DMAPV,VMAP) ! In any case, replace dummy array end if 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 XCmidV = (XCintV(N)+XCintV(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region of interest on map 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 if(nTG .eq. nTV .and. XCbegG(1,1).eq. XCbeV(1,1)) then ! Do this if smoother D maps needed call CopyDtoDVN(2,XCbegG,XCtbegG,XCtendG,XCtbegG,XCtendG,nLng,nLat, & nTG,nTG,ConsTG,VMAPD,DMAP,DDT,DDD) else call ArrR4Copy(nLngLatnTV,VMAPD,DMAP) ! In any case, replace dummy array end if 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 XCmidG = (XCintG(N)+XCintG(N+1))/2.0 - ((1.0 - RRS)*((400./Speed)*4.)/25.38) ! Mid region of interest on Map call FillMapOLS(XCbegG(1,N),XCbegG(2,N),XCmidG,nLng,nLat,ConstL,DMAP(1,1,N)) ! Remove 180 deg bad spots 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,VMAPD,DMAP,XCshift,DVfact,DDfact,NCoff,XCstrt,nCar,JDCar) if(bWrite3D) then 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(0,cPre,Nit,NiterT,NinterC,NCoff,XCintG,XCbegG,XCtbegG,XCtendG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D,V3DT,D3DT) call Extractd3D(bExtract,RRS,dRR,XCbeGG,XCintG,XCtbegG,XCtendG,nLng,nLat, & nMap,nTG,NinterC+1,nTmaxG,V3DT,D3DT,NCoff,XCstrt,nCar,JDCar) do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),1.0/RRS*RRS,DMAP(1,1,N)) end do 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 call FillMapL(2,XCbeV(1,N),XCbeV(2,N),nLng,nLat,ConstL,VMAP(1,1,N)) ! Remove seam end if end do 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 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 do N=1,nTG call ArrR4TimesConstant(-nLngLat,DMAP(1,1,N),RRS*RRS,DMAP(1,1,N)) end do cPre = 'nv3f' call write3D_infotd(1,cPre,Nit,NiterT,NinterD,NCoff,XCintG,XCbegG,XCtbegG,XCtendG,RRS,dRR, & nLng,nLat,nMapn,nTF,nTmaxG, & PWV,PWG,PWRV,PWRG,DEN1AU,CONRV,CONRD,CONSTV,CONSTG,CONVT,CONDT, & 90.0,Scale,VMAPD,DMAP,VM,DM,XCshift,DVfact,DDfact,V3D,D3D) end if call Say('IPSD2020','S','Stop','program successfully completed') end