clear;close all;clc; 

addpath('/Users/xudo627/donghui/CODE/Setup-E3SM-Mac/matlab-scripts-for-mosart/');
addpath('/Users/xudo627/donghui/mylib/m/');
addpath('functions/');

re = 6.37122e6;% Earth radius
datadir = '/Volumes/LaCie/PROJECTS/Pantanal/';
yrs = 1979;
yre = 2010;

show_domain      = 0;
show_sub_basin   = 0;
show_inund       = 0;
show_default_BFI = 0;
show_annual      = 0;

update = 0;

[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);
[gsim,gsim2] = search_GSIM_gauge(Sh.X,Sh.Y,-57.49,-18.04,show_domain);
[lon_gsim2,lat_gsim2,S2,yy,mm,da2,mu_gsim2,sd2,cv2] = get_GSIM_discharge(gsim,1);

forcing = 'CRUNCEPv7';%'GSWP3v1';%'CRUNCEPv7';
ex      = 'sub'; % sub, sur2
tag = ['outputs_' ex '_CLM' forcing];
area = areaint(Sh.Y,Sh.X)*4*pi.*(re^2); 
area2 = areaint(S2.Y,S2.X)*4*pi.*(re^2); 

xv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xv');
yv=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yv');
xc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','xc');
yc=ncread('inputdata_run/domain_lnd_Pantanal_all_template.nc','yc');

numc = length(xc);

yr1 = 1979;
yr2 = 2010;
if exist(['data/' tag '.mat'],'file') && ~update
    load(['data/' tag '.mat']);
else
    % Read MOSART
    [ntot, q, fdfall, fdf, R2, NSE, KGE, ind1, ind2, imax]                    = ...
             read_mosart(datadir, tag, yr1, yr2, mu_gsim, lon_gsim, lat_gsim, yr, area);
    for i = 1 : ntot
        fdf(i,:) = nanmean(fdfall((i-1)*numc+1:i*numc,:),1);
    end
    % Read ELM
    [fh2osfc,qrunoff,qdrai,fsat,fh2osfcall,qrunoffall,qdraiall,zwt,zwtall,twsall,et,etall,rain] = ...
             read_elm(datadir, tag, ntot, numc, ind2);  
    save(['data/' tag '.mat'], '-regexp', '^(?!(show_domain|show_sub_basin||show_inund|show_default_BFI|show_annual)$).');
end

yr1 = 1979;
yr2 = 2009;
[tmon] = get_time_number(yr1,yr2);
[xtrain,nr_uq] = get_parameters(ex,numc,ntot);
[muyr,qyr,z,zmu,trend,trendmu,slope,slopemu] = convert_mon_to_annual(mu_gsim,q,ntot,yr,yr1,yr2);

q = q(:,(yr1-yrs)*12+1:(yr2-yrs+1)*12);
qrunoff = qrunoff(:,(yr1-yrs)*12+1:(yr2-yrs+1)*12);
if strcmp(ex,'man') || strcmp(ex,'man2')
    figure; set(gcf,'Position',[10 10 900 300]);
%     for i = 1 : 50
%         c = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
%         if i == 21
%             subplot_tight(2,2,1,[0.08 0.06]);
%             plot(c,'k--','LineWidth',3); hold on; grid on;
%             title('Subnetwork manning coefficient','FontSize',15,'FontWeight','bold');
%         elseif i == 41
%             subplot_tight(2,2,3,[0.08 0.06]);
%             plot(c,'k--','LineWidth',3); hold on; grid on;
%             title('hillslope manning coefficient','FontSize',15,'FontWeight','bold');
%         end
%         tmp = nanmean(reshape(q(i,:)',[12,yr2-yr1+1]),2);
%         if i >= 21
%             plot(tmp,'-','LineWidth',0.5); 
%         end
%         
%     end
%     
    subplot_tight(1,2,1,[0.08 0.08]);
    plot(0.01:0.01:0.2,sqrt(R2(21:40)),'ko','MarkerSize',10,'LineWidth',1.5); hold on;
    plot([0.01 0.2],[sqrt(R2(25)) sqrt(R2(25))],'r-','LineWidth',2);grid on; ylim([0.1 0.25]);
    add_title(gca,'(a). Subnetwork Manning coefficient');
    set(gca,'FontSize',13);
    ylabel('Correlation coefficient','FontSize',15,'FontWeight','bold');
    legend('Simulations','Default','FontSize',14,'FontWeight','bold');
    
    subplot_tight(1,2,2,[0.08 0.08]);
    plot(0.05:0.05:0.5,sqrt(R2(41:50)),'ko','MarkerSize',10,'LineWidth',1.5); hold on;
    plot([0.05 0.5],[sqrt(R2(25)) sqrt(R2(25))],'r-','LineWidth',2);grid on; ylim([0.1 0.25]);
    add_title(gca,'(b). Hillslope Manning coefficient');
    set(gca,'FontSize',13);
    
end
ind2 = ind2(1);
% PLOT 
figure; set(gcf,'Position',[10 10 1600 400]);
p(1) = plot(tmon,mu_gsim(yr >= yr1 & yr <= yr2),'k-','LineWidth',3); hold on; grid on;
%p(2) = plot(tmon,q(ind1,:),'b--','LineWidth',3);
p(2) = plot(tmon,q(ind2,:),'b--','LineWidth',3);
xlim([tmon(1) tmon(end)]);
datetick('x','mmmyy','keeplimits');
set(gca,'FontSize',14);
ylabel('Discharge [m^{3}/s]','FontSize',16,'FontWeight','bold');
leg = legend(p,{'Observation','Simulation'});
leg.FontSize = 18; leg.FontWeight = 'bold';
ax = gca;
ax.XAxis.TickValues = tmon(1:12:end);
ax.XAxis.TickLabelRotation = 45;
[metric(1).R2,metric(1).RMSE,metric(1).NSE,metric(1).PBIAS,metric(1).MSE,metric(1).NSE1] = ...
      estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2),q(ind2,:)');
[metric(2).R2,metric(2).RMSE,metric(2).NSE,metric(2).PBIAS,metric(2).MSE,metric(2).NSE1] = ...
      estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2-1),q(ind2,1:end-12)');
strs{1} = {['PBIAS = ' num2str(metric(1).PBIAS) '%'],     ...
           ['\rho = ' num2str(round(sqrt(metric(1).R2),2))], ...
           ['NSE  = ' num2str(round(metric(1).NSE,2))]};
t = annotation('textbox',[0.65 0.7 0.3 0.2],'String',strs{1},       ... 
               'FitBoxToText','on','EdgeColor','none','Color','b',     ...
               'FontSize',16,'FontWeight','bold'); 

cmaps = getPyPlot_cMap('rainbow',10);
cmaps2 = getPyPlot_cMap('rainbow',256);
if strcmp(ex,'geo')
    w = 0.10;
    d = 0.15;
%     w = 0.08;
%     d = 0.09;
    man2 = load(['data/outputs_man2_CLM' forcing '.mat']);
    figure; set(gcf,'Position',[10 10 800 800]);
    %figure; set(gcf,'Position',[10 10 1200 800]);
    ax(1) = subplot_tight(2,1,1,[w d]);
    ssize = [50 75 100 125 150 175 200];
    rSD   = nanstd(q,[],2)./nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
    cstd  = [0.85 : 0.1 : 1.45];
    for i = 1 : length(cstd)
        if i == 1
            ii = find(rSD <= cstd(i));
        else
            ii  = find(rSD <= cstd(i) & rSD > cstd(i-1));
        end
            
        scatter(xtrain(ii,5),xtrain(ii,6),ssize(i),sqrt(R2(ii)),'filled', ...
                'MarkerEdgeColor','k','MarkerFaceAlpha',0.7,'LineWidth',1.5); hold on; 
    end
    colormap(gca,cmaps); 
    %plot(xtrain(:,5),xtrain(:,6),'k.'); hold on;
    set(gca,'FontSize',13); grid on;
    cb1 = colorbar('Location','west'); cb1.FontSize = 13;
    caxis([0.075 0.575]);
    xlim([2.6 20.2]); ylim([0.12 0.63]);
    xlabel('a_{w}','FontSize',15,'FontWeight','bold');
    ylabel('a_{d}','FontSize',15,'FontWeight','bold');
    
    ax(2) = subplot_tight(2,1,2,[w d]);
    tmp = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
    p(1) = plot(tmp,'k--','LineWidth',3); hold on; grid on;
    nr_prior = 0.01 : 0.01 : 0.2;
    for i = 2 : 2 : 20
        tmp = nanmean(reshape(man2.q(i,1:end-12)',[12,yr2-yr1+1]),2);
        plot(tmp,'-','Color',cmaps(i/2,:),'LineWidth',2);
    end
    xlim([1 12]);
    xticks([1 : 12]);
    xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
    cb2 = colorbar('XTickLabel',num2cell(nr_prior(2:2:20)),'XTick', 0.05 : 0.1 : 0.95,'Location','west');
    set(gca,'FontSize',13); 
    colormap(gca,cmaps);cb2.FontSize = 13;
    ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');
    legend(p,{'Observation'},'FontSize',14,'FontWeight','bold');
    
    add_title(ax(1),'(a). River geometry',18);
    add_title(ax(2),'(b). Main channel Manning coefficient',18);
    
    ax(1).Position(1) = ax(1).Position(1) - 0.025;
    ax(2).Position(1) = ax(2).Position(1) - 0.025;
    ax(1).Position(3) = ax(1).Position(3) - 0.2;
    ax(2).Position(3) = ax(2).Position(3) - 0.2;
    
    cb1.Position = [ax(1).Position(1) + ax(1).Position(3) + 0.01 ax(1).Position(2) ...
                    0.02 ax(1).Position(4)];
    cb2.Position = [ax(2).Position(1) + ax(2).Position(3) + 0.01 ax(2).Position(2) ...
                    0.02 ax(2).Position(4)];
    cb1.Label.String = 'Correlation coefficient';
    cb1.Label.FontSize = 15; cb1.Label.FontWeight = 'bold';
    cb2.Label.String = 'Manning coefficient';
    cb2.Label.FontSize = 15; cb2.Label.FontWeight = 'bold';
    
    ax(3) = axes('Position',[ax(1).Position(1) + ax(1).Position(3) + 0.1 ax(1).Position(2) 0.2 ax(1).Position(4)]);
    for i = 1 : length(ssize)
        scatter(0.15, 0.1 + 0.1*i,ssize(i),'MarkerFaceColor','none','MarkerEdgeColor','k','LineWidth',1.5); hold on;
        if i == 1
            str = ['0.80 < rSD <= ' num2str(cstd(i))];
            text(0.25, 0.1 + 0.1*i,str,'FontSize',14);
        else
            str = [num2str(cstd(i-1)) ' < rSD <= ' num2str(cstd(i))];
            text(0.25, 0.1 + 0.1*i,str,'FontSize',14);
        end
    end
    ylim([0 1]); xlim([0 1]);
    set(ax(3),'YTick',[],'Xtick',[]);ax(3).XColor = 'w';ax(3).YColor = 'w';
    
    ax(4) = axes('Position',[ax(2).Position(1) + ax(2).Position(3) + 0.1 ax(2).Position(2) 0.2 ax(2).Position(4)]);
    rSD = nanstd(man2.q(1:20,1:end-12),[],2)./nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
    plot(rSD,0.01 : 0.01 : 0.2,'ko','MarkerSize',10,'LineWidth',1.5); grid on; hold on;
    plot([1 1],[0.01 0.21],'-','Color',[0.5 0.5 0.5],'LineWidth',2);
    ylim([0.01 0.21]);
    yticklabels(""); set(ax(4),'FontSize',13);
    xticks([0.5 : 0.25 : 1.5]);
    xlabel('rSD','FontSize',15,'FontWeight','bold');
elseif strcmp(ex,'sub')
    w = 0.08;
    d = 0.09;
    figure; set(gcf,'Position',[10 10 1200 800]);
    ax(1) = subplot_tight(2,2,1,[w d]);
    fd = [0.25 0.5 0.75 1 1.5 2 2.5 3 4 5];
    tmp = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
    p1(1) = plot(tmp,'k--','LineWidth',3);hold on; grid on;
    for i = 1 : 10
        dist = (xtrain(:,1) - fd(i)).^2';
        ind = find(dist == min(dist));
        tmp = nanmean(reshape(q(ind,:)',[12,yr2-yr1+1]),2);
        plot(1:12,tmp,'-','Color',cmaps(i,:),'LineWidth',2); 
    end
    xlim([1 12]);
    xticks([1 : 12]);
    xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
    cb = colorbar('XTickLabel',num2cell(fd),'XTick', 0.05 : 0.1 : 0.95,'Location','west');
    colormap(gca,cmaps);
    cb.Position = [ax(1).Position(1) + ax(1).Position(3) + 0.01 ...
                   ax(1).Position(2)                            ...
                   0.01 ax(1).Position(4)];
    cb.FontSize = 13;
    cb.Label.String = 'f_{drain}';
    cb.Label.FontSize = 15; cb.Label.FontWeight = 'bold';
    set(gca,'FontSize',13);
    ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');
    leg1 = legend(p1,{'Observation'}); leg1.FontSize = 14;
    
    ax(2) = subplot_tight(2,2,2,[w d]);
    BFI_sims = NaN(numc,ntot);
    for i = 1 : ntot
        BFI_sims(:,i) = nanmean(qdraiall((i-1)*numc+1:i*numc,:),2)./nanmean(qrunoffall((i-1)*numc+1:i*numc,:),2);
    end
    BFI_sims(BFI_sims > 1) = 1;
    BFI_sims(BFI_sims < 0) = NaN;
    plot(xtrain(:,1),nanmean(BFI_sims,1),'ko','MarkerSize',10,'LineWidth',1.5); hold on; grid on;
    p2(1) = plot([0.01 5],[0.69 0.69],'k--','LineWidth',2);
    set(gca,'FontSize',13,'YAxisLocation','right');
    xlabel('f_{drain}','FontSize',15,'FontWeight','bold');
    ylabel('Base flow index','FontSize',15,'FontWeight','bold');
    leg2 = legend(p2,{'GSCD'}); leg2.FontSize = 14;
    
    ax(3) = subplot_tight(2,2,3,[w d]);
    yyaxis left;
    plot(xtrain(:,1),nanmean(q,2),'bo','MarkerSize',10,'LineWidth',1); hold on; grid on;
    mu_mu = nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
    mu_cv = nanstd(mu_gsim(yr >= yr1 & yr <= yr2))/nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
    mu_sd = nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
    p3(1) = plot([0 5],[mu_mu mu_mu],'b--','LineWidth',1.5);
    set(gca,'YScale','log');
    ylabel('Mean [m^{3}/s]','FontSize',15,'FontWeight','bold');

    yyaxis right;
    plot(xtrain(:,1),nanstd(q,[],2)./mu_sd,'rs','MarkerSize',10,'LineWidth',1); hold on
    p3(2) = plot([0 5],[1 1],'r--','LineWidth',1.5);
    set(gca,'FontSize',13);
    xlabel('f_{drain}','FontSize',15,'FontWeight','bold');
    ylabel('rSD [-]','FontSize',15,'FontWeight','bold');
    leg3 = legend(p3,{'Observed mean','rSD = 1'}); leg3.FontSize = 14;
    
    ax(4) = subplot_tight(2,2,4,[w d]);
    isig = find(trend ~= 0);
    inosig = find(trend == 0);
    plot(xtrain(inosig,1),slope(inosig)'./nanmean(q(inosig,:),2).*100,'ks','MarkerSize',10,'LineWidth',1.5); hold on; grid on;
    plot(xtrain(isig,1),slope(isig)'./nanmean(q(isig,:),2).*100,'gd','MarkerSize',10,'LineWidth',1.5); hold on; grid on;
    plot([0.01 5],[slopemu slopemu]./mu_mu.*100,'k--','LineWidth',2);
    legend('p-value > 0.05','p-value < 0.05','Observation');
    set(gca,'FontSize',13,'YAxisLocation','right');
    xlabel('f_{drain}','FontSize',15,'FontWeight','bold');
    ylabel('Sen slope [% yr^{-1}]','FontSize',15,'FontWeight','bold');
    ylim([-2.5 0.5]);
    
    add_title(ax(1),'(a)',18);
    add_title(ax(2),'(b)',18);
    add_title(ax(3),'(c)',18);
    add_title(ax(4),'(d)',18);

elseif strcmp(ex,'sur')
    w = 0.08;
    d = 0.09;
    figure; set(gcf,'Position',[10 10 1200 800]);
    ax(1) = subplot_tight(2,2,1,[w d]);
    fo = 0 : 0.02 : 0.16;
    cmaps = jet(length(fo));
    tmp = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
    p1 = plot(tmp,'k--','LineWidth',3);hold on; grid on;
    for i = 1 : length(fo)
        %dist = (xtrain(:,2) - fo(i)).^2';
        dist = (nanmean(fh2osfc,2) - fo(i)).^2;
        ind = find(dist == min(dist));
        disp(['fh2osfc = ' num2str(nanmean(fh2osfc(ind,:),2))]);
        tmp = nanmean(reshape(q(ind,:)',[12,yr2-yr1+1]),2);
        plot(1:12,tmp,'-','Color',cmaps(i,:),'LineWidth',2); 
    end
    xlim([1 12]);
    xticks([1 : 12]);
    xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
    cb = colorbar('XTickLabel',num2cell(fo),'XTick', 0+1/18 : 1/9 : 1 - 1/18,'Location','west');
    colormap(gca,cmaps);
    cb.Position = [ax(1).Position(1) + ax(1).Position(3) ...
                   ax(1).Position(2) ...
                   0.01 ax(1).Position(4)];
    cb.FontSize = 13;
    cb.Label.String = 'Wetland inundation';
    cb.Label.FontSize = 15; cb.Label.FontWeight = 'bold';
    leg1 = legend(p1,{'Observation'}); leg1.FontSize = 14;
    
    set(gca,'FontSize',13);
    ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');

    ax(2) = subplot_tight(2,2,2,[w d]);
    BFI_sims = NaN(numc,ntot);
    for i = 1 : ntot
        BFI_sims(:,i) = nanmean(qdraiall((i-1)*numc+1:i*numc,:),2)./nanmean(qrunoffall((i-1)*numc+1:i*numc,:),2);
    end
    BFI_sims(BFI_sims > 1) = 1;
    BFI_sims(BFI_sims < 0) = NaN;
    plot(nanmean(fh2osfc,2),nanmean(BFI_sims,1),'ko','MarkerSize',10,'LineWidth',1); hold on; grid on;
    p2 = plot([min(nanmean(fh2osfc,2)) 0.18],[0.69 0.69],'k--','LineWidth',2);
    set(gca,'FontSize',13,'YAxisLocation','right');
    xlabel('SWF','FontSize',15,'FontWeight','bold');
    ylabel('Base flow index','FontSize',15,'FontWeight','bold');
    leg2 = legend(p2,{'GSCD'}); leg2.FontSize = 14;
    
    ax(3) = subplot_tight(2,2,3,[w d]);
    yyaxis left;
    plot(nanmean(fh2osfc,2),nanmean(q,2),'bo','MarkerSize',10,'LineWidth',1); hold on; grid on;
    mu_mu = nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
    mu_cv = nanstd(mu_gsim(yr >= yr1 & yr <= yr2))/nanmean(mu_gsim(yr >= yr1 & yr <= yr2));
    p3(1) = plot([min(nanmean(fh2osfc,2)) 0.18],[mu_mu mu_mu],'b--','LineWidth',1.5);
    set(gca,'YScale','log');
    ylabel('Mean [m^{3}/s]','FontSize',15,'FontWeight','bold');

    yyaxis right;
    plot(nanmean(fh2osfc,2),nanstd(q,[],2)./nanmean(q,2),'rs','MarkerSize',10,'LineWidth',1); hold on
    p3(2) = plot([min(nanmean(fh2osfc,2)) 0.18],[mu_cv mu_cv],'r--','LineWidth',1.5);
    set(gca,'FontSize',13);
    xlabel('SWF','FontSize',15,'FontWeight','bold');
    ylabel('Coefficient of variation [-]','FontSize',15,'FontWeight','bold');
    leg3 = legend(p3,{'Observed mean','Observed coefficient of variation'}); leg3.FontSize = 14;
    
    ax(4) = subplot_tight(2,2,4,[w d]);
    isig = find(trend ~= 0);
    inosig = find(trend == 0);
%     plot(xtrain(inosig,4),slope(inosig),'ks','MarkerSize',10,'LineWidth',1); hold on; grid on;
%     plot(xtrain(isig,4),slope(isig),'gd','MarkerSize',10,'LineWidth',1.5); hold on; grid on;
    plot(nanmean(fh2osfc(inosig,:),2),slope(inosig),'ks','MarkerSize',10,'LineWidth',1); hold on; grid on;
    plot(nanmean(fh2osfc(isig,:),2),slope(isig),'gd','MarkerSize',10,'LineWidth',1.5); hold on; grid on;
    plot([min(nanmean(fh2osfc,2)) 0.18],[slopemu slopemu],'k--','LineWidth',2);
    legend('p-value > 0.05','p-value < 0.05','Observation');
    set(gca,'FontSize',13,'YAxisLocation','right');
    xlabel('SWF','FontSize',15,'FontWeight','bold');
    ylabel('Sen slope [m^{3} s^{-1} yr^{-1}x]','FontSize',15,'FontWeight','bold');

    add_title(ax(1),'(a)',18);
    add_title(ax(2),'(b)',18);
    add_title(ax(3),'(c)',18);
    add_title(ax(4),'(d)',18);

elseif strcmp(ex,'sur2')
    
    cmaps = getPyPlot_cMap('rainbow',9);
    
    inund = nanmean(fh2osfc,2);
    
    w = 0.08;
    d = 0.12;
    figure; set(gcf,'Position',[10 10 1200 800]);
    ax(1) = subplot_tight(2,2,1,[w d]);
    scatter(xtrain(:,2),xtrain(:,4),100,inund,'filled', ...
            'MarkerEdgeColor','k','MarkerFaceAlpha',0.7,'LineWidth',1.5); colormap(gca,cmaps); 
    caxis([0.01 0.19]);hold on;
    set(gca,'FontSize',13); grid on;
    cb1 = colorbar('Location','west'); cb1.FontSize = 13;
    xlabel('f_{over}','FontSize',15,'FontWeight','bold');
    ylabel('f_{c}','FontSize',15,'FontWeight','bold');

    ax(2) = subplot_tight(2,2,2,[w d]);
    fo = 0.02 : 0.02 : 0.18;

    tmp = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
    p1 = plot(tmp,'k--','LineWidth',3);hold on; grid on;
    for i = 1 : length(fo)
        %dist = (xtrain(:,2) - fo(i)).^2';
        dist = (nanmean(fh2osfc,2) - fo(i)).^2;
        ind = find(dist == min(dist));
        disp(['fh2osfc = ' num2str(nanmean(fh2osfc(ind,:),2))]);
        tmp = nanmean(reshape(q(ind,:)',[12,yr2-yr1+1]),2);
        plot(1:12,tmp,'-','Color',cmaps(i,:),'LineWidth',2); 
    end
    xlim([1 12]);
    xticks([1 : 12]);
    xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
    cb2 = colorbar('XTickLabel',num2cell(fo),'XTick', 0+1/18 : 1/9 : 1 - 1/18,'Location','west');
    colormap(gca,cmaps);
    leg1 = legend(p1,{'Observation'}); leg1.FontSize = 14;
    
    set(gca,'FontSize',13);
    ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');
    
    ax(3) = subplot_tight(2,2,3,[w d]);
    rSD = nanstd(q,[],2)./nanstd(mu_gsim(yr >= yr1 & yr <= yr2));
    scatter(xtrain(:,2),xtrain(:,4),100,nanmean(q,2),'filled', ...
            'MarkerEdgeColor','k','MarkerFaceAlpha',0.7,'LineWidth',1.5); colormap(gca,cmaps); 
    caxis([1250 5750]); hold on;
    set(gca,'FontSize',13); grid on;
    cb3 = colorbar('Location','west'); cb3.FontSize = 13;
    xlabel('f_{over}','FontSize',15,'FontWeight','bold');
    ylabel('f_{c}','FontSize',15,'FontWeight','bold');
    

    ax(4) = subplot_tight(2,2,4,[w d]);
    scatter(xtrain(:,2),xtrain(:,4),100,slope'./nanmean(q,2).*100,'filled', ...
            'MarkerEdgeColor','k','MarkerFaceAlpha',0.7,'LineWidth',1.5); colormap(gca,cmaps); 
    caxis([-1.25 -0.35]);hold on;
    set(gca,'FontSize',13); grid on;
    cb4 = colorbar('Location','west'); cb4.FontSize = 13;
    cb4.Ticks = [-1.2 : 0.1 : -0.4];
    xlabel('f_{over}','FontSize',15,'FontWeight','bold');
    ylabel('f_{c}','FontSize',15,'FontWeight','bold');
  
    cb1.Position = [ax(1).Position(1) + ax(1).Position(3) + 0.01 ax(1).Position(2) ...
                    0.01 ax(1).Position(4)];
    cb2.Position = [ax(2).Position(1) + ax(2).Position(3) + 0.01 ax(2).Position(2) ...
                    0.01 ax(2).Position(4)];
    cb3.Position = [ax(3).Position(1) + ax(3).Position(3) + 0.01 ax(3).Position(2) ...
                    0.01 ax(3).Position(4)];
    cb4.Position = [ax(4).Position(1) + ax(4).Position(3) + 0.01 ax(4).Position(2) ...
                    0.01 ax(2).Position(4)];
    cb1.Label.String = 'Wetland inundation [-]'; set(cb1,'YAxisLocation','right');
    cb1.Label.FontSize = 15; %cb1.Label.FontWeight = 'bold';
    cb2.Label.String = 'Wetland inundation [-]';
    cb2.Label.FontSize = 15; %cb2.Label.FontWeight = 'bold';
    cb3.Label.String = 'Annual streamflow [m^{3}/s]';
    cb3.Label.FontSize = 15; %cb3.Label.FontWeight = 'bold';
    cb4.Label.String = 'Sen slope [% per year]';
    cb4.Label.FontSize = 15; %cb4.Label.FontWeight = 'bold';
    
    add_title(ax(1),'(a)',18);
    add_title(ax(2),'(b)',18);
    add_title(ax(3),'(c)',18);
    add_title(ax(4),'(d)',18);
elseif strcmp(ex,'default')
    data_inund   = load(['data/outputs_geo_CLM' forcing]);
    xtrain_inund = get_parameters('geo',numc,100);
    dist = (xtrain_inund(:,5) - 7.2).^2 + (xtrain_inund(:,6) - 0.26).^2;
    ind_ex2      = find(dist == min(dist));
    disp(['aw = ' num2str(xtrain_inund(ind_ex2,5)) ', ad = ' num2str(xtrain_inund(ind_ex2,6))]);
    
    qdef = nanmean(reshape(q(1,:)',[12,yr2-yr1+1]),2);
    qobs = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,yr2-yr1+1]),2);
    qinu = nanmean(reshape(data_inund.q(ind_ex2,1:end-12)',[12,yr2-yr1+1]),2);
    
    cmap = getPyPlot_cMap('Set1',9);
    
    figure;
    plot(1:12,qobs,'k--','LineWidth',3); hold on; grid on;
    plot(1:12,qdef,'-','Color',cmap(1,:),'LineWidth',3);
    plot(1:12,qinu,'-','Color',cmap(3,:),'LineWidth',3);
    xlim([1 12]);
    xticks([1:12]);
    xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
    set(gca,'FontSize',14);
    ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');
    leg = legend('Observation','Default no inundation','Default inundation');
    leg.FontSize = 14; leg.FontWeight = 'bold';
    
elseif strcmp(tag,'outputs_geo_CLMGSWP3v1') || strcmp(tag,'outputs_geo_CLMCRUNCEPv7')
    load('data/Pantanal_geometry.mat');
    rwid = rwid./7.2.*xtrain(ind2(1),5);
    rdep = rdep./0.26.*xtrain(ind2(1),6);
    figure; set(gcf,'Position',[10 10 1200 400]);
    subplot(1,2,1);
    patch(xv,yv,rwid,'LineStyle','none'); colorbar; axis equal;
    title('River width [m]','FontSize',16,'FontWeight','bold');
    subplot(1,2,2);
    patch(xv,yv,rdep,'LineStyle','none'); colorbar; axis equal;
    title('River depth [m]','FontSize',16,'FontWeight','bold');
    
end

% ind2 = find(KGE == max(KGE));

[et_in,lon_in2,lat_in2,tmon_in2] = get_GLEAM_ET(Sh,1980,yr2);
[et_in2,~,~,~]                   = get_FLUXCOM_ET(Sh,1980,yr2);

days_of_month = [31;28;31;30;31;30;31;31;30;31;30;31];
et_in = et_in.*86400;
et_in2= et_in2.*86400;
et    = E1toE2(et,2);
rain  = rain.*86400;

figure; set(gcf,'Position',[10 10 1200 800]);
ax1 = axes('Position',[0.075 0.535 0.525 0.415]); 
C = 3600*24*30/area*1000;
a2 = nanmean(reshape(q(ind2,13:end)',[12,yr2-yr1]),2).*C;
%a1 = nanmean(reshape(q(ind1,13:end)',[12,yr2-yr1]),2).*C;
b = nanmean(reshape(qrunoff(ind2,13:end)',[12,yr2-yr1]),2).*86400.*30;
c = nanmean(reshape(mu_gsim(yr >= yr1+1 & yr <= yr2),[12,yr2-yr1]),2).*C;
% plot((a1-nanmean(a1))./nanmean(a1),'b--','LineWidth',3);hold on; grid on;
% plot((a2-nanmean(a2))./nanmean(a2),'r--','LineWidth',3);hold on; grid on;
% plot((b-nanmean(b))./nanmean(b),'r-','LineWidth',3);hold on; grid on;
% plot((c-nanmean(c))./nanmean(c),'k-','LineWidth',3); 
% plot(a1,'b--','LineWidth',3);hold on; grid on;
plot(a2,'r--','LineWidth',3);hold on; grid on;
plot(b,'r-','LineWidth',3);hold on; grid on;
plot(c,'k-','LineWidth',3); 

if 0
if strcmp(tag,'outputs_sur_CLMGSWP3v1') || strcmp(tag,'outputs_sur_CLMCRUNCEPv7') 
%     plot(nanmean(qrunoff_cal,1).*86400.*30,'k--','LineWidth',3);
leg1 = legend('Simulated Q','Simulated R','Observed Q'); 
elseif strcmp(tag,'outputs_geo_CLMGSWP3v1') ||   ...
       strcmp(tag,'outputs_geo_CLMCRUNCEPv7') || ...
       strcmp(tag,'outputs_geo_noinund_CLMGSWP3v1')
    dist = sqrt((xtrain(:,5) - 7.2).^2 + (xtrain(:,6) - 0.26).^2);
    ind3 = find(dist == min(dist));
    a3 = nanmean(reshape(q(ind3,13:end)',[12,yr2-yr1]),2).*C;
    plot(a3,'k--','LineWidth',3);hold on; grid on;
    leg1 = legend('Calibrated Q','Simulated R','Observed Q','Default geometry'); 
else
    leg1 = legend('Simulated Q','Simulated R','Observed Q'); 
end
set(gca,'xtick',[],'FontSize',15); xlim([0.5 12.5]);
ylabel('[mm/month]','FontSize',16,'FontWeight','bold');

leg1.FontSize = 15; leg1.FontWeight = 'bold';

ax2 = axes('Position',[0.075 0.05 0.525 0.415]);
bar(nanmean(reshape(nanmean(rain(:,13:end),1)',[12 yr2-yr1]),2),'FaceColor',[0.5 0.5 0.5]);hold on; grid on;
%plot(nanmean(reshape(nanmean(et(ind1,13:end),1)',[12 yr2-yr1]),2),'b--','LineWidth',2);
plot(nanmean(reshape(nanmean(et_in,1)',[12 yr2-yr1]),2),'b-','LineWidth',3); 
plot(nanmean(reshape(nanmean(et_in2,1)',[12 yr2-yr1]),2),'g-','LineWidth',3); 
plot(nanmean(reshape(nanmean(et(ind2,13:end),1)',[12 yr2-yr1]),2),'r--','LineWidth',3);
set(gca,'xtick',[1:12], ...
        'xticklabel',{'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'},...
        'FontSize',15);xlim([0.5 12.5]);
ylabel('[mm/month]','FontSize',16,'FontWeight','bold');
leg2 = legend({[forcing ' P'],'GLEAM','FLUXCOM','Simulated E'}); 
leg2.FontSize = 15; leg2.FontWeight = 'bold';

ax3 = axes('Position',[0.675 0.05 0.30 0.9]);
et_annual = nansum(nanmean(reshape(et(ind2,:)',[12 yr2-yr1+1]),2).*days_of_month);
C = 3600*24/area*1000;
ro_annual = nansum(nanmean(reshape(q(ind2,:)',[12,yr2-yr1+1]),2).*C.*days_of_month);
et_annual2= nansum(nanmean(reshape(nanmean(et_in2,1)',[12 yr2-yr1]),2).*days_of_month);
ro_annual2= nansum(nanmean(reshape(mu_gsim(yr >= yr1+1 & yr <= yr2).*C,[12,yr2-yr1]),2).*days_of_month);
pr_annual = nansum(nanmean(reshape(nanmean(rain,1)',[12 yr2-yr1+1]),2).*days_of_month);
X = categorical({'Sim','Obs'});
Y = [et_annual ro_annual; et_annual2 ro_annual2]; 
X = reordercats(X,{'Sim','Obs'});
bar(categorical({forcing}),pr_annual); hold on;
bar(X,Y,'stacked');
set(gca,'FontSize',15);
ylabel('[mm/year]','FontSize',16,'FontWeight','bold');
leg3 = legend('P','E','Q');
leg3.FontSize = 15;
leg3.FontWeight = 'bold';
fs = 18;
t1 = add_title(ax1,'(a). Streamflow (Q)',fs);
t2 = add_title(ax2,'(b). Evapotranspiration (E) and Precipitation (P)',fs);
t3 = add_title(ax3,'(c). Annual water balance',fs);
end
% Compare annual runoff trend
if show_annual
    figure; set(gcf,'Position',[10 10 600 1000]);
    subplot(3,1,1);
    %plot(qyr(ind1,:),'b--','LineWidth',2); hold on;
    plot(yr1:yr2,qyr(ind2,:),'r-','LineWidth',2); hold on; grid on;
    %plot(qyr(ind4,:),'k--','LineWidth',2);
    plot(yr1:yr2,muyr,'k-','LineWidth',3);
    leg1 = legend('Calbrated fdrai','Obs'); 
    leg1.FontSize = 15; leg1.FontWeight = 'bold';
    xlim([yr1  yr2]);
    ylabel('Streamflow [m^{3}/s]','FontSize',15,'FontWeight','bold');

    subplot(3,1,2);
    et1 = nanmean(reshape(nanmean(et_in,1)',[12, yr2-yr1]),1); et1 = [NaN et1];
    et2 = nanmean(reshape(nanmean(et_in2,1)',[12, yr2-yr1]),1);et2 = [NaN et2];
    etyr= nanmean(reshape(et(ind2,:)',[12, yr2-yr1+1]),1);
    plot(yr1:yr2,et1,'b-','LineWidth',2); hold on;grid on;
    plot(yr1:yr2,et2,'g-','LineWidth',2); 
    plot(yr1:yr2,etyr,'r-','LineWidth',2); 
    xlim([yr1  yr2]);
    ylabel('ET [W/s]','FontSize',15,'FontWeight','bold');
    leg1 = legend('GLEAM','FLUXCOM','Calbrated fdrai'); 
    leg1.FontSize = 15; leg1.FontWeight = 'bold';

    subplot(3,1,3);
    a=nanmean(reshape(nanmean(rain,1)',[12, yr2-yr1+1]),1);
    plot(yr1:yr2,a,'k-','LineWidth',2); grid on;
    xlim([yr1 yr2]);
    ylabel('Precipitation [mm/month]','FontSize',15,'FontWeight','bold');
end
if strcmp(ex,'geo') && show_default_BFI
    figure;
    patch(xv,yv,nanmean(qdraiall(1:numc,:),2)./nanmean(qrunoffall(1:numc,:),2),'LineStyle','none'); colorbar;
    hold on; plot(Sh.X,Sh.Y,'k-','LineWidth',3);
    caxis([0 1]); colormap(jet(10));
end
if show_inund
    
    figure;set(gcf,'Position',[10 10 1200 500]);
    f = fdfall((ind2-1)*numc+1:ind2*numc,13:end);
    f1 = nanmean(f,1);
    f1 = nanmean(reshape(f1',[12,yr2-yr1]),2);
    f2 = nanmean(reshape(fh2osfc(ind2,13:end)',[12,yr2-yr1]),2);
    subplot(1,2,1);
    plot(f1,'r--','LineWidth',2); hold on; grid on;
    plot(f2,'r:','LineWidth',2); hold on; grid on;
    plot(f1+f2,'r-','LineWidth',2);
    subplot(1,2,2);
    patch(xv,yv,nanmean(f,2),'LineStyle','none'); axis equal; caxis([0 0.25]); colorbar;

    glad = load('/Users/xudo627/DATA/GLAD/GLAD_Pantanal_half.mat');

    glad_swf = NaN(length(xc),12);
    for i = 1 : length(xc)
        [row,col] = find(glad.lon == xc(i) & glad.lat == yc(i));
        glad_swf(i,:) = glad.swf(row,col,:)./100;
    end
    fdf = fdfall((ind2(1)-1)*numc+1:ind2(1)*numc,:);

    figure; set(gcf,'Position',[10 10 1200 400]);
    subplot(1,2,1);
    patch(xv,yv,nanmean(fdf,2),'LineStyle','none');  hold on; colorbar; axis equal;
    plot(Sh.X,Sh.Y,'r-','LineWidth',2);
    caxis([0 0.25]);
    title('Floodplain inudnation','FontSize',16,'FontWeight','bold');
    subplot(1,2,2);
    patch(xv,yv,nanmean(glad_swf(:,3),2),'LineStyle','none'); hold on; colorbar; axis equal;
    plot(Sh.X,Sh.Y,'r-','LineWidth',2);
    caxis([0 0.25]);
    title('GLAD','FontSize',16,'FontWeight','bold');
end

if show_domain
    figure;
    plot(lon_gsim,lat_gsim,'go','LineWidth',2); hold on;
    plot(lon_gsim2,lat_gsim2,'bo','LineWidth',2);
    plot(xc(ioutlet),yc(ioutlet),'gd','LineWidth',2); 
    plot(xc(ioutlet2),yc(ioutlet2),'bd','LineWidth',2); 
    plot(Sh.X,Sh.Y,'k-','LineWidth',2);
    plot(S2.X,S2.Y,'g-','LineWidth',2);
    plot(xc(icontribution2),yc(icontribution2),'g.','LineWidth',2); 
end
if show_sub_basin
    figure; set(gcf,'Position',[10 10 1400 600]);
    plot(tmon,mu_gsim2(yy >= yr1 & yy <= yr2),'k-','LineWidth',3); hold on;
    plot(tmon,q2(ind1,:),'b--','LineWidth',3);
    plot(tmon,q2(ind2,:),'r--','LineWidth',3);
    xlim([tmon(1) tmon(end)]);
    datetick('x','mmmyy','keeplimits');
end