pro get_liq_data2, file_path, site, date, long_day, lwp, iwp

; IMPORTANT NOTE:
;
; The MWR is in MILLIMETERS even though the cdf file's units state
; centimeters... According to ARM, the correct units are millimeters.
; -Ben Foreback, July 6, 2009

multiplication_factor = double(90.)

;----------------------------------------------

if site eq 'sgp' then begin
  filename = file_path+'sgp.average.5min.'+date+'.000000.temp.cdf'
endif

spawn, 'gunzip '+filename+'*gz'

cdf_id = ncdf_open(filename)

;mwr_id = ncdf_varid(cdf_id,'vap')
mwr_id = ncdf_varid(cdf_id,'liq'); we want liq, not vap
lwc_id = ncdf_varid(cdf_id,'lwc_best_estimate')
iwc_id = ncdf_varid(cdf_id,'IWC_best_estimate')
base_time_id = ncdf_varid(cdf_id,'base_time')
time_offset_id = ncdf_varid(cdf_id,'time_offset')

ncdf_varget, cdf_id, mwr_id, mwr
;convert mwr from cm to g/m^2
;stop
;mwr = mwr * 1000. ; conversion from millimeters to grams/m^2
mwr = mwr * 10. ; conversion from centimeters to kg/m^2
;No conversion necessary for conversion from millimeters to kg/m^2

if lwc_id eq -1 or iwc_id eq -1 then return
ncdf_varget, cdf_id, lwc_id, lwc
ncdf_varget, cdf_id, iwc_id, iwc
;perform conversions
lwc = (10.^(float(lwc)/1000.))
iwc = (10.^(float(iwc)/1000.))
;stop
;result1 = where((lwc ge 8.e-8 and lwc lt 8.9e-8) or (lwc ge 9.e-9 and lwc lt 9.9e-9),count1)
result1 = where(lwc ge 7.e-8 and lwc lt 8.9e-8,count1)
result1b = where(lwc ge 8.e-9 and lwc lt 9.9e-9,count1b)
;lwc = lwc * multiplication_factor
if count1 gt 0 then lwc[result1] = -9999.
if count1b gt 0 then lwc[result1b] = 0.
;result2 = where((iwc ge 8.e-8 and iwc lt 8.9e-8) or (iwc ge 9.e-9 and iwc lt 9.9e-9),count2)
result2 = where(iwc ge 7.e-8 and iwc lt 8.9e-8,count2)
result2b = where(iwc ge 8.e-9 and iwc lt 9.9e-9,count2b)
;iwc = iwc * multiplication_factor
if count2 gt 0 then iwc[result2] = -9999.
if count2b gt 0 then iwc[result2b] = 0.
;stop
ncdf_varget, cdf_id, base_time_id, base_time
ncdf_varget, cdf_id, time_offset_id, time_offset


ncdf_close,cdf_id
epoch_zero = double(julday(1,1,1970,0,0,0))
;double_base_time = double(base_time)
long_day = (time_offset + base_time)/86400. + epoch_zero

;result_vap = where(mwr ge 0., count_vap)
;if count_ap le 0 then message, 'No valid values for MWR found!'

numtimes = n_elements(long_day)
lwp = make_array(numtimes,/double)
iwp = lwp
compare_flag = make_array(numtimes,/integer,value=1)

for i=0,numtimes-1 do begin
  liquid_column = lwc[*,i]
  ice_column = iwc[*,i]
  result1 = where(liquid_column ge 0., count1)
  result2 = where(ice_column ge 0., count2)
  if count1 gt 0 then begin
    valid_liquid_column = liquid_column[result1]
  endif else begin
    valid_liquid_column = 0.
  endelse
  if count2 gt 0 then begin
    valid_ice_column = ice_column[result2]
  endif else begin
    valid_ice_column = 0.
  endelse
  count2 = -1
  if count1 le 0 and count2 le 0 then begin
    compare_flag[i] = 0
    lwp[i] = -9999.
  endif else begin
    liq_total = total(valid_liquid_column)
    ice_total = total(valid_ice_column)
    ;lwp[i] = 90.*(liq_total + ice_total)
    lwp[i] = 90.*liq_total
    iwp[i] = 90.*ice_total
  endelse
endfor
;stop
result = where(mwr lt 0, count)
if count gt 0 then compare_flag[result] = 0


;stop
return
end