;----------------------------------------------------------------------------
;LIDAR CODE
;----------------------------------------------------------------------------
pro rd_clds_ascii, done
;----------------------------------------------------------------------------
;Purpose:	Reads the lidar CLD files (GUI to select which one).
;
;Returns:	depols - depolarization ratio
;		dtstrt - the date of the file
;		tmstrt - the times of each 10-min averaged beam
;----------------------------------------------------------------------------

common cld_data, beamnum, depols, icldbs, icldtp, dtstrt, tmstrt, max_layers,$
		 n_layers, num_recs, diffval

cd,'d:/lidar/CLD',current=odir

;-----Script to call up widget for selecting the input file------

done = 0
datafile = ' '
datafile = DIALOG_PICKFILE(/READ)
sl = strlen(datafile)

;-----Initialize the read variables

year = 0
month = 0
day = 0
hour = 0
minute = 0
second = 0
nlayers = 0
max_layers = 0

;-----Initialize the counter variable

j=0

;-----Initialize the array variables

dtstrt = strarr(25000)
tmstrt = strarr(25000)
icldbs = fltarr(25000,20)
icldtp = fltarr(25000,20)
diffval = fltarr(25000,20)
depols = fltarr(25000,20)
beamnum = intarr(25000)
n_layers = intarr(25000)

;-----Open the data file and read the data

openr,lun,datafile,/get_lun

while (not eof(lun)) do begin

;-----Read the date, time, and the number of layers for each time

readf, lun, year, month, day, hour, minute, second, nlayers
  if (year ge 2000) then begin
    year = year - 2000
  endif else begin
    year = year - 1900
  endelse

  beamnum(j) = j

  years = strtrim(string(year),2)
  months = strtrim(string(month),2)
  if (strlen(months) eq 1) then months='0'+months
  days = strtrim(string(day),2)
  if (strlen(days) eq 1) then days='0'+days
  ltime = months+'/'+days+'/'+years
  dtstrt(j) = ltime

  hours = strtrim(string(hour),2)
  minutes = strtrim(string(minute),2)
  if (strlen(minutes) eq 1) then minutes='0'+minutes
  seconds = strtrim(string(second),2)
  if (strlen(seconds) eq 1) then seconds='0'+seconds
  ltime = hours+':'+minutes+':'+seconds
  tmstrt(j) = ltime

  n_layers(j) = nlayers
  if (nlayers gt max_layers) then max_layers = nlayers

  for i=0, 20-1 do begin
      icldbs(j, i) = -999
      icldtp(j, i) = -999
      depols(j, i) = -999
      diffval(j,i) = -999
  endfor

;----Read each layer data (if there are no clouds then set these values)
  if (nlayers eq 0) then begin
      icldbs(j,0) = -500.00
      icldtp(j,0) = -999
      diffval(j,0) = -500.00
      depols(j,0) = -500.00
  endif else begin
    for i=0, nlayers-1 do begin
      readf, lun, cldbase, cldtop, intgrt_val, avg_depol, depol_int
      icldbs(j,i) = cldbase
      icldtp(j,i) = cldtop
      diffval(j,i) = cldtop - cldbase
      depols(j,i) = depol_int
    endfor
  endelse

  j = j + 1

endwhile

close, lun
free_lun, lun

;-----Down counter by one for correct sizing
num_recs = j-1

;-----Reform the Arrays for correct sizing
beamnum = beamnum(0:num_recs)
tmstrt = tmstrt(0:num_recs)
dtstrt = dtstrt(0:num_recs)
icldbs=icldbs(0:num_recs, 0:max_layers-1)
icldtp=icldtp(0:num_recs, 0:max_layers-1)
diffval=diffval(0:num_recs, 0:max_layers-1)
depols=depols(0:num_recs, 0:max_layers-1)

cd,odir

end
;---------------------------------------------------------------------
;END LIDAR CODE
;---------------------------------------------------------------------


;---------------------------------------------------------------------
;BEGIN THE MICROWAVE RADIOMETER CODE
;---------------------------------------------------------------------
pro rd_mwr,ldate,ltime
;---------------------------------------------------------------------
;Purpose:	Read in mwr data and interpolate to the lidar times
;
;Returns:	water - total liquid water column (cm)
;		vapor - total water vapor column (cm)
;
;I/O format:	rd_mwr,dtstrt[0],tmstrt
;---------------------------------------------------------------------

common mwr_data,water,vapor,ldtime,mwfile,mwvap,mwwat,mwtime

mwrawdir='d:/shebaar/microwave/raw'

;create the search datestring  (*980428*)
cd,mwrawdir,current=odir
mwarg='*'+strmid(ldate,6,2)+strmid(ldate,0,2)+strmid(ldate,3,2)+'*'
file=findfile(mwarg)
mwfile=file[0]
;if mwfile eq ' ' then print,'No MWR data' & return

;convert the lidar time strings into numeric values (ltime)
nlt=n_elements(ltime)
ldtime=fltarr(nlt)
for i=0,nlt-1 do begin
  if strlen(ltime[i]) eq 7 then begin
    hr=float('0'+strmid(ltime[i],0,1))
    min=float(strmid(ltime[i],2,2))/60.
  endif else begin
    hr=float(strmid(ltime[i],0,2))
    min=float(strmid(ltime[i],3,2))/60.
  endelse
  ldtime[i]=hr+min
endfor

;get data from the netCDF file

;mwfid=ncdf_open(mwfile)
;ncdf_varget,mwfid,0,mwbase_time		;sec since 1970-01-01 00:00:00
;ncdf_varget,mwfid,1,mwtime_off		;sec since base_time
;ncdf_varget,mwfid,4,mwvap		;Total water vapor (cm)
;ncdf_varget,mwfid,5,mwwat		;Total liquid water (cm)
;ncdf_close,mwfid

;Convert the time data into a time axis
;----------------------------------------
;base_date=julday(1,1,1970,0,0,0)
;base=jul_to_dt(base_date)
;!dt_base=base
;starttime=sec_to_dt(mwbase_time)
;mwtime=(starttime.second + mwtime_off)/3600.
;mwjtime=float(floor((starttime.julian - julday(1,1,starttime.year,0,0)+1)))$
;		+ mwtime/24
restore,mwfile
mwvap=vapor
mwwat=water
mwtime=(mwjtime-float(fix(mwjtime)))*24

;Interpolate the water & vapor data to lidar times (assure non-negative)
;-----------------------------------------------------------------------
vapor=interpol(mwvap,mwtime,ldtime)
water=interpol(mwwat,mwtime,ldtime)
water=water>0
vapor=vapor>0

cd,odir
end
;---------------------------------------------------------------------
;END THE MICROWAVE RADIOMETER CODE
;---------------------------------------------------------------------


;---------------------------------------------------------------------
;This is the MAIN LEVEL program
;---------------------------------------------------------------------

;(1) Run rd_clds_ascii.pro
;(2) Run rd_mwr.pro using the dtstrt[0] variable from lidar
;(3) Do a comparison of the two files

;These are the lidar common variables
;--------------------------------------
common cld_data, beamnum, depols, icldbs, icldtp, dtstrt, tmstrt, max_layers,$
		 n_layers, num_recs, diffval
common mwr_data,water,vapor,ldtime,mwfile,mwvap,mwwat,mwtime

rd_clds_ascii,0

rd_mwr,dtstrt[0],tmstrt

;set_plot,'ps'
;window,0
sd=size(depols,/dimensions)
!p.multi=[0,1,2]
plot,ldtime,depols[*,0],min_value=0,psym=1,ytitle='Lidar Depolarization Ratio',$
	title='Lidar Depol - MWRadiometer Comparison, '+dtstrt[0],xrange=[0,24],$
	xstyle=1
if n_elements(sd) ge 2 then begin
  for i=1,sd[1]-1 do begin
    ii=i+1
    if ii ge 3 then ii=ii+1
    if ii ge 8 then ii=7
    oplot,ldtime,depols[*,i],min_value=0,psym=ii
  endfor
endif
;plot,ldtime,water,xtitle='Time',ytitle='Column Liquid Water [cm]',xrange=[0,24]
plot,mwtime,mwwat,xtitle='Time',ytitle='Column Liquid Water [cm]',xrange=[0,24],$
	xstyle=1,min_value=0
;device,/close
;spawn,'lpr idl.ps'
;set_plot,'x'

!p.multi=0
;window,1
;t_depol=depols
;t_depol=t_depol>0
;plot,water,t_depol[*,0],psym=1,xtitle='Column Liquid Water [cm]',$
;	ytitle='Lidar Depolarization Ratio',$
;	title='Liquid Water vs. Depol Ratio, '+dtstrt[0]
;if n_elements(sd) ge 2 then begin
;  for i=1,sd[1]-1 do begin
;    ii=i+1
;    if ii ge 3 then ii=ii+1
;    if ii ge 8 then ii=7
;    oplot,water,t_depol[*,i],psym=ii
;  endfor
;endif

;print,'1st Layer Correlation: '+string(correlate(depols[*,0],water))
;print,'Raw Multiple Correlation: '+string(m_correlate(transpose(depols),water))

;iwhere=where(t_depol gt 0.16)
;if iwhere[0] ne -1 then t_depol[iwhere]=0
;nbeams=n_elements(t_depol[*,0])
;num=fltarr(nbeams)
;for i=0,nbeams-1 do num[i]=n_elements(where(t_depol[i,*] gt 0))
;if n_elements(sd) ge 2 then s_depol=total(t_depol,2)/num else s_depol=t_depol
;print,'Combined Layer Correlation: '+string(correlate(s_depol,water))



end
;--------------------------------------------------------------------