C+ C NAME: C ReadVIPSn C PURPOSE: C Read IPS velocity data into arrays C CATEGORY: C I/O C CALLING SEQUENCE: C call ReadVIPSn(iReadVIPSn,iProcessVIPSn,cWildVIPS, C & nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref, C & Radius, C & XElow,XEhigh,iEorW,XRlim, C & NLmax,NL,SRC,SRCsav, C & IYRF,IREC,IYRS,DOYS,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG, C & NSmax,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, C & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) C INPUTS: C ReadVIPS: C iReadVIPSn external integer function C iProcessVIPSn 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 DOYS(NLmax) real 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 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 CALLS: 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: iReadNagoyan, iProcessNagoyan C iReadUCSDn , iProcessUCSDn C iReadOotyn , iProcessOotyn C These are programmed in the form: C C function iReadNagoyan(iU,iFirst) C ... ( read a record from the data file ) ... C return C entry iProcessNagoyan(iRec,nCar,JDCar,NCoff,Radius, C cSrc,nYr,Doy,XCsc,XCobs,xLat,xVV,xGG,dRA, C rLngSun,rLngP,rLatP,rEloP ) C ... ( process the record to produce the output values ... C ... nYr,Doy,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 ReadVIPSn(iReadVIPSn,iProcessVIPSn,cWildVIPS, & nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref, & Radius, & XElow,XEhigh,iEorW,XRlim, & NLmax,NL,SRC,SRCsav, & IYRF,IREC,IYRS,DOYS,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG, & NSmax,IYRFsav,IRECsav,IYRSsav,DOYSsav,XDSsav,XLSsav,XLLsav, & XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav) real*8 JDCar(nCar), & MJDref, & xTT,JEpoch integer IYRF(NLmax), IYRFsav(NSmax), & IREC(NLmax), IRECsav(NSmax), & IYRS(NLmax), IYRSsav(NSmax) real DOYS(NLmax), DOYSsav(NSmax), ! Time: doy of year & XDS (NLmax), XDSsav (NSmax), ! Sun-Earth distance & XLS (NLmax), XLSsav (NSmax), ! Ecliptic longitude Sun & XLL (NLmax), XLLsav (NSmax), & XDL (NLmax), XDLsav (NSmax), & XCE (NLmax), XCEsav (NSmax), & XE (NLmax), XEsav (NSmax), ! Elongation point P & XC (NLmax), XCsav (NSmax), ! Carrington variable of Point P & YL (NLmax), YLsav (NSmax), ! Heliographic latitude of Point P & VV (NLmax), VVsav (NSmax), ! Solar wind velocity (km/s) & GG (NLmax), GGsav (NSmax) ! G-factor character cWildVIPS*(*), & SRC(NLmax)*9,SRCsav(NSmax)*9,cSrc*9, & cSay*9 /'ReadVIPSn'/, & cFile*256 integer iFrstYr /0/, & NS /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 include 'dirspec.h' include 'filparts.h' include 'openfile.h' save bOpen, bGetRecord, bSaveOverrun, & iU, iYr, NS, & iFrstYr, iFile logical bLOSCheckLoc integer NLOS /0/ save XCbeg,XCend,NLOS,NINSIDE,DLOS C------- C On first call check whether data files are present if (iFrstYr .eq. 0) then I = iSetFileSpec(cWildVIPS) I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile) iFile = index(cFile,cWildChar(:nYrLen)) if (iFile .eq. 0) call Say(cSay,'E','INVTMPLTE', & 'wildcard template must contain '//cWildChar(:nYrLen)) if (iSearch(1,cWildVIPS,cFile) .ne. 1) & call Say(cSay,'E','FILNOTFOUND', C & 'no files matching wildcard found :'//cWildVIPS) & 'no files matching wildcard found :') I = iSetFileSpec(cFile) I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile) I = 1 call Str2Flt_Int(.TRUE.) ! Save year of first file call Str2Flt(cFile(iFile:iFile+nYrLen-1),I,iFrstYr) iFile = index(cWildVIPS,cWildChar(:nYrLen)) ! Save position of wildcard %%%% in cWildVIPS end if C------- C Check the save scratch arrays for valid data points. C NOTE: NS can only be unequal zero in AUTO mode. In AUTO mode only XCtest1 C and XCtest2 change (both increase). All other data selection C criteria remain the same. NL = 0 ! Initialize data point counter LOS = 0 ! Counter for rejected lines of sight by bLOSCheckLoc if (NS .ne. 0) then ! Only in AUTO mode C------- C Scan the `save' arrays. Only points with XC >= XCtest1 are needed for C subsequent maps. Discard points with XC < XCtest1. N = 0 do I=1,NS if (XCsav(I) .ge. XCtest1) then N = N+1 SRCsav(N) = SRCsav(I) IYRFsav(N) = IYRFsav(I) IRECsav(N) = IRECsav(I) IYRSsav(N) = IYRSsav(I) DOYSsav(N) = DOYSsav(I) XDSsav (N) = XDSsav (I) XLSsav (N) = XLSsav (I) XLLsav (N) = XLLsav (I) XDLsav (N) = XDLsav (I) XCEsav (N) = XCEsav (I) XEsav (N) = XEsav (I) XCsav (N) = XCsav (I) YLsav (N) = YLsav (I) VVsav (N) = VVsav (I) GGsav (N) = GGsav (I) end if end do NS = N ! Valid # point in save arrays C-------- C Pick up all points in the save arrays which are needed for the current C map. Since all points with XC < XCtest1 have just been discarded, the C points are only tested for XC <= XCtest2. do I=1,NS if (XCsav(I) .le. XCtest2) then NL = NL+1 SRC(N) = SRCsav(I) IYRF(NL) = IYRFsav(I) IREC(NL) = IRECsav(I) IYRS(NL) = IYRSsav(I) DOYS(NL) = DOYSsav(I) XDS (NL) = XDSsav (I) XLS (NL) = XLSsav (I) XLL (NL) = XLLsav (I) XDL (NL) = XDLsav (I) XCE (NL) = XCEsav (I) XE (NL) = XEsav (I) XC (NL) = XCsav (I) YL (NL) = YLsav (I) VV (NL) = VVsav (I) GG (NL) = GGsav (I) endif end do end if C------- C If no file open then check which file is to be read. C (bOpen=.TRUE. only happens in AUTO mode. Then also bGetRecord=.FALSE., C indicating an unprocessed record is still in memory). C The data files are assumed to be organized in yearly files: one file C for each year of data. XCtest1-XCmargin is used to determine which yearly C file must be opened first. if (.not. bOpen) then XCt = XCtest1-XCmargin I = XCt if (I .lt. 1 .or. I .ge. nCar) call Say(cSay,'E','InvXCtest1', & 'start value outside range covered by JDCar') xTT = JDCar(I)+(XCt-I)*(JDCar(I+1)-JDCar(I)) call Julian(1,iYr,Doy,xTT,JEpoch) iYr = max(iYr,iFrstYr) ! Year of file to be opened iYr = iYr-1 ! Will be incremented again in a second bGetRecord = .TRUE. ! Make sure to read record end if C------- C Read data until C - end point XCtest2 is reached, or C - EOF is reached, or C - arrays are full C C bOpen and bGetRecord are used in subsequent calls. C If a record has been read succesfully, but is rejected because it's out C of range (iProc=2) or because the arrays are full (NL>=NLmax), then C the file is kept open (bOpen=.TRUE.) and bGetRecord is set .FALSE. C On the next call to ReadVIPSn, the procedure will skip straight to the C iProcessVIPS call to try processing the last record read in the previous C call before reading the next record from the file that is still open. bMore = .TRUE. bReading = .FALSE. ! Only used to control messages to screen do while (bMore) if (.not. bOpen) then ! Open next year's file iYr = iYr+1 iFirst = 1 ! Force ReadVIPSn to initialize record counter C------- C Complete the data file name by changing the file extension to the C proper year. Open data file. If bOpenFile returns 0 the file probably C didn't exist (exit do while loop). bOpenFile assigns iU if successfully C opening the file. cFile = cWildVIPS ! Replace wildcard by year I = Int2Str(iYr,cFile(iFile:iFile+nYrLen-1)) bOpen = bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY,iU,cFile,iRecl) bMore = bOpen ! Open error will exit while loop end if bReadErr = .FALSE. if (bOpen .and. bGetRecord) then if (.not. bReading) then call Say(cSay,'I','Reading',cFile) bReading = .TRUE. end if bMore = iReadVIPSn(iU,iFirst) .eq. 0 if (bMore) then iFirst = 0 ! Record counter has been initialized else ! Read error: EOF assumed iU = iFreeLun(iU) bOpen = .FALSE. bReadErr = .TRUE. end if end if if (bReadErr) then ! After read error, open new file bMore = .TRUE. else if (bMore) then iProc = iProcessVIPSn(iR,nCar,JDCar,NCoff,XCtest2,Radius, & cSrc,nYr,Doy,XCsc,XCobs,xLat,xVV,xGG,dRA, & rLngSun,rDisSun,rLngP,rLatP,rEloP) bMore = iProc .ne. 2 ! Check for XCsc past XCtest2 if (iProc .eq. 1) then ! Record contains valid data with XCsc <= XCtest2 call Julian(10,nYr,Doy,xTT,JEpoch) ! All XE (East+West), only positive XE (East) or only negative XE (West) Rdum = (1-abs(iEorW))*abs(rEloP)+iEorW*rEloP if ((XElow .le. Rdum .and. Rdum .le. XEhigh) .and. & ( abs(dRA) .le. XRlim ) .and. & (MJDref .eq. 0 .or. xTT .le. MJDref) & ) then if (XCtest1 .le. XCobs) then if (bAuto) then bSaveOverrun = NS .ge. NSmax if (.not. bSaveOverrun) then NS = NS+1 SRCsav(NS) = cSrc IYRFsav(NS) = iYr IRECsav(NS) = iR IYRSsav(NS) = nYr DOYSsav(NS) = Doy XDSsav (NS) = rDisSun XLSsav (NS) = rLngSun XLLsav (NS) = rLngP XDLsav (NS) = rLatP XCEsav (NS) = XCsc XEsav (NS) = rEloP XCsav (NS) = XCobs YLsav (NS) = xLat VVsav (NS) = xVV GGsav (NS) = xGG end if end if C-------- C Pick up all points needed for the current map.Test only for XC <= XCtest2 if (XCobs .le. XCtest2) then C-------- C Check what part of line of sight is inside the actual boundaries [XCbeg,XCend] if (bLOSCheckLoc(XCbeg,XCend,NLOS,NINSIDE,DLOS, & nYr,Doy,rLngP,rLatP,XCsc,rLngSun,rDisSun)) then if (NL .lt. NLmax) then NL = NL+1 SRC(NL) = cSrc IYRF(NL) = iYr IREC(NL) = iR IYRS(NL) = nYr ! Time of observation DOYS(NL) = Doy 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 call Say(cSay,'I','NLmax','Data arrays full: not all data processed') bMore = .FALSE. end if else LOS = LOS+1 end if end if end if end if end if C------- C At this point bMore = .FALSE. if iProc=2 (XCsc past XCtest2) or the data arrays are full bGetRecord = bMore .or. .not. bGetRecord end if end do if (.not. bAuto .and. bOpen) then ! If not in AUTO mode, make sure to close open file iU = iFreeLun(iU) bOpen = .FALSE. end if if (bSaveOverrun) then ! Can only happen in AUTO mode call Say(cSay,'W','NSmax','Save arrays overrun: contents discarded') NS = 0 bSaveOverrun = .FALSE. end if I = 0 I = I+Int2Str(NL,cFile(I+1:))+1 I = I+Str2Str('sources in',cFile(I+1:))+1 I = I+Flt2Str(NCoff+XCtest1,2,cFile(I+1:)) I = I+Str2Str('-',cFile(I+1:)) I = I+Flt2Str(NCoff+XCtest2,2,cFile(I+1:)) if (NLOS .ne. 0) then I = I+2 I = I+Int2Str(LOS,cFile(I+1:))+1 I = I+Str2Str('; sources rejected',cFile(I+1:)) end if call Say(cSay,'I','# obs',cFile) return entry ReadVipsLOSCheckn(XCbegIn,XCendIn,NLOSIn,NINSIDEIn,DLOSIn) XCbeg = XCbegIn XCend = XCendIn NLOS = NLOSIn NINSIDE = NINSIDEIn DLOS = DLOSIn return end