function [w,jd,depth]=ecomvel2(cdf,i,j,lev,iind)
% ECOMVEL extracts a velocity time series from a 4D netCDF ECOM output file.
%
%   USAGE:  [u,jd,depth]=ecomvel(cdf,i,j,lev,[tind]);
%
% cdf = cdf file e.g. 'ecomsi.cdf'
% (i,j) = station location in grid coordinates
% lev = sigma level from which to extract time series (0 is surface)
% w = complex vector of velocity with dimension time
% jd = time vector  (digital Julian day)
% depth = depth of time series (meters)

% Rich Signell (rsignell@usgs.gov)

% 3-27-96 fixed bug that was returning a slightly bad depth value 
%         due to averaging the sigma levels twice

nc=netcdf(cdf);
nz=length(nc('zpos'));
nt=length(nc('time'));
if(~exist('iind')),
  iind=1:nt;
end
sigma=nc{'sigma'}(:);
depth=nc{'depth'}(j,i);
depth=.5*(sigma(lev)+sigma(lev+1))*depth;

base_date=[0 0 0 0 0 0];
base_date(1:3)=nc.base_date(:);
t=nc{'time'}(iind);
jd0=julian(base_date);
jd=jd0+t;

ang=nc{'ang'}(j,i);
u=nc{'u',1}(iind,lev,j,i:i+1);
v=nc{'v',1}(iind,lev,j:j+1,i);
u=reshape(u,[length(iind) 1 1 2]);
v=reshape(v,[length(iind) 1 2 1]);
ucenter=squeeze(0.5*(u(:,:,:,1)+u(:,:,:,2))); %average u and v to center of grid
vcenter=squeeze(0.5*(v(:,:,1,:)+v(:,:,2,:)));
w=ucenter+sqrt(-1)*vcenter;      
w=w(:);
w=w.*exp(sqrt(-1)*ang);   %rotate vectors from grid coordinates to E,N