Documentation of tau_to_strf

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[strf, lat_out, lon_out] = tau_to_strf(taux, tauy, lat, lon, bc, H);

Help text

  [strf, lat_out, lon_out] = tau_to_strf(taux, tauy, lat, lon, bc, H);

Cross-Reference Information

Listing of function tau_to_strf

function [strf, lat_out, lon_out] = tau_to_strf(taux, tauy, lat, lon, bc, H);

global RADUS RADIAN DEGREE

[ny, nx] = size(taux);

%  Get f, beta and rho

beta = (2*7.292e-5)*cos(lat*RADIAN)/RADUS;
rho = 1e3;

%  First, get dtaux/dy

dtxdy = sph_grady(taux, RADIAN*lat', RADIAN*lon', 1);
%lat = DEGREE*ry2;

%  Find eastern and western boundaries

wb = repmat(NaN, [ny 1]);
eb = wb;
for i = 1:ny;
  land = find(isnan(dtxdy(i,:)));
  if ~isempty(land) & (length(land) < nx-2);
    ind = min(find(diff(land) > 2));
    wb(i) = land(ind)+1;
    eb(i) = land(ind+1)-1;
  else
    wb(i) = NaN;
    eb(i) = NaN;
  end
end

%  Interpolate tauy to new grid

dx = mean(diff(lon));
lon_out = [(lon(1) - dx/2); (lon + dx/2)];

ty2 = NaN * ones(ny, nx+1);
for i = 1:ny;
  ty2(i,:) = interp1(lon', tauy(i,:), lon_out');
end

%  Make sure eastern and western boundaries are kosher

for i = 1:ny;
  if ~isnan(wb(i));
    y1 = ty2(i, wb(i)+1); y2 = ty2(i, wb(i)+2); 
    if isnan(ty2(i, wb(i)));
      ty2(i, wb(i)) = 2*y1 - y2;
    end
    y1 = ty2(i, eb(i)-2); y2 = ty2(i, eb(i)-1);
    if isnan(ty2(i, eb(i)));
      ty2(i, eb(i)) = 2*y2 - y1;
    end
    if isnan(ty2(i, eb(i)+1));
      ty2(i, eb(i)+1) = 3*y2 - 2*y1;
    end
  end
end

%  Get tauy terms

ty = NaN * ones(size(ty2));
for i = 1:ny;
  if ~isnan(wb(i));
    ty(i, wb(i):(eb(i)+1)) = ty2(i, wb(i):(eb(i)+1)) - ty2(i, eb(i)+1);
  end
end

%  Get wind stress integrand

txint = dtxdy .* (RADUS * cos(RADIAN * lat) * diff(RADIAN * lon_out)'); %

%  Integrate westward

strf = NaN * ones(ny, nx+1);
for i = 1:ny;
  if ~isnan(wb(i));
    strf(i, eb(i)+1) = 0;
    tem = txint(i, wb(i):eb(i));
    tem = fliplr(cumsum(fliplr(tem)));
    strf(i,wb(i):eb(i)) = tem;
  end
end

%  Get constant

const = (1/(rho*H)) * (1 ./ beta) * ones(1, nx+1);

%  Get boundary condition

bc2 = NaN * ones(ny, 1);
for i = 1:ny;
  if ~isnan(eb(i));
    bc2(i) = bc(i, eb(i)+1);
  end
end
bc2 = bc2 * ones(1, nx+1);

%  Calculate strf

strf = const .* (ty - strf) + bc2;

lat_out = lat;