PRO get_ascat_daily, filename, mingmt, windspd, winddir, scatflag, radrain, sosmap
;
; This routine reads RSS ASCAT scatterometer daily bytemaps
; This code is similar to the quikscat bytemap read routine.  It has one additional parameter (SOS maps) and a longer time period for radiometer collocation (3 hours).
; The greater collocation time with radiometers has little effect since ASCAT retrievals are of good quality in rain except for at winds below 5 m/s.
; The SOS can be used to flag bad data and may perform better than the radiometer flag.  Data removal is up to the user.
;
; filename with full path is the only input argument.
; This routine assumes the file is *zipped*.  Remove compression from openr command if your files are unzipped.
;
; The routine returns:
;   mingmt, wspd, wdir, scatflag, radrain, sosmap
;   mingmt is the GMT time in minutes of day (mingmt = binaryvalue * 6)
;   windspd is the 10 meter surface wind speed in m/s  (windspd = binaryvalue * 0.2)
;   winddir is the wind direction in degrees the wind is blowing to, North = 0 degrees  (winddir = binaryvalue * 1.5)
;   scatflag is the scatterometer rain flag (0=no rain, 1=rain) derived from the rain bitflag variable
;   radrain is the collocated radiometer rain rate in mm/hr
;   sosmap is the sum of squares value which is a measure of the departure of the actual wind retrievals backscatter ratio
;          from the correspondent theoretical ratio described by the geophysical model function C-2015
;          (ratio of sum (sigma0_obs-sigma0_gmf)^2 scaled by the variance of observed sigma0, non-dimensional quantity, scale by 0.02)
;          The SOS is a "goodness of fit" measure. Where the values are high, the observation is different from what we expect and likely has greater error
;
;		SSMI, SSMIS, WindSat,  observations within 180 minutes of ASCAT
;		radrain contains -999. where there is no collocated radiometer data
;		                   -1. where there is radiometer data, no measureable rain, but nearby cells have rain
;						            0. where there is radiometer data, but there is no measureable rain
; 					            0.2 to 12.5  mm/hr radiometer rain rate for that cell
;
;  The center of the first cell of the 1440 column and 720 row map is at 0.125 E longitude and -89.875 latitude.
;  The center of the second cell is 0.375 E longitude, -89.875 latitude.
; 		XLAT=0.25*ILAT-90.125
;		  XLON=0.25*ILON-0.125
;
; Data format description : www.remss.com/missions/ascat
; RSS support available at: www.remss.com/support
;
; updated April 2016, D.Smith for new V02.1 data
; fixed typo in rain scaling, D.Smith,  April 2016

;binary data in file
binarydata= bytarr(1440,720,5,2)

;output products (lon,lat,asc/dsc)
mingmt  =fltarr(1440,720,2)
windspd =fltarr(1440,720,2)
winddir =fltarr(1440,720,2)
scatflag=intarr(1440,720,2)
radrain =fltarr(1440,720,2)
sosmap  =fltarr(1440,720,2)
radflag =intarr(1440,720)
krain=radflag


;determine if file exists
exist=findfile(filename,COUNT=cnt)
if (cnt ne 1) then begin
  print, 'FILE DOES NOT EXIST!!'
  mingmt[*,*,*]=1530

endif else begin

  ;open file, read binary data, close file
  close,2
  openr,2,filename, error=err, /compress	;compress keyword allows reading of gzip file
  if (err gt 0) then begin
  	print, 'ERROR 1 WITH FILE: ', filename
  endif else begin
  	readu,2,binarydata
  	close,2
  endelse

  ; multipliers to change binary data to real data
  xscale=[6.,.2,1.5, 1, 0.02]
      land=where((binarydata[*,*,1,*] eq 255B), ict_land)

  ; loop through asc/dsc  and all 4 variables
  for iasc=0,1 do begin
    for ivar=0,4 do begin

		; extract 1 variable, scale and assign to real array
        dat=binarydata[*,*,ivar,iasc]

        case ivar of
        	0: begin		;gmt time
        	    mingmt[*,*,iasc]=dat*xscale[ivar]
             end

			    1: begin		;wind speed
              windspd[*,*,iasc]=dat*xscale[ivar]
             end

			    2: begin		;wind direction
            	winddir[*,*,iasc]=dat*xscale[ivar]
             end

			    3: begin		;scatflag and radrain
			        scatflag[*,*,iasc] =dat-2*floor(dat/2)	; bit 1 is scat flag of 0=no rain, 1=rain
			        radflag =floor((dat-4*floor(dat/4))/2)  ; bit 2 is indicator of radiometer collocation, 0= no rad data, 1= 180 min rad data collocation
			        krain=floor(dat/4)						; bits 3-8 are rad rain rate with a value of 1 indicating adjacent rain
                temp=fltarr(1440,720)-999.
                radok=where(radflag eq 1, ict )

				      if (ict gt 0) then temp[radok]=0
                radok=where(radflag eq 1 and krain eq 1, ict) ; if there is adjacent rain, set radrain value to -1.0
              if (ict gt 0) then temp[radok]=-1
                radok=where(radflag eq 1 and krain ge 2, ict)
              if (ict gt 0) then temp[radok]=krain[radok]/5.-0.2
                radrain[*,*,iasc]=temp
             end

            4: begin	;sos map
            	  sosmap[*,*,iasc]=dat*xscale[ivar]
               end
         endcase

    endfor	;ivar
  endfor  ;iasc

endelse

bad=where(mingmt gt 1440)
windspd[bad]  = -999.
winddir[bad]  = -999.
scatflag[bad] = -999.
radrain[bad]  = -999.
sosmap[bad]   = -999.

land=where(mingmt eq 1530)
windspd[land]  = -555.
winddir[land]  = -555.
scatflag[land] = -555.
radrain[land]  = -555.
sosmap[land]   = -555.

return
END