% 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

maindir  = 'd:\checkouts\openearthmodels\riverlab\testcases\f27_mor1d_tabulated_crosssections\';
caseids  = {'c21_es'; ...
            'c22_es'; ...
            'c23_es'; ...
            'c24_es'; ...
            'c25_es'; ...
            'c26_es'; ...
            'c27_es'; ...
            'c28_es'  };
datadirs = {'c21_equilibrium_slope_ds_IBedCond1_us_IBedCond4\'; ...
            'c22_equilibrium_slope_ds_IBedCond1_us_IBedCond5\'; ...
            'c23_equilibrium_slope_ds_IBedCond0_us_IBedCond4\'; ...
            'c24_equilibrium_slope_ds_IBedCond0_us_IBedCond2\'; ...
            'c25_equilibrium_slope_ds_IBedCond0_us_IBedCond3\'; ...
            'c26_equilibrium_slope_ds_IBedCond1_us_IBedCond4_morfac10\'; ...
            'c27_equilibrium_slope_ds_IBedCond1_us_IBedCond1\'; ...
            'c28_equilibrium_slope_ds_IBedCond0_us_IBedCond0\' };

for i=5
        
    datadir = [maindir datadirs{i} 'matlab'];
    caseid  = caseids{i};
    
    cd(datadir);

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

    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;

    info=ncinfo(nc_dflowfm);
    ncdisp(nc_dflowfm);

    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_mor_bl              
    %            Size:       41x121
    %            Dimensions: nmesh1d_node,time
    %            Datatype:   double
    %            Attributes:
    %                        mesh          = 'mesh1d'
    %                        location      = 'node'
    %                        coordinates   = 'mesh1d_node_x mesh1d_node_y'
    %                        standard_name = ''
    %                        long_name     = 'Time-varying bottom level in flow cell center'
    %                        units         = 'm'
    %                        grid_mapping  = 'projected_coordinate_system'
    %                        _FillValue    = -999

    %% bovenstrooms al tijdserie

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

    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')

    close (1);

    %% bed level (x) @ T=0 and T=end

    figure (2);

    subplot(2,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,1),'-xr'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,1),'-ob');
    xlabel('x [m]')
    ylabel('bed level [m]')
    %ylim([-0.05 1.05])
    title('T = 0 days')
    %legend('SOBEK prof locs','SOBEK','FM','location','best')
    legend('SOBEK','FM','location','best')

    subplot(2,1,2);
    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,end),'-xr'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,end),'-ob');
    xlabel('x [m]')
    ylabel('bed level [m]')
    %ylim([-0.05 1.05])
    title('T = END')
    %legend('SOBEK prof locs','SOBEK','FM','location','best')
    legend('SOBEK','FM','location','best')

    figname = [caseid '_fig02_bedlevel_T0_TEND'];
    x0=11; y0=0; width=20; height=25;
    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')

    close (2)
    
    %% langetermijnontwikkeling

    F40 = find(timeday_sob3==40);

    figure (3);
    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,1  ),'-xr'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,1  ),'-or');
    p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F40),'-xb'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F40),'-ob');
    p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,end),'-xg'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,end),'-og');
    xlabel('x [m]')
    ylabel('bed level [m]')
    %ylim([-0.05 1.05])
    %title('tijdontwikkeling bodemligging')
    title({'tijdontwikkeling bodemligging',datadirs{i}})
    legend('SOBEK3 initial','FM initial',...
           'SOBEK3 40 days','FM 40 days',...
           'SOBEK3 4 years','FM 4 years',...
        'location','best')
    figname = [caseid '_fig03_bedlevel_T'];
    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')

    close (3)
    
    %% initiele ontwikkeling

    F01 = find(timeday_sob3==1);
    F02 = find(timeday_sob3==2);

    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,F01),'-xr'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F01),'-or');
    p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,F02),'-xb'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,F02),'-ob');
    p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,end),'-xg'); set(p,'markersize',8);
    p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,end),'-og');
    xlabel('x [m]')
    ylabel('bed level [m]')
    %ylim([-0.05 1.05])
    title('tijdontwikkeling bodemligging')
    legend('SOBEK3 initial','FM initial',...
           'SOBEK3 1 days','FM 1 days',...
           'SOBEK3 2 days','FM 2 days',...
        'location','best')
    figname = [caseid '_fig04_bedlevel_T'];
    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')

    close (4)

    %% tijdontwikkeling bovenrand

    figure (5);

    subplot(2,1,1);
    hold on; box on; grid on;
    p=plot( timeday_sob3 , sobek3_bl    ( 1   ,:  ),'-xr'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 1   ,:  ),'-or'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    ( 2   ,:  ),'-xb'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 2   ,:  ),'-ob'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    (end-1,:  ),'-xg'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl(end-1,:  ),'-og'); set(p,'markersize',8);
    xlabel('time [days]')
    ylabel('bed level [m]')
    %ylim([-0.05 1.05])
    title('tijdontwikkeling bodemligging')
     legend('SOBEK3 node 01','FM node 01',...
            'SOBEK3 node 02','FM node 02',...
            'SOBEK3 node -1','FM node -1',...
         'location','best')

    subplot(2,1,2);
    hold on; box on; grid on;
    p=plot( timeday_sob3 , sobek3_bl    ( 1   ,:  ),'-xr'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 1   ,:  ),'-or'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    ( 2   ,:  ),'-xb'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 2   ,:  ),'-ob'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    (end-1,:  ),'-xg'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl(end-1,:  ),'-og'); set(p,'markersize',8);
    xlabel('time [days]')
    ylabel('bed level [m]')
    xlim([0 5])
    title('tijdontwikkeling bodemligging')
     legend('SOBEK3 node 01','FM node 01',...
            'SOBEK3 node 02','FM node 02',...
            'SOBEK3 node -1','FM node -1',...
         'location','best')

    figname = [caseid '_fig05_bedlevel_Ft'];
    x0=11; y0=0; width=20; height=25;
    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')

    close (5)

    %% tijdontwikkeling bovenrand

    figure (6);

    %subplot(2,1,2);
    hold on; box on; grid on;
    p=plot( timeday_sob3 , sobek3_bl    ( 1   ,:  ),'-xr'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 1   ,:  ),'-or'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    ( 2   ,:  ),'-xb'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl( 2   ,:  ),'-ob'); set(p,'markersize',8);
    p=plot( timeday_sob3 , sobek3_bl    (end-1,:  ),'-xg'); set(p,'markersize',8);
    p=plot( timeday_fm1d , mesh1d_mor_bl(end-1,:  ),'-og'); set(p,'markersize',8);
    xlabel('time [days]')
    ylabel('bed level [m]')
    xlim([0.8 1.2])
    ylim([0.4 0.6])
    title('tijdontwikkeling bodemligging')
     legend('SOBEK3 node 01','FM node 01',...
            'SOBEK3 node 02','FM node 02',...
            'SOBEK3 node -1','FM node -1',...
         'location','best')

    figname = [caseid '_fig06_bedlevel_MorStt'];
    x0=11; y0=0; width=20; height=25;
    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')

    close (6)
    
    %  =================
    %return
    %  =================

end