C+ C NAME: C ReadVIPSn8 C ReadVIPSLOSCheckn C PURPOSE: C Read IPS velocity data into arrays C CATEGORY: C I/O C CALLING SEQUENCE: C call ReadVIPSn8_n(iReadVIPSn8,iProcessVIPSn8,cWildVIPS, C & nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref, C & Radius, C & XElow,XEhigh,iEorW,XRlim, C & NLmax,NL,cbvfiles,SRCVG,SRC,SRCVGsav,SRCsav, C & IYRF,IREC,IYRS,DOYS8,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG, C & NSmax,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) C C call ReadVIPSLOSCheck(XCbeg,XCend,NLOS,NINSIDE,DLOS) C INPUTS: C ReadVIPS: C iReadVIPSn8 external integer function C iProcessVIPSn8 external integer function C external functions control data file access C cWildVIPS character*(*) wildcard used to locate data files; C the wildcard must contain the substring C '%%' C nCar integer dimension of JDCar C JDCar(nCar) integer start times (Julian days) of Carrington rotations C bAuto logical C XCtest1 real modified Carrington variable for start of search C XCtest2 real modified Carrington variable for end of search C NCoff integer JDCar(I) is the start of rotation NCoff+I C MJDref real*8 if not equal 0, only data prior to C modified Julian day MJDref are used C Radius real reference distance (AU) C C XElow real C XEhigh real C iEorW integer C XRlim real C C NLmax integer max. # points read into arrays C XC,YL,VV,GG, XCsav,YLsav,VVsav,GGsav C NSmax integer C C ReadVIPSLOSCheck: C XCbeg real modified Carrington variable for start of map C XCend real modified Carrington variable for end of map C NLOS integer C NINSIDE integer C DLOS C OUTPUTS: C NL integer # valid data points in arrays XC,YL,VV,GG C SRC(NLmax) character*8 Character Source identifier C IYRF(NLmax) integer year of file from which observation was extracted (should be same as IYRS) C IREC(NLmax) integer record number on file IYRF C DOYS8(NLmax)real*8 time of observation: day of year (incl. fraction for time of day) C XDS (NLmax) real Sun-Earth distance C XLS (NLmax) real geocentric ecliptic longitude Sun C XLL (NLmax) real geocentric ecliptic lng(P)-lng(Sun) (deg) C XDL (NLmax) real geocentric ecliptic lat(P) (deg) C XCE (NLmax) real modified Carrington variable of sub-Earth point on Sun C XE (NLmax) real elongation (deg) (>0: East of Sun; <0: West of Sun) C XC (NLmax) real modified Carrington variable point P after traceback to C heliocentric distance Radius at speed VV C YL(NLmax) real heliographic latitude point P C VV(NLmax) real IPS velocity (km/s) C GG(NLmax) real IPS disturbance factor C NSmax integer # data points saved in XCsav,YLsav,VVsav C SRCsav(NSmax) character scratch arrays (used internally only), C IYRFsav(NSmax) integer C IRECsav(NSmax) integer C IYRSsav(NSmax) integer C DOYSsav(NSmax) real*8 C XDSsav (NSmax) real C XLLsav (NSmax) real C XDLsav (NSmax) real C XCEsav (NSmax) real C XEsav (NSmax) real C XCsav (NSmax) real C YLsav (NSmax) real C VVsav (NSmax) real C GGsav (NSmax) real C RESTRICTIONS: C Data files are read using a logical unit number assigned by iGetLun. C Since a data file may be left open on return, there is a potential C danger if the same unit is used somewhere. Assign unit numbers C with iGetLun to make sure this doesn't happen. C FUNCTIONS/SUBROUTINES: C bOpenFile, itrim, SunNewcomb, iSetFileSpec, iGetFileSpec, Julian C EXTERNAL: C INCLUDE: C include 'filparts.h' C include 'sun.h' C COMMON BLOCKS: C SIDE EFFECTS: C PROCEDURE: C > All points in the range [XCtest1,XCtest2] are extracted from the data files C >!!!! The modified Carrington variables are in units offset by NCoff from C the regular Carrington rotation number (i.e. in units which can be C used as index to JDCar) C > Two external user-defined are needed to read the data files. C External functions exist for reading the Nagoya and the UCSD velocity C IPS data: iReadNagoyan8, iProcessNagoyan8 C iReadUCSDn8 , iProcessUCSDn8 C iReadOotyn8 , iProcessOotyn8 C These are programmed in the form: C C function iReadNagoyan8(iU,iFirst) C ... ( read a record from the data file ) ... C return C entry iProcessNagoyan8(iRec,nCar,JDCar,NCoff,Radius, C cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA, C rLngSun,rLngP,rLatP,rEloP ) C ... ( process the record to produce the output values ... C ... nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA, ... C ... rLngSun,rLngP,rLatP,rEloP ) C return C end C C Input values to iReadVIPSn C iU integer logical unit number assigned by bOpenFile C iFirst integer 1 : initializes record counter C (used when reading 1st record of file) C 0 : adds one to record counter C Output values of iProcessVIPSn: C iRec integer record # on file C cSrc character Source name C nYr integer year and .. C Doy real .. day of year of observation C XCsc real modified Carington variable for sub-SC point C (NCoff already subtracted) C XCobs real modified Carrington variable of point P after C trace-back to distance Radius (NCoff already subtracted) C xLat real heliographic latitude point P (deg) C xVV real IPS velocity (km/s) C xGG real IPS disturbance factor C dRA real RA(point P)-RA(Sun) (deg) C rLngSun real geocentric ecliptic longitude Sun (deg) C rLngP real geocentric ecliptic longitude of point P (deg) C (relative to Sun) C rLatP real geocentric ecliptic latitude of point P (deg) C rEloP real geocentric source elongation (deg) C NOTE: the sign is used to indicate East/West of Sun C rEloP>0: East of Sun; rEloP<0: West of Sun C MODIFICATION HISTORY: C JUN-1994, Paul Hick (UCSD) C OCT-1996, Paul Hick (UCSD); added Ooty routines C- subroutine ReadVIPSn8_n(iReadVIPSn,iProcessVIPSn,cWildVIPS, & nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref, & Radius,XElow,XEhigh,iEorW,XRlim, & NLmax,NL,cvgfiles,SRCVG,SRC,SRCVGsav,SRCsav, & IYRF,IREC,IYRS,DOYS8,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG) real*8 JDCar(nCar), & dLngSun, dLatSun, dDisSun, & MJDref, & xTT,JEpoch, & Doy8, & DOYS8(NLmax) ! Time: doy of year integer IYRF(NLmax), & IREC(NLmax), & IYRS(NLmax) real XDS (NLmax), ! Sun-Earth distance & XLS (NLmax), ! Ecliptic longitude Sun & XLL (NLmax), & XDL (NLmax), & XCE (NLmax), & XE (NLmax), ! Elongation point P & XC (NLmax), ! Carrington variable of Point P & YL (NLmax), ! Heliographic latitude of Point P & VV (NLmax), ! Solar wind velocity (km/s) & GG (NLmax) ! G-factor character cWildVIPS*(*), & cWildIPS*80, & cvgfiles(*)*11, & SRC(NLmax)*9,cSrc*9,SRCVG(NLmax)*66,cSRCVG*66, & cSay*9 /'ReadVIPSn8'/, & cFile*80,cFiles*11 integer iFrstYr /0/ real XCmargin /1./ logical bAuto, & bOpen /.FALSE./, & bGetRecord, & bReading, & bOpenFile, & bMore, & bReadErr, & bSaveOverrun /.FALSE./ integer nYrLen /4/ ! # digits used to identify year in file name integer Str2Str, Flt2Str external EARTH8 include 'dirspec.h' include 'filparts.h' include 'openfile.h' include 'sun.h' logical bLOSCheckLocn8 logical bLOSCheck integer NLOS /0/ character cFmt*40 /'(A9,3I2,F6.2,F5.2,I4,I5,I4,3I5,I4,F8.5)'/ character cMon*3, & cSrce*9,cSRCVGe*66 print*,' ' print*,'In readvipsn8_n.f' C------- C On first call check whether data files are present call Julian(1,iYr1,Doy,JDCar(int(XCtest1)),JEpoch) write (*,'(A,F9.4,A,I5)') 'In readvipsn8_n. The IPS V analysis begins after Doy',Doy,' of',iYr1 call Julian(1,iYr2,Doy,JDCar(int(XCtest2+1.0)),JEpoch) write (*,'(A,F9.4,A,I5)') 'In readvipsn8_n. The IPS V analysis ends before Doy',Doy,' of',iYr2 iFirst = 1 Iyrno = 1 if(iYr1.ne.iYr2) Iyrno = 2 NL = 0 Ifst = 1 do N=1,Iyrno if(Ifst.eq.1) iYr = iYr1 if(Ifst.eq.2) iYr = iYr2 write (*,*) cWildVIPS cWildIPS = cWildVIPS III = iSetFileSpec(cWildIPS) III = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile) iFile = index(cFile,cWildChar(:nYrLen)) III = 1 call Str2Flt_Int(.TRUE.) call Str2Flt(cFile(iFile:iFile+4-1),III,0) iFile = index(cWildIPS,cWildChar(:4)) cFile = cWildIPS III = Int2Str(iYr,cFile(iFile:iFile+4-1)) C write (*,*) cFile IL1 = iGetFileSpec(FIL__NAME,FIL__TYPE,cFiles) open (13, file=cFile,status='old',recl=120,access='sequential',form='formatted',iostat=iReadIPS) if(iReadIPS.ne.0.and.Ifst.eq.2) then close(13) print*, 'No 2nd year IPS file found.' go to 1111 ! assume if no file exists on second year, the second year has not yet happened end if IL1 = iSetFileSpec(cFile) IL1 = iGetFileSpec(FIL__NAME,FIL__TYPE,cFiles) print *, cFiles do II=1,8000 read (13,cFmt,iostat=iReadNagoyan8) cSrce,iY,iM,iD,UT,DP,iHLA,iHLO,iGLA,iGLO,iROT,iVEL,iDVEL,GVAL write(cSRCVGe(1:66),cFmt) cSrce,iY,iM,iD,UT,DP,iHLA,iHLO,iGLA,iGLO,iROT,iVEL,iDVEL,GVAL C print*, cSRCVGe if(iY.lt.50) iY = iY + 100 C was iprocessnagoyan begin C entry iProcessNagoyan8(iR,nCar,JDCar,NCoff,XCend,Radius, C & cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA, C & rLngSun,rDisSun,rLngP,rLatP,rEloP) if (iVEL .eq. 0) then ! No velocity available: reject iProcessNagoyan8 = 0 go to 1111 end if iR = iRecc cSRCVG = cSRCVGe cSrc = cSrce nYr = 1900+iY cMon = ' ' ! Make sure iM is used call DATE_DOY(0,nYr,cMon,iM,iD,iDoy) Doy8 = iDoy+dble(UT/24.) XCsc = XMAP_SC_POS8(EARTH8,nYr,Doy8,-nCar,JDCar) C print*,'to here 0',iDoy,UT,Doy8,XCsc,XCtest2 if (XCsc .eq. BadR4()) then ! Time outside JDCar range iProcessNagoyan8 = 0 go to 1111 end if if (XCsc .lt. XCtest1) then ! Sub-Earth before XCtest1 iProcessNagoyan8 = 1 iRecc = II + 1 go to 1111 end if ifirst = 0 iRecc = (1-iFirst)*iRecc+1 ! iFirst=1 initializes record counter. ifirst turns zero at the time XCsc > XCtest1 ! Sub-Earth position if (XCsc .gt. XCtest2) then ! Sub-Earth past XCend iProcessNagoyan8 = 2 go to 1111 end if C Print*,'to here 1a',xCsc,XCtest2 XCobs = iROT+(360.-iGLO)/360. ! Point-P at DP XCobs = XCobs+Sun__Spiral*(Radius-DP)/iVEL ! Point-P traced to Radius XCobs = XCobs-NCoff ! Subtract offset C Print*,'to here 1b',iRot,-iGLO,XCoff,XCobs xLat = iGLA ! Heliographic lat(P) xVV = iVEL ! IPS velocity xGG = GVAL ! Set to 0 if GVAL not available rLng = iHLO ! Heliocentric ecliptic lng(P)-lng(Earth) rLat = iHLA ! Heliocentric ecliptic lat(P) rP = DP ! Distance Sun-P (AU) call POINT_ON_LOS(rLng,rLat,rP,rEloSun,rEloP,i) c call PointOnLos(rLng,rLat,rP,rEloSun,rEloP,i) rLngP = rLng ! Geocentric ecliptic lng(P)-lng(Sun) rLatP = rLat ! Geocentric ecliptic lat(P) rEloP = -i*rEloP ! rEloP>0 : East of Sun ! <0 : West of Sun call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun) C print*,'to here 2', XCobs rDisSun = dDisSun rLngSun = dLngSun ! Geocentric ecliptic lng(Sun) rSun = rLngSun rLng = rSun+rLng ! Geocentric ecliptic lng(P) Doy = sngl(Doy8) call ECLIPTIC_EQUATOR(0,nYr,Doy,rLng,rLat) ! rLng = RA of P rLat = dLatSun ! Geocentric ecliptic lat(Sun) call ECLIPTIC_EQUATOR(0,nYr,Doy,rSun,rLat) ! rSun = RA of Sun dRA = rLng-rSun ! RA(P)-RA(Sun) if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360. C print*,'to here 3',dRA,XCobs C was iprocessnagoyan end if (XCobs .le. XCtest2) then C-------- C Check what part of line of sight is inside the actual boundaries [XCbeg,XCend] C if (bLOSCheckLocn8(XCbeg,XCend,NLOS,NINSIDE,DLOS, C & nYr,Doy8,rLngP,rLatP,XCsc,rLngSun,rDisSun)) then C if(NL.gt.0) print*, 'before bLOSCheckLocn8',NL,XCbeg,XCend,LOS,NINSIDE,DLOS,nYr,Doy8,rLngP,rLatP,XCsc,rLngSun,rDisSun C bLOSCheck = bLOSCheckLocn8(XCbeg,XCend,NLOS,NINSIDE,DLOS,iYr,Doy8,rLngLOS,rLatLOS,ObsXC,rLngSun,rDisSun) C if(NL.gt.0) print*, ' after bLOSCheckLocn8',NL,XCbeg,XCend,LOS,NINSIDE,DLOS,nYr,Doy8,rLngP,rLatP,XCsc,rLngSun,rDisSun C print*,'bLOSCheck =',bLOSCheck,XCbeg,XCend C bLOSCheckLocn8 = NLOS .eq. 0 C if (.not.bLOSCheck) go to 1111 ObsLng = rLngSun-180. Inside = 0 do J=1,NLOS ! Loop over all segments along LOS P1 = rLngLOS ! Topocentric ecliptic longitude diff. with Earth-Sun line P2 = rLatLOS ! Topocentric ecliptic latitude P3 = (J-.5)*DLOS ! Distance along LOS (units of ObsR) ! Convert to heliocentric coordinates call POINT_ON_LOS(P1,P2,P3,ELO,ELOS,iEorW) P3 = P3*rDisSun ! Convert to AU P1 = ObsLng+P1 ! Heliocentric longitude of P call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,P1,P2) ! Convert to heliographic coordinates P1 = XMAP_OBS_POS(ObsXC,P1) ! Modified Carrington variable of P if (XCbeg .le. P1 .and. P1 .le. XCend) Inside = Inside+1 ! ok end do C print*,'to here 4', P1,P2,P3 bLOSCheck = Inside .ge. NINSIDE if (XCobs .le. XCtest2) then if(bLOSCheck) then C-------- C Check what part of line of sight is inside the actual boundaries [XCbeg,XCend] if (NL .lt. NLmax) then NL = NL+1 cvgfiles(NL) = cFiles SRCVG(NL) = cSRCVG SRC(NL) = cSrc IYRF(NL) = iYr IREC(NL) = iR IYRS(NL) = nYr ! Time of observation DOYS8(NL)= Doy8 XDS (NL) = rDisSun ! Sun-Earth distance (AU) XLS (NL) = rLngSun ! Geocentric ecliptic longitude of Sun (deg) XLL (NL) = rLngP ! Geocentric ecliptic lng(P)-lng(Sun) XDL (NL) = rLatP ! Geocentric ecliptic lat(P) XCE (NL) = XCsc ! Carrington variable sub-Earth point XE (NL) = rEloP XC (NL) = XCobs YL (NL) = xLat VV (NL) = xVV GG (NL) = xGG c write (*,'(4F10.3)') XC(NL),YL(NL),VV(NL),GG(NL),TTsav(I) c write (*,'(I4,F6.2,4F10.2)') ID,UT,XCtest1,XCobs,XCtest2 else go to 1111 end if else LOS = LOS+1 end if end if if(NL.gt.0.and.NL.le.100.and.NL.ne.NLlast) then if(NL.eq.1) then print*,' ' print*,'The first check 100 in readvipsn8_n (check)' end if print*,SRCVG(NL) write(*,'(I4,1X,A,3I7,3F10.4)') NL,SRC(NL),IYRF(NL),IREC(NL),IYRS(NL),DOYS8(NL),XDS(NL),XLS(NL) write(*,'(3X,9F8.3)') XLL(NL),XDL(NL),XCE(NL),XE(NL),XC(NL),YL(NL),VV(NL),GG(NL) end if NLlast = NL end if 1111 continue end do Ifst = Ifst + 1 close(13) end do write(*,'(A,I6,A)'),'There were ',NL,' sources processed' return end