%function [cor,pg,taus,taub,jd]=sterms(cdf,i,j);

cdf='ecomsi.cdf';
ista=4;
i=40;
j=26;    

tau=ecomtau(cdf);
    [d,x,y]=kslice(cdf,'depth');
    [h1,x,y]=kslice(cdf,'depth');
    [h2,x,y]=kslice(cdf,'depth');
    lat=30.
    f= 2. * 7.292E-5 * sin(lat*pi/180);
    for k=1:9,
      [w(:,k),jd,depth(k)]=ecomvel(cdf,i,j,k);
    end
%
% coriolis
    wbar=mean(w.').';
    corx=-f*imag(wbar);
    cory= f*real(wbar);
    cor=corx+sqrt(-1)*cory;
% pressure gradient
    g=9.81;
    e1=ecomelev(cdf,i+1,j);
    e2=ecomelev(cdf,i-1,j);
    e3=ecomelev(cdf,i,j+1);
    e4=ecomelev(cdf,i,j-1);
    dx=2/3*(h1(i,j)+h1(i+1,j)+h1(i-1,j));
    detax= (e1-e2)/dx;
    dy=2/3*(h2(i,j)+h2(i,j+1)+h2(i,j-1));
    detay= (e3-e4)/dy;
    pg=(detax+sqrt(-1)*detay)/g;
    taus=-tau/1000/d(i,j);
    Cd=z0tocd(.003,d(i,j)-abs(depth(9))) 
    taub=Cd*abs(w(:,9)).*w(:,9)/d(i,j);
%
% du/dt (use tsepic data for station ista)
%
     [wtf,jdts]=tsvel('tsepic.cdf',ista);
     wts=mean(wtf.').';
     njd=length(jd);
     jdm=.5*(jd(1:njd-1)+jd(2:njd));  % mid points
     dudt(1)=0; 
     for k=1:length(jdm)-1;
        ind=find(jdts>jdm(k) & jdts<jdm(k+1));
        dudt(k+1)=mean(diff(wts))/(jdts(2)-jdts(1))/3600/24;
        tautmp=Cd*abs(wtf(ind,9)).*wtf(ind,9)/d(i,j);
        taub(k+1)=mean(tautmp(:));
     end
     dudt(njd)=0;
     terms=[dudt(:) cor(:) pg(:) taus(:) taub(:)];
