function LNew=RefineBCN(LOld,SplitArray,NOld,code)
% LNew=RefineBCN(LOld,SplitArray,NOld,code)
% Refines FDM boundary condtion lis according to SplitArray,
% code is 'x', 'y' or 'z' and defined the direction to split
% NOld is the number of cells in that direction
% SplitArray has the following form [i1 i2 i3 i4; N1 N2 N3 N4; I1 I2 I3 I4]
% The first line are the cell indices of the current model to be
% subdivided (splitted up). Only the cells to be splitted need to be
% mentioned.
% The second line indicates into how many new cells this old cell will be
% splitted. Any positive integer is acceptable.
% The third line is optoinal.
% If it is ommitted, then the conductances of the boundary list
% or the flow of a well list will be equally divided over the new refined
% grid cells.
% Else it will indicate over which new cells the conductances and the flows
% of the well list will be (equally) distributed.
% 0  = the center of the new subcells gets the complete concuctance or flow
% >0 = The conductance or flow is equally divided over the first new cells
%      up to the number given by the value. If this number is larger than the
%      the number of subcells indicated in the second line it is the same as
%      division equally over all cells.
% <0 = The conductances or flow is equally divided over the last new cells
%      as indicated by the number given on the third line. Use -1 to
%      attribute all weight on the last cell, -4 to divie the total
%      conductance or flow over the last 4 of thenew cells. If abs(value)
%      is larger than the total number of new cells, then the conductance
%      or flow will be equally distributed over all new cells. The result
%      will then be the same as a large positive number > NNew or omitting
%      the third line of Splitarray altogether.
% TO 100604

code=lower(code(1));

switch lower(code(1))
    case 'x', LCol=4;
    case 'y', LCol=3;
    case 'z', LCol=2;
    otherwise
        error('code must be ''x'',''y'' or ''z'' to indicate splitting direction');
end

IOld=ones(2,NOld);
IOld(1,:)=1:NOld;
IOld(2,SplitArray(1,:))=SplitArray(2,:);

NNew=sum(IOld(2,:));

if size(SplitArray,1)<3
    SplitArray(3,:)=SplitArray(2,:);
else
    I=find(SplitArray(3,:)<-abs(SplitArray(2,:)));
    if ~isempty(I), SplitArray(3,I)=-SplitArray(2,I); end
    I=find(SplitArray(3,:)>+abs(SplitArray(2,:)));
    if ~isempty(I), SplitArray(3,I)=+SplitArrah(2,I); end
end

INew=zeros(2,NNew);
INew(1,:)=1:NNew;
INew(2,end)=sum(IOld(2,:)); 
for iOld=NOld:-1:1,
   INew(2,1:sum(IOld(2,1:iOld)))=iOld;
end

FracNew=ones(1,NNew);  % the fractions of the old conductance and flow in the new cells
for i=1:size(SplitArray,2)
    iOld=SplitArray(1,i);
    I=find(INew(2,:)==iOld);
    FracNew(I)=0;
    frac=abs(SplitArray(3,i))/SplitArray(2,i);
    if SplitArray(3,i)>0        % divide over the first so many new cells
        FracNew(I(1:SplitArray(3,i)))=frac;
    elseif SplitArray(3,i)<0    % divide over the last so many new cells
        FracNew(I(end+SplitArray(3,i)+1:end))=frac;
    else
        FracNew(fix(mean(I)))=1; % middle cell gets all conductance or flow
    end
end

% Last
%% Divide the conductance the the NSplit (skip wells)
LNew=[];
for iNew=1:NNew
    I=find(LOld(:,LCol)==INew(2,iNew)); % This is the corresponding iOld
    if ~isempty(I) && FracNew(iNew)>0
        L=LOld(I,:);
        L(:,LCol)=iNew;
        if size(L,2)==6 || size(L,2)==7,   % not well but DRN, CHD, RIV, GHB
            L(:,6)=L(:,6)*FracNew(iNew);
        else
            L(:,5)=L(:,5)*FracNew(iNew);   % WELLS
        end
        LNew=[LNew;L];
    end
end
