PRO rd_sbuv_pmf,file,asc,$                       ;inputs
                v8tot,v8umk,v8vmr,npts, $        ;outputs
                lat,lon,sza,hour,refl,doy,oz, $
		nofilter=nofilter, no_saafilter=no_saafilter, $
		filter_oz=filter_oz, v8ak=v8ak
		
; This is a program to read in SBUV V8.7 PMF files
; file is the name of the file to be read
; asc=1 to read ascending data; asc=0 to read descending data 
; (descending data periods: N17 all; N11 1997-2001, N14 2002-2005, N16 2010-2014)

; v8tot= Profile Total Ozone
; v8umk= Ozone Profile (DU) in 21 Layers
; v8vmr= Ozone Profile (VMR) on 15 Pressure Levels
; npts= Number of Profiles after filtering
; lat= Profile Latitude
; lon= Profile Longitude
; sza= Profile Solar Zenith Angle
; hour= Profile Time of Measurement
; refl= Profile Reflectivity
; doy= Day of Year
; oz=full data array before filtering (dimensions: [2000,total number profiles])
; set no_filter=1 to not apply filters
; filter_oz=full data array after filtering (dimensions: [2000,npts])
; v8ak is profile averaging kernel

; The V8.7 layer data are defined in 21 pressure intervals that are equally 
; spaced in log base 10 pressure coordinates.
; a=0.-indgen(21)*0.2
; pr=1013.25*10^a, where pr is the pressure of the lower boundary.
; There is no upper limit to the top level. 

; Profiles of ozone volume mixing ratio are provided on the following pressure levels
; Volume Mixing Ratio Levels 
;[  0.5,  0.7,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  7.0, $
;   10.0, 15.0, 20.0, 30.0, 40.0, 50.0]

; If nofilter is set, then ignore all algorithm filters.
if (keyword_set(nofilter)) then filter=0 else filter=1

; If no_saafilter is set, then do not filter OMPS NP for South Atlantic Anomaly
if (keyword_set(no_saafilter)) then saa_filter=0 else saa_filter=1

npar=2000 ;number of parameters per file

v8tot=9999.
lat=9999.
lon=9999.
sza=9999.
refl=9999.
hour=9999.
v8umk=9999.
v8vmr=9999.
v8ak=9999.
doy=9999.
npts=0
print,file

; FIRST READ FILES AND COUNT NUMBER OF RECORDS
close, 2
openr, 2, file,error=err, /swap_if_little_endian
temp=fstat(2)
size=temp.size
IF (err ne 0 or size eq 0) then begin
   print,  'This file is not available'
   print, 'error= '+string(err), string(size)
   close,2
 goto, over
ENDIF
close, 2

; reads f77 file as a unix file,later accounts for f77 4 byte
;  definitions at beginning and end of record
;  reduce rec to exclude 2 standard headers and the trailer record
; (# parameters + 2) * 4 bytes = # bytes
nbytes=(npar+2)*4
rec=size/nbytes-3

if (rec le 0) then begin
   print,  'This file is empty'
 goto, over
ENDIF

hdr=bytarr(nbytes)
oz=fltarr(npar+2,rec)

openr, 2, file,error=err, /swap_if_little_endian

readu, 2, hdr
readu, 2, hdr
readu,2,oz
oz=oz(1:npar,*)
close,2

; CHECK FOR NON FINITE VALUES
ff=finite(oz)
q=where(ff eq 0,nq)
if (nq ne 0) then oz(q)=9999.

; PICK OFF SPECIFIC VARIABLES
tozcod=reform(oz(36,*))
procod=reform(oz(481,*))
lsn=reform(oz(2,*))
v8tot=reform(oz(183,*))
surf_pr=reform(oz(67,*))
v8umk=reform(oz(142:162,*))
v8vmr=reform(oz(185:199,*))
v8ak=reform(oz(500:899,*))
lat=reform(oz(98,*))
lon=reform(oz(99,*))
sza=reform(oz(8,*))
refl=reform(oz(37,*))
hour=reform(oz(1,*)/3600.)
doy=reform(oz(4,*))
finres=reform(oz(225:234,*))
lngwav=fix(reform(oz(482,*)))
n=n_elements(procod)
filter_oz=oz

;  avgres is calculated as sum(abs(finres(1:lngwav)))/nwaves
avgres=fltarr(n)
for ii=0, n-1 do begin
   lw=lngwav(ii)-1
   if (lw gt 0) then begin
      tmp=reform(finres(1:lw,ii))
      q=where(abs(tmp) lt 99, nn)
      avgres(ii)=total(abs(tmp(q))/nn)
   endif else begin
      avgres(ii)=9999.
   endelse
endfor

; Filter the data
if (filter) then begin
   errst=strarr(n)
   for ii=0, n-1 do errst(ii)=string(procod(ii), format='(I3.3)')
   prof_code=fix(strmid(errst,2,1))
    asn_code=fix(strmid(errst,1,1))
    val_code=fix(strmid(errst,0,1))

   for ii=0, n-1 do errst(ii)=string(tozcod(ii), format='(I3.3)')
   toz_code=fix(strmid(errst,2,1))

   ; validation code=0 or 1
   ; If n14 data, include profile error=3 and total error=3
   r1=strpos(file, '.n14')
   if (r1 eq -1) then begin ; NOT N14
      if (asc eq 0) then $ ;descending orbit
          oo=where(avgres lt 0.3 and prof_code le 2 and (toz_code le 2 or toz_code eq 5) and asn_code eq 1, npts) else $
          oo=where(avgres lt 0.3 and prof_code le 2 and (toz_code le 2 or toz_code eq 5) and asn_code eq 0, npts)
      endif else begin   ; N14 include prof_code=3
      if (asc eq 0) then $ ;descending orbit
         oo=where(avgres lt 0.35 and prof_code le 3 and toz_code ne 4 and asn_code eq '1', npts) else $
         oo=where(avgres lt 0.35 and prof_code le 3 and toz_code ne 4 and asn_code eq '0', npts)
   endelse
endif else begin
   errst=strarr(n)
   for ii=0, n-1 do errst(ii)=string(procod(ii), format='(I3.3)')
   prof_code=fix(strmid(errst,2,1))
    asn_code=fix(strmid(errst,1,1))
    val_code=fix(strmid(errst,0,1))

   ; No Filter, check for Ascending/Descending and out of bounds Total Ozone 
   if (asc eq 0) then $ ;descending orbit
   oo=where(asn_code eq 1 and v8tot gt 50. and v8tot lt 700., npts) else $
   oo=where(asn_code eq 0 and v8tot gt 50. and v8tot lt 700., npts)
endelse

if (npts gt 1) then begin
   v8tot=v8tot(oo)
   v8umk=v8umk(*,oo)
   v8vmr=v8vmr(*,oo)
   v8ak=v8ak(*,oo)
   lat=lat(oo)
   lon=lon(oo)
   sza=sza(oo)
   refl=refl(oo)
   hour=hour(oo)
   doy=doy(oo)
   filter_oz=reform(filter_oz(*,oo))   
endif else begin
   v8tot=9999.
   lat=9999.
   lon=9999.
   sza=9999.
   refl=9999.
   hour=9999.
   v8umk=9999.
   v8vmr=9999.
   v8ak=9999.
   doy=9999.
   npts=0
   filter_oz=9999.
   goto, over
endelse

r1=strpos(file, '.npp')
if (r1 and saa_filter) then begin
   restore, '/misc/smh02/data/sbuvdat/progs/v8/saa_filter_omps.sav'
   for ii=44, 90 do begin
      saa=where(lat ge lolat(ii) and lat lt hilat(ii) and $
                lon ge saa_filter_lolon(ii) and lon le saa_filter_hilon(ii), nbad)
      if (nbad ne 0) then v8tot(saa)=9999.
   endfor

   oo=where(v8tot ne 9999., npts)
   print, n, npts
   if (npts gt 1) then begin
      v8tot=v8tot(oo)
      v8umk=v8umk(*,oo)
      v8vmr=v8vmr(*,oo)
      v8ak=v8ak(*,oo)
      lat=lat(oo)
      lon=lon(oo)
      sza=sza(oo)
      refl=refl(oo)
      hour=hour(oo)
      doy=doy(oo)
      filter_oz=reform(oz(*,oo))
   endif else begin
      v8tot=9999.
      lat=9999.
      lon=9999.
      sza=9999.
      refl=9999.
      hour=9999.
      v8umk=9999.
      v8vmr=9999.
      v8ak=9999.
      doy=9999.
      filter_oz=9999.
      goto, over
   endelse
endif

over:
end