pro read_newoutput_alldays_no2,monthofinterest,aot440m_min,filename,ddfittolerance,sza_min,sza_max,$
day_plot,month_plot,year_plot, $
hour,min,sec,day,month,dayofyear,year,$
  cv,cvf,cvc,rv,rvf,rvc,s,sf,sc,$
aot1020m,aot870m,aot670m,aot500m,aot440m,aot380m,aot340m,$
aot440,aot670,aot870,aot1020,w440,w670,ri440,ri670,ri870,ri1020,rr440,rr670,rr870,rr1020,$
sza_mfrsr,hour_mfrsr,min_mfrsr,sec_mfrsr, ozone,pres_md12, m_tau368,m_w368,a_tau368,dd_calc368,dd_mfrsr368,niter368, i1,count,m_total,m_dif,ct_meas,n_real,k_real,no2,asym368


;		01_18.13: process new format outout: add asym and Jacobian columns to the end
;		09_07.12: process new format outout                   
;		07_28.05: process measurements with NO2               
;		04_15.04: turn-off day selection: get all days        
;		04_09.04: read new output from the fitting program

openr,1,filename                    
nrows = file_lines(filename)-1
x='' & a=dblarr(70,nrows)
for i=0,0 do readf,1,x	& print,x ; read header
readf,1,a				; read aeronet dump
close,1

     day =reform(a(0,*))
    month=reform(a(1,*))
     year=reform(a(2,*))
     hour=reform(a(3,*))
     min =reform(a(4,*))
     sec =reform(a(5,*))
dayofyear=reform(a(6,*))

aot1020m=reform(a(7,*))
aot870m =reform(a(8,*))
aot670m =reform(a(9,*))
aot500m =reform(a(10,*))
aot440m =reform(a(11,*))
aot380m =reform(a(12,*))
aot340m =reform(a(13,*))

aot440  =reform(a(14,*))
aot670  =reform(a(15,*))
aot870  =reform(a(16,*))
aot1020 =reform(a(17,*))

w440 	=reform(a(18,*))
w670 	=reform(a(19,*))
w870 	=reform(a(20,*))
w1020	=reform(a(21,*))

rr440	=reform(a(22,*))
rr670	=reform(a(23,*))
rr870	=reform(a(24,*))
rr1020	=reform(a(25,*))

ri440	=reform(a(26,*))
ri670	=reform(a(27,*))
ri870	=reform(a(28,*))
ri1020	=reform(a(29,*))


cvf	=reform(a(30,*))
cvc	=reform(a(31,*))
rvf	=reform(a(32,*))
rvc	=reform(a(33,*))
sf	=reform(a(34,*))
sc	=reform(a(35,*))

hour_mfrsr =reform(a(36,*))
min_mfrsr  =reform(a(37,*))
sec_mfrsr  =reform(a(38,*))
timeofday_mfrsr=double(hour_mfrsr + min_mfrsr/60. + sec_mfrsr/3600.)

sza_mfrsr  =reform(a(39,*))
ozone      =reform(a(40,*))	; True measured ozone
no2        =reform(a(41,*))	; True measured No2  
pres_md12  =reform(a(42,*))

n_real     =reform(a(44,*))
k_real     =reform(a(45,*))

m_tau368   =reform(a(46,*))	; MFRSR AOT
m_w368     =reform(a(47,*))
a_tau368   =reform(a(48,*))	; aeronet extrapolated AOT

tau_ray    =reform(a(49,*))
tau_oz     =reform(a(50,*))
l_rad      =reform(a(51,*))
m_total    =reform(a(52,*))		;  total transmittance
m_dif      =reform(a(53,*))		;  diffuse transmittance
m_dirn     =reform(a(54,*))		;  direct normal transmittance
m_v0raw    =reform(a(55,*))		;  V0 raw (measured, not adjusted

e_dir_ray368	=reform(a(58,*))
e_dir_rat368	=reform(a(59,*))  &  e_dir_aer=e_dir_ray368*e_dir_rat368

e_dif_ray368	=reform(a(60,*))
e_dif_rat368	=reform(a(61,*))  &  e_dif_aer=e_dif_ray368*e_dif_rat368

e_tot_ray=e_dir_ray368 + e_dif_ray368
e_tot_aer=e_dir_aer    + e_dif_aer
ct_calc = e_tot_aer/e_tot_ray

ct_meas=m_total/e_tot_ray	; ct_meas_d=(m_dirn*cos(!pi/180.*sza_mfrsr)+m_dif)/e_tot_ray

dd_calc368	=reform(a(62,*))
dd_mfrsr368	=reform(a(63,*))
niter368  	=reform(a(64,*))

asym368  	=reform(a(67,*))

;       FILTERING THE DATA


i1=where(                      $
                (sza_mfrsr    ge sza_min and sza_mfrsr    le sza_max) $    
        and     (aot440m  ge  aot440m_min) $    ;  select only days with AOT440m > aot440m_min
;       and     (month    eq monthofinterest) $    ;  select only data for specific month       
;       and     (year     eq 2005           ) $    ;  select only data for specific year        
;         and     (day eq day_plot and month eq month_plot and year eq year_plot) $     
;       and     (m_tau368 gt  0.0       )$                          
;       and     (e_dir_ray368 gt 0.0 ) $
;       and     (niter368 le 6 ) $
        and     (abs(dd_calc368 - dd_mfrsr368) lt  ddfittolerance       )$                          
        and     (k_real ge 0.0  and k_real le 1.0     )$    ;  This filter removes non retrievals
                ,count )

print,'read_newoutput:  Number of points i1= ',count


return
end