%--------------------------------------------------------------------------
%--------
% Purpose:
% Calculate the rate of surface wind stress working of surface geostrophic
% flow
% Specific region of interest via find statement.

% NbLongitudes        longitude in degrees east, from Aviso grid
% NbLatitudes         latitude in degrees north, from Aviso grid
% [xi yi] = meshgrid(NbLongitudes,NbLatitudes);
% dA                  area of grid cells on Aviso grid in m^2
% dA = (1/3 cosd(yi) mperdeg).^2
% ngood               number of weeks of data that contributed to mean
% utaux_bar               <u_g tau_x >, where < > is 6-year average, W/m^2
% vtauy_bar               <v_g tau_y >, 
%     where u_g is the geostrophic velocity from Aviso,
%     updated and tau_x is the surface stress from QuikSCAT.
%     see Scott and Xu, 2009: 
% @Article{ScottXu2009,
%   author =   "Robert B. Scott  and Yongsheng Xu",
%   title =    "An update on the wind power input to the surface geostrophic flow of the World Ocean",
%   journal =  DSRI,
%   year =     "2009",
%   number = {3},
%   volume =   "56",
%   pages =   {295--304}}
%--------------------------------------------------------------------------
%--------

load Pin_u_tau_20000105_20051228_stress3_aviso_upd_smooth0.mat
% find index of region of interest: e.g. the ACC
idx_acc = find(yi>=-60 & yi<=-40 & ngood>26);

% find the area averaged WPI = <u_g tau_x + v_g tau_y>, where <> is a
% 6-year average

kurts_no = sum((utaux_bar(idx_acc)+vtauy_bar(idx_acc)).*dA(idx_acc))./sum(dA(idx_acc));

disp(['area averaged wind power input from QuikSCAT and updated AVISO 40S to 60S = ' ...
    num2str(kurts_no*1000) 'mW/m^2'])