;
;  Plot a monthly summary
;
pro histogram_average_1d_rad;,site,searchdate,avg_int
site='sgp'
searchdate='200004'
avg_int='1hr'

;
;  common containing the monthly averaged data
;
common monthly_averaged_data, base_time,time_offset,height,reflectivity,$
  reflectp, reflectc, cbase, ctop, velocity, clutter,maxclutter,spectrum,$
  lwp,vceil,numtimes,numheights,allclutter,goodclutter,temp,press,$
  vap,dir,dsol,sfpres,numlayers,twsflux,zenith,solarratio,sfcalbedo,$
  column_cfrac

;
;  common containing the monthly averaged data
;
common monthly_averaged_data_rad, base_time_rad,time_offset_rad,height_rad,$
  numtimes_rad,numheights_rad,$
  IWC_best_estimate,re_ice_best_estimate,ice_micro_error_estimate, $
  iwc_source,ice_effective_radius_source,ice_tau_solar, $
  ice_omega_solar,ice_g_solar,ice_tau_IR,ice_omega_IR,ice_g_IR, $
  ice_RadiationParameterization_source_solar,ice_RadiationParameterization_source_ir, $
  lwc_best_estimate,re_liquid_best_estimate,liquid_micro_error_estimate, $
  lwc_source,liquid_effective_radius_source,liquid_tau_solar,liquid_omega_solar, $
  liquid_g_solar,liquid_tau_IR,liquid_omega_IR,liquid_g_IR, $
  liquid_RadiationParameterization_source_solar,$
  liquid_RadiationParameterization_source_ir,$   
  ice_fraction_CCM3,solar_flux_down_clear,solar_flux_up_clear, $
  IR_flux_up_clear,IR_flux_down_clear,downwelling_solar_diffuse_clear, $
  downwelling_solar_direct_clear,TOA_solar_up_clear,TOA_solar_down_clear, $    
  TOA_IR_up_clear,solar_flux_down,solar_flux_up,IR_flux_up,IR_flux_down, $
  downwelling_solar_diffuse,downwelling_solar_direct,TOA_solar_up,TOA_solar_down, $
  TOA_IR_up
;
;  Read in the monthly averaged data
;
get_average_data,site,searchdate,avg_int
get_average_data_rad,site,searchdate,avg_int
;
;  Convert heights to km
;
height=height/1000.0
;
;  Number of histograms
;
num_histograms=4
;
;  Select which data to plot
;
dataname=strarr(num_histograms)
dataname=['Ice Water Path',$
          'Liquid Water Path',$
          'Optical Depth',$
          'Column Cloud Fraction']

datatag=strarr(num_histograms)
datatag=['iwp',$
         'lwp',$
         'optical_depth',$
         'column_cfrac']

cb_title=strarr(num_histograms)
cb_title=['g/m^2',$
          'g/m^2',$
          'Optical depth',$
          'Percent']
;
;  Create some histograms
;
for datanum=0,num_histograms-1 do begin
;for datanum=3,3 do begin

;  Get the data

;***************
;  Ice Water Path
;***************
  if datanum eq 0 then begin

;  Array to hold the ice water path data
    data=make_array(numtimes,/float,value=0)

    for i=0,numtimes-1 do begin
      ;print,i
      temp_IWC=IWC_best_estimate[i,*]
      result=where(temp_IWC gt -8888,count)
      if count gt 0 then begin
        data[i]=total(temp_IWC[result]*90.0)
        ;print,temp_IWC[result],'IWC',data[i]
      endif

;  Determine if there is missing data
      result=where(temp_IWC ne -8888,count)
      if count eq 0 then data[i]=-8888

    endfor

;  Log it for ease of plotting
;    result=where(data gt -8888,count)
;    if count gt 0 then data[result]=alog10(data[result])

;  Histogram values 
;    result=where(data gt -8888,count)
    result=where(data gt 0,count)
    binsize=0.1
  endif

;***************
;  Liquid Water Path
;***************
  if datanum eq 1 then begin

;  Array to hold the liquid water path data
    data=make_array(numtimes,/float,value=0)

    for i=0,numtimes-1 do begin
      ;print,i
      temp_lwc=lwc_best_estimate[i,*]
      result=where(temp_lwc gt -8888,count)
      if count gt 0 then begin
        data[i]=total(temp_lwc[result]*90.0)
        ;print,temp_lwc[result],'lwp',data[i]
      endif

;  Determine if there is missing data
      result=where(temp_lwc ne -8888,count)
      if count eq 0 then data[i]=-8888

    endfor

;  Log it for ease of plotting
;    result=where(data gt -8888,count)
;    if count gt 0 then data[result]=alog10(data[result])

;  Histogram values 
;    result=where(data gt -8888,count)
    result=where(data gt 0,count)
    binsize=0.1

  endif 

;***************
;  Optical Depth
;***************
  if datanum eq 2 then begin

;  Array to hold the optical depth data
    data=make_array(numtimes,/float,value=0)

    for i=0,numtimes-1 do begin
      ;print,i
      temp_ice=ice_tau_solar[i,*]
      result=where(temp_ice gt -8888,count_ice)
      if count_ice gt 0 then begin
        data[i]=total(temp_ice[result])
        ;print,temp_ice[result],'ice',data[i]
      endif
 
      temp_liq=liquid_tau_solar[i,*]
      result=where(temp_liq gt -8888,count_liq)
      if count_liq gt 0 then begin
        data[i]=data[i]+total(temp_liq[result])
        ;print,temp_liq[result],'liquid',data[i]
      endif

;  Determine if there is missing data
      result=where(temp_ice ne -8888,count_ice)
      result=where(temp_liq ne -8888,count_liq)
      if count_ice eq 0 and count_liq eq 0 then data[i]=-8888

    endfor

;  Log it for ease of plotting
;    result=where(data gt -8888,count)
;    if count gt 0 then data[result]=alog10(data[result])

;  Histogram values 
;    result=where(data gt -8888,count)
    result=where(data gt 0,count)
    binsize=0.1

  endif

  if datanum eq 3 then begin
    data=column_cfrac

;  Histogram values 
;    result=where(data gt -8888,count)
    result=where(data ge 0,count)
    binsize=1.0

  endif

;********************
;  Histogram the data
;********************

  hist=histogram(data[result],binsize=binsize,omax=dmax,omin=dmin)
  bins=(findgen(n_elements(hist))*binsize)+dmin
  print,dmax,dmin

;  Create the plot
  window,0,/pixmap
  loadct,39
  !p.background=!d.n_colors-1
  !p.charsize=1.5
  !p.multi=[1,1,1]
  !x.style=1
  !y.style=1
  pos_vect=[0.1,0.1,0.9,0.84]
  xaxis='log'  ;lin,log

  if xaxis eq 'lin' then begin
  plot,bins,hist,$
    yrange=[min(hist)-1,max(hist)+1],$
    xrange=[0,dmax],$
    psym=10,$
    xtitle=cb_title[datanum],ytitle='Density per Bin',$
    color=0,position=pos_vect
  endif else if xaxis eq 'log' then begin

  plot,bins,hist,$
    yrange=[min(hist)-1,max(hist)+1],$
    psym=10,$
    xtitle=cb_title[datanum],ytitle='Density per Bin',$
    color=0,position=pos_vect,/xlog
  endif
;
;  Add a title
;
  smonth=strmid(searchdate,4,2)
  syear=strmid(searchdate,0,4)
  sday='01'
  if smonth eq 1 then smon='Jan'
  if smonth eq 2 then smon='Feb'
  if smonth eq 3 then smon='Mar'
  if smonth eq 4 then smon='Apr'
  if smonth eq 5 then smon='May'
  if smonth eq 6 then smon='Jun'
  if smonth eq 7 then smon='Jul'
  if smonth eq 8 then smon='Aug'
  if smonth eq 9 then smon='Sep'
  if smonth eq 10 then smon='Oct'
  if smonth eq 11 then smon='Nov'
  if smonth eq 12 then smon='Dec'
  if site eq 'sgp' then begin
    capsite='SGP'
    datastream='average.'+avg_int
  endif else if site eq 'nsa' then begin
    capsite='NSA'
    datastream='average.'+avg_int
  endif else if site eq 'twpc1' then begin
    capsite='TWP  Manus'
    datastream='average.'+avg_int
  endif else if site eq 'twpc2' then begin
    capsite='TWP  Nauru'
    datastream='average.'+avg_int
  endif
  syr=strcompress(string(syear),/remove_all)
  xyouts,0.5,0.95,align=0.5,charsize=1.7,/normal,color=0,$
  'Histogram of Hourly Averages!C'+ $
  capsite+' - '+smon+' '+syr + '!C' + $
  dataname[datanum]
;
;  Create the gif file
;
  strmon=string(smonth)
  if smonth lt 10 then strmon='0'+strmon
  strday=string(sday)
  if sday lt 10 then strday='0'+strday
  curdate=strcompress(syr+strmon+strday,/remove_all)

  image_dir='/data/mace4/arm_data/temp_image/'
  image_name=site+'.'+datastream+'.'+curdate+'.'+ $
    datatag[datanum]+'.hist.'+xaxis
  print,image_dir+image_name
  write_gif,image_dir+image_name+'.gif',tvrd( )
;
;  Remove the .png file
;
  unix_command='rm '+image_dir+image_name+'.png'
  spawn,unix_command

endfor  ;loop through the data to histogram

end