function pool_check10(Z,NX,NY,VOcean)

% input  Z(current elevation profile)
% output ZS, drain, pool, drainage, ZBeachLevel, AOcean
global ZS drain pool drainage ocean ZBeachLevel AOcean;

pool     = zeros(NX,NY);    %matrix of pooled areas
drain    = zeros(NX,NY);    %matrix of draining points
drainage = zeros(NX,NY);    %matrix of drainage points(points connecting drains and pools)
ocean    = zeros(NX,NY);

ZS = zeros(NX*NY,3);            %Matrix of Z to be sorted.

%copy all NX*NY cell in 1D array ZS 
for j = 1:NY
    for i = 1:NX
        ZS(i+NX*(j-1),1) = Z(i,j);
        ZS(i+NX*(j-1),2) = i;
        ZS(i+NX*(j-1),3) = j;
    end
end

ZS = sortrows(ZS);               %sort the matrix with ascending elevations

%Set the main draining point, lowest points of all
x = ZS(1,2);
y = ZS(1,3);
drain(x,y) = 1;

P = 1;  %Pool number to identify pools and associated drainages (pool index or pool ID)

sumVOcean = 0;      %sum of ocean volume elements, is zero to begin
AOcean = 1;         %number of ocean cells

xLow = ZS(1,2);
yLow = ZS(1,3);
ZBeachLevel      = Z(xLow,yLow);    %elevation at next lowest ocean point (will be updated in the ocean loop)
ocean(xLow,yLow) = 1;               %begin with lowest elevation as an ocean

%Loop over all cell starting the one above the lowest
for I = 2:NX*NY
    x = ZS(I,2);
    y = ZS(I,3);

    if (sumVOcean < VOcean)
        ocean(x,y) = 1;                 %set point as part of ocean
        sumVOcean = sumVOcean + (Z(x,y) - ZBeachLevel)*AOcean;  %add new ocean layer to all existing ocean points
        AOcean = AOcean + 1;                    %increase ocean area
        ZBeachLevel = Z(x,y);           %move the beach up!
    else

        xU = mod(x-1-1,NX)+1;       % normally x-1 but observe p.b.c.
        xD = mod(x+1-1,NX)+1;       % normally x+1 but observe p.b.c.
        yL = mod(y-1-1,NY)+1;       % normally y-1 but observe p.b.c.
        yR = mod(y+1-1,NY)+1;       % normally y+1 but observe p.b.c.

        pL = pool(x,yL);            %Store IDs for surrounding pools (and current maximum elevation of that pool)
        pR = pool(x,yR);
        pU = pool(xU,y);
        pD = pool(xD,y);

        %zero out those pools that have a drainage point already
        %assign values to capitalized PL,PR,PU,PD
        if pL && all(all(drainage~=pL));  %%%Establish if pool is present and  
            PL = pL;
        else
            PL = 0;
        end

        if pR && all(all(drainage~=pR));  %%%no drainage point already
            PR = pR;
        else
            PR = 0;
        end

        if pU && all(all(drainage~=pU));  %%%exists for that pool (ie. 
            PU = pU;
        else
            PU = 0;
        end

        if pD && all(all(drainage~=pD));  %%%non-draining pool)
            PD = pD;
        else
            PD = 0;
        end

        %set DL,DR,DU,DD to pool ID if they have drainage point already
        DL = pL && any(any(drainage==pL));  %%%Establish if pool is present and  
        DR = pR && any(any(drainage==pR));  %%%drainage point does already
        DU = pU && any(any(drainage==pU));  %%%exist for that pool (ie. 
        DD = pD && any(any(drainage==pD));  %%%draining pool)

        %set dL,dR,dU,dD to 1 to indicate a nearby drain or drainage point that
        %is lower in elevantion
        dL = drain(x,yL) || drainage(x,yL) || ocean(x,yL); %%%Establish if point contacts any
        dR = drain(x,yR) || drainage(x,yR) || ocean(x,yR); %%%draining areas or drainage points
        dU = drain(xU,y) || drainage(xU,y) || ocean(xU,y);
        dD = drain(xD,y) || drainage(xD,y) || ocean(xD,y);

        %%%%If connected to a pool that does not yet have a drainage point, the
        %%%%point must either be a pool or a point draining that pool (ie.
        %%%%drainage)

        if PL||PR||PU||PD
            PHigh = max([PL,PR,PU,PD]); % highest nearby pool
            if dL||dR||dU||dD||DL||DR||DU||DD
                drainage(x,y) = PHigh; %drainage is set to ID of pool it drains
            else
                pool(x,y) = PHigh;
            end

            if PL
                pool(pool==PL) = PHigh;
            end
            if PR
                pool(pool==PR) = PHigh;
            end
            if PU
                pool(pool==PU) = PHigh;
            end
            if PD
                pool(pool==PD) = PHigh;
            end

        %%%%if connected to a pool with a drainage point it drains freely
        %%%%to that pool as the drainage point is at a lower elevation.
        %%%%Also, if the point is touching a drain or drainage point it
        %%%%also drains.
        elseif pL||pR||pU||pD||dL||dR||dU||dD
            drain(x,y) = 1;
    %         checked(x,y) = 1;

        %%%%If the point doesn't contact a pool, drain, or drainage point,
        %%%%then it must be a new pool.
        else
            pool(x,y) = P;
    %         checked(x,y) = 1;
            P = P+1; % increase pool index
        end
    end
end 