pro get_ciret2_5min2_mininput, ciret2_fname, ciret2_iwc5, $
                   ciret2_dge5, ciret2_error_frac, ciret2_extinction

common five_min_data_input2, hrfrac5,jday5,height5,dbz5,nppts5,npts5,cbase5,ctop5, $
    vel5, width5,temp5, pressure5, n_layers5,bzl15,tzl15,btl15,ttl15,bzl25,tzl25,btl25,ttl25, $
    bzl35,tzl35,btl35,ttl35,bzl45,tzl45,btl45,ttl45,fluxclear5,solrat5,lwp5,vceil5,vap5


common minnis_tau_input, tau_minnis,tau_ir_minnis,jday_minnis

common iter_tau_results, iter_tau_iter, minnis_tau_iter, ba_iter, aa_iter, am_iter, bm_iter, jday_iter

; puts the iwc and re onto the time-height grid.

print, ciret2_fname
print, n_elements(ciret2_fname)



for jjjkkk=0,n_elements(ciret2_fname)-1 do begin

cdfid=ncdf_open(ciret2_fname[jjjkkk])   ; open the netcdf file

print, 'opened ',ciret2_fname[jjjkkk]

time_did=ncdf_dimid(cdfid,'Time')
height_did=ncdf_dimid(cdfid,'Height')

ncdf_diminq, cdfid, time_did, char_strng, num_times
ncdf_diminq, cdfid, height_did, char_strng, num_heights

; get the id's of the variables to be read

num_sec_id=ncdf_varid(cdfid,'base_time')
time_offset_id=ncdf_varid(cdfid,'time_offset')
N0s_id=ncdf_varid(cdfid,'N0s')
mlambda_id=ncdf_varid(cdfid,'mlambda')
iwc_id=ncdf_varid(cdfid,'Ice_Water_Content')
height_id=ncdf_varid(cdfid,'Height')
dbz_id=ncdf_varid(cdfid,'RadarReflectivity')
temp_id=ncdf_varid(cdfid,'Temperature')
dge_id=ncdf_varid(cdfid,'Dge')
vtau_id=ncdf_varid(cdfid, 'optical_depth_path')
vdq_id=ncdf_varid(cdfid, 'DopplerVelocity')



; get the data



ncdf_varget, cdfid, num_sec_id, num_sec
ncdf_varget, cdfid, time_offset_id, time_offset
ncdf_varget, cdfid, height_id, height

ijday_ihrfrac_fm_numsec, num_sec, time_offset, hrfrac, jday, yy, mm, dd, hh, mi, ss


ncdf_varget, cdfid, iwc_id, iwc_in
ncdf_varget, cdfid, N0s_id, N0s_in
ncdf_varget, cdfid, mlambda_id, mlambda_in
ncdf_varget, cdfid, dbz_id, dbz_in
ncdf_varget, cdfid, temp_id, temp
ncdf_varget, cdfid, dge_id, Dge_in
ncdf_varget, cdfid, vtau_id, vtau
ncdf_varget, cdfid, vdq_id, vdq_array
ncdf_varget, cdfid, temp_id, temperature



;if max(radrat) gt 0.15 then iwc_in[where(radrat gt 0.1)]=-9999.

ncdf_close, cdfid

       dge=fltarr(num_times,num_heights)
       extinction=fltarr(num_times,num_heights)
       dge[*,*]=-9999.
       extinction[*,*]=-9999.

; following loops attempts to solve equation 2.3 of Fu 1996.
count_it=0

; run through the time loop twice.  The first time iterate on the coeefecients if minnis tau exists.  If a mean set of
; coeffecients can be calculated, calcualte them and redo the retrieval with the tuned coeffecients.  If not, then
; use the mean coeffecients.

; check the iteration vectors and see if a good set of coeffecients can be derived for this date before iterating

aa_coef=-9999.

time_dif_index1=where(abs(jday_iter-jday[0]) lt 0.5)
if max(time_dif_index1) gt 0 then begin
    iter_tau_dif=(abs(iter_tau_iter[time_dif_index1]-minnis_tau_iter[time_dif_index1]))
    if n_elements(where(iter_tau_dif lt 0.1)) gt 1 then begin
    num1=n_elements(where(iter_tau_dif lt 0.1))
    aa_vector=aa_iter[time_dif_index1]
    am_vector=am_iter[time_dif_index1]
    ba_vector=ba_iter[time_dif_index1]
    bm_vector=bm_iter[time_dif_index1]

    aa_coef1=mean(aa_vector[where(iter_tau_dif lt 0.1)])
    am_coef1=mean(am_vector[where(iter_tau_dif lt 0.1)])
    ba_coef1=mean(ba_vector[where(iter_tau_dif lt 0.1)])
    bm_coef1=mean(bm_vector[where(iter_tau_dif lt 0.1)])
    endif else begin
       num1=0
    endelse

endif
time_dif_index2=where(abs(jday_iter-jday[n_elements(jday)-1]) lt 0.5)
if max(time_dif_index2) gt 0 then begin
    iter_tau_dif=(abs(iter_tau_iter[time_dif_index2]-minnis_tau_iter[time_dif_index2]))
    if n_elements(where(iter_tau_dif lt 0.1)) gt 1 then begin
    num2=n_elements(where(iter_tau_dif lt 0.1))
    aa_vector=aa_iter[time_dif_index2]
    am_vector=am_iter[time_dif_index2]
    ba_vector=ba_iter[time_dif_index2]
    bm_vector=bm_iter[time_dif_index2]

    aa_coef2=mean(aa_vector[where(iter_tau_dif lt 0.1)])
    am_coef2=mean(am_vector[where(iter_tau_dif lt 0.1)])
    ba_coef2=mean(ba_vector[where(iter_tau_dif lt 0.1)])
    bm_coef2=mean(bm_vector[where(iter_tau_dif lt 0.1)])
    endif else begin
       num2=0
    endelse

endif

if n_elements(time_dif_index1) gt 1 and n_elements(time_dif_index2) gt 1 and n_elements(aa_coef1) gt 0 and $
n_elements(aa_coef2) gt 0 then begin
    aa_coef=(1./(float(num1)+float(num2)))*((float(num1)*aa_coef1)+(float(num2)*aa_coef2))
    am_coef=(1./(float(num1)+float(num2)))*((float(num1)*am_coef1)+(float(num2)*am_coef2))
    ba_coef=(1./(float(num1)+float(num2)))*((float(num1)*ba_coef1)+(float(num2)*ba_coef2))
    bm_coef=(1./(float(num1)+float(num2)))*((float(num1)*bm_coef1)+(float(num2)*bm_coef2))
endif
if n_elements(time_dif_index1) gt 1 and n_elements(time_dif_index2) le 1 and n_elements(aa_coef1) gt 0 then begin
    aa_coef=(1./(float(num1)))*((float(num1)*aa_coef1))
    am_coef=(1./(float(num1)))*((float(num1)*am_coef1))
    ba_coef=(1./(float(num1)))*((float(num1)*ba_coef1))
    bm_coef=(1./(float(num1)))*((float(num1)*bm_coef1))
endif
if n_elements(time_dif_index1) le 1 and n_elements(time_dif_index2) gt 1 and n_elements(aa_coef2) gt 0 then begin
    aa_coef=(1./(float(num2)))*((float(num2)*aa_coef2))
    am_coef=(1./(float(num2)))*((float(num2)*am_coef2))
    ba_coef=(1./(float(num2)))*((float(num2)*ba_coef2))
    bm_coef=(1./(float(num2)))*((float(num2)*bm_coef2))
endif

if aa_coef lt 0 then begin

aa_vector=0. & ba_vector=0. & am_vector=0. & bm_vector=0.
for jj=0,num_times-1 do begin
;for jj=505,550 do begin

time_dif_index=where(abs(jday_minnis-jday[jj]) lt 20.d/1440.d)

if max(time_dif_index) gt 0 then begin

if tau_minnis[time_dif_index[0]] le 0. and (1.-exp(-tau_ir_minnis[time_dif_index[0]])) gt 0. and (1.-exp(-tau_ir_minnis[time_dif_index[0]])) lt 0.95 then tau_minnis[time_dif_index[0]]=2.*tau_ir_minnis[time_dif_index[0]]

if tau_minnis[time_dif_index[0]] gt 0. and tau_minnis[time_dif_index[0]] lt 10. then begin

ciret2_coef_iterate_minnis, tau_minnis[time_dif_index[0]], ((float(dbz_in[jj,*]))/100.), ((float(vdq_array[jj,*]))/10.), $
temperature[jj,*], height, min_tau_dif, min_aa, min_ba, min_am, min_bm



    if abs(min_tau_dif/tau_minnis[time_dif_index[0]]) lt 0.5 then begin

         if max(aa_vector) eq 0. then begin
          aa_vector=min_aa
          ba_vector=min_ba
          bm_vector=min_bm
          am_vector=min_am
         endif else begin
          aa_vector=[aa_vector,min_aa]
          ba_vector=[ba_vector,min_ba]
          bm_vector=[bm_vector,min_bm]
          am_vector=[am_vector,min_am]
         endelse

    endif  ; if min_tau_dif/tau_minnis lt 0.1 then begin

endif   ; if tau_minnis gt 0. and tau_minnis lt 100. then begin

endif   ; if max(time_dif_index) gt 0 then begin

if n_elements(min_tau_dif) gt 0 and max(time_dif_index) gt 0 then $
print, jj, num_times, min_ba, min_bm,abs(min_tau_dif)/tau_minnis[time_dif_index[0]]

endfor  ; for jj=0,num_times-1 do begin

if n_elements(where(ba_vector lt 2. and bm_vector lt 3.)) gt 1  then begin
    aa_coef=mean(aa_vector[where(ba_vector lt 2. and bm_vector lt 3.)])
    ba_coef=mean(ba_vector[where(ba_vector lt 2. and bm_vector lt 3.)])
    am_coef=mean(am_vector[where(ba_vector lt 2. and bm_vector lt 3.)])
    bm_coef=mean(bm_vector[where(ba_vector lt 2. and bm_vector lt 3.)])
endif else begin
    bm_coef=2.2515733
    am_coef=0.014726900
    aa_coef=0.20
    ba_coef=1.887
endelse

endif   ; if aa_coef lt 0 then begin

print, 'Using coeffecients ',am_coef, bm_coef, aa_coef, ba_coef

; now run through the time loop again and do the retrieval with the coeffecients

for jj=0,num_times-1 do begin
; run through the vertical loop and calculate the layer mean dbz, vdq, temperature
dbz_sum=0. & count_dbz=0 & vdq_sum=0. & temp_sum=0.
for k=0,num_heights-1 do begin
if mlambda_in[jj,k] gt -9000 then begin
    dbz_sum=dbz_sum+((float(dbz_in[jj,k]))/100.) & count_dbz=count_dbz+1
    vdq_sum=vdq_sum+((float(vdq_array[jj,k]))/10.)
    temp_sum=temp_sum+((float(temperature[jj,k])))
endif   ; if N0s_in[jj,k] ne 0 and N0s_in[jj,k] gt -10000 then begin
endfor  ; for k=0,num_heights-1 do begin

    if count_dbz gt 0. then begin
    dbz_layer=dbz_sum/float(count_dbz) & vdq_layer=vdq_sum/float(count_dbz) & temp_layer=temp_sum/float(count_dbz)
    endif
; run through the vertical loop again and calculate the extinction and microphysics

for k=0,num_heights-1 do begin

if mlambda_in[jj,k] gt -9000 and count_dbz gt 0 then begin

ciret2_itercoef_retrieval2, dbz_layer, vdq_layer, temp_layer, ((float(dbz_in[jj,k]))/100.), ((float(vdq_array[jj,k]))/10.), $
extinc, iwc, lambda, N0,am_coef,bm_coef,aa_coef,ba_coef

    extinction[jj,k]=extinc*100.     ; 100 converts from 1/cm to 1/m
    iwc_in[jj,k]=fix(1000.*(alog10(iwc*1.e6)))   ; 1.e6 converts from g/cm3 to g/m3
    mlambda_in[jj,k]=fix(1000.*(alog10(lambda*1.e-4)))    ; 1.e4 converts from 1/cm to 1/micron
    N0s_in[jj,k]=fix(1000.*(alog10(N0/10.)))  ; 10 converts to 1/liter/micron from 1/cm3/cm

endif   ; if N0s_in[jj,k] ne 0 and N0s_in[jj,k] gt -10000 then begin
endfor  ; for k=0,num_heights-1 do begin

time_dif_index=where(abs(jday_minnis-jday[jj]) lt 20.d/1440.d)
    if max(time_dif_index) gt 0 then begin
    
 ;  x=extinction[jj,*]
;if tau_minnis[time_dif_index[0]] gt 0. and tau_minnis[time_dif_index[0]] lt 10. and max(x) gt 0. then begin
;     print, total(x[where(x gt 0.)]*90.),tau_minnis[time_dif_index[0]],jj
;endif
   endif

    for k=0,num_heights-1 do begin

    if mlambda_in[jj,k] gt -9000 and count_dbz gt 0 then begin

         N0s=10.^((float(N0s_in[jj,k]))/1000.) & mlambda=10.^((float(mlambda_in[jj,k]))/1000.)

; first integral from size of 0 to 30 microns has aspect ratio of 1.
    asp=1.0 & lh=0. & rh=29.
    term3_30=(exp_integral_soln(3., -mlambda, lh, rh))
    term2_30=(exp_integral_soln(2., -mlambda, lh, rh))
    denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_30
    numer=asp*term3_30
    fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
    dge_30=fraction*(numer/denom)

    asp=0.6 & lh=30. & rh=79.
    term3_80=(exp_integral_soln(3., -mlambda, lh, rh))
    term2_80=(exp_integral_soln(2., -mlambda, lh, rh))
    denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_80
    numer=asp*term3_80
    fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
    dge_80=fraction*(numer/denom)


    asp=0.3 & lh=80. & rh=199.
    term3_200=(exp_integral_soln(3., -mlambda, lh, rh))
    term2_200=(exp_integral_soln(2., -mlambda, lh, rh))
    denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_200
    numer=asp*term3_200
    fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
    dge_200=fraction*(numer/denom)

    asp=0.2 & lh=200. & rh=499.
    term3_500=(exp_integral_soln(3., -mlambda, lh, rh))
    term2_500=(exp_integral_soln(2., -mlambda, lh, rh))
    denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_500
    numer=asp*term3_500
    fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
    dge_500=fraction*(numer/denom)

    asp=0.1 & lh=500. & rh=5000.
    term3_5000=(exp_integral_soln(3., -mlambda, lh, rh))
    term2_5000=(exp_integral_soln(2., -mlambda, lh, rh))
    denom=(1.+ ((asp)*((sqrt(3.))/4.)))*term2_5000
    numer=asp*term3_5000
    fraction=(exp_integral_soln(0., -mlambda, lh, rh))/(exp_integral_soln(0., -mlambda, 0., 5000.))
    dge_5000=fraction*(numer/denom)

    dge[jj,k]=dge_30+dge_80+dge_200+dge_500+dge_5000



;     dge[jj,k]=term1/(term2+term3)

    endif else begin

       if mlambda_in[jj,k] gt -9000 then begin
         dge[jj,k]=50.
       endif

    endelse

    endfor
    ;print, float(jj)/float(num_times)
       if count_it eq 100 then begin
         print, jj, num_times
         count_it=0
       endif
       count_it=count_it+1

endfor
print, 'done calculating dge'

; go through and calcuate dge


; loop through the observations, match times, then go through heights

total_ciret2_tau_vector=0.

for j=0,n_elements(hrfrac5)-1 do begin
    for k=0,n_elements(height5)-1 do begin

    time_index=where( (abs(jday-jday5[j])) le 0.0035 )
    z_index=where( (abs(height-height5[k])) lt 90. )

       max_val=-9999.
    for jjj=0,n_elements(time_index)-1 do begin
    for kkk=0,n_elements(z_index)-1 do begin
       if max(time_index) gt 0 and max(z_index) gt 0 then $
        max_val=max([max_val,dge[time_index[jjj],z_index[kkk]]])
    endfor
    endfor

    if max_val eq -9999. then z_index=where( (abs(height-height5[k])) le 200. )

       max_val=-9999.
    for jjj=0,n_elements(time_index)-1 do begin
    for kkk=0,n_elements(z_index)-1 do begin
       if max(time_index) gt 0 and max(z_index) gt 0 then $
         max_val=max([max_val,dge[time_index[jjj],z_index[kkk]]])
    endfor
    endfor


if max_val gt -1 then begin

    count=0. & mean_iwc=0. & mean_dge=0. & mean_extinc=0.
for jjj=0,n_elements(time_index)-1 do begin
    for kkk=0,n_elements(z_index)-1 do begin

    if dge[time_index[jjj],z_index[kkk]] gt -9999. and extinction[time_index[jjj],z_index[kkk]] gt -9999. then begin
       count=count+1.
;print, hrfrac[time_index[jjj]],height[z_index[kkk]],' in ',hrfrac5[j],height5[k]
       mean_iwc=mean_iwc+(10.^((float(iwc_in[time_index[jjj],z_index[kkk]]))/1000.))
       mean_dge=mean_dge+dge[time_index[jjj],z_index[kkk]]
       mean_extinc=mean_extinc+extinction[time_index[jjj],z_index[kkk]]
    endif


    endfor
endfor

       if count gt 0. then begin
    ciret2_iwc5[k,j]=mean_iwc/count & ciret2_error_frac[k,j]=0.04
    ciret2_dge5[k,j]=(mean_dge/count)*(1.e-6)
    ciret2_extinction[k,j]=mean_extinc/count
    ;print, ciret2_dge5[k,j], mean_dge, count
       endif else begin
    ciret2_iwc5[k,j]=-9999.
    ciret2_dge5[k,j]=-9999.
    ciret2_extinction[k,j]=-9999.
       endelse

endif else begin

    if ciret2_iwc5[k,j] le 0. then ciret2_iwc5[k,j]=-9999.
    if ciret2_iwc5[k,j] le 0. then ciret2_dge5[k,j]=-9999.
    if ciret2_iwc5[k,j] le 0. then ciret2_extinction[k,j]=-9999.

endelse

    endfor


; check this column against minnis tau.  If bad, don't use.
time_dif_index=where(abs(jday_minnis-jday5[j]) lt 20.d/1440.d)

if max(time_dif_index) ge 0 then begin

if tau_minnis[time_dif_index[0]] gt 0. and tau_minnis[time_dif_index[0]] lt 100. then begin

ciret2_extinction_vector=ciret2_extinction[*,j]
    if max(ciret2_extinction[*,j]) gt 0. then begin
total_ciret2_tau=total(90.*ciret2_extinction_vector[where(ciret2_extinction_vector gt 0.)])

if max(total_ciret2_tau_vector) eq 0. then begin
	if (total_ciret2_tau-tau_minnis[time_dif_index[0]])/tau_minnis[time_dif_index[0]] lt 10. then begin
    total_ciret2_tau_vector=total_ciret2_tau
    tau_minnis_vector=tau_minnis[time_dif_index[0]]
    	endif else begin
ciret2_iwc5[*,j]=-9999.
ciret2_dge5[*,j]=-9999.
ciret2_extinction[*,j]=-9999.
ciret2_error_frac[*,j]=-9999.
	endelse
endif else begin
	if (total_ciret2_tau-tau_minnis[time_dif_index[0]])/tau_minnis[time_dif_index[0]] lt 10. then begin
    total_ciret2_tau_vector=[total_ciret2_tau_vector,total_ciret2_tau]
    tau_minnis_vector=[tau_minnis_vector,tau_minnis[time_dif_index[0]]]
    	endif else begin
ciret2_iwc5[*,j]=-9999.
ciret2_dge5[*,j]=-9999.
ciret2_extinction[*,j]=-9999.
ciret2_error_frac[*,j]=-9999.
	endelse
endelse

endif  ; if max(ciret2_extinction[*,j]) gt 0. then begin

endif   ; if tau_minnis[time_dif_index[0]] gt 0. and tau_minnis[time_dif_index[0]] lt 100. then begin

endif   ; if max(time_dif_index) ge 0 then begin



endfor





endfor

; finally, compute comparison stats between minnis and ciret2 if possible

if max(total_ciret2_tau_vector) gt 0. then begin

r=correlate(total_ciret2_tau_vector,tau_minnis_vector)
bias=(abs((mean(total_ciret2_tau_vector))-(mean(tau_minnis_vector))))/(mean(tau_minnis_vector))

if bias gt 1.0 then begin
ciret2_iwc5[*,*]=-9999.
ciret2_dge5[*,*]=-9999.
ciret2_extinction[*,*]=-9999.
ciret2_error_frac[*,*]=-9999.

print, 'removing ciret2 due to bad minnis tau comparison ',r, bias

endif else begin    ; if r lt 0.5 and bias gt 0.5 then begin

print, 'keeping ciret2 due to good minnis tau comparison ',r, bias

endelse ;if r lt 0.5 and bias gt 0.5 then begin

endif

return
end