function [w,jd,depth]=ecomvel(cdf,i,j,zi,iind)
% ECOMVEL extracts a velocity time series from a 4D netCDF ECOM output file.
%
%   USAGE:  [u,jd,depth]=ecomvel(cdf,i,j,[zi],[tind]);
%
% cdf = cdf file e.g. 'ecomsi.cdf'
% (i,j) = station location in grid coordinates
% zi = z levels from which to extract time series 
%    If zi=nan or missing, the default is to extract data at all sigma levels
% iind= time indices 
%
% w = complex vector of velocity 
% jd = time vector  (digital Julian days)
% depth = depth(s) of velocity time series (meters)
%
% examples:  To extract velocity at all levels and all time steps
%            for cell(i=30,j=47) from the file 'ecomsi.cdf': 
%
%              [w,jd,depths]=ecomvel('ecomsi.cdf',30,47);
%
%            To extract velocity at 5 and 25 m for time steps 1:10,
%            for cell(i=30,j=47) from the file 'ecomsi.cdf': 
%
%              [w,jd,depths]=ecomvel('ecomsi.cdf',30,47,[-5 -25],1:10);
%
%            To extract velocity at 5 and 25 m for all time steps
%            for cell(i=30,j=47) from the file 'ecomsi.cdf': 
%
%              [w,jd,depths]=ecomvel('ecomsi.cdf',30,47,nan,1:10);

% Rich Signell (rsignell@usgs.gov)

ncid=mexcdf('open',cdf,'nowrite');
[nam,nz]=mexcdf('diminq',ncid,'zpos');
[nam,nt]=mexcdf('diminq',ncid,'time');
if(exist('iind')==1),
  istart=min(iind);
  icount=min(max(iind),nt)-istart+1;
else
  istart=0;
  icount=nt;
end
sigma=mexcdf('varget',ncid,'sigma',0,nz);
sigma=.5*(sigma(1:(nz-1))+sigma(2:nz));
depth=mexcdf('varget',ncid,'depth',[j-1 i-1],[1 1]);

base_date=[0 0 0 0 0 0];
base_date(1:3)=mexcdf('attget',ncid,'global','base_date');
t=mexcdf('varget',ncid,'time',istart,icount);
jd0=julian(base_date);
jd=jd0+t;

ang=mexcdf('varget',ncid,'ang',[j-1 i-1],[1 1],1);
u=mexcdf('varget',ncid,'u',[istart 0 j-1 i-1],[icount -1 1 2],1);
v=mexcdf('varget',ncid,'v',[istart 0 j-1 i-1],[icount -1 2 1],1);
u=reshape(mean(u),nz,icount)';
v=reshape(mean(v),nz,icount)';

depth=sigma*depth;
w=u+sqrt(-1)*v;       %average u and v to center of grid
w=w.*exp(sqrt(-1)*ang);   %rotate vectors from grid coordinates to E,N  

if(exist('zi')==1),
  if(~isnan(zi(1))),
  m=length(zi);
  if(min(zi)<min(depth)), disp('requested level below data!'),return,end
  if(max(zi)>max(depth)), disp('requested level above data!'),return,end
  for k=1:m,
    lev2=max(find(depth>zi(k)));
    lev1=lev2+1;
    frac=(zi(k)-depth(lev1))/(depth(lev2)-depth(lev1));
    wmod(:,k)=w(:,lev1)+frac*(w(:,lev2)-w(:,lev1));
  end
  w=wmod;
  depth=zi;
  end
end
