;
; Makes a graph (png file) of the BVMB model output. 
;
; INPUT: 
;
;   filename		string naming the file
;
; OPTIONS:
;   ntime=ntime		inteer number of time slice from output file
;   /usa		flag for map projection over USA with states
;   /polar		flag for polar map projection
;   /web		copy output file to ~/public_html
;   /x			use X windows device
;
; OUTPUT:
;  plot.png		output graphics file created
;
; USES: cbar.pro	makes a simple colorbar
;
;
; David Noone <dcn@colorado.edu> - Mon Nov 30 23:56:07 MST 2009
;
;
pro bvmbplot,filename,ntime=ntime,usa=usa,polar=polar,web=web,x=x

  nn = 0
  if keyword_set(ntime) then nn = ntime

;
; Read the data
;
  print,'Reading data...'
  ncid = ncdf_open(filename,/nowrite)
  ncdf_varget,ncid,'lon',lons
  ncdf_varget,ncid,'lat',lats
  ncdf_varget,ncid,'VOR',vor 
  ncdf_varget,ncid,'Z',hgt
  ncdf_varget,ncid,'U',ug
  ncdf_varget,ncid,'V',vg
  ncdf_close,ncid
;
; Reform for just the time we want
;
  vor = reform(vor[*,*,nn])
  ug  = reform(ug [*,*,nn])
  vg  = reform(vg [*,*,nn])
  hgt = reform(hgt[*,*,nn])
;
; Wrap the data at grenwich, so the plot is pretty
;
  lons = [lons,359.99]
  vor = [vor[*,*],vor[0,*]]
  ug  = [ug [*,*],ug [0,*]]
  vg  = [vg [*,*],vg [0,*]]
  hgt = [hgt[*,*],hgt[0,*]]

  nlons = n_elements(lons)
  nlats = n_elements(lats)

;
; Set up the plot
;

  if keyword_set(x) then begin
    set_plot,'X',/copy
  endif else begin
    set_plot,'Z',/copy
    device,set_resolution=[800,800]
    device,set_pixel_depth=24,decomposed=0
    device,z_buffering=0
  endelse

;  loadct,23
  loadct,39

  !p.font=1
  !p.color=0 
  !p.background=255

;
; Set up the map projection: either USA or whole NH (polar)
;

  length = 1
  stride = 2
  if (keyword_set(usa)) then begin	
      print,'Using USA projection...'
      map_set,41,255,/azim,limit=[15, -135, 70, -75],/iso

  endif else begin

     if (keyword_set(polar)) then begin
       print,'Using polar projection...'
       map_set,90,0,/STEREOGRAPHIC,/isotropic, $
            limit=[90,0,15,0]
     endif else begin
       print,'Using global projection...'
       length = 3.
       map_set,0,180,/isotropic, $
            limit=[89,0,1,0]
     endelse 

  endelse

;
; Do a color shading for the relative vorticity
;
   nc = 41
   clevvor = 1.0e-4*(2.*findgen(nc)/(nc-1) - 1)
   cmax = max(clevvor)
   cmin = min(clevvor)
   i = where(vor gt cmax,n)
   if (n gt 0) then vor[i] = cmax - 1e-7
   i = where(vor lt cmin,n)
   if (n gt 0) then vor[i] = cmin + 1e-7
   contour,vor,lons,lats,/cell_fill,/overplot, $
           levels=clevvor
;
; Continental outlines and grid
;
  map_continents,thick=3,color=70
  map_continents,thick=2,color=70,/usa
  map_grid,linestyle=2,color=20,thick=2

;
; Overlay the geopotential height contours
;
   nc = 50
   cint = 50.
   clev =  cint*2*floor(min(hgt)/(2*cint)) + cint*findgen(nc)
   contour,hgt,lons,lats,/follow,/overplot, $
           levels=clev,thick=5,c_charsize=2

;
; Overlay the vector field for u and v
;
   ipts = where((findgen(nlons) mod stride) eq 0)
   u2 = ug[ipts,*]
   v2 = vg[ipts,*]
   x2 = lons[ipts]

   jpts = where((findgen(nlats) mod stride) eq 0)
   u2 = u2[*,jpts]
   v2 = v2[*,jpts]
   y2 = lats[jpts]
   

   velovect,u2,v2,x2,y2,/overplot,length=length,thick=2

;
; Add a color bar
;
  cbar,clevvor[0:*:4]*1e5,position=[.1,.09,0.9,.12],format='(f5.1)', $
    charsize=2
  xyouts,0.7,0.13,'Vorticity x1e5',/norm,charsize=3

;
; Write out the graphics to a file
;
  if (not keyword_set(x)) then begin
    image = tvrd(/true)
    write_png,'plot.png',image
    if (keyword_set(web)) then spawn,'cp plot.png ~/public_html/plot.png'
  endif
;
  print,'All done. See plot.png'
end