function [njoin,meshnod1,meshele1,bathy1,meshnod2,meshele2,bathy2] =...  
         submesh(lnod,meshnod,meshele,bathy)
% submesh.m   Nov 3, 1995 alterations for speed
%  [njoin,meshnod1,meshele1,bathy1,meshnodsub,meshelesub,bathysub] =...  
%         submesh(lnod,meshnod,meshele,bathy)  ;
% submesh extraction routine
% lnod is the list of nodes in the submesh
% meshnod1, meshele1  are the old mesh
% meshnod2,meshele2 are the sub mesh
%
% e.g.
% load mesh1all
% L = newaxis(meshnod);
%lnod = inaxis(meshnod(:,2),meshnod(:,3),L);
% submesh2(lnod,meshnod,meshele,bathy)  ;
disp('Submesh Dec. 1995')
tic

i1 = meshele(:,4);
j1 = meshele(:,2);
i2 = meshele(:,2);
j2 = meshele(:,3);
i3 = meshele(:,3);
j3 = meshele(:,4);

l = [1:length(meshele(:,1))]';
S = sparse([i1;i2;i3],[j1;j2;j3],[l;l;l]);
X = S(:,lnod);
[l,m,lele] = find(X);
lele = union(lele,[]);

%disp([' first loop lele  ',num2str(toc)]);tic
 
% Now find the nodes which are in the elements, but not
% in the node list.  These are just outside the region!    lout
lout = nan*ones(3*length(lele),1); louti = 0;
for i = 1:length(lele)
n1 = find(meshnod(lnod,1)==meshele(lele(i),2));
n2 = find(meshnod(lnod,1)==meshele(lele(i),3));
n3 = find(meshnod(lnod,1)==meshele(lele(i),4));

if isempty(n1)
louti = louti+1; lout(louti) =  meshele(lele(i),2) ;
end
if isempty(n2)
louti = louti+1; lout(louti) = meshele(lele(i),3) ;
end
if isempty(n3)
louti = louti+1; lout(louti) = meshele(lele(i),4) ;
end 
end
lout = lout(find(~isnan(lout)));
lout = union(lout,[]);
%disp([' second loop lout ',num2str(toc)]);tic


% Renumber the nodes so that subnodes are at end of list
nn = length(meshnod);
ne = length(meshele);
lnew = [lout' ;lnod];
lold = 1:nn;
lold(lnew) = NaN*lold(lnew);
lold = find(~isnan(lold));
lrenum = [lold'; lnew];
lchange = lrenum;
for i = 1:length(lrenum)
lchange(lrenum(i)) = i;
end

%disp([' renumbering old  ',num2str(toc)]);tic

meshnodnew = meshnod;
meshnodnew(:,2) = meshnod(lrenum,2);
meshnodnew(:,3) = meshnod(lrenum,3);
bathynew = bathy;
bathynew(:,2) = bathy(lrenum,2);

meshelenew = meshele;
meshelenew(:,2) = lchange(meshele(:,2));
meshelenew(:,3) = lchange(meshele(:,3));
meshelenew(:,4) = lchange(meshele(:,4));

%disp([' renumbering new  ',num2str(toc)]);tic

nodsold = 1:(nn-length(lnew));
nodsbnd = (nn-length(lnew)+1):(nn-length(lnod)) ;
nodsnew = (nn-length(lnod)+1):nn;

%disp([' Find elements nend, nbegin  ',num2str(toc)]);tic

% New method to find elements nend, nbegin 

i1 = meshelenew(:,4);
j1 = meshelenew(:,2);
i2 = meshelenew(:,2);
j2 = meshelenew(:,3);
i3 = meshelenew(:,3);
j3 = meshelenew(:,4);

l = [1:length(meshelenew(:,1))]';
S = sparse([i1;i2;i3],[j1;j2;j3],[l;l;l]);
X = S(:,nodsnew);
[l,m,nend] = find(X);
nend = union(nend,[]);

X = S(:,nodsold);
[l,m,nbegin] = find(X);
nbegin = union(nbegin,[]);

% disp([' New method for Find elements nend, nbegin  ',num2str(toc)]);tic

meshele1 = meshelenew(nbegin,:);
meshele2 = meshelenew(nend,:);
meshele2(:,2:4) = meshele2(:,2:4)-length(lold);
meshnod1 = meshnodnew(1:(nn-length(lnod)),:);
bathy1 = bathynew(1:(nn-length(lnod)),:);
meshnod2 = meshnodnew([(length(lold)+1):nn],:);
bathy2 = bathynew([(length(lold)+1):nn],:);
njoin = length(lout);

 
disp([ 'Done with Submesh ',num2str(toc)])

disp(['Main  Mesh  ',int2str(length(meshele(:,1))),' elements  '...
,int2str(length(meshnod(:,1))),' nodes'])
disp(['Outer Mesh  ',int2str(length(meshele1(:,1))),' elements  '...
,int2str(length(meshnod1(:,1))),' nodes'])
disp(['Sub   Mesh  ',int2str(length(meshele2(:,1))),' elements  '...
,int2str(length(meshnod2(:,1))),' nodes'])
