;
;  Plot a monthly summary
;
pro histogram_average_2d;,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

;**********************************************************************
;  Hour of the day vs column cloud fraction
;**********************************************************************

;
;  Strings
;
dataname='Hour of the Day vs. Column Cloud Fraction'
datatag='hr_column_cfrac'
cb_title='%'
;
;  Calculate the file starting and ending date and time
;
standard_date_time,base_time+time_offset(0),syear,smonth,sday,shour,smin,ssec
print,'file start',syear,smonth,sday,shour,smin,ssec
standard_date_time,base_time+time_offset(numtimes-1),eyear,emonth,eday,ehour,emin,esec
print,'file end',eyear,emonth,eday,ehour,emin,esec
;
;  Calculate the standard date arrays
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Set up hour bins with mid points 0.5,1.5,2.5,3.5,4.5 hrs
;
hour_bins=0.5
dhour_bins=0.5
end_hour_bins=23.5
while hour_bins[n_elements(hour_bins)-1] lt end_hour_bins $
 do hour_bins=[hour_bins,hour_bins[n_elements(hour_bins)-1]+(2.*dhour_bins)]
;
;  Set up cfrac bins with mid points 5,15,25,35,45,55,65, %
;
cfrac_bins=5.0
dcfrac_bins=5.0
end_cfrac_bins=105.0
while cfrac_bins[n_elements(cfrac_bins)-1] lt end_cfrac_bins $
  do cfrac_bins=[cfrac_bins,cfrac_bins[n_elements(cfrac_bins)-1]+(2.*dcfrac_bins)]
;
;  Set up the 2D histogram arrays
;
hour_cfrac_freq= $ 
  make_array(n_elements(hour_bins),n_elements(cfrac_bins),/float,value=0) 
;  
; Build the 2D histogram for temp-thick
;
for j=0,n_elements(hour_bins)-1 do begin
  for k=0,n_elements(cfrac_bins)-1 do begin
    result=where(hh ge hour_bins[j]-dhour_bins and $
                 hh lt hour_bins[j]+dhour_bins and $
		 column_cfrac ge cfrac_bins[k]-dcfrac_bins and $
                 column_cfrac lt cfrac_bins[k]+dcfrac_bins,count) 

    hour_cfrac_freq[j,k]=count

    if hour_bins[j]-dhour_bins eq 5 then begin
    print,hour_bins[j]-dhour_bins,hour_bins[j]+dhour_bins,$
          cfrac_bins[k]-dcfrac_bins,cfrac_bins[k]+dcfrac_bins,hour_cfrac_freq[j,k]

    endif

  endfor  ;end of loop through thick bins
endfor  ;end of loop through temp bins

;  Change the counts array to a frequency array
;hour_cfrac_freq=float(hour_cfrac_freq)/float(total(hour_cfrac_freq))

;  Frequency array per hour
;for j=0,n_elements(hour_bins)-1 do begin
;  hour_cfrac_freq[j,*]=float(hour_cfrac_freq[j,*])/float(total(hour_cfrac_freq[j,*]))
;endfor



;******************
;  Set up 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
!x.ticklen=-0.02
!y.ticklen=-0.02
pos_vect=[0.15,0.15,0.7,0.8]
cb1=[0.9,0.15,0.95,0.8]

dmax=max(hour_cfrac_freq)
dmin=min(hour_cfrac_freq)

contour,hour_cfrac_freq,hour_bins,cfrac_bins,/nodata,$
  xstyle=4,ystyle=4,position=pos_vect

px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Scale and create the image
image=bytscl(hour_cfrac_freq,max=dmax,min=dmin,top=254)
result=where(hour_cfrac_freq eq 0)
image[result]=255
tv,congrid(image,sx,sy),px[0],py[0]
 
;  Set up the axes
contour,hour_cfrac_freq,hour_bins,cfrac_bins,/noerase,/nodata,$
  color=0,xtitle='hours',ytitle='percent cloudy',$
  xstyle=1,ystyle=1,position=pos_vect,yticks=11,$
  xrange=[hour_bins[0]-dhour_bins,hour_bins[n_elements(hour_bins)-1]+dhour_bins],$
  yrange=[cfrac_bins[0]-dcfrac_bins,cfrac_bins[n_elements(cfrac_bins)-1]+dcfrac_bins]

colorbar_david_fanning,maxrange=dmax,minrange=dmin,$
  title='counts',$
  ncolors=254,format='(F6.2)',$
  vertical=1,color=0,position=cb1,charsize=1.5,$
  divisions=8
;
;  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
;
;  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+'.2dhist'
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

end