%% Analyzing the simulation reusult of the shallow Dtuch top system.

clear variables; close all;

load name
load(basename)

load underneath
 
%% Build the model grid

[xGr,yGr,Z,xm,ym,Zm,DX,DY,DZ,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);


%% Show the heads at the center of the parcels 


H=readDat(['-' basename '.hds'],'','',1); % get heads for these dates only and only layer 1

totim=tne(1,1)+cumsum([H.pertim]);  % compute realworkd time using totim and tne(1,1)
for i=1:length(H), H(i).totim=totim(i); end

Nt=length(H);

%% ix of water divide of all cross sections
%  get the right-most acitive cell of each cross section
% this is the cell corresponding to the center of the parcel

IXend=NaN(Ny,1);
for iy=1:Ny, IXend(iy)=find(squeeze(IBOUND(iy,:,1))~=0,1,'last'); end

IDIV=cellIndex(IXend,(1:Ny)',ones(Ny,1),Nx,Ny,Nz); % global indices of water divide Ny*1

%% Get the heads in the center of the parcels
Phi=NaN(Ny,Nt); % The maximum head in all cross sections and all times
for it=1:length(H)
    Phi(:,it)=H(it).values(IDIV);
end

%% Plot the head for all cross section as a function of time
figure;

plot([H.totim],Phi);  % this plots all sections together
xlabel('t [d]'); ylabel('head [NAP'); title(sprintf('max head in %d sections over time',Ny));
datetick; grid on


%% ===== GHG ==== GVG ==== GLG =====

%% Periods to select for GHG and GLG

[Y M D]=datevec(tne(:,1));  % all dates in tne ---> Datevec: year month day

y=unique(Y); y=y(2:end);    % all years in tne except first, which is skipped

%% Indices in the tne array for GHX computation
Iper=find(Y>=y(1) & (D==14 |D==28));  % all 14 an 28th month day, skip first year

H=readDat([basename '.hds'],Iper,'',1); % get heads for these dates only and only layer 1
H=maskHC(H,1000);
%% contatenate heads and fold so that time is third dimension

h=[H.values]; % concatenate all heads horizontally
NROW=H(1).NROW;
NCOL=H(1).NCOL;
NT  =size(h,2)/NCOL;

h=reshape(h,NROW,NCOL,NT); % make time third dimension

ghg=zeros(NROW,1);  % initialize ghg
glg=zeros(NROW,1);  % initialize glg

for iy=1:length(y);  % for all (8) years
    
    % select year to get min and max 3 values
    I=find(Y(Iper)==y(iy));
    
    % sort per row, to get max and min value in each cross section (ignoring NaNs)
    % this assumes that the lowest heads in a year correspond to the center
    % of the cross section !!! This is generally true for shallow grw in
    % parcels surrounded by ditches
    hM=sort(squeeze(max(h,[],2)),2);  % max head in all cross sections sorted over all measurements in year
    hm=sort(squeeze(min(h,[],2)),2);  % min head in all cross sections sorted over all measurements in year
    
    ghg=ghg+mean(hM(:,3:end),2);      % add average of higest three in year
    glg=glg+mean(hm(:,1:3)  ,2);      % add average of lowest three in year
end
ghg=ghg/length(y);
glg=glg/length(y);
gvg=ghg-(1/6)*(ghg-glg); % Cultuurtech Vadamecum 2000

figure; hold on;
plot(1:NROW,ghg,'b'); plot(1:NROW,glg,'r'); plot(1:NROW,gvg,'g'); legend('GHG','GLG','GVG');
xlabel('Section Nr'); ylabel('GXG [mNAP]'); title(sprintf('GXG over period %d - %d',y(1),y(end)));

for iP=1:length(P);
    P(iP).GHG=ghg(iP);
    P(iP).GVG=gvg(iP);
    P(iP).GLG=glg(iP);
end

%% === water budgets ===== water budgets ===== water budgets =====

%% Compute budegets

ISec=[1 2 3 10]; % only for these cross sections

B=readBud([basename '.bgt'],'','','','',ISec); % selected cross sections for water balance
B=maskHC(B,IBOUND);

%% Terms for the water budget

% Allocate memory. The column is for the sections and the rows are for the
% times or rather time steps during which these flows are present
QDRN=NaN(Ny,Nt);   % flow throug the drains at ground surface
QRCH=NaN(Ny,Nt);   % flow form recharge
QGHB=NaN(Ny,Nt);   % flow through the general head boundaries (to and form the ditches)
QSTO=NaN(Ny,Nt);   % storage change during the time step
QWEL=NaN(Ny,Nt);   % prescribe flow for the seeapage from or to the regional aquifer

%% Fill the arrays using the flow terms from the budget fiile
% first sum along the second (x) dimension and than along the third (z)
% dimension to obtain the total for all Ny separate cross sections
for i=1:length(B)
    try QDRN(:,i)=sum(sum(B(i).term{strmatchi('DRAINS'       ,Label)},2),3); catch end
    try QRCH(:,i)=sum(sum(B(i).term{strmatchi('RECHARGE'     ,Label)},2),3); catch end
    try QGHB(:,i)=sum(sum(B(i).term{strmatchi('HEADDEPBOUNDS',Label)},2),3); catch end
    try QSTO(:,i)=sum(sum(B(i).term{strmatchi('STORAGE'      ,Label)},2),3); catch end
    try QWEL(:,i)=sum(sum(B(i).term{strmatchi('WELLS'        ,Label)},2),3); catch end
end

%% Total these flows to see of the water balance is zero (or almost)
QTOT=sum(QDRN,1)+sum(QRCH,1)+sum(QGHB,1)+sum(QSTO,1)+sum(QWEL,1);
QABS=abs(sum(QDRN,1))+abs(sum(QRCH,1))+abs(sum(QGHB,1))+abs(sum(QSTO,1))+abs(sum(QWEL,1));


figure;
plot([H.totim],QTOT,'b'); xlabel('t [d]'); ylabel('Q [m3/d]'); title('Running budget of entire model'); grid on; hold on
plot([H.totim],QABS,'r');
legend('Qtot','QAbs');
datetick; grid on;

%% print these flows for all cross sections for any time
% of for one cross section of all times
% or the sum
fprintf('         QWEL     QSTO     QDRN     QGHB     QRCH');
display([sum(QWEL,2)  sum(QSTO,2) sum(QDRN,2) sum(QGHB,2) sum(QRCH,2)]);

%% Show running water budget
j=1;
POS=[max(0,QWEL(j,:));max(0,QSTO(j,:));max(0,QGHB(j,:));max(0,QDRN(j,:));max(0,QRCH(j,:))];
NEG=[min(0,QWEL(j,:));min(0,QSTO(j,:));min(0,QGHB(j,:));min(0,QDRN(j,:));min(0,QRCH(j,:))];

figure; hold on
area([B.period],POS');
area([B.period],NEG');
legend('QWEL','QSTO','QGHB','DRN','QRCH','QWEL','QSTO','QGHB','DRN','QRCH');
xlabel('t [d]');
ylabel('Flow m2/d');
title(sprintf('Water budget over %d cross sections as a function of time',Ny));
datetick; grid on;