function TS=read_ncom_ts(datadir,dlat,dlon,ef,vf,tf,sf)
%  TS=read_ncom_ts(DATADIR, LAT, LON, sshflag, velflag, tempflag, salflag) returns time series of
%  model variables from the model archive grid point nearest LAT/LON. Data are returned as a 
%  structured array TS with fields:
%    U/V - velocity eastward (U) and northward (V) components (m/s), 
%    T - potential temperature (deg C), 
%    S - salinity, 
%    SSH - surface elevation (meters),
%    NCOMLAT/NCOMLON - model grid point latitude and longitude, 
%    Z - the vector of model grid depths (meters)
%    DATE - model date in numeric YYYYMMDDHH format
%  
%  DATADIR is the the path of the directory of the model archive.
%  LAT and LON are the northward latitude and eastward longitude of the desired location.
%  sshflag, velflag, tempflag, and salflag are flags (set to 1 for true, 0 for false) for reading
%  the SSH, velocity, temperature, and salinity variables.
%
%  The model was run for 20 years -- 3 years corresponding to 1 Jan 1993 through 31 Dec 2012.

if (datadir(length(datadir)) ~= filesep())
  datadir=[datadir filesep()];
end;

[lon,lat,z]=read_ncom_grid([datadir 'grid.bin']);
[npts,nlevs]=size(z);

ntimes=58440;

dst=sqrt((dlat-lat).^2+(dlon-lon).^2);
idx=find(min(dst)==dst);
idx=idx(1);
TS.NCOMLAT=lat(idx);
TS.NCOMLON=lon(idx);

nd=max(find(z(idx,:)<0));
TS.Z=squeeze(z(idx,1:nd));

fnum=floor((idx-1)/100);
i1=idx-fnum*100;

fname=[datadir 'ncom_Col_' sprintf('%3.3i',fnum) '.nc'];
tlon=ncread(fname,'longitude',[i1],[1]);
tlat=ncread(fname,'latitude',[i1],[1]);
dpth=ncread(fname,'depth',[i1,1],[1,nlevs]);
if max(abs([tlon tlat dpth]-[TS.NCOMLON TS.NCOMLAT z(idx,:)]))>0
  error(['Error reading ' fname])
end

if ef
 e=squeeze(ncread(fname,'e',[i1,1],[1,ntimes]));
 e(e>100)=nan;
 TS.SSH=e;
end

if vf
 u=squeeze(ncread(fname,'u',[i1,1,1],[1,nd,ntimes]))';
 v=squeeze(ncread(fname,'v',[i1,1,1],[1,nd,ntimes]))';
 u(u>900)=nan;
 v(v>900)=nan;
 [u,v]=rot20(u,v); %Rotate velocity from rotated NCOM grid to E/W
 TS.U=u;
 TS.V=v;
end

if tf
 t=squeeze(ncread(fname,'t',[i1,1,1],[1,nd,ntimes]))';
 t(t>900)=nan;
 TS.T=t;
end

if sf
 s=squeeze(ncread(fname,'s',[i1,1,1],[1,nd,ntimes]))';
 s(s>900)=nan;
 TS.S=s;
end

date=ncread(fname,'date');
date(date<0)=2999010100;
TS.DATE=date;

