;read_ascat_averaged_bytemaps

PRO get_ascat_averaged, filename, windspd, winddir, scatflag, radrain, sosmap, ierr

; this routine will read the RSS ASCAT scatterometer time_averaged bytemap files
;
; The only input argument is the filename with full path.
; This routine assumes the file is zipped.  Remove compression from the openr command if your files are unzipped.
;
; The routine returns:
;   wspd, wdir, scatflag, radrain, sosmap
;   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)  (extracted from bit rain data)
;   radrain is the collocated radiometer rain rate in mm/hr
;		derived from an SSMIS, SSMIS, or WindSat, observation collocated within 180 minutes of the scatterometer observation
;		radrain contains -999. where there is no collocated radiometer data
;						    0. where there is radiometer data, but there is no measureable rain
;						   -1. where there is radiometer data, no measureable rain, but nearby cells have rain
;						    0.2 to 12.5  mm/hr radiometer rain rate for that cell
;   sosmap is the sum-of-squares map, which is a measure of the departure of the actual wind retrievals backscatter ratio
;          from the 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)
;          (sosmap = binaryvalue * 0.02)
;
;  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
;
; The data contained in the daily files are scalar (wind speed, sosmap) and vector (wind direction) averaged
; to produce the values in the time-averaged files.  A data value for a given cell is only provided in a
; time-averaged file if a minimum number of data exist within the time period being produced
; (3-day maps, 2 obs;  week maps, 5 obs;  month maps, 20 obs)
;     3-day   =(average of 3 days ending on file date)
;     weekly  =(average of 7 days ending on Saturday of file date)
;     monthly =(average of all days in month)
;
; 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


ierr=0

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

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


;determine if file exists
exist=findfile(filename,COUNT=cnt)
if (cnt ne 1) then begin
  print, 'FILE DOES NOT EXIST!!'
  ierr=1
  return

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
  	ierr=1
  	return

  endif else begin
  	readu,2,binarydata
  	close,2
  endelse

  ; multipliers to change binary data to real data
  xscale=[.2,1.5,1, 0.02]

  ; loop through 4 variables
    for ivar=0,3 do begin

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

        nodata=where((dat gt 250B) and (dat le 254B), ict_bad)
        land=where((dat eq 255B), ict_land)

        case ivar of
         0: begin		;wind speed
            windspd[*,*]=dat*xscale[ivar]
    		 end

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

	     	 2: begin		;scatflag and radrain
				    scatflag[*,*] =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= 60 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)
            temp[radok]=0
            radok=where(radflag eq 1 and krain eq 1) ; if there is adjacent rain, set radrain value to -1.0
            temp[radok]=-1
            radok=where(radflag eq 1 and krain ge 2)
            temp[radok]=krain[radok]/5.-0.2
            radrain[*,*]=temp
         end

         3: begin		;sos map
            sosmap[*,*]=dat*xscale[ivar]
         end

         endcase

    endfor	;ivar
endelse

land=where(windspd eq 51.0, ict_land)
if (ict_land gt 0) then begin
  windspd[land]  = -555.
  winddir[land]  = -555.
  scatflag[land] = -555.
  radrain[land]  = -555.
  sosmap[land]   = -555.
endif

bad=where(windspd gt 50.0, ict_bad)
if (ict_bad gt 0) then begin
windspd[bad]  = -999.
winddir[bad]  = -999.
scatflag[bad] = -999.
radrain[bad]  = -999.
sosmap[bad]   = -999.
endif

return
END