function [LAYnams,LAYvals,LAYCBD,NLAY]=getLayers(XLSF,zGr)
% [LAYnams,LAYvals,LAYCBD,NLAY]=getLayers(XLSF,zGr)
% XLSF is workbook for this model with sheet LAY
% nLay is number of aquifers obtained from worksheet LAY in XLSF
% nZ   is number of layer interfaces is size(Z,3)
% nZ and nLay are desired to allow checking whether nZ matches
% nLay given LAYCBD in the layer definition.
% Make sure that the list of layers is consistent, apply rules
% for flexibility as defined in the manual
% The highest layer number must be>= size(IBOUND,3)=nZ
% Missing Layers get the values of the next defined layer
% Layers numbered<1 are ignored
%
% TO 091212 091224

if nargin<2,
    help getLayers
    error('Not enough input arguments in getLayers');
end

[LAYnams,layvals]=getExcelData(XLSF,'LAY','Horizontal');

nLay=max(layvals(:,strmatchi('LAYER',LAYnams)));

JP=strmatchi('LAYER',LAYnams);  % column which holds stress period number

% Only the layer that have Layer>0 are valid
layvals=layvals(layvals(:,JP)>0,:);
layvals=sortrows(layvals,JP);

% Fill specified layers in the overall layer matrix, gaps are permittede

LAYvals=NaN(nLay,size(layvals,2)); % Allocate memory to hold all layers
for iP=size(layvals,1):-1:1       % fill in the specified layers, backward
    LayNr=min(nLay,layvals(iP,JP)); % if LayNr>NLAY, put it in Layer(NLAY,...)
    LAYvals(LayNr,:)=layvals(iP,:); % insert given layers
end
% filling in backward means, that the first specified layers wins
% if layer numbers are  5 6 6 35 38 41, while NLAY=30, layers 5, the first
% 6 and 35 will be used in the end (not 41).

% fill in the gaps
for iLay=size(LAYvals,JP):-1:1
    if isnan(LAYvals(iLay,JP))
        LAYvals(iLay,:)=LAYvals(iLay+1,:); % use next one up
    end
end

if any(isnan(LAYvals(:,JP)))
    error(['You must at least specify the last Layer =%d\n',...
        'see worksheet LAY. It must match size(IBOUND,3)']);
end

%% evaluate input of zGr should be (:,:,nZ) or (1:nZ) or nZ
if length(size(zGr))>1
    zGr=zGr(1,1,:);
    nZ=length(zGr);
elseif length(zGr)>1
    nZ=length(zGr);
else
    nZ=zGr(1);
end

%% Check for missing layers at the end, allow users to just specify the
% first or so layer and let mfLab fill in the rest

iCBD =strmatchi('LAYCBD',LAYnams);

LAYCBD = LAYvals(:,iCBD);  % confining beds

% perhaps maximum laynumber < layers in model, then add missing layers at
% the end using the properties of the last layer defined
NLAY=nLay+(nZ-1-nLay-sum(LAYCBD)+LAYCBD(end))/(1+LAYCBD(end));

LAYvals(NLAY,1)=NLAY; % allocate memory to store the rest
for iL=nLay+1:NLAY
    LAYvals(iL,:)=LAYvals(nLay,:);
end

%% update LAYCBD
LAYvals(end,iCBD)=0; % no confining bed below lowest layer
LAYCBD = LAYvals(:,iCBD);  % confining beds

%% Check number of layers nlay with zGr and LAYCBD

NLAY=nZ-1-sum(LAYCBD(1:end-1)); %subtract the number of confining beds

if NLAY~=size(LAYvals,1)
    error(['The number of layers deduced from Z-array (%d)\n',...
        'does not match the highest layer number from worksheet LAY (%d).\n',...
        'Check that the highest layer number in worksheet LAY corresponds\n',...
        'to NLAY (size (IBOUND,3)) of your model!\n\n',...
        'You may also check LAYCBD in the LAY workhseet\n',...
        'NLAY=max(LayrNr)=size(Z,3)-1+sum(LAYCBD(1:NLAY-1)\n',...
        '%3d =      %3d     =  %3d -1 - %3d  ... Check if this is the case !!\n'],...
        NLAY,nLay,nLay,nZ-1-sum(LAYCBD(1:end-1)),nZ,sum(LAYCBD(1:end-1)));
end

