;
;  create the plot of the 5 min averaged data
;

pro plot_daily_5min,site,curdate

site=strcompress(site,/remove_all)

;
;  Get the file using the date as an input
;

curdate=strcompress(curdate,/remove_all)
ncdf_file=strcompress('Output/'+site+'.average.5min.'+curdate+'.000000.cdf',/remove_all)

;
;  Open the netcdf file
;

cdf_id=ncdf_open(ncdf_file)

;
;  Get the dimension id's
;

time_did=ncdf_dimid(cdf_id, 'time')
height_did=ncdf_dimid(cdf_id, 'height')

;
;  Get the dimensions
;

ncdf_diminq, cdf_id, time_did, char_string, num_times
ncdf_diminq, cdf_id, height_did, char_string, num_heights

;
;  Get the variable id's
;

baset_id=ncdf_varid(cdf_id, 'base_time')
timeo_id=ncdf_varid(cdf_id, 'time_offset')
height_id=ncdf_varid(cdf_id, 'height')
reflect_id=ncdf_varid(cdf_id, 'ReflectivityBestEstimate')
reflectp_id=ncdf_varid(cdf_id, 'ReflectivityPossiblePoints')
reflectc_id=ncdf_varid(cdf_id, 'ReflectivityCloudyPoints')
clutter_id=ncdf_varid(cdf_id, 'qc_ReflectivityClutterFlagMode')
cbase_id=ncdf_varid(cdf_id, 'CloudBaseBestEstimate')
ctop_id=ncdf_varid(cdf_id, 'RadarFirstTop')
velocit_id=ncdf_varid(cdf_id, 'MeanDopplerVelocity')
spectr_id=ncdf_varid(cdf_id, 'SpectralWidth')

down_short_id=ncdf_varid(cdf_id,'down_short')
down_long_id=ncdf_varid(cdf_id,'down_long')
up_short_id=ncdf_varid(cdf_id,'up_short')
up_long_id=ncdf_varid(cdf_id,'up_long')

fluxcl_id=ncdf_varid(cdf_id, 'twsflux_clearskyfit')
solarra_id=ncdf_varid(cdf_id, 'SolarRatio')
lwp_id=ncdf_varid(cdf_id, 'liq')
if lwp_id eq -1 then lwp_id=ncdf_varid(cdf_id, 'liq2')
vceil_id=ncdf_varid(cdf_id,'cloud1')
if vceil_id eq -1 then vceil_id=ncdf_varid(cdf_id,'first_cbh')

ciret_iwc_id=ncdf_varid(cdf_id,'IWC_mean')
ciret_emiss_id=ncdf_varid(cdf_id,'emissivity_mean')
ciret_numtot_id=ncdf_varid(cdf_id,'Number_total_mean')
ciret_re_id=ncdf_varid(cdf_id,'radius_effective_mean')
ciret_opt_id=ncdf_varid(cdf_id,'fl_vis_opt_depth_mean')
;
;  Get the variables
;

if baset_id ne -1 then ncdf_varget, cdf_id, baset_id, base_time
if timeo_id ne -1 then ncdf_varget, cdf_id, timeo_id, time_offset
if height_id ne -1 then ncdf_varget, cdf_id, height_id, height

if reflect_id ne -1 then ncdf_varget, cdf_id, reflect_id, reflectivity
if reflectp_id ne -1 then ncdf_varget, cdf_id, reflectp_id, reflectp
if reflectc_id ne -1 then ncdf_varget, cdf_id, reflectc_id, reflectc
if cbase_id ne -1 then ncdf_varget, cdf_id, cbase_id, cbase
if ctop_id ne -1 then ncdf_varget, cdf_id, ctop_id, ctop
if velocit_id ne -1 then ncdf_varget, cdf_id, velocit_id, velocity
if spectr_id ne -1 then ncdf_varget, cdf_id, spectr_id, spectrum
if clutter_id ne -1 then ncdf_varget, cdf_id, clutter_id,clutter

if down_long_id ne -1 then ncdf_varget, cdf_id, down_long_id,down_long
if down_short_id ne -1 then ncdf_varget, cdf_id, down_short_id,down_short
if up_long_id ne -1 then ncdf_varget, cdf_id, up_long_id, up_long
if up_short_id ne -1 then ncdf_varget, cdf_id, up_short_id, up_short

if fluxcl_id ne -1 then ncdf_varget, cdf_id, fluxcl_id, fluxclear
if solarra_id ne -1 then ncdf_varget, cdf_id, solarra_id, solarratio
if lwp_id ne -1 then ncdf_varget, cdf_id, lwp_id, lwp
if vceil_id ne -1 then ncdf_varget, cdf_id, vceil_id, vceil

;
;  Close the netcdf file
;

ncdf_close,cdf_id

;
;  Calculate Julian Date (jday) and Fractional Hour (hrfrac)
;

ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss

;
;  Change height to km 
;
height=height/1000.0
;
;  Get and filter clutter
;
if clutter_id ne -1 then begin
  result=where(clutter eq 0,count)
  if count gt 0 then reflectivity[result]=-9999
endif

;
;  Change the format of reflectivity
;
if reflect_id ne -1 then begin
  reflectivity=transpose(reflectivity)
;  reflectivity=float(reflectivity)
;  result=where(reflectivity gt -9999,count)
;  if count ne 0 then reflectivity[result]=reflectivity[result]/100.
endif

;
;  Start the plot
;

!x.margin=[10,3]
!y.margin=[3,2]
!y.omargin=[0,3]
!p.multi=[0,2,3]  ;eightplot
;!p.multi=[0,1,1]
!p.charsize=2  ;make letters bigger
loadct, 5  ;colortable
!p.background=!d.n_colors-1  ;change background color to white
blanks=[' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ']
window,0, xsize=900,ysize=700,/pixmap
stime=hrfrac[0]
etime=hrfrac[num_times-1]
sheight=height[0]
eheight=height[num_heights-1]
print,' heights',sheight,eheight
print, 'times',stime,etime
;
;  Set up place holder plot
;
holdx=[0,1]
holdy=[0,1]

;*********
;  Plot 1
;*********
if reflect_id ne -1 then begin

contour, reflectivity, hrfrac, height, /nodata, $
  title='Radar Reflectivity (DbZe)',$
  xstyle=1,ystyle=1,$
  xmargin=[5,8],$
  xrange=[stime,etime],yrange=[sheight,eheight],color=0

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

image_scale, 'ReflectivityBestEstimate', reflectivity, image,$
  stime,etime,hrfrac,sheight,eheight,height,dmax,dmin

tv,congrid(image,sx,sy),px[0],py[0]

axis,stime,sheight,xaxis=0,color=0,xrange=[stime,etime],xstyle=1
axis,stime,eheight,xaxis=1,color=0,xrange=[stime,etime],xstyle=1,$
  xtickformat='(A1)'

axis,stime,sheight,yaxis=0,color=0,yrange=[sheight,eheight],ystyle=1
axis,etime,sheight,yaxis=1,color=0,yrange=[sheight,eheight],ystyle=1,$
  ytickformat='(A1)'

if cbase_id ge 0 then begin
  oplot,hrfrac,cbase/1000.,color=0,psym=4,symsize=0.5
  plots,[(etime-stime)*0.85,(etime-stime)*0.933],$
    [(eheight-sheight)*0.75,(eheight-sheight)*0.75],$
    /data,linestyle=0,color=0,psym=4,symsize=0.5
  xyouts,(etime-stime)*0.70,(eheight-sheight)*0.75,$
    /data,'Cloud Base',color=0,charsize=1
endif

if ctop_id ge 0 then begin
  oplot,hrfrac,ctop/1000.,color=50,psym=4,symsize=0.5
  plots,[(etime-stime)*0.85,(etime-stime)*0.95],$
    [(eheight-sheight)*0.85,(eheight-sheight)*0.85],$
    /data,linestyle=2,color=50,psym=4,symsize=0.5
  xyouts,(etime-stime)*0.70,(eheight-sheight)*0.85,$
    /data,'Cloud Top',color=50,charsize=1
endif

colorbar,'ReflectivityBestEstimate',image,dmax,dmin,cb_title,px,py,sx,sy

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Merged Moments',/data,color=0,alignment=0.5
endelse ;end of plot 1

;*********
;  Plot 2
;*********
 
if fluxcl_id ne -1 or down_short_id ne -1 then begin

; find the max and min value for fluxclear
  if fluxcl_id ne -1 then begin
    fmax=max(fluxclear)
    result=where(fluxclear eq -9999,count)
    if count ne 0 then fluxclear[result]=1e6
    fmin=min(fluxclear)
  endif else begin
    fmax=0
    fmin=0
  endelse

; find the max and min value for down short
  if down_short_id ne -1 then begin
    pmax=max(down_short)
    result=where(down_short eq -9999,count)
    if count ne 0 then down_short[result]=1e6
    pmin=min(down_short)
  endif else begin
    pmax=0
    pmin=0
  endelse

; find the max and min of both maxes and mins
  dmax=max(fmax,pmax)
  dmin=min(fmin,pmin)

  if down_short_id ne -1 then begin
    plot, hrfrac,down_short,$
      yrange=[dmin,dmax],xrange=[stime,etime],$
      xstyle=1,ystyle=1,psym=1,$
      ytitle='w/m!E2',title='Flux',color=0
    if fluxcl_id ne -1 then oplot,hrfrac,fluxclear,psym=6,color=60
  endif else begin
    plot,hrfrac,fluxclear,$
      yrange=[dmin,dmax],xrange=[stime,etime],$
      xstyle=1,ystyle=1,psym=1,$
      ytitle='w/m!E2',title='Flux',color=0
  endelse

  ypt=(dmax-dmin)/100.
  xpt=(hrfrac[n_elements(hrfrac)-1]-hrfrac[0])/20.

  plots,[xpt*9,xpt*9],[72*ypt,72*ypt],/data,psym=1,color=0
  xyouts,xpt*10,ypt*70,/data,'Solar Down',color=0,charsize=1
  plots,[xpt*9,xpt*9],[62*ypt,62*ypt],/data,psym=6,color=60
  xyouts,xpt*10,ypt*60,/data,'Flux Clear - Long',color=60,charsize=1

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Skyrad or Sirs',/data,color=0,alignment=0.5
endelse

;*********
;  Plot 3
;*********

if velocit_id ne -1 then begin

  result=where(velocity gt -9999,count)
  if count gt 0 then velocity[result]=velocity[result]/1000.

  velocity=transpose(velocity)
  contour, velocity, hrfrac, height, /nodata,$
    title='Velocity (m/s)',xstyle=1,ystyle=1,$
    xmargin=[5,8],xrange=[stime,etime],$
    yrange=[sheight,eheight],color=0

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

  image_scale, 'MeanDopplerVelocity', velocity, image,$
    stime,etime,hrfrac,sheight,eheight,height,dmax,dmin
  tv,congrid(image,sx,sy),px[0],py[0]

  axis,stime,sheight,xaxis=0,color=0,xrange=[stime,etime],xstyle=1
  axis,stime,eheight,xaxis=1,color=0,xrange=[stime,etime],xstyle=1,$
    xtickformat='(A1)'

  axis,stime,sheight,yaxis=0,color=0,yrange=[sheight,eheight],ystyle=1
  axis,etime,sheight,yaxis=1,color=0,yrange=[sheight,eheight],ystyle=1,$
    ytickformat='(A1)'

  if vceil_id ge 0 then begin
    oplot,hrfrac,vceil/1000.,color=40,psym=4,symsize=0.5
    plots,[(etime-stime)*0.85,(etime-stime)*0.95],$
      [(eheight-sheight)*0.85,(eheight-sheight)*0.85],$
      /data,color=40,psym=4,symsize=0.5
    xyouts,(etime-stime)*0.63,(eheight-sheight)*0.85,$
      /data,'Laser Ceilometer',color=40,charsize=1
  endif

  colorbar,'MeanDopplerVelocity',image,dmax,dmin,cb_title,px,py,sx,sy
endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Merged Moments',/data,color=0,alignment=0.5
endelse

;********
;  Plot 4
;********

if solarra_id ne -1 then begin
plot, hrfrac, solarratio, yrange=[0,1],xrange=[stime,etime],psym=6,$
  title='Solar Transmisivity Ratio = Solar Down/Long',$
  xstyle=1,ystyle=1,color=0
endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Skyrad or Sirs',/data,color=0,alignment=0.5
endelse

;********
;  Plot 5
;********

if spectr_id ne -1 then begin

  result=where(spectrum gt -9999,count)
  if count gt 0 then spectrum[result]=spectrum[result]/1000.

  spectrum=transpose(spectrum)
  contour, spectrum, hrfrac, height, /nodata,$
    title='Spectrum Width (m/s)',$
    xtitle='Time (UTC)',$
    xstyle=1,ystyle=1,$
    xmargin=[5,8],$
    xrange=[stime,etime],yrange=[sheight,eheight],color=0

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

  image_scale, 'SpectralWidth', spectrum, image, $
    stime,etime,hrfrac,sheight,eheight,height,dmax,dmin
  tv,congrid(image,sx,sy),px[0],py[0]

  axis,stime,sheight,xaxis=0,color=0,xrange=[stime,etime],xstyle=1
  axis,stime,eheight,xaxis=1,color=0,xrange=[stime,etime],xstyle=1,$
    xtickformat='(A1)'

  axis,stime,sheight,yaxis=0,color=0,yrange=[sheight,eheight],ystyle=1
  axis,etime,sheight,yaxis=1,color=0,yrange=[sheight,eheight],ystyle=1,$
    ytickformat='(A1)'

  colorbar,'SpectralWidth',image,dmax,dmin,cb_title,px,py,sx,sy
endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Merged Moments',/data,color=0,alignment=0.5
endelse

;********
;  Plot 6
;********

if lwp_id ne -1  then begin
  result=where(lwp gt -9999,count)
  if count gt 0 then lwp[result]=lwp[result]*10000.0

  dmax=max(lwp)
  result=where(lwp eq -9999,count)
  if count ne 0 then lwp[result]=1e6
  dmin=min(lwp)

  if dmax eq -9999 and dmin eq 1e6 then begin
     plot,holdx,holdy,/nodata
     xyouts,0.5,0.5,'LWP was all -9999',/data,color=0,alignment=0.5
  endif else begin
     plot, hrfrac, lwp,title='LWP - Microwave Radiometer',$
       ytitle='g/m!E2', $
       xtitle='Time (UTC)',$
       xstyle=1,ystyle=1,psym=4,$
       yrange=[dmin,dmax],xrange=[stime,etime],color=0
  endelse
endif
;
;  Add title to page
;
if site eq 'nsa' then sitetitle='NSA '+curdate
if site eq 'sgp' then sitetitle='SGP '+curdate
if site eq 'twpc2' then sitetitle='TWP C2 '+curdate
xyouts,0.5,0.95,sitetitle,/normal,color=0,charthick=2,charsize=3,alignment=0.5
;
;  Create gif file
;
dep=strpos(ncdf_file,'cdf')
print, dep
giffile=strmid(ncdf_file,0, dep)

print, giffile
write_gif, giffile+'gif', tvrd( )
;
;  Remove png file
;
unix_command='rm '+giffile+'png'
spawn,unix_command

end