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

pro plot_daily_5min_stret_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')
cbase_id=ncdf_varid(cdf_id,'CloudBaseBestEstimate')
ctop_id=ncdf_varid(cdf_id,'RadarFirstTop')

re_dong_id=ncdf_varid(cdf_id,'re_dong')   ;1d
n_mean_id=ncdf_varid(cdf_id,'N_mean')     ;1d
tau_dong_id=ncdf_varid(cdf_id,'tau_dong') ;1d
lwc_id=ncdf_varid(cdf_id,'lwc')           ;2d
re_id=ncdf_varid(cdf_id,'re')             ;2d
tau_id=ncdf_varid(cdf_id,'tau')           ;2d
;
;  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 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 re_dong_id ne -1 then ncdf_varget,cdf_id,re_dong_id,re_dong    ;1d
if n_mean_id ne -1 then ncdf_varget,cdf_id,n_mean_id,n_mean       ;1d
if tau_dong_id ne -1 then ncdf_varget,cdf_id,tau_dong_id,tau_dong ;1d
if lwc_id ne -1 then ncdf_varget,cdf_id,lwc_id,lwc                ;2d
if re_id ne -1 then ncdf_varget,cdf_id,re_id,re                   ;2d
if tau_id ne -1 then ncdf_varget,cdf_id,tau_id,tau                ;2d
;
;  Close the netcdf file
;
ncdf_close,cdf_id
;
;  Only continue if there are stratus retrievals
;
if re_dong_id ne -1 and n_mean_id ne -1 and $
   tau_dong_id ne -1 and lwc_id ne -1 and $
   re_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
;
;  Change the format of reflectivity
;
reflectivity=float(reflectivity)
result=where(reflectivity gt -8888,count)
if count gt 0 then reflectivity[result]=reflectivity[result]/100.0
;
;  Filter reflectivity with clutter
;
if reflect_id ne -1 and clutter_id ne -1 and site eq 'sgp' then begin
  result=where(clutter eq 0,count)
  if count gt 0 then reflectivity[result]=-9999
endif
;
;  Transpose 2d vars
;
reflectivity=transpose(reflectivity)
lwc=transpose(lwc)
re=transpose(re)
tau=transpose(tau)
;
;  Determine the subsetting indicies
;
;***************************
;  Enter start time below
;***************************
stime=0.0
closetime=hrfrac-stime
closest=min(abs(closetime),indst)
stime=hrfrac[indst]
;***************************
;  Enter end time below
;**************************
etime=24.0
closetime=hrfrac-etime
closest=min(abs(closetime),indet)
etime=hrfrac[indet]
;***************************
;  Enter start height below
;***************************
sheight=0.0
closeheight=height-sheight
closest=min(abs(closeheight),indsh)
sheight=height[indsh]
;***************************
;  Enter end height below
;**************************
eheight=3.0
closeheight=height-eheight
closest=min(abs(closeheight),indeh)
eheight=height[indeh]
;
;  Subset the arrays
;
reflectivity=reflectivity[indst:indet,indsh:indeh]
height=height[indsh:indeh]
cbase=cbase[indst:indet]
ctop=ctop[indst:indet]
lwc=lwc[indst:indet,indsh:indeh]
tau=tau[indst:indet,indsh:indeh]
re=re[indst:indet,indsh:indeh]

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

;!x.margin=[11,6]
!x.margin=[9,9]
!y.margin=[3,2]
!y.omargin=[0,3]
!p.multi=[0,2,3]  
!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

dx=((etime-stime)/20.0) + stime

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

result=where(reflectivity le -8888,count)
if count gt 0 then reflectivity[result]=1e6

dmin=-65.0
dmax=15.0
image=bytscl(reflectivity,min=dmin,max=dmax)

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 cbase to km
result=where(cbase gt -9999,count)
if count gt 0 then cbase[result]=cbase[result]/1000.0

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

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

oplot,hrfrac,ctop,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 - Effective Radius and Number Concentration
;******************************************************

if re_dong_id ne -1 then begin

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

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

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

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

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

; plot ciret_numtot
  oplot,hrfrac,n_mean,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 3 lwc
;************************************
 
if lwc_id ne -1 then begin

contour, lwc, hrfrac, height, /nodata, $
  title='Liquid Water Content (g/m!E3!N)',$
  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

result=where(lwc le -8888,count)
if count gt 0 then lwc[result]=1e6
dmin=0.0
dmax=1.0
image=bytscl(lwc,min=dmin,max=dmax)

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,'lwc',image,dmax,dmin,cb_title,px,py,sx,sy

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

;******************************************************
;  Plot 4 - Optical Depth - 1d
;******************************************************

if tau_dong_id ne -1 then begin

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

; plot tau_dong
  plot,hrfrac,tau_dong,$
    yrange=[dmin,dmax],xrange=[stime,etime],$
    xstyle=1,ystyle=1,psym=4,$
    title='Optical Depth',color=0

; key
  dy=(dmax-dmin)/10.0
  xyouts,dx,9.0*dy+dmin,'opt depth',color=0,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 6 - effective radius - 2d
;******************************************************

if re_id ne -1 then begin

contour, re, hrfrac, height, /nodata, $
  title='Effective Radius (um)',$
  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

dmax=max(re)
result=where(re le -8888,count)
if count gt 0 then re[result]=1e6
dmin=min(re)

image=bytscl(re,min=dmin,max=dmax)

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,'re',image,dmax,dmin,cb_title,px,py,sx,sy

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No 2D Effective Radius',/data,color=0,alignment=0.5
endelse

;******************************************************
;  Plot 5 - tau - 2d
;******************************************************

if tau_id ne -1 then begin

contour, tau, hrfrac, height, /nodata, $
  title='Extinction Coefficient (1/km)',$
  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

dmax=max(tau)
result=where(tau le -8888,count)
if count gt 0 then tau[result]=1e6
dmin=min(tau)
dmin=0.0
dmax=100.0
image=bytscl(tau,min=dmin,max=dmax)

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,'tau',image,dmax,dmin,cb_title,px,py,sx,sy

endif else begin
  plot,holdx,holdy,/nodata
  xyouts,0.5,0.5,'No 2D Optical Depth',/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,'Stratus 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+'stret.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