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) = NX*NY;
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