%
%Victor Chavarrias (victor.chavarrias@deltares.nl)
%
%$Revision$
%$Date$
%$Author$
%$Id$
%$HeadURL$
%

%% PREAMBLE

clear
clc
fclose all;

%% ADD OET

%change to your path to vtools/general in OET. 
path_add_fcn='p:\dflowfm\projects\2020_d-morphology\modellen\checkout\openearthtools_matlab\applications\vtools\general\';
addpath(path_add_fcn)
addOET(path_add_fcn)

%% INPUT

in_read.branch={'rhine','lek'}; %branch
% simdef.D3D.dire_sim=fullfile(pwd,'../','\01_simulations\r004\dflowfm\');
simdef.D3D.dire_sim=fullfile(pwd,'../','\01_simulations\r006\dflowfm\');
simdef.flg.which_v=1; 
%   1=etab
%   2=h
%   3=dm Fak
%   4=dm fIk
%   5=fIk
%   6=I
%   7=elliptic
%	8=Fak
%   9=detrended etab based on etab_0
%   10=depth averaged velocity
%   11=velocity
%   12=water level
%   13=face indices
%   14=active layer thickness
%   15=bed shear stress
%   16=specific water discharge
%   17=cumulative bed elevation
%   18=water discharge 
%   19=bed load transport in streamwise direction (at nodes)
%   20=velocity at the main channel
%   21=discharge at main channel
%   22=cumulative nourished volume of sediment
%   23=suspended transport in streamwise direction
%   24=cumulative bed load transport
%   25=total sediment mass (summation of all substrate layers)
%   26=dg Fak
%   27=total sediment thickness (summation of all substrate layers)
%   28=main channel averaged bed level
%   29=sediment transport magnitude at edges m^2/s
%   30=sediment transport magnitude at edges m^3/s
%   31=morphodynamic width [m]
%   32=Chezy 
%   33=cell area [m^2]

%% do not change

simdef.flg.which_p=3; %plot type 

%% CALL

simdef=D3D_comp(simdef);
simdef=D3D_simpath(simdef);
in_read.kt=0; %output time index
out_read=D3D_read(simdef,in_read);
nTt=out_read.nTt;
nn=out_read.nn;
etab=NaN(nn,nTt);
tim=NaN(1,nTt);
for kT=1:nTt
    in_read.kt=kT; %output time index
    out_read=D3D_read(simdef,in_read);
    np=numel(out_read.z);
    etab(1:np,kT)=out_read.z;
    tim(kT)=out_read.time_r;
    fprintf('data read %4.2f %% \n',kT/nTt*100)
end
etab=etab(1:np,:);
s=out_read.SZ;

%% PLOT

cmap=brewermap(nTt,'YlOrRd');
figure
hold on
for kT=1:nTt
plot(s,etab(:,kT),'color',cmap(kT,:))
end
ylabel(out_read.zlabel)
xlabel('streamwise coordinate [m]')
han.cbar=colorbar;
han.cbar.Label.String='time [days]';
caxis([tim(1),tim(end)]/3600/24)
colormap(cmap)

