function [wobs,wmod,wind,jdmod,gind]=tscomp(cdfobs,cdfmod,sta,z);
%function [wobs,wmod,wind,jdmod,gind]=tscomp(cdfobs,cdfmod,sta,z);
disp(['station: ' cdfobs]);
cdfmod2='tsepic.jan.cdf';
depth=mcvgt(cdfmod,'depth');
sigma=mcvgt(cdfmod,'sigma');
depth=depth(sta)*sigma;
% which levels bracket depth? 
lev2=max(find(depth>z));
lev1=lev2+1;
% get model time series 
[w2,jdmod]=tsvel(cdfmod,lev2);
w2=w2(:,sta);
[w1,jdmod]=tsvel(cdfmod,lev1);
w1=w1(:,sta);
frac=(z-depth(lev1))/(depth(lev2)-depth(lev1));
wmod=w1+frac*(w2-w1);
%
% get model time series2 
[w2,jdmod2]=tsvel(cdfmod2,lev2);
w2=w2(:,sta);
[w1,jdmod2]=tsvel(cdfmod2,lev1);
w1=w1(:,sta);
frac=(z-depth(lev1))/(depth(lev2)-depth(lev1));
wmod2=w1+frac*(w2-w1);
jd1=jdmod(1);
jd2=jdmod2(length(jdmod2));
dt=jdmod(2)-jdmod(1);
nt=(jd2-jd1)/dt;
time=[0:nt-1]*fix(dt*24*3600+.5);
jd=jd1+time/(24*3600);
ind1=tind(jd,gregorian(jdmod(length(jdmod))));
ind2=tind(jdmod2,gregorian(jdmod(length(jdmod))))-1;
jdmod=jd;
wmod([ind1:length(jdmod)])=wmod2([ind2:length(wmod2)-2]);

% lowpass model data
[wmod,jdmod]=plfilt(wmod,jdmod);
startstop=ss(jdmod);
%
%loc=mcvgt(cdfmod,'loc');
%ii=loc(1,sta);
%jj=loc(2,sta);
%
% get observed time series
%
% add one day on both ends so lowpass covers it
start=startstop(1,:);
jd=julian(start);
start=gregorian(jd-1);
stop=startstop(2,:);
jd=julian(stop);
stop=gregorian(jd+1);
%
chdir /usr3/rps/massbay/data/hour
[wobs,jdobs]=getvel(cdfobs,start,stop);
chdir /usr3/rps/massbay/heatflux
[wind,jdw]=getvel('bulk9010.cdf',start,stop,'east_str','north_str');
chdir /usr3/rps/models/ecomsi/mbay/nas
[wobs,jdobs]=plfilt(wobs,jdobs);      % lowpass observed hourly data
[wind,jdw]=plfilt(wind,jdw);      % lowpass observed hourly data
wobs=interpm(jdobs(:),wobs(:),jdmod);  % interpolate onto model time base
wind=interpm(jdw(:),wind(:),jdmod);  % interpolate onto model time base
wmod = wmod*100;   %convert to cgs

%fill gaps in observed data with model data
ind=find(wobs==NaN);
gind=find(wobs~=NaN);
frac=100*length(ind)/length(wobs);
disp([num2str(frac) ' percent gaps'])
%wobs(ind)=wmod(ind);
  

