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

;
;  Plot a monthly summary
;
pro plot_average_1d_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=5
;
;  Select which data to plot
;
dataname=strarr(num_plots)
dataname=['Optical Depth',$
          'Liquid Water Path',$
	  'Ice Water Path',$
	  'TOA Albedo',$
	  'OLR']
datatag=strarr(num_plots)
datatag=['optical_depth',$
         'lwp',$
	 'iwp',$
	 'toa_albedo',$
	 'olr']
cb_title=strarr(num_plots)
cb_title=['log(optical depth)',$
          'g/m^2',$
	  'log(g/m^2)',$
	  'albedo',$
	  'W/m^2']
;
;  Loop through the variables to plot
;
for datanum=0,num_plots-1 do begin  ;loop through the different variables
;for datanum=2,2 do begin

print, datatag[datanum]

;
;  Pick up the data array
;

;***************
;  Optical Depth
;***************
if datanum eq 0 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])

endif

;***************
;  Liquid Water Path
;***************
if datanum eq 1 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_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
endif 

;***************
;  Ice Water Path
;***************
if datanum eq 2 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])

endif 

;***************
;  TOA Albedo
;***************
if datanum eq 3 then begin

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

;  Calculate TOA Albedo
  result=where(TOA_solar_up ge 0 and TOA_solar_down ge 0,count)
  if count gt 0 then data[result]=TOA_solar_up[result]/TOA_solar_down[result]

endif 

;**************
;  OLR
;**************
if datanum eq 4 then data=TOA_IR_up

;
;  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
;
;  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
;
loadct,39 ;5  ;load the color table
!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
!p.background=!d.n_colors-1  ;make the background white
!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=data[week1:week2]
weektime=julian_day[week1:week2]
plot,weektime,week,color=0,$
  xstyle=1,ystyle=1,$
  yrange=[dmin,dmax],ytitle=cb_title[datanum],psym=4,symsize=0.5,$
  xtickunits='time',xtickformat='label_date'
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=data[week1:week2]
weektime=julian_day[week1:week2]
plot,weektime,week,color=0,$
  xstyle=1,ystyle=1,$
  yrange=[dmin,dmax],ytitle=cb_title[datanum],psym=4,symsize=0.5,$
  xtickunits='time',xtickformat='label_date'
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=data[week1:week2]
weektime=julian_day[week1:week2]
plot,weektime,week,color=0,$
  xstyle=1,ystyle=1,$
  yrange=[dmin,dmax],ytitle=cb_title[datanum],psym=4,symsize=0.5,$
  xtickunits='time',xtickformat='label_date'
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=data[week1:week2]
weektime=julian_day[week1:week2]
plot,weektime,week,color=0,$
  xstyle=1,ystyle=1,$
  yrange=[dmin,dmax],ytitle=cb_title[datanum],psym=4,symsize=0.5,$
  xtickunits='time',xtickformat='label_date'
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=data[week1:week2]
weektime=julian_day[week1:week2]
if eday eq 31 then begin 
  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.
  numticks=1
endif

plot,weektime,week,/nodata  ;plot to get size of device
;  Get position of the plot window in device pixels
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
;  Sixe of the plot window in pixels
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
;  cut the size of the plot window down according to adjust variable
sx=sx*adjust
plot,weektime,week,color=0,$
  xstyle=1,ystyle=1,$
  yrange=[dmin,dmax],ytitle=cb_title[datanum],psym=4,symsize=0.5,$
  xtickunits='time',xtickformat='label_date',$
  position=[px[0],py[0],px[0]+sx,py[1]],/noerase,/device
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 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