function substring (instr[1]:string, first[1]:numeric, last[1]:numeric)
local instr, first, last, main, p2
begin
   main = stringtochar (instr)          ; convert input string to char array
   if (last .ge. first) then            ; check requested end position
      p2 = last                         ; go to last position specified
   else                                 ; but if less than first:
      p2 = dimsizes (main) - 2          ; go to last avail char in main string
   end if
   return (chartostring (main(first:p2)))       ; extract substring
end

; wind_2.ncl
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
   diri = "./"
   fili = systemfunc ("cd "+diri+" ; ls uv_1975_08_02_00.nc")
   print(fili)
   nfil = dimsizes(fili)

FILENAME1  = fili(0)
FIELDNAME1 = "U10"
FIELDNAME2 = "V10"
; FIELDNAME3 = "tas"
TITLE1     = ""+FILENAME1+" "+FIELDNAME1+""
TITLE2     = ""+FILENAME1+" "+FIELDNAME2+""
; TITLE3     = ""+FILENAME1+" "+FIELDNAME3+""
  wks  = gsn_open_wks("x11","w2_wind")             ; send graphics to PNG file
  res                         = True
  res@vcRefMagnitudeF         = 5.              ; make vectors larger
  res@vcRefLengthF            = 0.050           ; reference vector
  res@vcGlyphStyle            = "CurlyVector"   ; turn on curly vectors
  res@vcMinDistanceF          = 0.022           ; thin the vectors
  res@vcRefAnnoOrthogonalPosF = -1.0            ; move ref vector up
  res@gsnLeftString           = "Rotational Wind" 

  do nf=0,nfil-1
  filename   = fili(nf)                 ;  filename  = "T2_2000_01_01.nc"
  f  = addfile(filename,"r")    ; Open NetCDF file and read variable.
  u = (f->$FIELDNAME1$(0,:,:))
  v = (f->$FIELDNAME2$(0,:,:))
  copy_VarAtts(f->$FIELDNAME1$,u)
  copy_VarAtts(f->$FIELDNAME2$,v)

;*************************************************
; calculate vorticity: Wrap version maintains meta data
;*************************************************
  vrt = uv2vrG_Wrap(u,v)         

;*************************************************
; calculate rotational wind component 
; --
; note: the calculation uses a procedure, so memory 
; must be preallocated.
;*************************************************    
  ur    = new ( dimsizes(u), typeof(u), u@_FillValue )
  vr    = new ( dimsizes(v), typeof(v), v@_FillValue )

  vr2uvg(vrt, ur,vr)          
  ur@long_name  = "Zonal Rotational Wind"
  ur@units      = u@units
  vr@long_name  = "Meridional Rotational Wind"
  vr@units      = v@units
  copy_VarCoords(u, ur ) 
  copy_VarCoords(u, vr ) 
;*************************************************
; plot results
;*************************************************    
                                                ; plot 1st time step
  res@tiMainString =  substring(fili(nf),0,15)          ; put file name on each plot
  plot= gsn_csm_vector_map(wks,ur(:,:),vr(:,:),res)
;  plot= gsn_csm_vector_map(wks,ur(0,:,:),vr(0,:,:),res)
 end do
end