;
;------------------------------------------------------------------
;
pro mkx,x,y
if x eq 1.e12 or y eq 1.e12 then return
width = .015 & hw = .5*width
plots,[x-hw,x+hw,x,x,x],[y,y,y,y-hw,y+hw],/norm
return
end
;
;------------------------------------------------------------------
;
pro gridvec,un,vn,glon,glat,proj,spval=spv,position=pos,toplon=tlon,$
  perimlat=plat,scalemag=scmag,scalelen=sclen,xscale=xsc,yscale=ysc,$
  color=clr,ps=ips
;
; Plot un+vn vectors on pre-existing map (projection proj) at glat 
;   latitudes and glon longitudes. 
; Do not plot where either un=spval or vn=spval (default spval=1.e36)
; Arrows will be scaled such that a vector of magnitude "scalemag" (m/s) 
;   will have length "scalelen" (normalized 0.->1.)
;   (default scalemag=500. and default scalelat=.05, so the default 
;    scaling will show a 500 m/s arrow .05 in length 
; A horizontal scale arrow (labeled "scalemag" (m/s)) will be centered 
;   over the point (xscale,yscale) (normalized coords)
;
if proj ne "GLOBAL CYLINDRICAL EQUIDISTANT" and $
   proj ne "POLAR STEREOGRAPHIC" and $
   proj ne "SATELLITE VIEW" and $
   proj ne "MOLLWEIDE" then begin
  print,'>>> gridvec: unknown proj=',proj
  return
endif
if n_elements(spv) le 0 then spv = spval
if n_elements(pos) eq 0 then pos = [0.,0.,1.,1.]
if n_elements(tlon) eq 0 then tlon = 0.
if n_elements(plat) eq 0 then plat = 0.
if n_elements(scmag) eq 0 then scmag = 100.
if n_elements(sclen) eq 0 then sclen = .05
if n_elements(clr) eq 0 then clr = !d.table_size-1 ; default white
if n_elements(ips) eq 0 then ips = 0
jmx = n_elements(glat)
imx = n_elements(glon)
xll=pos(0) & yll=pos(1) & xur=pos(2) & yur=pos(3)
;
; If scale magnitude is 0., then use the maximum magnitude rounded
; to the nearest 50 m/s:
;
if scmag eq 0. then begin
  good = where(un ne spv and vn ne spv)
  mag = sqrt(un(good)^2 + vn(good)^2)
  maxmag = max(mag,min=minmag)
  scmag = rnd(maxmag,50.)
  if scmag eq 0. then scmag = float(maxmag)
  print,'gridvec: using scale magnitude of ',scmag,' m/s'
endif
;
; Plot arrows over lat,lon grid:
;
ascale = sclen/scmag
case proj of
  "GLOBAL CYLINDRICAL EQUIDISTANT": begin
    for j=0,jmx-1 do begin
      for i=0,imx-1 do begin
        if un(i,j) ne spv and vn(i,j) ne spv then begin
          tail = convert_coord([glon(i)],[glat(j)],/data,/to_normal)
          xh = tail(0)+ascale*un(i,j) 
          yh = tail(1)+ascale*vn(i,j) 
          mkarrow,tail(0),tail(1),xh,yh,color=clr
        endif
      endfor
    endfor
  end
  "POLAR STEREOGRAPHIC": begin
    for j=0,jmx-1 do begin
      for i=0,imx-1 do begin
        if un(i,j) ne spv and vn(i,j) ne spv then begin
          polcoords,glon(i),glat(j),xt,yt,perimlat=plat,toplon=tlon,$
            position=pos
          rlon = !dtor*(glon(i)-tlon)
          if plat ge 0. then begin          
            du = vn(i,j)*sin(rlon)-un(i,j)*cos(rlon)
            dv =-vn(i,j)*cos(rlon)-un(i,j)*sin(rlon)
          endif else begin
            du = un(i,j)*cos(rlon)+vn(i,j)*sin(rlon) 
            dv =-un(i,j)*sin(rlon)+vn(i,j)*cos(rlon)
          endelse
          xh = xt+ascale*du & yh = yt+ascale*dv      
          hlon = 1.e12 & hlat = 1.e12
          polcoords,hlon,hlat,xh,yh,perimlat=plat,toplon=tlon,position=pos,$
            /todata
          if (plat ge 0. and hlat ge plat) or $
             (plat lt 0. and hlat le plat) then mkarrow,xt,yt,xh,yh,color=clr
        endif
      endfor
    endfor
  end
  "SATELLITE VIEW": begin
    for j=0,jmx-1 do begin
      for i=0,imx-1 do begin
        if un(i,j) ne spv and vn(i,j) ne spv then begin
          tail = convert_coord([glon(i)],[glat(j)],/data,/to_normal)
          xt = tail(0) & yt = tail(1)
          rlon = !dtor*(glon(i)-tlon)
          if glat(j) ge 0. then begin          
            du = vn(i,j)*sin(rlon)-un(i,j)*cos(rlon)
            dv =-vn(i,j)*cos(rlon)-un(i,j)*sin(rlon)
          endif else begin
            du = un(i,j)*cos(rlon)+vn(i,j)*sin(rlon) 
            dv =-un(i,j)*sin(rlon)+vn(i,j)*cos(rlon)
          endelse
          xh = xt+ascale*du
          yh = yt+ascale*dv      
;
; Do not plot if arrow head is outside projection:
;
          head = convert_coord([xh],[yh],/normal,/to_data)
          if head(0) ne 1.e12 and head(1) ne 1.e12 then $
            mkarrow,xt,yt,xh,yh,color=clr
        endif
      endfor
    endfor
  end
  "MOLLWEIDE": begin
    for j=0,jmx-1 do begin
      for i=0,imx-1 do begin
        if un(i,j) ne spv and vn(i,j) ne spv then begin
          tail = convert_coord([glon(i)],[glat(j)],/data,/to_normal)
          xt = tail(0) & yt = tail(1)
          rlon = !dtor*(glon(i)-tlon)
          if glat(j) ge 0. then begin          
            du = vn(i,j)*sin(rlon)-un(i,j)*cos(rlon)
            dv =-vn(i,j)*cos(rlon)-un(i,j)*sin(rlon)
          endif else begin
            du = un(i,j)*cos(rlon)+vn(i,j)*sin(rlon) 
            dv =-un(i,j)*sin(rlon)+vn(i,j)*cos(rlon)
          endelse
          xh = xt+ascale*du
          yh = yt+ascale*dv      
          head = convert_coord([xh],[yh],/normal,/to_data)
          if head(0) ne 1.e12 and head(1) ne 1.e12 then $
            mkarrow,xt,yt,xh,yh,color=clr
        endif
      endfor
    endfor
  end
  else: print,'gridvec: unknown projection: ',proj
endcase
;
; Add scale arrow:
;
if n_elements(xsc) and n_elements(ysc) then begin
  lab = strcompress(string(scmag,format="(f6.1)"),/remove_all) + $
        ' (M/S)'
  xlen = scmag*ascale
  sclr = !d.table_size-1    ; default white scale arrow
  if ips gt 0 then sclr = 0 ; black scale arrow if doing postscript
  mkarrow, xsc-.5*xlen,ysc, xsc+.5*xlen,ysc,color=sclr
  xyouts,xsc,ysc-.03,lab,/norm,alignment=.5
endif
return
end