function [A,X,Y] = binningSO_sum(xin,yin,zin,xres,yres,westl,eastl,southl,northl)

% Bins data into a grid of equal resolution gres in x and y
% xin, yin and zin are column vector all with the same size
% global northl southl westl eastl xres yres
% Aug 2, 2011. Adapted for the entire Southern Ocean. Gets rid of profiles
% or data points outside the grid (such that we have nonzero values in
% xbin, ybin and hence accumarray doesn't give error)
% Jan 27, 2016: Added grid as output

% xin=WINDEVE(:,1);
% yin=WINDEVE(:,2);
% zin=WINDEVE(:,3);

% xin=longitude;
% yin=latitude;
% zin=sa_mld;

% Define Grid
xedges=westl:xres:eastl;
yedges=southl:yres:northl;
lnx=length(xedges); lny=length(yedges);
%xborder=westl*ones(lny,1);
%yborder=southl*ones(lnx,1);
gsX=abs(westl-eastl)/xres; % Grid Size will be 151 x 76
gsY=abs(northl-southl)/yres;


% Get rid of all profiles outside the grid boundaries
        % Define polygon
        RLat=[northl northl southl southl northl];
        RLon=[westl eastl eastl westl westl]; % W long

        % Extract region
        in=inpolygon(xin,yin,RLon,RLat);
        
         xin=xin(in);
         yin=yin(in);
         zin=zin(in);


% Create a border with nan such that in every bin there's a value (or nan),
% otherwise accumarray doesn't work.
% maxgsX=360/0.25;
% maxgsY=180/0.25;
% fxedges=-180:xres:180;
% fyedges=-90:yres:90;
if gsX==gsY
    x = [xin; xedges'];
    y = [yin; yedges']; % yborder; (yborder(1:length(xborder)-length(yborder)))
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif gsY>gsX   
    x = [xin; xedges'; (xedges(1:gsY-gsX))']; %xborder];%
    y = [yin; yedges']; % yborder; (yborder(1:length(xborder)-length(yborder)))
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif gsX>gsY
    nrep=gsX/gsY; % 6 for the sO grid 2x2, 3 for grid 2x4
    x = [xin; xedges']; % xborder
    % y = [yin; yedges'; yedges'; (yedges(1:gsX-gsY))']; % yborder; (yborder(1:length(xborder)-length(yborder))) % yedges'; yedges'; yedges'; yedges'   
    yfaked=[yedges';yedges'; yedges(1:end-nrep+1)']; %; yedges';yedges'; yedges'
    y = [yin; yfaked]; % for entire SO 240 x 1440
    z = [zin; nan(lnx,1)]; % need to modify the faked brder accordingly if grid size changes!
end


% Use histc to get the index for each bin.
[nx, xbin]=histc(x,xedges);
[ny, ybin]=histc(y,yedges);

% accumarray computes nanmean on values corresponding to each bin.
A = accumarray([xbin ybin],z,[],@nansum); % gives the accumulated sum by default
A=A';
A(A==0)=NaN;

% Remove faked edge
A=A(1:end-1,1:end-1);
% edges{1} = xedges;
% edges{2} = yedges;
% N = hist3([x y],'Edges',edges); %hist3 ignores NaN
% N=N';
% meanZ = A./N;

% % need to modify if grid is not square!!
% x = [xin; xedges'; (xedges(1:end-1))']; %xborder];%
% y = [yin; yedges']; % yborder; (yborder(1:length(xborder)-length(yborder)))
% z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
% % x,y,z have to be all the same size
% % need to to add a border of nan values to the grid so that at list there's
% % 1 value (or nan) in each bin

% Jan 27, 2016: Added grid as output:
% Grid to plot in a map
[X,Y]=meshgrid(westl+xres/2:xres:eastl-xres/2,southl+yres/2:yres:northl-yres/2);
X=X-xres/2;
Y=Y-yres/2;
