PRO rd_npp_lp_prof_daily,yr,mn,dy,first,inst,intoz,nprofs, $
                   lat,lon,hour

; This is a program to read in the OMPS_LP profiles and convert
; to partial column ozone in SBUV layers. 
; Read OMPS files collected into daily format

; Check for year format
if (yr lt 1900.) then yy=yr+1900 else yy=yr
syr=strmid(strcompress(fix(yy),/remove_all),2,2)
lyr=strcompress(fix(yy),/remove_all)

; The V8 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. 
nlevs=21
a=0.-indgen(nlevs)*0.2
hipress=1013.25*10^a
lopress=[hipress(1:20),0.0]

if (first eq 1) then begin
   inst=0
   oops1:
   PRINT, 'Choose Instrument:   (1) OMPS_LP v26 01/2012-present (Feb 2023)'
   READ,inst
   if (inst gt 1) then begin
      print, 'Try Again'
      goto, oops1
   endif
   first=0
endif

mth=string(mn,format='(I2.2)')
day=string(dy,format='(I2.2)')

CASE 1 of
   (inst eq 1): begin
      name='/science/npp/omps_lp/data/L2_v2.6/Y'
       filename=strcompress(name+lyr+'/M'+mth+ $ +
              '/OMPS-NPP_LP-L2-O3-DAILY_v2.6_'+lyr+'m'+mth+day+'.h5',/remove_all)
   end
ENDCASE

thalt=indgen(61)+0.5
intoz=9999.
nprofs=0
lat=9999.
lon=9999.
hour=9999.

print,filename
exist=FILE_TEST(filename)
if (exist eq 0) then goto, over

; OPEN FILE AND READ DATA
r=h5_parse(filename,/read_data)

hour=r.GeolocationFields.SecondsInDay._DATA
lat= r.GeolocationFields.Latitude._DATA
lon= r.GeolocationFields.Longitude._DATA
asc= r.GeolocationFields.AscendingDescendingFlag._DATA
saa= r.GeolocationFields.SwathLevelQualityFlags._DATA

pmcflag=  r.DataFields.ASI_PMCFlag._DATA 
qmvflag=  r.DataFields.QMV._DATA
o3quality=r.DataFields.O3Quality._DATA
o3status= r.DataFields.O3Status._DATA
o3conv=   r.DataFields.O3Convergence._DATA
o3value=  r.DataFields.O3Value._DATA

pressure=r.AncillaryData.Pressure._DATA
temp=   r.AncillaryData.Temperature._DATA
tropalt=r.AncillaryData.TropopauseAltitude._DATA

nprofs=n_elements(hour)

; Filter the data (recommended by Natalya, email 2/27/2023)
; I also included ascending only, saa=0
q=where(o3value lt 0., nq)
if (nq ne 0) then o3value(q)=9999.
q=where(o3conv lt 0 or o3conv gt 10 or $
        o3status lt 2 or o3status gt 7 or $
        o3quality ne 0 or $
        pmcflag ne 0 or $
        qmvflag ne 0 or $
        asc ne 0 or $
        saa ne 0, nq)
if (nq ne 0) then begin
   o3value(*,q)=9999.
   lat(q)=9999.
   lon(q)=9999.
   hour(q)=9999.
endif

; Convert Number Density to VMR, then integrate to get layer data
intoz=fltarr(nlevs,nprofs)+9999.
for iprof=0, nprofs-1 do begin
   tmpnd=reform(o3value(*,iprof))
   q=where(tmpnd ne 9999., nq)
   if (nq eq 0) then goto, next_prof

   tmpnd=reform(tmpnd(q))
   tmppr=reform(pressure(q,iprof))
   tmpt=reform(temp(q,iprof))   

   tmpvmr=1.e6*(1.38e-19*tmpt)/tmppr*tmpnd

; REVERSE PROFILES AND COMPUTE AVERAGE OZONE BETWEEN LEVELS
   press_bounds=reverse(tmppr)
   logpr=alog(press_bounds)
   midpr=(logpr(0:nq-2)+logpr(1:nq-1))/2.
   aoz=interpol(reverse(tmpvmr), logpr, midpr)
;stop
; INTEGRATE TO SBUV LAYERS 
   bpr=hipress
   tpr=lopress
   w=where(bpr lt press_bounds(nq-1) and $
           tpr gt press_bounds(0), nw)

   vmrtodu_int_cumoz, press_bounds, bpr(w), tpr(w), aoz, tot
      
   intoz(w,iprof)=tot
next_prof:
endfor

over:	   
end