
% Set these for single or centered differences
j1 = 2:nsta;
j0 = 1:(nsta-1);

% When dynamic height INCREASES along the cruise track, set
% sign= -1 for velocity to right of cruise track in NH, to left in SH.
% otherwise, sign = 1
sign = 1;

% Example: crossing ACC going south, dynamic height decreases along the
% cruise track (surface relative to deep), velocity is to left; to make it
% positive, set sign to 1.  On northbound leg, set sign to -1 to keep the
% ACC as a positive cross-track component.


pp = p(ones(nsta,1),:)';

[sv sig] = svan(S,T,pp);

% D - Dref = - integral delta dp from pref to p.
% Integrating from surface down (pref = 0), dp is positive:
% Units come out in 1e-4 m2 s-2, hence the scaling factor:

dynht = -cumsum(sv)*2 * 1.e-4;   % assuming 2-dbar spacing

yy =  0.5 * (y(j0) + y(j1));
xx =  0.5 * (x(j0) + x(j1));
% f is 2 omega sin theta; omega is 2 pi radians per 86637 sec (sidereal day)
% one sidereal day shoud be 23h56m4s=86164s    --- Shuiming CHEN,June 22,95
f = 4*pi/86164 * sin(yy*pi/180);
dydeg = (y(j1) - y(j0));
dxdeg = (x(j1) - x(j0));
dy = lat_to_m(dydeg, yy);
dx = lon_to_m(dxdeg, yy);
ds = sqrt(dx.^2 + dy.^2);
fds = f .* ds;
ffds = fds(:,ones(1,ndepth))';

% fu = - dD/dy; fv = dD/dx
u = sign * (dynht(:,j1) - dynht(:,j0))./ffds;

%  Given the definition above of dynht, this is u at each depth referenced
%  to u at the surface; so add the true surface velocity to get the
% absolute velocity at each depth.  (Recall that sign = 1 gives velocity to
% left of track in the NH, to right in SH.)

      eval(['save ' sprintf(['ug%.0f_%.0f.mat'], min(jsta), max(jsta))...
             ' u ds ndepth nsta jsta sign x y xx yy -mat']);

