;
;  This program plots the 2d variables from the 
;  averaged file .comb.  The variables are mostly
;  retrievals and radiation data
;


;
;  Plot a monthly summary
;
pro plot_average_2d_rad,site,searchdate,avg_int
;site='sgp'
;searchdate='200004'
;avg_int='1hr'
;
;  common containing the monthly averaged data
;
common monthly_averaged_data_rad, base_time,time_offset,height,numtimes,numheights,$
  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_rad,site,searchdate,avg_int
;
;  Convert heights to km
;
height=height/1000.0
;
;  Select how many plots to make
;
num_plots=8
;
;  Select which data to plot
;
dataname=strarr(num_plots)
dataname=['Best Estimate of Ice Water Content',$
          'Best Estimate of Liquid Water Content',$
	  'Best Estimate of Liquid Effective Radius',$
	  'Best Estimate of Ice Effective Radius',$
	  'Source Index of the Ice Water Content Algorithm',$
	  'Source Index of the Liquid Water Content Algorithm',$
          'Net Solar Radiation',$
          'Net IR Radiation']
datatag=strarr(num_plots)
datatag=['IWC_best_estimate',$
          'lwc_best_estimate',$
	  're_liquid_best_estimate',$
	  're_ice_best_estimate',$
	  'iwc_source',$
	  'lwc_source',$
          'net_solar_radiation',$
          'net_ir_radiation']
cb_title=strarr(num_plots)
cb_title=['log(g/m^3)',$
          'log(g/m^3)',$
	  'um',$
	  'um',$
	  'iwc_source',$
	  'lwc_source',$
          'W/m^2',$
          'W/m^2']
;
;  Loop through the variables to plot
;
for datanum=0,num_plots-1 do begin  ;loop through the different variables
;for datanum=7,7 do begin

print, datatag[datanum]

;
;  Pick up the data array
;
if datanum eq 0 then data=IWC_best_estimate  
if datanum eq 1 then data=lwc_best_estimate  
if datanum eq 2 then data=re_liquid_best_estimate
if datanum eq 3 then data=re_ice_best_estimate
if datanum eq 4 then data=iwc_source
if datanum eq 5 then data=lwc_source
if datanum eq 6 then begin
  data=make_array(numtimes,numheights,/float,value=-9999)
  result=where(solar_flux_down ge 0 and $
               solar_flux_up ge 0,count)
  if count gt 0 then data[result]=solar_flux_down[result]-solar_flux_up[result]
;  Change the missing flag from -9999 to -8888
  result=where(solar_flux_down eq -9999 and $
               solar_flux_up eq -9999,count)
  if count gt 0 then data[result]=-8888
endif
if datanum eq 7 then begin
  data=make_array(numtimes,numheights,/float,value=-9999)
  result=where(IR_flux_down ge 0 and $
               IR_flux_up ge 0,count)
  if count gt 0 then data[result]=IR_flux_up[result]-IR_flux_down[result]
  result=where(IR_flux_down eq -9999 and $
               IR_flux_up eq -9999,count)
  if count gt 0 then data[result]=-8888
endif
;
;  Do data conversions based on the data array
;
result=where(data gt -8888,count)
if (datanum eq 0 or datanum eq 1) and count gt 0 then begin
  data[result]=alog10(data[result]) ;take the log
endif
if datanum eq 2 and count gt 0 then begin
  data[result]=data[result]*10.0e6  ;convert from m to um
endif
if datanum eq 3 and count gt 0 then begin
  data[result]=data[result]*10.0e6  ;convert from m to um
endif
;
;  Find the max and min data values
;
result=where(data gt -8888,count)
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(data[result])
  dmin=min(data[result])
endelse
print, 'dmax',dmax,'dmin',dmin
;
;  Color the data
;
if datanum ne 4 and datanum ne 5 then begin
;  Bytscl
  image=bytscl(data,max=dmax,min=dmin,top=253)
endif else if datanum eq 4 or datanum eq 5 then begin  ;source
  dmax=7
  dmin=1
  image=data
  source_colors=[0,30,60,90,150,190,217,253]
  for source=dmin,dmax do begin
    result=where(data eq source,count)
    if count gt 0 then image[result]=source_colors[source]
  endfor
endif
;
;  Put all no cloud values to white
;
result=where(data eq -9999,count)
if count gt 0 then image[result]=255
;
;  Put all missing data to gray
;
result=where(data eq -8888,count)
if count gt 0 then image[result]=254

;******************************
;if dmax eq dmin then dmin=0
;if datanum eq 4 then stop
;
;  skip this data num if all the data is missing
;
;if dmax eq -9999 and dmin eq 1e6 then begin
;  print,'no data to plot'
;  goto,no_data_to_plot
;endif
;*******************************

;
;  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
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)

;********************************
;  Set up the plot
;********************************

;  Set up colors
loadct,39 ;5  ;load the color table
!p.background=!d.n_colors-1  ;make the background white
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr
;  make gray for missing data in color 253
r_curr[254]=200
g_curr[254]=200
b_curr[254]=200

;  Set up the rest of the plot
!y.omargin=[2,5]  ;set aside some room at the top of the page for the title
!x.margin=[8,15]  ;change to accomodate the colorbar
!y.margin=[2,2]  ;change to move the images closer together
!p.multi=[0,1,5]  ;put 5 plots on a page
!p.charsize=2  ;increase the character size
!x.ticklen=-0.02
xwin=800  ;xsize of the plot window
ywin=800  ;ysize of the plot window
window,0,xsize=xwin,ysize=ywin,/pixmap
dummy=label_date(date_format=[' %D %M !C %H:%I '])

;********
;  Plot 1
;********
plot1='yes'
if plot1 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,1,0,0,0)
week2=num_sec_ep(syear,smonth,8,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_1
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
if datanum ne 4 and datanum ne 5 then $
  colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_1:
endif
;********
;  Plot 2
;********
plot2='yes'
if plot2 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,8,0,0,0)
week2=num_sec_ep(syear,smonth,15,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_2
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
if datanum ne 4 and datanum ne 5 then $
  colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_2:
endif
;********
;  Plot 3
;********
plot3='yes'
if plot3 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,15,0,0,0)
week2=num_sec_ep(syear,smonth,22,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_3
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
if datanum ne 4 and datanum ne 5 then  $
  colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_3:
endif
;********
;  Plot 4
;********
plot4='yes'
if plot4 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,22,0,0,0)
week2=num_sec_ep(syear,smonth,29,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_4
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
if datanum ne 4 and datanum ne 5 then  $
  colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_4:
endif
;********
;  Plot 5
;********
plot5='yes'
if plot5 eq 'yes' then begin  

week1=num_sec_ep(syear,smonth,29,0,0,0)
week2=num_sec_ep(syear,smonth,31,23,59,60)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_5
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
if eday eq 31 then begin ;determine if 30 or 31 days
  adjust=3./7.   ;31 days
  numticks=3
endif else if eday eq 30 then begin 
  adjust=2./7.  ;30 days
  numticks=2
endif else if eday eq 29 then begin
  adjust=1./7.  ;29 days
  numticks=1
endif

tv,congrid(week,sx*adjust,sy),px[0],py[0]
sz=size(congrid(week,sx*adjust,sy))

contour,week,weektime,height,xstyle=1,ystyle=1,/nodata,/noerase,$
	position=[px[0],py[0],px[0]+sz[1]-1,py[0]+sz[2]-1],/device
tv,congrid(week,sx*adjust,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xrange=[weektime[0],weektime[n_elements(weektime)-1]],$
  xtickformat='(A1)'
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
if datanum ne 4 and datanum ne 5 then  $
  colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_5:
endif
;
;  Add a title
;
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,$
;'Monthly Summary of Hourly Averages for '+
dataname[datanum]+'!C'+capsite+' - '+smon+' '+syr
;
;  Create a code for the sources
;
if datanum eq 4 or datanum eq 5 then begin
  xyouts,0.5,0.05,'1=mace et al 1998',/normal,color=source_colors[1],charsize=1.2 
  xyouts,0.5,0.07,'2=mace et al 2001',/normal,color=source_colors[2],charsize=1.2 
  xyouts,0.5,0.09,'3=ccm3 Parameterization (Kiehl et al 1998)',/normal,color=source_colors[3],charsize=1.2
  xyouts,0.5,0.11,'4=Liu and Illingworth, 2000',/normal,color=source_colors[4],charsize=1.2
  xyouts,0.5,0.13,'5=Dong et al 1998 parameterization',/normal,color=source_colors[5],charsize=1.2
  xyouts,0.5,0.15,'6=Dong and Mace 2001 retrieval',/normal,color=source_colors[6],charsize=1.2
  xyouts,0.5,0.17,'7=frisch et al 1998 lwc retrieval',/normal,color=source_colors[7],charsize=1.2
endif


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

gifdir='/data/mace4/arm_data/temp_image/'
giffile=site+'.'+datastream+'.'+curdate+'.'+datatag[datanum]+'.gif'
print, gifdir+giffile
write_gif,gifdir+giffile,tvrd( )
;
;  Remove the .png file
;
png_file=site+'.'+datastream+'.'+curdate+'.'+datatag[datanum]+'.png'
unix_command='rm '+gifdir+png_file
spawn,unix_command
;
;  Skip to here to skip plotting if data is all missing
;
no_data_to_plot:

endfor  ;end loop through variables

end