function [meshnodjoin,meshelejoin,bathyjoin]=...
    joinmesh(njoin,meshnod1,meshele1,bathy1,meshnod2,meshele2,bathy2)          
%  [lnod,meshnodjoin,meshelejoin,bathyjoin]=...
%    joinmesh(njoin,meshnod1,meshele1,bathy1,meshnod2,meshele2,bathy2);   
% Join the submeshes back together 
% njoin is the number of overlapping nodes
% 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);
%  [njoin,meshnod1,meshele1,meshnodsub,meshelesub] =...  
%         submesh(lnod,meshnod,meshele,bathy)  ;
%   Mess around with meshnodsub and meshelesub,
%   but don't touch the overlap nodes  1:njoin
%   (Elements using them are OK to change, delete, add)
%   Splice the two sets back together:
%  [meshnod,meshele,bathy]=...
%      joinmesh(njoin,meshnod1,meshele1,bathy1,meshnod2,meshele2,bathy2);       


% Join the two sets back together:
meshelejunk=meshele2;
meshelejunk(:,2:4) = meshele2(:,2:4)+length(meshnod1)-njoin;
meshelejoin = [meshele1; meshelejunk];
meshnodjoin = [meshnod1; meshnod2((njoin+1):length(meshnod2),:) ];
bathyjoin = [bathy1; bathy2((njoin+1):length(meshnod2),:) ];
%figure(3)
%ploteles(meshnodjoin,meshelejoin);
%title('Rejoined')


meshnodjoin(:,1) = [1:length(meshnodjoin)]';
meshelejoin(:,1) = [1:length(meshelejoin)]';
bathyjoin(:,1) = [1:length(bathyjoin)]';


