pro read_rss_mwir_sst, filename, sst, error, mask


;	This routine reads version-4.0 RSS MW_IR Optimal Interpolated SST fusion daily files
;	You do not need to unzip files before reading them
;
;   INPUT
;   filename  with path in form : mwir.fusion.yyyy.doy.v04.0
;	  yyyy = year
;	  doy  = day of year

;	  OUTPUT
;	  sst  (a 4096 x 2048 real*4 array of data)
;	  error(a 4096 x 2048 real*4 array of the interpolation error estimate)
;	  mask (a 4096 x 2048 integer*1 array of data masking information)
;			mask values are
;       bit 0 = 1 if land is present      
;       bit 1 = 1 if sea ice is present
;       bit 2 = 1 if IR satellite data were utilized for this pixel
;       bit 3 = 1 if MW data are utilized for this pixel (if no MW data are available, only previous day data exist)  
;       bit 4 = 1 if the pixel is near land and data is suspect (bad data)  
;       bit 5 = 1 if the pixel is in an unclassified sea ice region.  Could be sea ice or not.
;       bit 6 = not used
;       bit 7 = not used
;
; Latitude/Longitude
;   xcell = grid cell values between 0 and xmax (where xmax = 4095)  IDL is zero-based
;   ycell = grid cell values between 0 and ymax (where ymax = 2047)
;   dx=360./xmax
;   dy=180./ymax
;   Center of grid cell Longitude  = dx*(xcell+1)-dx/2.    (degrees east)
;   Center of grid cell Latitude   = dy*(ycell+1)-(90+dy/2.)  (-90 to 90)
;
;   Please read the data description on www.remss.com
;   See the verify.txt file to determine if you are using this routine correctly
;   Contact http://www.remss.com/support with questions about the data or program.
;
;   updated 8/2010 d.smith
;   updated to version 4.0, 11/2014 d.smith 


binarydata= bytarr(4096,2048)
sst       = fltarr(4096,2048)
error     = fltarr(4096,2048)
mask      = bytarr(4096,2048)
sst_scale = 0.15
sst_offset= -3.0
err_scale = 0.005
err_offset= 0.0


;determine if file exists
exist=File_Test(filename)
if (exist ne 1) then begin
  print, 'file does not exist ', filename
  return
endif


openr,lun,filename, error=err,  /get_lun, /compress   ;compress keyword allows reading of gzip file, remove if data already unzipped
if (err gt 0) then begin
  print, 'error opening file ', filename
  return
endif

readu,lun,binarydata
sst =float(binarydata)
good=where(binarydata le 250)  ;scale good data
sst[good]=binarydata[good] * sst_scale + sst_offset

readu,lun,binarydata
error=float(binarydata)
good=where(binarydata le 250)  ;scale good data
error[good]=binarydata[good] * err_scale + err_offset

readu,lun,binarydata
mask=binarydata
free_lun, lun

mask_land=where(mask and 1)
if mask_land[0] ne -1 then begin
  sst[mask_land]=255
  error[mask_land]=255
endif

mask_ice=where(mask and 2)
if mask_ice[0]  ne -1 then begin
  sst[mask_ice]=252
  error[mask_ice]=252
endif

mask_bad=where(mask and 16)
if mask_bad[0]  ne -1 then begin
  sst[mask_bad]=254
  error[mask_bad]=254
endif

;we advise users to decide if they need to use this flag sea ice flag
mask_bad=where(mask and 32)
if mask_bad[0]  ne -1 then begin
  sst[mask_bad]=252  ;no good ice mask to know for sure, but could be sea ice 
  error[mask_bad]=252
endif
return
end