% The following is a modification of geost.m, customized 
%  for the equatorial part of P17c, stations 74 to 123.
% From 74 to 80, indices 1 to 7, the spacing is 0.5 degree.
% From 80 to 116, indices 7 to 43, the spacing is 1/10 degree.
% From 116 to 118, indices 43 to 45, the spacing is 1/4 degree.
% From 118 to 123, indices 45 to 50, the spacing is 0.5 degree.


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


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

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

% Integrate from surface down, pressure increasing,
% to get dyn ht at depth minus dynht at surface; I have put in
% a minus sign to get D at surface minus D at depth.
% Units come out in 1e-4 m2 s-2, hence the scaling factor:

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

% Allocate space for a dynamic height array gridded at 1/6 degree:
dynht_fi = zeros(ndepth,73);
y_f = (6:-(1/6):-6)';

% Fill in the parts of this array that coincide with the station positions:
j_station = [1 4 7 10 13 16 19:55 58 61 64 67 70 73];
dynht_fi(:,j_station) = dynht(:,[1:43 45:50]); % Note: 3.25S station not on grid 

% Interpolate the remainder:
d_dh = diff(dynht')';
d_y = diff(y)';
deriv_dy = d_dh ./ d_y(ones(ndepth,1),:);

j_int1 = [1 4 7 10 13 16 58 61 64 67 70];
j_int01 = [1:6 45:49];
j_int2 = [56 57];
j_int02 = [43 44];

dy = (y_f(j_int1+1) - y(j_int01))';
dynht_fi(:,j_int1+1) = dynht(:,j_int01) + ...
                deriv_dy(:,j_int01) .* dy(ones(ndepth,1),:);

dy = (y_f(j_int1+2) - y(j_int01))';
dynht_fi(:,j_int1+2) = dynht(:,j_int01) + ...
                deriv_dy(:,j_int01) .* dy(ones(ndepth,1),:);

dy = (y_f(j_int2) - y(j_int02))';
dynht_fi(:,j_int2) = dynht(:,j_int02) + ...
                deriv_dy(:,j_int02) .* dy(ones(ndepth,1),:);

dynht_f = dynht_fi;

% Filter the interpolated dynamic height.
for ii = 1:6
   dynht_f(:,2:(nsta_f-1)) = ...
          0.25 * (dynht_f(:,1:(nsta_f-2)) + dynht_f(:,3:nsta_f)) + ...
          0.5  * (dynht_f(:,2:(nsta_f-1)));
end


yy =  0.5 * (y_f(j0) + y_f(j1));

% f is 2 omega sin theta; omega is 2 pi radians per 86164 sec (sidereal day)
f = 4*pi/86164 * sin(yy*pi/180);
dy = (y_f(j1) - y_f(j0)) * 1.1e5;  % last number converts deg to meters
fdy = f .* dy;
ffdy = fdy(:,ones(1,ndepth))';

% fu = - dD/dy;
u = -(dynht_f(:,j1) - dynht_f(:,j0))./ffdy;

%  Given the definition above of dynht, this is u at the surface 
%  minus u at each depth.
%  From the whole array, subtract the row for a reference depth to
%  get u relative to that depth.

%  Use the equatorial form of geostrophy within 0.5 degrees of the
%  equator.  The equator is index 37 in the dynamic height array,
%  so indices 34 through 40 are the range from -0.5 to 0.5 inclusive.
%  In the u array, the grid points bounded by this range are 34-39.
%  To estimate the second derivative, use only the (smoothed) dynamic
%  heights at +-0.5 and the equator.

%  Alternative: estimate second deriv from a larger range, say +-1,
%  indices 31-43.

beta = 2*pi/(86637 * 6323000);
Dyy = (dynht_f(:,31) + dynht_f(:,43) - 2 * dynht_f(:,37)) ...
                                     /( (y_f(43)-y_f(37))^2 * 110000^2);
u(:,34:39) = -Dyy(:,ones(34:39))/beta;
