C+ C NAME: C BRead_WSO C PURPOSE: C Read file with synoptic map data for the magnetic field at a C source surface (in muT) into a 2D array C This routine applies to the regular WSO files downloaded from C the WSO web site. C CALLING SEQUENCE: integer function BRead_WSO(cFile, nMax, ZFile, nLng, nLat) C INPUTS: C cFile*(*) character fully-qualified name of source surface file C nMax integer Max. # allowed grid points in synoptic map C OUTPUTS: C BRead_WSO integer return status (0=Failure,1=Success) C C ZFile(nLng,nLat) real 2D synoptic map with radial magnetic field C in muT (micro Tesla) C nLng integer # grid longitudes C nLat integer # grid latitudes C INCLUDE: include 'openfile.h' include 'filparts.h' include 'ftspar.h' C CALLS: C ArrR4Bad, Say, nrInterpol, Str2Flt, ArrR4Step, ArrR4Copy, ArrR4Reverse C bOpenFile, iFreeLun, iOSgunzip, iOSDeleteFile C FTNOPN, FTGISZ, FTGKYJ, FTGKYJ, FTGPVE, FTCLOS C SEE ALSO: C BListAll, MapReadSingle, MapReadTimes, BList_WSO, BRead_WSO_NOAA C PROCEDURE: C > The BRead* functions are called as external functions by MapReadSingle C and MapReadTimes. C > The magnetic source surfaces files contain the radial component of C the magnetic field at some source surface. Currently they come in C two flavors: C - containing 73x29 arrays at a 5 deg resolution covering 360 degree C in longitude and [-70,70] in latitude, or C - containing 73x30 arrays at a 5 deg resolution in longitude covering C 360 degree, and covering -14.5/15 to +14.5/15 in 30 equal steps in C sine latitude. C > The array ZFile(nLng,nLat) is assumed to refer to a single rotation C covering 360 degrees in longitude, and [-90,+90] in latitude). C The part of the array above and below the latitudes available in the C file are flagged as 'bad'. C > The start longitude (i.e. the heliographic longitude of column C Z[nLng,*) is determined from the file name by function BList_WSO. C MODIFICATION HISTORY: C JUL-1993, Paul Hick (UCSD/CASS) C JUL-2001, Paul Hick (UCSD/CASS) C Removed a redundant call to iGetLun C MAY-2002, Paul Hick (UCSD/CASS) C Revamped to work with the tomography program. C Support for GRPACK has been dropped. C JUL-2002, Tamsen Dunn (UCSD/CASS) C Check for empty record at end of file added. C AUG-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Bug fix. For incomplete files the missing columns were not C explicitly flagged as bad, so they contained garbage. C DEC-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added reading of FITS files (as created by IDL routine C wso_write.pro). C- character cFile*(*) integer nMax real ZFile(*) integer nLng integer nLat character cSay*9 /'BRead_WSO'/ character cSeverity /'E'/ !------- ! Maps are returned on a 73x37 grid with a resolution of 5 degrees in longitude ! and latitude. parameter (nLngFile = 73) ! (73-1)*5= 360; # longitudes parameter (nLatFile = 37) ! (37-1)*5=2*90 parameter (dLng = 360.0/(nLngFile-1)) parameter (dLat = 180.0/(nLatFile-1)) !------- ! The source surface magnetic field is stored on a grid of 73x29 ! points, or 73x30 points. ! ! The resolution in longitude is 5 degrees, i.e. 73 longitudes cover ! 360 degrees. The latitudes either cover [-70,+70] degrees in steps ! of 5 degrees (29 points), or cover sine(latitude) in [-14.5,+14.5]/15 ! in equal steps in sine latitude (30 points). parameter (nLatSteps = 29) ! (29-1)*5=2*70 parameter (TopLat = 70.0) ! Maximum latitude parameter (nSineSteps= 30) parameter (TopSine = 14.5/15.) ! Maximum sine(lat) parameter (dSine = 2.*TopSine/(nSineSteps-1))! Step size in sine(lat) real XLat(nLatSteps) real ZLat(nLatSteps) real XSine(nSineSteps) ! covers 30 sine(lat) real Ztmp (nSineSteps) real Ytmp (nSineSteps) character cRec*80 character cEOF*7 /'_______'/ character cUnzippedFile*256 character cInputFile*256 character cType*5 integer nDim(2) real ZBuf(nLngFile,nSineSteps) logical anyv logical bOpenFile logical bSinLat logical bInit logical bUnzip logical bFits locMap (I,J) = I + nLngFile*(J-1) BRead_WSO = 0 ! Assume something will go wrong bInit = .TRUE. LatOff = (nLatFile-nLatSteps)/2 ! Latitudinal margin; =4 bins (20 deg) north and south !------- ! In theory it would only be necessary to flag top and bottom LatOff ! rows as bad. For incomplete files (does happen!) some columns would ! not be filled in, so to be safe we flag the whole map as bad. call ArrR4Bad(nLngFile*nLatFile, ZFile) !------- ! Check first whether cFile is a gzipped file (with extension .gz) ! If it is, the name of the unzipped file is returned in cUnzippedFile. bUnzip = iOSgunzip(cFile,cUnzippedFile) .eq. 1 if (bUnzip) then cInputFile = cUnzippedFile else cInputFile = cFile end if I = iSetFileSpec(cInputFile) I = iGetFileSpec(FIL__TYPE,FIL__TYPE,cType) bFits = cType .eq. '.fts' .or. cType .eq. '.fit' .or. cType .eq. '.fits' if (bFits) then iU = iGetLun(cInputFile) call Say(cSay,'I','#'//cInputFile,'reading') I = 0 ! Open fits file call FTNOPN(iU, cInputFile, FTS__READONLY, I) if (I .eq. 0) then call FTGISZ(iU, 2, nDim, I) ! Get array dimensions if (I .ne. 0) call say_fts(cSay,cSeverity,I) ! nDim should be 73x30 or 73x29 if (nDim(1) .ne. nLngFile .or. & (nDim(2) .ne. nLatSteps .and. nDim(2) .ne. nSineSteps)) & call Say(cSay,'E','Invalid','dimensions: '//cInputFile) call FTGKYJ(iU, 'CARROT' , iRotStart, cRec, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) call FTGKYJ(iU, 'CARRLONG', iLngStart, cRec, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) call FTGKYJ(iU, 'STEPLONG', iLngStep , cRec, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) ! iRotStart, iLngStart are CR and heliographic longitude for ! the first column in the map. Add one to iRotStart to switch ! to the last column (this assumes that the map is exactly one ! rotation with the heliographic longitude in first and last ! column the same. iRotStart = iRotStart-1 nullv = 0 anyv = .FALSE. nBuf = nDim(1)*nDim(2) ! Get data call FTGPVE(iU, 1, 1, nBuf, nullv, ZBuf, anyv, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) call FTCLOS(iU, I) if (I .ne. 0) call say_fts(cSay,cSeverity,I) iU = iFreeLun(iU) bSinLat = nDim(2) .eq. nSineSteps iLng = iLngStep do N=nDim(1),1,-1 iRot = 0 iLng = iLng-iLngStep do J=1,nDim(2) Ztmp(J) = ZBuf(N,J) end do if (bSinLat) then ! Ztmp contains equal steps in sin(latitude) !------- ! Convert from equal steps in sine(latitude) to equal steps ! in latitude by a spline interpolation if (bInit) then ! Initialize latitude arrays needed for interpolation bInit = .FALSE. !call Say(cSay,'I','Converting','from equal steps in sine(latitude)') call ArrR4Step(-TopLat ,dLat ,nLatSteps ,XLat ) ! Equal steps in latitude (degrees) call ArrR4Step(-TopSine,dSine,nSineSteps,XSine) ! Equal steps in sine(latitude) do I=1,nSineSteps ! Convert to degrees XSine(I) = asind(XSine(I)) end do end if call nrInterpol(nSineSteps,XSine,Ztmp, nLatSteps,XLat,ZLat, Ytmp,4) else ! Ztmp already in equal steps in latitude call ArrR4Copy(nLatSteps,Ztmp,ZLat) end if I = nLngFile - nint( (iRot-iLng/360.0)*(nLngFile-1) ) do J=LatOff+1,LatOff+nLatSteps ZFile( locMap(I,J) ) = ZLat(J-LatOff) end do end do else call Say(cSay,'E','#'//cInputFile,'read error') end if else iAct = OPN__REOPEN+OPN__STOP+OPN__READONLY+OPN__TEXT iRecl = 0 if (bOpenFile(iAct,iU,cInputFile,iRecl)) then !------- ! Look for first record starting with 'CT' read (iU,'(A)',iostat=iEOF) cRec do while (iEOF .eq. 0 .and. cRec(:2) .ne. 'CT') ! Read past header records read (iU,'(A)',iostat=iEOF) cRec end do if (iEOF .ne. 0 .or. cRec(:7) .eq. cEOF) & call Say(cSay,'E',cFile,'read error or empty file') !------- ! Each longitude is stored in four records. The first of these ! starts with CT0000:000 (containing rotation number and ! heliographic longitude). The differences between rotation and ! longitude in each record and the values in the first record are ! used to determine the location in the output array: the column ! with values iRotStart, iLngStart is written to column nLngFile; ! subsequent records are put in column nLngFile-1, nLngFile-2,..,1. read (cRec(3: 6),'(I4)') iRotStart ! First rotation number read (cRec(8:10),'(I3)') iLngStart ! First longitude do while ( iEOF .eq. 0 .and. cRec(:7) .ne. cEOF .and. itrim(cRec) .ne. 0 ) !------- ! First line: get Carrington rotation and longitude read (cRec(3: 6),'(I4)') iRot ! Rotation number read (cRec(8:10),'(I3)') iLng ! Longitude iRot = iRot-iRotStart iLng = iLng-iLngStart !------- ! First line contains 6 numbers following CT0000:000 N = 8 call Str2Flt(cRec(11:),N,Ztmp) !------- ! Read the following three lines. Lines 2 and 3 each contain ! 8 numbers. The third contains 7 numbers (equal steps in latitude) ! or 8 numbers (equal steps in sine latitude). I = N do J=1,3 read (iU,'(A)') cRec N = 8 call Str2Flt(cRec,N,Ztmp(I+1)) I = I+N end do !-------- ! Latitudes are stored from high to low latitude. We need them in reverse order. call ArrR4Reverse(I,Ztmp,Ztmp) ! Reverse order !-------- ! I is total # numbers read for each longitude (29 for equal ! steps in latitude, or 30 for equal steps in sine(latitude)) bSinLat = I .eq. nSineSteps if (bSinLat) then ! Ztmp contains equal steps in sin(latitude) !------- ! Convert from equal steps in sine(latitude) to equal steps ! in latitude by a spline interpolation if (bInit) then ! Initialize latitude arrays needed for interpolation bInit = .FALSE. !call Say(cSay,'I','Converting','from equal steps in sine(latitude)') call ArrR4Step(-TopLat ,dLat ,nLatSteps ,XLat ) ! Equal steps in latitude (degrees) call ArrR4Step(-TopSine,dSine,nSineSteps,XSine) ! Equal steps in sine(latitude) do I=1,nSineSteps ! Convert to degrees XSine(I) = asind(XSine(I)) end do end if call nrInterpol(nSineSteps,XSine,Ztmp, nLatSteps,XLat,ZLat, Ytmp,4) else ! Ztmp already in equal steps in latitude call ArrR4Copy(nLatSteps,Ztmp,ZLat) end if I = nLngFile - nint( (iRot-iLng/360.0)*(nLngFile-1) ) do J=LatOff+1,LatOff+nLatSteps ZFile( locMap(I,J) ) = ZLat(J-LatOff) end do read (iU,'(A)',iostat=iEOF) cRec end do iU = iFreeLun(iU) if ( (iRot .eq. 0 .and. iLng .ne. -360) .or. & (iRot .eq. 1 .and. iLng .ne. 0) ) & call Say(cSay,'W','incomplete','rotation does not cover 360 deg') else call Say(cSay,'E','#'//cInputFile,'read error') end if end if nLng = nLngFile nLat = nLatFile if (bUnzip) I = iOSDeleteFile(cUnzippedFile) BRead_WSO = 1 return end