C+ C NAME: C ReadVIPSn C PURPOSE: C Reads IPS velocity data from yearly data files into arrays C CATEGORY: C I/O C CALLING SEQUENCE: subroutine ReadVIPSn(iReadVIPS,iProcessVIPS,cWildVIPS, & nCar,JDCar,bAuto,XCtst1,XCtst2,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) C INPUTS: C iReadVIPS external integer function C iProcessVIPS external integer function C external functions controlling 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 XCtst1 real modified Carrington variable for start of search C XCtst2 real modified Carrington variable for end of search C NCoff integer JDCar(I) is the start of rotation NCoff+I C MJDref double precision 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 lower limit on elongation from Sun C XEhigh real upper limit on elongation from Sun C iEorW integer 0: Both east and west data C 1: East data only C -1: West data only C XRlim real limit on difference RA of point-P and C RA of Sun. C C NLmax integer max. # points read into arrays C XC,YL,VV,GG, XCsav,YLsav,VVsav,GGsav C NSmax integer C OUTPUTS: C NL integer # valid data points in arrays XC,YL,VV,GG C SRC (NLmax) character*9 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 IYRS(NLmax) integer time of observation: year 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 C NSmax integer # data points saved in XCsav,YLsav,VVsav C C SRCsav (NSmax) character*9 C IYRFsav(NSmax) integer scratch arrays (used internally only) 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 href=iGetLun= to make sure this does not happen. C CALLS: C iSetFileSpec, iGetFileSpec, Say, iSearch, Str2Flt, Int2Str, bOpenFile, iFreeLun C itrim, SunNewcomb, Julian, bLOSCheckLoc C INCLUDE: include 'dirspec.h' include 'filparts.h' include 'openfile.h' C PROCEDURE: C > All points in the range [XCtst1,XCtst2] 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: iReadNagoya, iProcessNagoya C iReadUCSD , iProcessUCSD C iReadOoty , iProcessOoty C These are programmed in the form: C C function iReadNagoya(iU,iFirst) C ... ( read a record from the data file ) ... C return C entry iProcessNagoya(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 iReadVIPS C 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 nCar integer # of Carrington rotation start times C JDCar double precision Start times of Carrington rotation in C modified Julian days. The first entry C JDCar(1) is for rotation NCoff C NCoff integer Carrington rotation offset C XCend real maximum permitted Carrington variable. C When Earth moves passed XCend then all C data are assumed collected. C Radius real reference distance (AU) C C Output values of iProcessVIPS: C C iProcessVIPS 0: bad IPS observation C 1: good IPS observation C 2: all data collected (Earth passed XCend) C 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 C (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 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/CASS) C OCT-1996, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added Ooty routines C- integer iReadVIPS integer iProcessVIPS character cWildVIPS*(*) integer nCar double precision JDCar(nCar) logical bAuto real XCtst1 real XCtst2 integer NCoff double precision MJDref real Radius real XElow real XEhigh integer iEorW real XRlim integer NLmax integer NL character SRC (NLmax)*9 character SRCsav (NSmax)*9 integer IYRF (NLmax) integer IREC (NLmax) integer IYRS (NLmax) real DOYS (NLmax) real XDS (NLmax) real XLS (NLmax) real XLL (NLmax) real XDL (NLmax) real XCE (NLmax) real XE (NLmax) real XC (NLmax) real YL (NLmax) real VV (NLmax) real GG (NLmax) integer NSmax integer IYRFsav (NSmax) integer IRECsav (NSmax) integer IYRSsav (NSmax) real DOYSsav (NSmax) ! Time: doy of year real XDSsav (NSmax) ! Sun-Earth distance real XLSsav (NSmax) ! Ecliptic longitude Sun real XLLsav (NSmax) real XDLsav (NSmax) real XCEsav (NSmax) real XEsav (NSmax) ! Elongation point P real XCsav (NSmax) ! Carrington variable of Point P real YLsav (NSmax) ! Heliographic latitude of Point P real VVsav (NSmax) ! Solar wind velocity (km/s) real GGsav (NSmax) ! G-factor double precision xTT double precision JEpoch double precision FLINT8 character cSrc*9 character cSay*8 /'ReadVIPS'/ character cHideLogical*(FIL__LENGTH) character cFile*(FIL__LENGTH) integer iFrstYr /0/ integer NS /0/ real XCmargin /1.0/ logical bOpen /.FALSE./ logical bGetRecord logical bReading logical bOpenFile logical bMore logical bReadErr logical bSaveOverrun /.FALSE./ integer nYrLen /4/ ! # digits used to identify year in file name integer Str2Str integer Flt2Str save bOpen, bGetRecord, bSaveOverrun, & iU, iYr, NS, & iFrstYr, iFile logical bLOSCheckLoc integer NLOS /0/ ! By default bLOSCheckLoc accepts all lines of sight save XCbeg,XCend,NLOS,NINSIDE,DLOS !------- ! On first call check whether data files are present if (iFrstYr .eq. 0) then !------- ! Little fudge: replace first * by ????. I = index(cWildVIPS,'*') if (I .ne. 0) cWildVIPS(I:) = cWildChar(:nYrLen)//cWildVIPS(I+1:) I = iSetFileSpec(cWildVIPS) I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile) iFile = index(cFile,cWildChar(:nYrLen)) if (iFile .eq. 0) call Say(cSay,'E','Wildcard','template must contain '//cWildChar(:nYrLen)) if (iSearch(1,cWildVIPS,cFile) .ne. 1) call Say(cSay,'E',cWildVIPS,'no matching files found') I = iSetFileSpec(cFile) I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile) I = 1 call Str2Flt(cFile(iFile:iFile+nYrLen-1),I,Tmp) iFrstYr = nint(Tmp) ! Save year of first file iFile = index(cWildVIPS,cWildChar(:nYrLen)) ! Save position of wildcard %%%% in cWildVIPS end if !------- ! Check the save scratch arrays for valid data points. ! NS can only be unequal zero in AUTO mode. In AUTO mode only XCtst1 ! and XCtst2 change (both increase). All other data selection ! 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 !------- ! Scan the "save" arrays. Only points with XC >= XCtst1 are needed for ! subsequent maps. Discard points with XC < XCtst1. N = 0 do I=1,NS if (XCsav(I) .ge. XCtst1) 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 !-------- ! Pick up all points in the save arrays which are needed for the current ! map. Since all points with XC < XCtst1 have just been discarded, the ! points are only tested for XC <= XCtst2. do I=1,NS if (XCsav(I) .le. XCtst2) then NL = NL+1 SRC (NL) = 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 !------- ! If no file open then check which file is to be read. ! (bOpen=.TRUE. only happens in AUTO mode. Then also bGetRecord=.FALSE., ! indicating an unprocessed record is still in memory). ! The data files are assumed to be organized in yearly files: one file ! for each year of data. XCtst1-XCmargin is used to determine which yearly ! file must be opened first. The value is kept larger inside the range 1,nCar ! to avoid problems accessing the JDCar array. XCtst1 itself is already ! supposed to be set to a safe value in SetXCRange. if (.not. bOpen) then XCt = max(1.0,min(XCtst1-XCmargin,1.0*nCar)) xTT = FLINT8(1,nCar,JDCar,dble(XCt),0.0d0) 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 !------- ! Read data until ! - end point XCtst2 is reached, or ! - EOF is reached, or ! - arrays are full ! bOpen and bGetRecord are used in subsequent calls. ! If a record has been read succesfully, but is rejected because it is out ! of range (iProc=2) or because the arrays are full (NL>=NLmax), then ! the file is kept open (bOpen=.TRUE.) and bGetRecord is set .FALSE. ! On the next call to ReadVIPS, the procedure will skip straight to the ! iProcessVIPS call to try processing the last record read in the previous ! 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 ReadVIPS to initialize record counter !------- ! Complete the data file name by changing the file extension to the ! proper year. Open data file. If bOpenFile returns 0 the file probably ! did not exist (exit do while loop). bOpenFile assigns iU if successfully ! opening the file. cFile = cWildVIPS ! Replace wildcard by year I = Int2Str(iYr,cFile(iFile:iFile+nYrLen-1)) iRecl = 0 bOpen = bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY+OPN__NOMESSAGE,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','Read',cHideLogical(cFile)) bReading = .TRUE. end if bMore = iReadVIPS(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 = iProcessVIPS(iR,nCar,JDCar,NCoff,XCtst2,Radius, & cSrc,nYr,Doy,XCsc,XCobs,xLat,xVV,xGG,dRA, & rLngSun,rDisSun,rLngP,rLatP,rEloP) bMore = iProc .ne. 2 ! Check for XCsc past XCtst2 if (iProc .eq. 1) then ! Record contains valid data with XCsc <= XCtst2 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 (XCtst1 .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 !-------- ! Pick up all points needed for the current map.Test only for XC <= XCtst2 if (XCobs .le. XCtst2) then !-------- ! 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,XCtst1,XCobs,XCtst2 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 !------- ! At this point bMore = .FALSE. if iProc=2 (XCsc past XCtst2) 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+XCtst1,2,cFile(I+1:)) I = I+Str2Str('-',cFile(I+1:)) I = I+Flt2Str(NCoff+XCtst2,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','nobs',cHideLogical(cFile)) return C+ C NAME: C ReadVIPSLOSCheckn C PURPOSE: C CALLING SEQUENCE: entry ReadVipsLOSCheckn(XCbegIn,XCendIn,NLOSIn,NINSIDEIn,DLOSIn) C INPUTS: C XCbegIn real modified Carrington variable for start of map C XCendIn real modified Carrington variable for end of map C NLOSIn integer C NINSIDEIn integer C DLOSIn real C SEE ALSO: C ReadVIPS C PROCEDURE: C Input values are copied to internally saved variable for use by ReadVIPS C MODIFICATION HISTORY: C- XCbeg = XCbegIn XCend = XCendIn NLOS = NLOSIn NINSIDE = NINSIDEIn DLOS = DLOSIn return end