function zonebudget(B,Zonearray,zones,allflag)
% zonebudget(B [,Zonearray[,zones,[allflag]])
% compute budget for each zone in zones as spatially given by Zonearray
% use allflag=1 to get the budget for all zones in zones combined
% B is budget struct generated by readBud
% Zonearray is a 3D array (NROW,NCOL,NLAY) with zone numbers
% zones is the set of zones for which the budget is desired
%
% TO 091129


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

switch nargin
    case 1,
        Zonearray=ones(length(B(1).rows),length(B(1).cols),length(B(1).lays));
        zones=1;
        allflag=0;
    case 2,
        zones=1;
        allflag=0;
    case 3
        allflag=0;
    case 4
        allflag=1;
        ZA=zeros(size(Zonearray));
        for iz=1:length(zones)
            ZA=ZA | Zonearray==zones(iz);
        end
        zones=1;
end

for i=1:length(B)
    for iz=1:length(zones)
        Qtoti=0; Qtoto=0; 
        if ~allflag, ZA=(Zonearray==zones(iz)); end
        if allflag, what='all zones'; else what=sprintf('zone %d',zones(iz)); end

        fprintf('\nBudget [consistent model units] for period=%d, itsp=%d:\n\n',B(i).period,B(i).tstp);
        fprintf('BUDGET FOR %-20s  -----IN-----  ----OUT-----\n',upper(what));

        for iterm=1:length(B.label)
            Lbl=B(i).label{iterm};
            if strmatch('FLOWFRONTFACE',B(i).label{iterm})
                QFF=diff(ZA,1,1).*B(i).term{iterm}(1:end-1,:,:);
                QFFi= sum(QFF(QFF>0));  Qtoti=Qtoti+QFFi;
                QFFo=-sum(QFF(QFF<0));  Qtoto=Qtoto+QFFo;
                fprintf('total for %20s = %12.2f  %12.2f\n',Lbl,QFFi,QFFo);
                
            elseif strmatch('FLOWRIGHTFACE',B(i).label{iterm})
                QRF=diff(ZA,1,2).*B(i).term{iterm}(:,1:end-1,:);
                QRFi= sum(QRF(QRF>0));  Qtoti=Qtoti+QRFi;
                QRFo=-sum(QRF(QRF<0));  Qtoto=Qtoto+QRFo;
                fprintf('total for %20s = %12.2f  %12.2f\n',Lbl,QRFi,QRFo);
                
            elseif strmatch('FLOWLOWERFACE',B(i).label{iterm})
                QLF=diff(ZA,1,3).*B(i).term{iterm}(:,:,1:end-1);
                QLFi= sum(QLF(QLF>0));  Qtoti=Qtoti+QLFi;
                QLFo=-sum(QLF(QLF<0));  Qtoto=Qtoto+QLFo;
                fprintf('total for %20s = %12.2f  %12.2f\n',Lbl,QLFi,QLFo);
                
            else
                Q=sum(B(i).term{iterm}(ZA));
                Qi= sum(Q(Q>0));  Qtoti=Qtoti+Qi;
                Qo=-sum(Q(Q<0));  Qtoto=Qtoto+Qo;
                fprintf('total for %20s = %12.2f  %12.2f\n',Lbl,Qi,Qo);
                
            end
        end
        fprintf('          %20s   ------------  ------------\n',' ');
        fprintf('total for %20s = %12.2f  %12.2f\n',what,Qtoti,Qtoto);
        fprintf('total net %20s = %12.2f  \n',what,Qtoti-Qtoto);
    end
end
