function NewLayer=JoinLayers(OldLayer,JoinArray,DZ,code)
% NewLayer=JoinLayers(OldLayer,JoinArray,DZold,code)
% joins layer array OldLayer according to JoinArray
% JoinArray has form [OldLayers; NewLayers]
% DZold is the array with the DZ of the original model
% code is 'k', 'c' or 'z'
% k joins according to transmissivities
% c joins according to resistances
% z yields the new cell elevations, then OldLayer=Zold and NewLayer=Znew
% Z is the array with top and bottom elevations of cells
% TO 100601

if ~all(JoinArray(1,:)==sort(JoinArray(1,:)))
    error('JoinArray(1,:) must be 1: number of old layers');
end
NewLayers=unique(JoinArray(2,:));

if ~all((1:length(NewLayers)) == NewLayers)
    error('JoinArray(2,:) must contain all new layers 1:NNew<NNold');
end

NewLayer=NaN(size(OldLayer,1),size(OldLayer,2),length(NewLayers));

for iLay=1:length(NewLayers)
    IZ=find(JoinArray(2,:)==iLay);
    if code=='k'
        NewLayer(:,:,iLay)=sum(OldLayer( :,:,IZ).*DZ(:,:,IZ),3)./sum(DZ(:,:,IZ),3);
    elseif code=='c'
        NewLayer(:,:,iLay)=sum(DZ(:,:,IZ),3)./sum(DZ(:,:,IZ)./OldLayer( :,:,IZ),3);
    elseif code=='z' || size(OldLayer,3)==size(DZ,3)+1  % cell elevations
        if iLay==1,
            NewLayer(:,:,1)=OldLayer(:,:,1);
        end
        NewLayer(:,:,iLay+1)=OldLayer(:,:,IZ(end)+1);    
    else
        error('code must be ''k'' or ''c'' or ''z''');
    end
end
