Global Index (short | long) | Local contents | Local Index (short | long)
[strf, lat_out, lon_out] = tau_to_strf(taux, tauy, lat, lon, bc, H);
[strf, lat_out, lon_out] = tau_to_strf(taux, tauy, lat, lon, bc, H);
This function calls | |
---|---|
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, lat1] = sph_grady1(taux, RADIAN*lat', RADIAN*lon', 1); lat1 = lat1*DEGREE; tauy2 = repmat(NaN, size(dtxdy)); for i = 1:nx; tauy2(:,i) = interp1(lat, tauy(:,i), lat1); end [ny, nx] = size(dtxdy); strf = repmat(NaN, [ny nx]); %%% Start latitude iteration for i = 1:ny; % Find eastern and western boundaries wb = []; eb = []; land = find(isnan(dtxdy(i,:))); land = [0 land nx+1]; if ~isempty(land) & (length(land) < nx-2); ind = find(diff(land) > 2); wb = land(ind)+1; eb = land(ind+1)-1; wrap = iselement(wb, 1); for j = 1:length(wb); ty = tauy2(i, else wb(i) = NaN; eb(i) = NaN; 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;