function [xGr,yGr,ZGR,xm,ym,ZM,dx,dy,DZ,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr,LAYCBD)
% [xGr,yGr,ZGR,xm,ym,ZM,dx,dy,DZ,Nx,Ny,Nz]=modelsize3(xGr,yGr,xGr,LAYCBD)
% makes values in xGr,yGr,zGr unique sorted
% xGr runs upward
% yGr runs downward
% zGr runs downward and is a 3d vector (1,1,Nz)
% if zGr in the input list is a 3D array, then ZGR,ZM and DZ in the output
% will also by 3D arrays.
% A 3D zGr input is taken as is and so it must
% already be sorted such that the the last layer has the lowest elevation.
% LAYCBD is the vector NLAY long indicating whether a layer as a confining
% unit below it.
% then compute the other values
% TO 090104 091207

% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later


FLGZ3D=numel(zGr)>max(size(zGr));  % zGr is given as a vector

if ~FLGZ3D
    
    xGr=unique(xGr(:)');
    yGr=unique(yGr(:) ); yGr=flipud(yGr);
    zGr=flipud(unique(zGr(:)));

    if nargin<4, LAYCBD=zeros(size(zGr(1:end-1))); end

    xm=0.5*(xGr(1:end-1)+xGr(2:end));
    ym=0.5*(yGr(1:end-1)+yGr(2:end));
    zm=0.5*(zGr(1:end-1)+zGr(2:end));

    dx=abs(diff(xGr));
    dy=abs(diff(yGr));
    dz=abs(diff(zGr));

    Nx=length(dx);
    Ny=length(dy);
    Nz=length(LAYCBD);
    
    ZGR=zGr(:);  ZGR=reshape(ZGR,[1,1,length(ZGR)]);
    ZM=NaN(size(LAYCBD));
    DZ=NaN(size(LAYCBD));
    
    k=1;
    for iLay=1:length(LAYCBD)
        ZM(iLay)=0.5*(ZGR(k)+ZGR(k+1));
        DZ(iLay)= ZGR(k)-ZGR(k+1);
        k=k+1;
        if LAYCBD(iLay)
            k=k+1;
        end
    end

    DZ=abs(DZ);
    
else
    if nargin<4, LAYCBD=zeros(size(zGr,3)-1,1); end
    fprintf('modelsize3: numel>max(size(zGr)) --> full 3D Z array input is assumed.\n');
    
    if length(size(xGr))==length(size(zGr)) && all(size(xGr)==size(zGr))
        xm=0.5*(xGr(:,1:end-1,:)+xGr(:,2:end,:));
        dx=abs(diff(xGr,1,2));
        Nx=size(dx,2);
    else
        xGr=unique(xGr(:)');
        xm=0.5*(xGr(1:end-1)+xGr(2:end));
        dx=abs(diff(xGr));
        Nx=length(dx);
    end
    if length(size(yGr))==length(size(zGr)) && all(size(yGr)==size(zGr))
        ym=0.5*(yGr(1:end-1,:,:)+yGr(2:end,:,:));
        dy=abs(diff(yGr,1,1));
        Ny=size(dy,1);
    else
        yGr=unique(yGr(:) ); yGr=flipud(yGr);
        ym=0.5*(yGr(1:end-1)+yGr(2:end));
        dy=abs(diff(yGr));
        Ny=length(dy);
    end
   
    Nz=length(LAYCBD);

    if zGr(1,1,end)>zGr(1,1,end),   % sort if necessary
        zGr=zGr(:,:,end:-1:1);
        fprintf('modelsize3: zGr array now sorted\n');
    end
    ZM=NaN(Ny,Nx,Nz);
    DZ=NaN(Ny,Nx,Nz);
    k=1;
    for iLay=1:Nz
        ZM(:,:,iLay)=0.5*(zGr(:,:,k)+zGr(:,:,k+1));
        DZ(:,:,iLay)=zGr(:,:,k)-zGr(:,:,k+1);
        k=k+1;
        if LAYCBD(iLay)
            k=k+1;
        end
    end
    
    DZ=abs(DZ);
    
    % Fill in a full nodal grid array
    ZGR               =zeros(Ny+1,Nx+1,Nz+1);
    
    ZGR(2:end-1,[1 end],:)=0.5*(zGr(1:end-1,[1 end],:)+zGr(2:end,[1 end],:));
    ZGR([1 end],2:end-1,:)=0.5*(zGr([1 end],1:end-1,:)+zGr([1 end],2:end,:));

    ZGR([1 end],[1 end],:)=zGr([1 end],[1 end],:);
    
 
    if size(zGr,1)>1 && size(zGr,2)>1   
        ZGR(2:end-1,2:end-1,:)=(zGr(1:end-1,1:end-1,:)+...
                                zGr(2:end  ,1:end-1,:)+...
                                zGr(2:end  ,2:end  ,:)+...
                                zGr(1:end-1,2:end  ,:))/4;
    end

    
    
end 
