% Compare water depth for steady uniform flow with zero bed slope
% Andries Paarlberg
% 24-11-2017

clc;
clear all;
close all;

if isempty(which('vs_ini'))
    run('d:\checkouts\openearthtools\matlab\oetsettings.m')
end

closefigs = 0;

figdir  = '../doc/figures/';
caseid  = 'c81';

nc_sobek3   = ['..\sobek3_dimr\dflow1d\output\gridpoints.nc'];
nc_dflowfm  = ['..\DFM_OUTPUT_' caseid '\' caseid '_map.nc'];

info=ncinfo(nc_sobek3);
ncdisp(nc_sobek3);
time_sob3=ncread(nc_sobek3,'time');

sobek3_x    = ncread(nc_sobek3,'x_coordinate');
sobek3_y    = ncread(nc_sobek3,'y_coordinate');
sobek3_wd   = ncread(nc_sobek3,'water_depth');
sobek3_wl   = ncread(nc_sobek3,'water_level');
sobek3_bl   = sobek3_wl - sobek3_wd;

% inlezen op basis van data die Willem heeft gemaakt (ODSview rsegsub naar csv en vervolgens xlsx)
xlsfile = ['..\sobek3_dimr\dflow1d\output\rsegsub.xlsx'];
Tmain   = xlsread(xlsfile, 'rsegsub', 'A6:A10000');
Qmain   = xlsread(xlsfile, 'rsegsub', 'B6:AO10000');
Amain   = xlsread(xlsfile, 'rsegsub', 'IH6:JU10000');
Umain   = Qmain./Amain;
Umain   = Umain';

info=ncinfo(nc_dflowfm);
ncdisp(nc_dflowfm);
time=ncread(nc_dflowfm,'time');
time_fm1d   = ncread(nc_dflowfm,'time');

mesh1d_node_x  = ncread(nc_dflowfm,'mesh1d_node_x');
mesh1d_node_y  = ncread(nc_dflowfm,'mesh1d_node_y');
mesh1d_x       = mesh1d_node_x;
mesh1d_y       = mesh1d_node_y;
mesh1d_wd      = ncread(nc_dflowfm,'mesh1d_waterdepth');
mesh1d_wl      = ncread(nc_dflowfm,'mesh1d_s1');

mesh1d_mor_bl  = ncread(nc_dflowfm,'mesh1d_mor_bl');

mesh1d_sxtot   = squeeze(ncread(nc_dflowfm,'mesh1d_sxtot'));
mesh1d_sytot   = squeeze(ncread(nc_dflowfm,'mesh1d_sytot'));
%mesh1d_ucxq    = squeeze(ncread(nc_dflowfm,'mesh1d_ucxq'));
%mesh1d_ucyq    = squeeze(ncread(nc_dflowfm,'mesh1d_ucyq'));

mesh1d_sbcx    = squeeze(ncread(nc_dflowfm,'mesh1d_sbcx'));

mesh1d_uuu     = squeeze(ncread(nc_dflowfm,'mesh1d_uuu'));

mesh1d_umain   = mesh1d_uuu;

% cross-section data
mesh1d_mor_crs_z = squeeze(ncread(nc_dflowfm,'mesh1d_mor_crs_z'));
mesh1d_mor_crs_n = squeeze(ncread(nc_dflowfm,'mesh1d_mor_crs_n'));

% parameters needed to compute sediment transport (kg/s)
Bmain = 300;
mesh1d_czs = squeeze(ncread(nc_dflowfm,'mesh1d_czs')); 
C = mesh1d_czs(:,1);
g = 9.81;
vonkar = 0.41;
mesh1d_d50   = squeeze(ncread(nc_dflowfm,'mesh1d_DXX04'));
D50 = mesh1d_d50(:,1);

%% bovenstrooms als tijdserie

timeday_sob3 = time_sob3./(3600*24);
timeday_fm1d = time_fm1d./(3600*24);

F_01d = find(timeday_fm1d==1);
F_04d = find(timeday_fm1d==4);

figure (1);

subplot(2,2,1);
hold on; box on; grid on;
plot(timeday_sob3,sobek3_wl( 2,:),'-xr');
plot(timeday_fm1d,mesh1d_wl( 2,:),'-ob');
xlabel('time [days]')
ylabel('water level')
%ylim([5 6])
title('us')
legend('SOBEK','FM','location','best')

subplot(2,2,3);
hold on; box on; grid on;
plot(timeday_sob3,sobek3_wl(end-1,:),'-xr');
plot(timeday_fm1d,mesh1d_wl(end-1,:),'-ob');
xlabel('time [days]')
ylabel('water level')
%ylim([4 5])
title('ds')

subplot(2,2,2);
hold on; box on; grid on;
plot(timeday_sob3,sobek3_wl( 2,:),'-xr');
plot(timeday_fm1d,mesh1d_wl( 2,:),'-ob');
xlabel('time [days]')
ylabel('water level')
%ylim([5.6 5.7])
xlim([0 5])
title('us')
legend('SOBEK','FM','location','best')

subplot(2,2,4);
hold on; box on; grid on;
plot(timeday_sob3,sobek3_wl(end-1,:),'-xr');
plot(timeday_fm1d,mesh1d_wl(end-1,:),'-ob');
xlabel('time [days]')
ylabel('water level')
%ylim([4.6 4.7])
xlim([0 5])
title('ds')

figname = [caseid '_fig01_waterlevel_us_ds_Ft'];
x0=11; y0=0; width=16.5; height=16.5;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (1);
end

%% water level(x) after 1 day and @ end of simulation

figure (2);
subplot(2,1,1);
hold on; box on; grid on;
plot(sobek3_x( 2:10),sobek3_wl( 2:10,F_01d),'-xr');
plot(mesh1d_x( 1:10),mesh1d_wl( 1:10,F_01d),'-xb');
xlabel('x [m]')
ylabel('water level')
title('after 1 day')
legend('SOBEK','FM','location','best')

subplot(2,1,2);
hold on; box on; grid on;
plot(sobek3_x( 2:10),sobek3_wl( 2:10,end),'-xr');
%plot(mesh1d_x( 1:10),mesh1d_wl( 1:10,end),'-xb');
xlabel('x [m]')
ylabel('water level')
title('end of simulation')
%legend('SOBEK','FM','location','best')
legend('FM','location','best')

figname = [caseid '_fig02_waterlevel_x'];
x0=11; y0=0; width=16.5; height=20.0;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (2);
end

%% Umain

sobek_Umain_x = 50:50:2000;
ylims = [0.2 1];

% 25    49    73    97
lw=1.5;
figure (3);
subplot(2,1,1);
hold on; grid on; box on;
p=plot(mesh1d_node_x,-mesh1d_umain(:,F_01d),'-+'); set(p,'linewidth',lw)
p=plot(sobek_Umain_x(2:end),Umain(2:end,F_01d)       ,'-o'); set(p,'linewidth',lw)
xlabel('x (m)')
ylabel('Umain (m/s)')
ylim(ylims)
legend('fm1d','sobek3\_dimr','location','northwest')
title('after spin up')

subplot(2,1,2);
hold on; grid on; box on;
p=plot(mesh1d_node_x,-mesh1d_umain(:,end),'-+'); set(p,'linewidth',lw)
%p=plot(sobek_Umain_x(2:end),Umain(2:end,end)       ,'-o'); set(p,'linewidth',lw)
xlabel('x (m)')
ylabel('Umain (m/s)')
ylim(ylims)
%legend('fm1d','sobek3\_dimr','location','northwest')
legend('fm1d','location','northwest')
title('end simulation')

figname = [caseid '_fig03_Umain'];
x0=11; y0=0; width=16.5; height=20.0;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (3);
end

%% langetermijnontwikkeling

ylims = [-0.5 +0.5];

T1 = 30;
T2 = 200;
T3 = 730;
F1 = find(timeday_fm1d==T1);
F2 = find(timeday_fm1d==T2);
F3 = find(timeday_fm1d==T3);

figure (4);
subplot(1,1,1);
hold on; box on; grid on;
%p=plot( profdata(:,1) + sobek3_x(1) , profdata(:,2)   ,' sg'); set(p,'markerfacecolor','g','markeredgecolor','k','markersize',6);
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F1),'-xr'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F1),'-or');
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F2),'-xb'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F2),'-ob');
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F3),'-xg'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F3),'-og');
xlabel('x [m]')
ylabel('bed level [m]')
ylim(ylims)
title('tijdontwikkeling bodemligging lange termijn')
%legend('SOBEK3 initial','FM initial',...
%       ['SOBEK3 ' num2str(T2) ' days'],['FM ' num2str(T2) ' days'],...
%       ['SOBEK3 ' num2str(T3) ' days'],['FM ' num2str(T3) ' days'],...
legend(['FM ' num2str(T1) ' days'],...
       ['FM ' num2str(T2) ' days'],...
       ['FM ' num2str(T3) ' days'],...
    'location','best')

figname = [caseid '_fig04_bedlevel_T_lang'];
x0=11; y0=0; width=20; height=15;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (4);
end

%% initiele ontwikkeling

ylims = [-0.5 +0.5];

T1 = 0;
T2 = 1;
T3 = 5;
F1 = find(timeday_fm1d==T1);
F2 = find(timeday_fm1d==T2);
F3 = find(timeday_fm1d==T3);

figure (5);
subplot(1,1,1);
hold on; box on; grid on;
%p=plot( profdata(:,1) + sobek3_x(1) , profdata(:,2)   ,' sg'); set(p,'markerfacecolor','g','markeredgecolor','k','markersize',6);
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F1),'-xr'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F1),'-or');
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F2),'-xb'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F2),'-ob');
%p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F3),'-xg'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F3),'-og');
xlabel('x [m]')
ylabel('bed level [m]')
ylim(ylims)
title('tijdontwikkeling bodemligging initieel')
%legend('SOBEK3 initial','FM initial',...
%       ['SOBEK3 ' num2str(T2) ' days'],['FM ' num2str(T2) ' days'],...
%       ['SOBEK3 ' num2str(T3) ' days'],['FM ' num2str(T3) ' days'],...
legend(['FM ' num2str(T1) ' days'],...
       ['FM ' num2str(T2) ' days'],...
       ['FM ' num2str(T3) ' days'],...
    'location','best')

figname = [caseid '_fig05_bedlevel_T_initieel'];
x0=11; y0=0; width=20; height=15;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (5);
end

%% profielen

%mesh1d_mor_crs_z % 19x41x2401
% npoints in prof   x   grid points    x    time
% 1 en 2 zijn leeg!!
% zijn dat up en dwn?

profloc = [0 2000 50:50:1950];

dx = 50;

T1 = 0;
T2 = 10;
T3 = 100;
T4 = 730;
F1 = find(timeday_fm1d==T1);
F2 = find(timeday_fm1d==T2);
F3 = find(timeday_fm1d==T3);
F4 = find(timeday_fm1d==T4);

if ishandle(6); close (6); end
figure (6);

ylims = [-1 +5];

x=50;
n=3;
subplot(3,1,1);
hold on; grid on; box on;
profn = squeeze(mesh1d_mor_crs_n(:,n,:));
profz = squeeze(mesh1d_mor_crs_z(:,n,:));
plot( (profn(:,F1)),profz(:,F1),'-ok');
plot( (profn(:,F2)),profz(:,F2),'-+r');
plot( (profn(:,F3)),profz(:,F3),'-dg');
plot( (profn(:,F4)),profz(:,F4),'-^b');
plot(-(profn(:,F1)),profz(:,F1),'-ok');
plot(-(profn(:,F2)),profz(:,F2),'-+r');
plot(-(profn(:,F3)),profz(:,F3),'-dg');
plot(-(profn(:,F4)),profz(:,F4),'-^b');
p=plot([-Bmain/2 -Bmain/2],ylims,'k'); set(p,'LineWidth',2);
p=plot([+Bmain/2 +Bmain/2],ylims,'k'); set(p,'LineWidth',2);
xlabel('x [m]')
ylabel('bed level [m AD]')
title(['x = ' num2str(x) ' m'])
ylim(ylims)
legend(['FM ' num2str(T1) ' days'],...
       ['FM ' num2str(T2) ' days'],...
       ['FM ' num2str(T3) ' days'],...
       ['FM ' num2str(T4) ' days'],...
    'location','best')

x=500;
n=x/dx+2;
subplot(3,1,2);
hold on; grid on; box on;
profn = squeeze(mesh1d_mor_crs_n(:,n,:));
profz = squeeze(mesh1d_mor_crs_z(:,n,:));
plot( (profn(:,F1)),profz(:,F1),'-ok');
plot( (profn(:,F2)),profz(:,F2),'-+r');
plot( (profn(:,F3)),profz(:,F3),'-dg');
plot( (profn(:,F4)),profz(:,F4),'-^b');
plot(-(profn(:,F1)),profz(:,F1),'-ok');
plot(-(profn(:,F2)),profz(:,F2),'-+r');
plot(-(profn(:,F3)),profz(:,F3),'-dg');
plot(-(profn(:,F4)),profz(:,F4),'-^b');
p=plot([-Bmain/2 -Bmain/2],ylims,'k'); set(p,'LineWidth',2);
p=plot([+Bmain/2 +Bmain/2],ylims,'k'); set(p,'LineWidth',2);
xlabel('x [m]')
ylabel('bed level [m AD]')
title(['x = ' num2str(x) ' m'])
ylim(ylims)

x=1000;
n=x/dx+2;
subplot(3,1,3);
hold on; grid on; box on;
profn = squeeze(mesh1d_mor_crs_n(:,n,:));
profz = squeeze(mesh1d_mor_crs_z(:,n,:));
plot( (profn(:,F1)),profz(:,F1),'-ok');
plot( (profn(:,F2)),profz(:,F2),'-+r');
plot( (profn(:,F3)),profz(:,F3),'-dg');
plot( (profn(:,F4)),profz(:,F4),'-^b');
plot(-(profn(:,F1)),profz(:,F1),'-ok');
plot(-(profn(:,F2)),profz(:,F2),'-+r');
plot(-(profn(:,F3)),profz(:,F3),'-dg');
plot(-(profn(:,F4)),profz(:,F4),'-^b');
p=plot([-Bmain/2 -Bmain/2],ylims,'k'); set(p,'LineWidth',2);
p=plot([+Bmain/2 +Bmain/2],ylims,'k'); set(p,'LineWidth',2);
xlabel('x [m]')
ylabel('bed level [m AD]')
title(['x = ' num2str(x) ' m'])
ylim(ylims)

figname = [caseid '_fig06_profielen'];
x0=11; y0=0; width=16.5; height=16.5;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (6);
end

%% profielen

%mesh1d_mor_crs_z % 19x41x2401
% npoints in prof   x   grid points    x    time
% 1 en 2 zijn leeg!!
% zijn dat up en dwn?

profloc = [0 2000 50:50:1950];

dx = 50;

T1 = 0;
T2 = 50;
T3 = 100;
T4 = 730;
F1 = find(timeday_fm1d==T1);
F2 = find(timeday_fm1d==T2);
F3 = find(timeday_fm1d==T3);
F4 = find(timeday_fm1d==T4);

if ishandle(7); close (7); end
figure (7);

ylims = [-1 +5];

x=50;
n=x/dx+2;
hold on; grid on; box on;
profn = squeeze(mesh1d_mor_crs_n(:,n,:));
profz = squeeze(mesh1d_mor_crs_z(:,n,:));
plot( (profn(:,F1)),profz(:,F1),'-ok');
plot( (profn(:,F2)),profz(:,F2),'-+r');
plot( (profn(:,F3)),profz(:,F3),'-dg');
plot( (profn(:,F4)),profz(:,F4),'-^b');
plot(-(profn(:,F1)),profz(:,F1),'-ok');
plot(-(profn(:,F2)),profz(:,F2),'-+r');
plot(-(profn(:,F3)),profz(:,F3),'-dg');
plot(-(profn(:,F4)),profz(:,F4),'-^b');
p=plot([-Bmain/2 -Bmain/2],ylims,'k'); set(p,'LineWidth',2);
p=plot([+Bmain/2 +Bmain/2],ylims,'k'); set(p,'LineWidth',2);
xlabel('x [m]')
ylabel('bed level [m AD]')
title(['x = ' num2str(x) ' m'])
ylim(ylims)
legend(['FM ' num2str(T1) ' days'],...
       ['FM ' num2str(T2) ' days'],...
       ['FM ' num2str(T3) ' days'],...
       ['FM ' num2str(T4) ' days'],...
    'location','best')

figname = [caseid '_fig07_profiel'];
x0=11; y0=0; width=25; height=16.5;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
%print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

if closefigs
    close (7);
end
