function [bcsnew, bndelenew] = bcsriver(bcs,meshnod,bndele,L)
% [bcsnew, bndelenew] = bcsriver(bcs,meshnod,bndele);
%  This establishes which nodes of bcs are to be which
% river numbers.  
% Then it fills in the normal oriented length vector
% Messy stuff, examine bcs to identify the bcs nodes
i = sqrt(-1);
figure(1); clf
if nargin==4
% L is already given
else
% create a new L
disp('bcsriver.m  Needs an axis.  Do the rubberband box thing.')
L = newaxis(meshnod);
end

plot(meshnod(bndele(:,2),2),meshnod(bndele(:,2),3),'og')
axis(L)
title(' Click on title line to increment rivernumber RIVERNOW==0')
xlabel('click down here to end')

hold on
flag = 0; river = 0;
rivers = zeros(length(bcs),1);

bcsnew = [bcs(:,1:3) rivers];
bndelenew=bndele;
nodenew=[];
while flag~=-i
  [node,flag]=picknode(meshnod(bndele(:,2),:));
  if flag ==i
   river = river+1;
   title([' RIVER == ', int2str(river)])
  elseif (flag ~=-i & river~=0)
%   rivers(node) = river;
    lrepeat = find(bndele(node,2)==bcsnew(:,1))
    if isempty(lrepeat)
     bcsnew=[bcsnew ;[bndele(node,2), 0 , 0 , river] ];
     bndelenew(node,5)=5;
      nodenew = [nodenew node];
     plot(meshnod(bndele(node,2),2),meshnod(bndele(node,2),3),'*r')
    end
  end
end
hold off

%lr1 = [58 59 60];
%lr2 = [54 55 56 57];
%rivers = zeros(length(bcs),1);
%rivers(lr1) = 1+rivers(lr1);
%rivers(lr2) = 2+rivers(lr2);
%bcsnew = [bcs(:,1:3) rivers];

%  Find and eliminate any new "5" segments which don't have two river nodes
for node = nodenew
l = find(bndelenew(node,3)==bndelenew(nodenew,2) );
if isempty(l)
bndelenew(node,5) = 1;
end
end


l = find(bcsnew(:,4)>0);
l=l(:)';
for i = l
midnod=bcsnew(i,1);
lnb=bndele(find( bndele(:,3)==midnod),2);
lne=bndele(find(bndele(:,2)==midnod ),3);

ny =  -(meshnod(lne,2)-meshnod(lnb,2))/2;
nx = (meshnod(lne,3)-meshnod(lnb,3))/2;
bcsnew(i,2:3) = [nx ny];
end

figure(2);clf
plot(meshnod(bndele(:,2),2),meshnod(bndele(:,2),3),...
meshnod(bndele(:,2),2),meshnod(bndele(:,2),3),'or')
axis('normal')
axis(L)

hold on
arrowuv(meshnod(bcsnew(l,1),2),meshnod(bcsnew(l,1),3),...
    bcsnew(l,2),bcsnew(l,3),1,L)
hold off
