C+ C NAME: C BList_NSO_NOAA C PURPOSE: C Retrieve information about source surface file for magnetic field. C This routine applies to the WSO files we download from NOAA (courtesy C of Nick Arge). C CALLING SEQUENCE: integer function BList_NSO_NOAA(cFile, NCoff, TFound,XCFoundBeg,XCFoundEnd,NrFound,R0,Scale) C INPUTS: C cFile*(*) character fully-qualified name of file containing C source surface map C NCoff integer Carrington offset subtracted from C TFound, XCFoundBeg, XCFoundEnd C OUTPUTS: C BList_NSO_NOAA logical return status (0=Failure,1=Success) C (actually is always 1) C TFound real Time assigned to source surface map C XCFoundBeg real Start Carrington variable C XCFoundEnd real End Carrington variable C NrFound integer version number C R0 real source surface distance (in AU) C Scale real conversion factor from muT to nT) C The data read by BRead_NSO_NOAA will be C multiplied by this constant, i.e. C BRead_NSO_NOAA returns data in nT. C INCLUDE: include 'openfile.h' include 'filparts.h' include 'sun.h' include 'ftspar.h' C EXTERNAL: external EARTH C CALLS: C iSetFileSpec, iGetFileSpec, Str2Flt, Str2Str, Flt2Str, iFilePath, iFreeLun C bOpenFile, Julian, SunNewcomb, ECLIPTIC_HELIOGRAPHIC, N_CARRINGTON C Str2Dbl, Str2Flt_Format C SEE ALSO: C BListAll, MapReadSingle, MapReadTimes, BList_WSO, BRead_WSO_NOAA C PROCEDURE: C The WSO files from NOAA have file names of the type C wso1986_010_1.dat, i.e. containing a Carrington rotation number, C a heliographic longitude in degrees, and a file version number. C Records in the file cover latitudes from =70 to +70 degrees. C Each record is prefixed with a string of type CT1986:010. C The CR and longitude in the file name correspond to the LAST C record in the file (NOT the first as in the regular WSO files). C C For the example wso1986_010_1.dat, the output values would be: C C TFound = 1986+(1-10/360) - 0.152 = 1986.9722 - 0.1520 = 1986.8202 C XCFoundBeg = 1985.9722 C XCFoundEnd = 1986.9722 C NrFound = 1 C R0 = 0.069786090 (= 15 solar radii in AU) C Scale = 1000.0 (= is conversion from muT to nT) C C If the date file 'datefile.dat' is present in the same C directory as the magnetic field files than the time TFound C is updated using information from that file. The first date C file we got had 152 entries for which TFound was decreased by C 0.152 +/- 0.004 rotations (about -4 days). C MODIFICATION HISTORY: C JUL-2002, Paul Hick (UCSD/CASS) C AUG-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added automatic determination of data file format. C- character cFile*(*) integer NCoff real TFound real XCFoundBeg real XCFoundEnd integer NrFound real R0 character cSay *14 /'BList_NSO_NOAA'/ character cSeverity /'E'/ character cDir *(FIL__LENGTH) character cName*(FIL__LENGTH) character cType*7 character cDateName*12 /'datefile.dat'/ character cDateFile*(FIL__LENGTH) character cFmt*30 !/'(F11.4,I4,2I3,2I5,I4,I2)'/ character cStr*(FIL__LENGTH) character cTmp*20 character cEOF*7 /'_______'/ real rVec(3) real*8 dVec(8) real*8 JD real*8 JEpoch real*8 dLng real*8 dLat real*8 dDis real*8 dvar logical bOpenFile logical bFound integer tt(2) integer Str2Str integer Flt2Str BList_NSO_NOAA = 0 I = iSetFileSpec(cFile) I = iGetFileSpec(0,FIL__DIRECTORY,cDir) I = iGetFileSpec(FIL__NAME,FIL__NAME,cName) I = iGetFileSpec(FIL__TYPE,FIL__TYPE,cType) I = iFilePath(cDir,0,' ',cDateName,cDateFile) !------- ! Retrieve the time assigned to the synoptic map, ! and the start and end Carrington variable. ! For now we assume that this can be obtained from ! the file name, e.g. nso1889_360_1 I = 3 call Str2Flt(cName,I,rVec) if (I .eq. 2) then ! Missing version number: set to 1 rVec(3) = 1 I = 3 end if if (I .ne. 3) then cStr = cFile call Say(cSay,'E','#'//cStr,'invalid file name') end if !------- ! Here we assume that the file names contain ! a rotation number and a start longitude TFound = rVec(1)-NCoff+(1.0-rVec(2)/360.0) XCFoundBeg = TFound-1.0 XCFoundEnd = TFound NrFound = nint(rVec(3)) !------ ! Empirical corrections based on content of 150 ! entries in the date file (dT = 0.152 +/- 0.004) TFound = TFound-0.152 R0 = 15.0*sngl(SUN__RAU) ! Source surface at 15 Rsun Scale = 1000.0 ! Conversion factor from file units to nT !------- ! Check the header of input file for the time of observation ! For Fits files this should be in keyword OBSTIME. ! For Ascii files this should be in a record OBSTIME=YYYY_DOY_hhmmss bFound = .FALSE. if (cType(:4) .eq. '.fts') then ! Fits file iU = iGetLun(cFile) I = 0 ! Open fits file call FTNOPN(iU, cFile, FTS__READONLY, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) if (I .eq. 0) then call FTGKYS(iU, 'OBSTIME', cStr, cTmp, I) bFound = I .eq. 0 call FTGKYE(iU, 'RSS', rVec(1), cTmp, I) if (I .eq. 0) R0 = rVec(1)*sngl(SUN__RAU) call FTCLOS(iU, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) end if iU = iFreeLun(iU) else ! Ascii (.dat) file iAct = OPN__REOPEN+OPN__STOP+OPN__READONLY+OPN__TEXT I = 0 if (bOpenFile(iAct,iU,cFile,I)) then cTmp = 'OBSTIME=' I = itrim(cTmp) !------- ! Look for first record starting with 'CT' read (iU,'(A)',iostat=iEOF) cStr ! Read past header records do while (iEOF .eq. 0 .and. cStr(:2) .ne. 'CT' .and. .not. bFound) bFound = cStr(:I) .eq. cTmp ! Check for 'OBSTIME=' if (bFound) then ! OBSTIME= found cStr = cStr(I+1:) ! YYYY_DOY_hhmmss time string else read (iU,'(A)',iostat=iEOF) cStr end if end do if (iEOF .ne. 0 .or. cStr(:len(cEOF)) .eq. cEOF) then cStr = cFile call Say(cSay,'E','#'//cStr,'read error or empty file') end if iU = iFreeLun(iU) end if end if ! If bFound = .TRUE. then the OBSTIME keyword was found. ! If bFound = .FALSE. then try the date file. if (bFound) then call Time2Split('SMEI',cStr,tt) call Time2Carrington(0,tt,ivar,dvar) TFound = ivar-NCoff+sngl(dvar) else if (bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY+OPN__NOMESSAGE,iU,cDateFile,iRecl)) then !------- ! Look for date file to check for better time TFound read (iU,cStr,iostat=I) if (I .eq. 0) then n = 8 call Str2Dbl(cStr,n,dVec) call Str2Flt_Format(cFmt) call Say(cSay,'I','Format',cFmt) end if read (cStr,cFmt,iostat=I) JD,iDoy,iMon,iDay,iYr,iRot,iLng,iVer !read (iU,cFmt,iostat=I) JD,iDoy,iMon,iDay,iYr,iRot,iLng,iVer do while (I .eq. 0 .and. .not. bFound) ! Loop over all records in date file !------- ! This assumes that the actual file names never use longitude 0. if (iLng .eq. 0) then iRot = iRot+1 iLng = 360 end if !------- ! Construct a file name from the information found in the date ! file, and compare it with the actual file name. If there is ! a match use the matching record to update TFound ! If the file name does not have a version number then assume that the entry in ! the date file with version number one is the one needed. write (cStr,'(A3,I4.4,A1,I3.3,A1,I1)') 'nso',iRot,'_',iLng,'_',iVer bFound = cStr .eq. cName .or. (iVer .eq. 1 .and. cStr .eq. cName(:itrim(cName))//'_1') if (bFound) then call Julian(1,iYr,Doy,JD+2440000d0,JEpoch) ! Time of last magnetogram call SunNewcomb(0,iYr,Doy,dLng,dLat,dDis) ! Ecliptic coordinates Sun rLng = sngl(dLng+180.d0) ! Ecliptic coordinates Earth rLat = -sngl(dLat) call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,rLng,rLat) ! Heliographic coordinates Earth N0 = N_CARRINGTON(iYr,Doy) ! Integer Carrington rotation number TFile = N0-NCoff+(1.0-rLng/360.0) ! Carrington variable I = 0 I = I+Str2Str('T=' , cStr(I+1:)) I = I+Flt2Str(NCoff+TFile ,-2,cStr(I+1:)) I = I+Str2Str(',' , cStr(I+1:))+1 I = I+Flt2Str((TFile-TFound)*sngl(SUN__SYNODICP),-2,cStr(I+1:))+1 I = I+Str2Str('days from estimate (+=East)',cStr(I+1:)) !call Say(cSay,'I',cName,cStr) TFound = TFile ! Update TFound end if read (iU,cFmt,iostat=I) JD,iDoy,iMon,iDay,iYr,iRot,iLng,iVer end do iU = iFreeLun(iU) if (.not. bFound) then I = 0 I = I+Str2Str('no time info; estimate is T=',cStr(I+1:)) I = I+Flt2Str(NCoff+TFound,-2,cStr(I+1:)) !call Say(cSay,'I',cName,cStr) end if end if BList_NSO_NOAA = 1 return end