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

pro plot_daily_5min_ciret_six,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')
clutter_id=ncdf_varid(cdf_id, 'qc_ReflectivityClutterFlagMode')

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')

est_clr_id=ncdf_varid(cdf_id,'est_clr_flx_obs')
bsfd_id=ncdf_varid(cdf_id,'bsfd')
calc_sfd_id=ncdf_varid(cdf_id,'calc_sfd')
calc_clr_id=ncdf_varid(cdf_id,'clr_flx_calc')
tb_diff_id=ncdf_varid(cdf_id,'Tb_diff')

cld_base_id=ncdf_varid(cdf_id,'cld_base')
top_base_id=ncdf_varid(cdf_id,'top_base')
cld_base_temp_id=ncdf_varid(cdf_id,'cld_base_temp')

;
;  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 clutter_id ne -1 then ncdf_varget, cdf_id, clutter_id,clutter

if ciret_iwc_id ne -1 then ncdf_varget,cdf_id,ciret_iwc_id,ciret_iwc
if ciret_emiss_id ne -1 then ncdf_varget,cdf_id,ciret_emiss_id,ciret_emiss
if ciret_numtot_id ne -1 then ncdf_varget,cdf_id,ciret_numtot_id,ciret_numtot
if ciret_re_id ne -1 then ncdf_varget,cdf_id,ciret_re_id,ciret_re
if ciret_opt_id ne -1 then ncdf_varget,cdf_id,ciret_opt_id,ciret_opt

if est_clr_id ne -1 then ncdf_varget,cdf_id,est_clr_id,est_clr
if bsfd_id ne -1 then ncdf_varget,cdf_id,bsfd_id,bsfd
if calc_sfd_id ne -1 then ncdf_varget,cdf_id,calc_sfd_id,calc_sfd
if calc_clr_id ne -1 then ncdf_varget,cdf_id,calc_clr_id,calc_clr
if tb_diff_id ne -1 then ncdf_varget,cdf_id,tb_diff_id,tb_diff

if cld_base_id ne -1 then ncdf_varget,cdf_id,cld_base_id,cld_base
if top_base_id ne -1 then ncdf_varget,cdf_id,top_base_id,top_base
if cld_base_temp_id ne -1 then ncdf_varget,cdf_id,cld_base_temp_id,cld_base_temp
;
;  Close the netcdf file
;
ncdf_close,cdf_id

;
;  Only continue if there are retrievals
;
if ciret_iwc_id ne -1 and ciret_emiss_id ne -1 and $
   ciret_numtot_id ne -1 and ciret_re_id ne -1 and $
   ciret_opt_id ne -1 then begin
;
;  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)
;  done in image_scale.pro
;  reflectivity=float(reflectivity)
;  result=where(reflectivity gt -9999,count)
;  if count ne 0 then reflectivity[result]=reflectivity[result]/100.
endif

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

;!x.margin=[10,3]
!x.margin=[11,6]
!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
window,0, xsize=900,ysize=700,/pixmap
result=where(ciret_re gt -9999,count)
;stime=hrfrac[0]
stime=hrfrac[result[0]]
;etime=hrfrac[num_times-1]
etime=hrfrac[result[n_elements(result)-1]]
dx=((etime-stime)/20.0) + stime

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 - Reflectivity
;***************************
if reflect_id ne -1 then begin

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

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)'

; convert cld_base to km
result=where(cld_base gt -9999,count)
;if count gt 0 then cld_base[result]=(cld_base[result]-2*318.0)/1000.0
if count gt 0 then cld_base[result]=cld_base[result]/1000.0

oplot,hrfrac,cld_base,psym=4,symsize=0.5,color=0

; convert top_base to km
result=where(top_base gt -9999,count)
;if count gt 0 then top_base[result]=(top_base[result]-2*318.0)/1000.0
if count gt 0 then top_base[result]=top_base[result]/1000.0

oplot,hrfrac,top_base,psym=4,symsize=0.5,color=50

;key
xyouts,dx,13,'cloud top',color=50,charsize=1,/data
xyouts,dx,12,'cloud base',color=0,charsize=1,/data

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

endif  ;end of plot 1

;*********
;  Plot 2
;*********
 
if ciret_iwc_id ne -1 then begin
  dmax=max(ciret_iwc)
  result=where(ciret_iwc eq -9999,count)
  if count ne 0 then ciret_iwc[result]=1e6
  dmin=min(ciret_iwc)
  dmax=100
  dmin=0.1 
  print,'iwc',dmax,dmin

  plot,hrfrac,ciret_iwc*1000.0,/ylog,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=1,psym=4,title='Ice Water Content',$
    ytitle='IWC (mg/m!E3!N)',color=0

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Ciret IWC',/data,color=0,alignment=0.5
endelse

;******************************************************
;  Plot 3 - Effective Radius and Number Concentration
;******************************************************

if ciret_re_id ne -1 then begin

; convert ciret_re to microns
  result=where(ciret_re gt -9999,count)
  if count gt 0 then ciret_re[result]=ciret_re[result]*1000.0

; determine the data min max range for ciret_re
  dmax=max(ciret_re)
  result=where(ciret_re eq -9999,count)
  if count ne 0 then ciret_re[result]=1e6
  dmin=min(ciret_re)
  dmax=120
  dmin=0

; plot ciret_re
  plot,hrfrac,ciret_re,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=8,psym=4,title='Size & Concentration',$
    ytitle='R_eff (um)',color=0

; convert ciret_numtot to number/liter
  result=where(ciret_numtot gt -9999,count)
  if count gt 0 then ciret_numtot[result]=ciret_numtot[result]/10000.

; determine the data min max range for ciret_re
  dmax=max(ciret_numtot)
  result=where(ciret_numtot eq -9999,count)
  if count gt 0 then ciret_numtot[result]=1e9
  dmin=min(ciret_numtot)

; create the new axis for the plot
  axis,etime,sheight,yaxis=1,color=100,/save,$
    yrange=[dmin,dmax],ystyle=1,ytitle='Concen(/liter)'

; plot ciret_numtot
  oplot,hrfrac,ciret_numtot,color=100,psym=1

; key
  dy=(dmax-dmin)/10.0
  xyouts,dx,9.0*dy+dmin,'size',color=1,charsize=1,/data
  xyouts,dx,8.0*dy+dmin,'concen',color=100,charsize=1,/data

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Ciret RE',/data,color=0,alignment=0.5
endelse

;******************************************************
;  Plot 4 - Emittance and Optical Depth
;******************************************************

if ciret_emiss_id ne -1 then begin

; convert ciret_emiss to microns
;  result=where(ciret_emiss gt -9999,count)
;  if count gt 0 then ciret_emiss[result]=ciret_emiss[result]*1000.0

; determine the data min max range for ciret_re
  dmax=max(ciret_emiss)
  result=where(ciret_emiss eq -9999,count)
  if count ne 0 then ciret_emiss[result]=1e6
  dmin=min(ciret_emiss)
  dmax=1
  dmin=0

; plot ciret_emiss
  plot,hrfrac,ciret_emiss,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=8,psym=4,$
    title='Emittance and Vis. Optical Depth',$
    ytitle='Layer Emittance',color=0

; convert ciret_opt to number/liter
;  result=where(ciret_opt gt -9999,count)
;  if count gt 0 then ciret_opt[result]=ciret_opt[result]/10000.

; determine the data min max range for ciret_opt
  dmax=max(ciret_opt)
  result=where(ciret_opt eq -9999,count)
  if count gt 0 then ciret_opt[result]=1e9
  dmin=min(ciret_opt)

; create the new axis for the plot
  axis,etime,sheight,yaxis=1,color=100,/save,$
    yrange=[dmin,dmax],ystyle=1,ytitle='Vis. Optical Depth'

; plot ciret_opt
  oplot,hrfrac,ciret_opt,color=100,psym=1

; key
  dy=(dmax-dmin)/10.0
  xyouts,dx,9.0*dy+dmin,'emittance',color=0,charsize=1,/data
  xyouts,dx,8.0*dy+dmin,'opt depth',color=100,charsize=1,/data

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Ciret OPT',/data,color=0,alignment=0.5
endelse

;******************************************************
;  Plot 5 - Solar Forcing
;******************************************************

if est_clr_id ne -1 then begin

  solar1=make_array(n_elements(est_clr),/float,value=-9999)
  result=where(est_clr ne -9999 and bsfd ne -9999,count)
  if count gt 0 then solar1[result]=(est_clr[result]-bsfd[result])/est_clr[result]

  solar2=make_array(n_elements(est_clr),/float,value=-9999)
  result=where(calc_clr ne -9999 and calc_sfd ne -9999,count)
  if count gt 0 then solar2[result]=(calc_clr[result]-calc_sfd[result])/calc_clr[result]

; determine the data min max range for radiation
  max1=max(solar1)
  max2=max(solar2)
  dmax=max([max1,max2])
  result=where(solar1 eq -9999,count)
  if count ne 0 then solar1[result]=1e6
  result=where(solar2 eq -9999,count)
  if count ne 0 then solar2[result]=1e6
  min1=min(solar1)
  min2=min(solar2)
  dmin=min([min1,min2])


; plot radiation
  plot,hrfrac,solar1,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=1,psym=4,xtitle='Time (UTC)',$
    title='Solar Forcing',$
    ytitle='W/m!E2!N',color=0

; plot ciret_numtot
  oplot,hrfrac,solar2,color=100,psym=1

; key
  dy=(dmax-dmin)/10.0
  xyouts,dx,9.0*dy+dmin,'est_clr/bsfd',color=0,charsize=1,/data
  xyouts,dx,8.0*dy+dmin,'calc_sfd/calc_clr',color=100,charsize=1,/data

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Solar Forcing',/data,color=0,alignment=0.5
endelse

;******************************************************
;  Plot 5 - Radiation
;******************************************************

if cld_base_temp_id ne -1 then begin

; determine the data min max range for cld_base_temp
  dmax=max(cld_base_temp)
  result=where(cld_base_temp eq -9999,count)
  if count ne 0 then cld_base_temp[result]=1e6
  dmin=min(cld_base_temp)
  ;dmax=120
  ;dmin=0

; plot cld_base_temp
  plot,hrfrac,cld_base_temp,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=8,psym=4,title='Cloud Base Temp & Tb_diff',$
    ytitle='Kelvin',color=0,xtitle='Time (UTC)'

; plot tb_diff 
  if tb_diff_id ne -1 then begin

; determine the data min max range for tb_diff
    dmax=max(tb_diff)
    result=where(tb_diff eq -9999,count)
    if count ne 0 then tb_diff[result]=1e6
    dmin=min(tb_diff)

; create the new axis for the plot
    axis,etime,sheight,yaxis=1,color=100,/save,$
      yrange=[dmin,dmax],ystyle=1,ytitle=''

; plot tb_diff
    oplot,hrfrac,tb_diff,color=100,psym=1
  endif

; key
  dy=(dmax-dmin)/10.0
  xyouts,dx,9.0*dy+dmin,'cld base temp',color=0,charsize=1,/data
  xyouts,dx,8.0*dy+dmin,'tb_diff',color=100,charsize=1,/data

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No Cloud Base Temp',/data,color=0,alignment=0.5
endelse


;
;  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
if site eq 'twpc1' then sitetitle='TWP C1 '+curdate
xyouts,0.5,0.97,'Cirrus Cloud Properties '+sitetitle,$
  /normal,color=0,charthick=1,charsize=1.7,alignment=0.5
;
;  Create gif file
;
dep=strpos(ncdf_file,'cdf')
giffile=strmid(ncdf_file,0, dep)
giffile=giffile+'ciret.six.'
print, giffile+'gif'
write_gif, giffile+'gif', tvrd( )
;
;  Remove png file
;
unix_command='rm '+giffile+'png'
spawn,unix_command
;
;  end of retrieval exists
;
endif

end