;
;  Plot a monthly summary
;
pro histogram_average_2d_opt_press;,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,top_press_layer_1,top_press_layer_2,top_press_layer_3,$
  top_press_layer_4,top_press_layer_5

;
;  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

;
;  Filter reflectivity with clutter for sgp
;
if site eq 'sgp' then begin
  
;  Manipulate the data if it is Reflectivity
  reflectivity=float(reflectivity)
  result=where(reflectivity gt -8888, count)
  if count gt 0 then reflectivity[result]=reflectivity[result]/100.

;  Assign all bad clutter reflectivity data -9999
  if n_elements(clutter) ne 0 then begin
   result=where(clutter eq 0, count)
   if count gt 0 then reflectivity[result]=-9999
  endif

;  Turn missing data to gray
  result=where(reflectp eq 0,count)
  if count gt 0 then reflectivity[result]=-8888
endif

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

;
;  Strings
;
dataname='Total Column Optical Depth vs. Cloud Top Pressure of the Top Layer'
datatag='opt_dep_top_press'
cb_title='counts'
;
;  Put together the layer top pressure array
;
top_press_layer=make_array(numtimes,5,/float,value=-9999)
top_press_layer[*,0]=top_press_layer_1
top_press_layer[*,1]=top_press_layer_2
top_press_layer[*,2]=top_press_layer_3
top_press_layer[*,3]=top_press_layer_4
top_press_layer[*,4]=top_press_layer_5
;
;  Pull out the top pressure of the highest layer to
;  put in one array.  
;  top_press=-8888 when numlayers=0 and missing reflectivity
;  top_press=-9999 when numlayers>0 and there is reflectivity (no cloud).
;  top_press=-9999 when numlayers>0 and missing pressure data  
;
top_press=make_array(numtimes,/float,value=-9999)
for i=0,numtimes-1 do begin
;  If there is a layer than pull out the top one
  if numlayers[i] ge 1 then begin
    top_press[i]=top_press_layer[i,numlayers[i]-1]
;  If there are 0 layers then either give -8888 for no reflectivity data (missing)
;  or -9999 for no cloud
  endif else begin
    if reflectp[i,0] eq -8888 or reflectp[i,0] eq 0 then top_press[i]=-8888 else $
    if reflectp[i,0] gt 0  then top_press[i]=-9999 else stop
    ;print,top_press[i],reflectp[i,0]
    ;print,transpose(reflectivity[i,*]),format='(167(F9.2,1X))'
  endelse
  ;print,transpose(top_press_layer[i,*]),format='(5(F9.2,1x))'
  ;print,numlayers[i]
  ;print,i,top_press[i]
endfor

;***************
;  Calculate Optical Depth
;***************

;  Array to hold the optical depth data
opt_depth=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
    opt_depth[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
    opt_depth[i]=opt_depth[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 opt_depth[i]=-8888

endfor

;*****************************************
;  Histogram
;*****************************************

;
;  Set up optical depth bins with mid points 0.5,1.5,2.5,3.5,4.5 hrs
;
opt_bins=5.0
dopt_bins=5.0
end_opt_bins=160.0
while opt_bins[n_elements(opt_bins)-1] lt end_opt_bins $
 do opt_bins=[opt_bins,opt_bins[n_elements(opt_bins)-1]+(2.*dopt_bins)]
;
;  Set up pressure bins with mid points 5,15,25,35,45,55,65, %
;
press_bins=12000.0
dpress_bins=1000.0
end_press_bins=100000.0
while press_bins[n_elements(press_bins)-1] lt end_press_bins $
  do press_bins=[press_bins,press_bins[n_elements(press_bins)-1]+(2.*dpress_bins)]
press_bins=reverse(press_bins)
stop
;
;  Set up the 2D histogram arrays
;
opt_press_freq= $ 
  make_array(n_elements(opt_bins),n_elements(press_bins),/float,value=0) 
;  
; Build the 2D histogram for opt_press
;
for j=0,n_elements(opt_bins)-1 do begin
  for k=0,n_elements(press_bins)-1 do begin
    result=where(opt_depth ge opt_bins[j]-dopt_bins and $
                 opt_depth lt opt_bins[j]+dopt_bins and $
		 top_press ge press_bins[k]-dpress_bins and $
                 top_press lt press_bins[k]+dpress_bins,count) 

    opt_press_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(opt_press_freq)
dmin=min(opt_press_freq)

contour,opt_press_freq,opt_bins,press_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(opt_press_freq,max=dmax,min=dmin,top=254)
result=where(opt_press_freq eq 0)
image[result]=255
tv,congrid(image,sx,sy),px[0],py[0]
 
;  Set up the axes - convert press from pascals to millibars
contour,opt_press_freq,opt_bins,press_bins/100.0,/noerase,/nodata,$
  color=0,xtitle='optical depth',ytitle='pressure (mb)',$
  xstyle=1,ystyle=1,position=pos_vect,$
  xrange=[opt_bins[0]-dopt_bins,opt_bins[n_elements(opt_bins)-1]+dopt_bins],$
  yrange=[(press_bins[0]-dpress_bins)/100.0, $
          (press_bins[n_elements(press_bins)-1]+dpress_bins)/100.0]

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