clear;close all;clc;

% This script is to calibrate runoff generation processes

save_cal = 0;

addpath('/Users/xudo627/developments/getPanoply_cMap/');
addpath('/Users/xudo627/donghui/mylib/m/');
addpath('/Users/xudo627/donghui/CODE/Setup-E3SM-Mac/matlab-scripts-for-mosart/');
re = 6.37122e6;% Earth radius

cmap1 = getPanoply_cMap('NEO_ceres_net');
cmap2 = getPanoply_cMap('NEO_omi_ozone_to3');

tag = 'runall';
yr1 = 1979;
yr2 = 2010;

load('data/BFI3.mat');
if strcmp(tag,'runall')
    a = load(['data/outputs_run5_CLMCRUNCEPv7.mat'],'qrunoffall','qdraiall','fh2osfcall','fdfall','qrunoff','zwtall','etall');
    b = load(['data/outputs_run_CLMCRUNCEPv7.mat'],'qrunoffall','qdraiall','fh2osfcall','fdfall','qrunoff','zwtall','etall');
    qrunoffall = [a.qrunoffall; b.qrunoffall];
    qdraiall   = [a.qdraiall;   b.qdraiall  ];
    fh2osfcall = [a.fh2osfcall; b.fh2osfcall];
    fdfall     = [a.fdfall;     b.fdfall    ];
    zwtall     = [a.zwtall;     b.zwtall    ];
    etall      = [a.etall;      b.etall     ];
    qrunoff    = [a.qrunoff;    b.qrunoff   ];
    clear a b;
else
    load(['data/outputs_' tag '_CLMCRUNCEPv7.mat'],'qrunoffall','qdraiall','fh2osfcall','fdfall','qrunoff','zwtall','etall');
end
data_def = load('data/outputs_geo_CLMCRUNCEPv7.mat');
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');
area = ncread('inputdata_run/MOSART_Pantanal_all_template_c220126.nc','area');

% Process GLAD 
load('/Users/xudo627/DATA/GLAD/GLAD_Pantanal_half_annual.mat');

[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge('BR_0002835',1);

BFI_sims = nanmean(qdraiall./qrunoffall,2);
BFI_sims_def = nanmean(data_def.qdraiall./data_def.qrunoffall,2);
numc = length(BFI);
ntot = length(BFI_sims)/numc;
ntot_def = length(BFI_sims_def)/numc;

glad = NaN(numc,22);
for i = 1 : length(xc)
    [irow,icol] = find(lon == xc(i) & lat == yc(i));
    glad(i,:) = swf(irow,icol,:);
end
glad1= nanmean(glad(:,1:12),2)./100;

rng(1);
if strcmp(tag,'run')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1); % fdrain
    xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1); % fover
    xtrain(:,3)=ones(ntot,1).*0.1;
    xtrain(:,4)= 0 + (1-0)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26; 
elseif strcmp(tag,'run2')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1); % fdrain
    xtrain(:,2)= ones(ntot,1).*0.5; % fover
    xtrain(:,3)= ones(ntot,1).*0.1;
    xtrain(:,4)= 0 + (1-0)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 +(xtrain(:,2) - 0.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
elseif strcmp(tag,'run3')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,3)= 0.001 + (1-0.001)*rand(ntot,1); % micro_sigma
    xtrain(:,4)= 0.001 + (0.5-0.001)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 +(xtrain(:,2) - 0.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
elseif strcmp(tag,'run4')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,3)= 10.^(-3 + (0+3)*rand(1,ntot)); % micro_sigma
    xtrain(:,4)= 10.^(-3 + (-0.155+3)*rand(1,ntot)); %fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 +(xtrain(:,2) - 0.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
elseif strcmp(tag,'run5')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,3)= ones(ntot,1).*0.1; % micro_sigma
    xtrain(:,4)= 10.^(-3 + (-0.155+3)*rand(1,ntot)); %fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 +(xtrain(:,2) - 0.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
elseif strcmp(tag,'runall')
    rng(1); ntot = 1500;
    xtrain = NaN(ntot,6);
    xtrain(1:1000,1)= 0.1 + (5-0.1)*rand(1000,1);
    xtrain(1:1000,2)= 0.1 + (5-0.1)*rand(1000,1);
    xtrain(1:1000,3)= ones(1000,1).*0.1; % micro_sigma
    xtrain(1:1000,4)= 10.^(-3 + (-0.155+3)*rand(1,1000)); %fc   
    xtrain(1:1000,5)= ones(1000,1).*7.2;  
    xtrain(1:1000,6)= ones(1000,1).*0.26;  
    
    xtrain(1001:1500,1)= 0.1 + (5-0.1)*rand(500,1); % fdrain
    xtrain(1001:1500,2)= 0.1 + (5-0.1)*rand(500,1); % fover
    xtrain(1001:1500,3)=ones(500,1).*0.1;
    xtrain(1001:1500,4)= 0 + (1-0)*rand(500,1); % fc   
    xtrain(1001:1500,5)= ones(500,1).*7.2;  
    xtrain(1001:1500,6)= ones(500,1).*0.26; 
    
elseif strcmp(tag,'sub')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1); % fdrain
    xtrain(:,2)= ones(ntot,1).*0.5; % fover
    xtrain(:,3)= ones(ntot,1).*0.1;
    xtrain(:,4)= ones(ntot,1).*0.4; % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 +(xtrain(:,2) - 0.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
elseif strcmp(tag,'sur')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= ones(ntot,1).*2.5;
    xtrain(:,2)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,3)= 0.001 + (1-0.001)*rand(ntot,1); % micro_sigma
    xtrain(:,4)= 0.01 + (1-0.01)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;   
elseif strcmp(tag,'sur3')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= ones(ntot,1).*2.5;
    xtrain(:,2)= ones(ntot,1).*0.5;
    xtrain(:,3)=ones(ntot,1).*0.1;
    xtrain(:,4)= 0 + (0.5-0)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26; 
elseif strcmp(tag,'hydrology')
    xtrain = NaN(ntot,6);
    xtrain(:,1)= 0.1 + (5-0.1)*rand(ntot,1);
    xtrain(:,2)= ones(ntot,1).*0.5; % fover
    xtrain(:,3)= ones(ntot,1).*0.1; % micro_sigma
    xtrain(:,4)= 0.001 + (0.5-0.001)*rand(ntot,1); % fc   
    xtrain(:,5)= ones(ntot,1).*7.2;  
    xtrain(:,6)= ones(ntot,1).*0.26;  
    dist = (xtrain(:,1) - 2.5).^2 + (xtrain(:,4) - 0.4).^2;
    idef = find(dist == min(dist));
end



BFI_cal   = NaN(numc,1);
BFI_def   = NaN(numc,1);
qruno_cal = NaN(numc,size(qrunoffall,2));
qruno_def = NaN(numc,size(qrunoffall,2));
qdrai_cal = NaN(numc,size(qrunoffall,2));
qdrai_def = NaN(numc,size(qrunoffall,2));
zwt_cal   = NaN(numc,size(qrunoffall,2));
et_cal    = NaN(numc,size(qrunoffall,2));
et_def    = NaN(numc,size(qrunoffall,2));
inund_cal = NaN(numc,1);
inund_def = NaN(numc,1);
inund_frac_cal = NaN(numc,1);
fdrai   = NaN(numc,1);
fover   = NaN(numc,1);
fc      = NaN(numc,1);
sigma   = NaN(numc,1);

inund = fdfall(:,20*12+1:end) + fh2osfcall(:,20*12+1:end);
inund0 = data_def.fdfall(:,20*12+1:end) + data_def.fh2osfcall(:,20*12+1:end);
inund = nanmean(inund,2);
inund0 = nanmean(inund0,2);
inund_frac = nanmean(fh2osfcall(:,20*12+1:end),2)./inund;

w1 = 1;%1;
w2 = 5;%10;
w3 = 1;%1.5;
for i = 1 : numc
    disp(i);
    tmp1 = BFI_sims(i : numc : end);
    tmp2 = inund(i : numc :end);
    tmp3 = qrunoffall(i:numc:end,:);
    tmp4 = qdraiall(i:numc:end,:);
    tmp5 = inund_frac(i : numc :end);
    
    tmp6 = BFI_sims_def(i : numc : end);
    tmp7 = inund0(i : numc : end);
    tmp8 = data_def.qrunoffall(i : numc : end,:);
    tmp9 = data_def.qdraiall(i : numc : end,:);
    
    tmp10 = zwtall(i:numc:end,:);
    tmp11 = etall(i:numc:end,:);
    
    ind = find(xtrain(:,1) <= 0.2);
    tmp1(ind) = NaN;
    tmp2(ind) = NaN;
    tmp3(ind,:) = NaN;
    tmp4(ind,:) = NaN;
    tmp5(ind) = NaN;
    
    sens = NaN(ntot,1);
    for j = 1 : ntot
        tmp = nanmean(reshape(tmp3(j,:)',[12,32]),1)./1000 .* area(i);
        sens(j) = sens_slope(tmp(1:31)').*numc;
    end
    
    obj =  w1*abs(log(tmp2./glad1(i))) + ...
           w2*abs(log(tmp1./BFI(i)))   + ...
           w3*abs(log(sens./-59));%abs(log(tmp1./BFI(i)));%0.5 * abs(log(tmp2./glad1(i))) + 0.5 * abs(log(tmp1./BFI(i)));% + 0.5 * abs(log(tmp2./glad1(i))) +
    
%     obj =  abs(1 - tmp2./glad1(i)) + ...
%            abs(1 - tmp1./BFI(i))   + ...
%            abs(1 - sens./-59);
    obj =  abs((tmp2 - nanmean(tmp2))./nanstd(tmp2) - (glad1(i)- nanmean(tmp2))./nanstd(tmp2) ) + ...
           abs((tmp1 - nanmean(tmp1))./nanstd(tmp1) - (BFI(i)- nanmean(tmp1))./nanstd(tmp1))   + ...
           abs((sens - nanmean(sens))./nanstd(sens) - (-59- nanmean(sens))./nanstd(sens));
    ind = find(obj == min(obj)); ind = ind(1);
    
    BFI_cal(i) = tmp1(ind);
    BFI_def(i) = tmp6(54);
    inund_cal(i) = tmp2(ind);
    inund_def(i) = tmp7(54);
    inund_frac_cal(i) = tmp5(ind);
    qruno_cal(i,:) = tmp3(ind,:);
    qruno_def(i,:) = tmp8(54,:);
    qdrai_def(i,:) = tmp9(54,:);
    zwt_cal(i,:)   = tmp10(ind,:); 
    et_cal(i,:)    = tmp11(ind,:);
    
    fdrai(i)   = xtrain(ind,1);
    fover(i)   = xtrain(ind,2);
    sigma(i)   = xtrain(ind,3);
    fc(i)      = xtrain(ind,4); 
end

if save_cal
    save('data/Calibrated_hydrology.mat','fdrai','fover','sigma','fc');
end

figure; set(gcf,'Position',[10 10 750 800]);
w = 0.06;
d = 0.08;
ax1(1) = subplot_tight(3,2,1,[w d]);
patch(xv,yv,BFI,'LineStyle','none'); hold on; caxis([0 1]); axis equal; colormap(gca,cmap1);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'FontSize',13);
ax1(2) = subplot_tight(3,2,2,[w d]);
patch(xv,yv,BFI_def,'LineStyle','none'); hold on; caxis([0 1]); axis equal; colormap(gca,cmap1);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]); c1 = colorbar('east');
ax1(3) = subplot_tight(3,2,3,[w d]);
patch(xv,yv,glad1,'LineStyle','none'); hold on; caxis([0 0.2]); axis equal; colormap(gca,cmap2);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'FontSize',13);
ax1(4) = subplot_tight(3,2,4,[w d]);
patch(xv,yv,inund_def,'LineStyle','none'); hold on; caxis([0 0.2]);axis equal; colormap(gca,cmap2);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]); c2 = colorbar('east');

ax1(1).Position(1) = ax1(1).Position(1) + 0.02;
ax1(3).Position(1) = ax1(3).Position(1) + 0.02;
ax1(2).Position(1) = ax1(2).Position(1) - 0.02;
ax1(4).Position(1) = ax1(4).Position(1) - 0.02;
c1.Position = [ax1(2).Position(1) + ax1(2).Position(3) + 0.01 ax1(2).Position(2) 0.02 ax1(2).Position(4)];
c2.Position = [ax1(4).Position(1) + ax1(4).Position(3) + 0.01 ax1(4).Position(2) 0.02 ax1(4).Position(4)];
c1.FontSize = 14;
c2.FontSize = 14;

x0 = ax1(3).Position(1);
y0 = 0.05;
d0 = ax1(3).Position(2) - 0.05 - y0;
w0 = ax1(4).Position(1) + ax1(4).Position(3) - x0;
ax1(5) = axes('Position',[x0 y0 w0 d0]);
qruno_cal = qruno_cal ./1000 .* area; 
qruno_def = qruno_def ./1000 .* area;
yobs = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,32]),1)';
ydef = nanmean(reshape(nansum(qruno_def,1)',[12,32]),1)';
ycal = nanmean(reshape(nansum(qruno_cal,1)',[12,32]),1)';
yobs(end) = NaN;
ydef(end) = NaN;
ycal(end) = NaN;
[z(1), trend(1),~,pvalue(1)] = mk_test(yobs);
[z(2), trend(2),~,pvalue(2)] = mk_test(ydef);
[z(3), trend(3),~,pvalue(3)] = mk_test(ycal);
[sens(1), ~] = sens_slope(yobs); sens(1) = sens(1)/nanmean(yobs).*100;
[sens(2), ~] = sens_slope(ydef); sens(2) = sens(2)/nanmean(ydef).*100;
[sens(3), ~] = sens_slope(ycal); sens(3) = sens(3)/nanmean(ycal).*100;
plot(yr1:yr2,yobs,'k--','LineWidth',3); hold on; grid on;
plot(yr1:yr2,ydef,'r-','LineWidth',3);
xlim([yr1 yr2]);
set(gca,'FontSize',13);
ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');

KGE(1) = estimateKGE(BFI,BFI_def);
KGE(2) = estimateKGE(BFI,BFI_cal);
KGE(3) = estimateKGE(glad1,inund_def);
KGE(4) = estimateKGE(glad1,inund_cal);
[R2(1),RMSE(1),NSE(1),PBIAS(1)] = estimate_evaluation_metric(BFI,BFI_def);
disp(['Spatial correlation coefficient for default BFI is ' num2str(sqrt(R2(1)))]);
[R2(2),RMSE(2),NSE(2),PBIAS(2)] = estimate_evaluation_metric(BFI,BFI_cal);
disp(['Spatial correlation coefficient for calibrated BFI is ' num2str(sqrt(R2(2)))]);
[R2(3),RMSE(3),NSE(3),PBIAS(3)] = estimate_evaluation_metric(glad1,inund_def);
disp(['Spatial correlation coefficient for default Surface water is ' num2str(sqrt(R2(3)))]);
[R2(4),RMSE(4),NSE(4),PBIAS(4)] = estimate_evaluation_metric(glad1,inund_cal);
disp(['Spatial correlation coefficient for calibrated Surface water is ' num2str(sqrt(R2(4)))]);
sens   = round(sens,1);
PBIAS  = round(PBIAS);
pvalue = round(pvalue,2);

strs{1} = {['PBIAS = ' num2str(PBIAS(1)) ' [%]'], ...
           ['\rho = ' num2str(round(sqrt(R2(1)),2))]};
strs{2} = {['PBIAS = ' num2str(PBIAS(3)) ' [%]'], ...
           ['\rho = ' num2str(round(sqrt(R2(3)),2))]};
x0 = ax1(2).Position(1);
y0 = ax1(2).Position(2)+ax1(2).Position(4)-0.1;
w0 = 0.2;
d0 = 0.1;
t(1) = annotation('textbox',[x0 y0 w0 d0],'String',strs{1},       ... 
                      'FitBoxToText','on','EdgeColor','none','Color','k',     ...
                      'FontSize',14); 
x0 = ax1(4).Position(1);
y0 = ax1(4).Position(2)+ax1(4).Position(4)-0.1;
t(2) = annotation('textbox',[x0 y0 w0 d0],'String',strs{2},       ... 
                      'FitBoxToText','on','EdgeColor','none','Color','k',     ...
                      'FontSize',14); 
strs{3} = ['Observed:  Sen slope = ' num2str(sens(1)) ' [%/yr], p-value < 0.01'];
strs{4} = ['Default:      Sen slope = ' num2str(sens(2)) '.0 [%/yr], p-value = ' num2str(pvalue(2))];
t(3) = annotation('textbox',[0.425 0.225 0.3 0.1],'String',strs{3},       ... 
                'FitBoxToText','on','EdgeColor','none','Color','k',     ...
                'FontSize',13,'FontWeight','bold'); 
t(4) = annotation('textbox',[0.425 0.200 0.3 0.1],'String',strs{4},       ... 
                'FitBoxToText','on','EdgeColor','none','Color','r',     ...
                'FontSize',13,'FontWeight','bold');

            
add_title(ax1(1), '(a). GSCD BFI');
add_title(ax1(2), '(b). Default BFI');
add_title(ax1(3), '(c). GLAD SWF');
add_title(ax1(4), '(d). Default SWF');
add_title(ax1(5), '(e). Annual streamflow');

% Calibrated results
cols = getPyPlot_cMap('Set1',9);

figure; set(gcf,'Position',[10 10 750 800]);
w = 0.06;
d = 0.08;
ax2(1) = subplot_tight(3,2,3,[w d]);
patch(xv,yv,BFI_cal,'LineStyle','none'); hold on; caxis([0 1]); axis equal; colormap(gca,cmap1);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'FontSize',13); c3 = colorbar('west');
ax2(2) = subplot_tight(3,2,4,[w d]);
patch(xv,yv,inund_cal,'LineStyle','none'); hold on; caxis([0 0.2]);axis equal; colormap(gca,cmap2);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[],'FontSize',13); c4 = colorbar('west');
ax2(1).Position(1) = ax2(1).Position(1);
ax2(2).Position(1) = ax2(2).Position(1) + 0.01;
c3.Position = [ax2(1).Position(1) + ax2(1).Position(2) + 0.01 ax2(1).Position(2) 0.02 ax2(1).Position(4)];
c4.Position = [ax2(2).Position(1) + ax2(2).Position(2) + 0.01 ax2(2).Position(2) 0.02 ax2(2).Position(4)];
c3.FontSize = 14;
c4.FontSize = 14;

x0 = ax2(1).Position(1);
y0 = 0.05;
d0 = ax2(1).Position(2) - 0.05 - y0;
w0 = ax2(2).Position(1) + ax2(2).Position(3) - x0;
ax2(3) = axes('Position',[x0 y0 w0 d0]);
plot(yr1:yr2,yobs,'k--','LineWidth',3); hold on; grid on;
plot(yr1:yr2,ycal,'-','Color',cols(2,:),'LineWidth',3);
xlim([yr1 yr2]);
set(gca,'FontSize',13);
ylabel('Discharge [m^{3}/s]','FontSize',15,'FontWeight','bold');

strs{1} = {['PBIAS = ' num2str(PBIAS(2)) ' [%]'], ...
           ['\rho = ' num2str(round(sqrt(R2(2)),2))]};
strs{2} = {['PBIAS = ' num2str(PBIAS(4)) ' [%]'], ...
           ['\rho = ' num2str(round(sqrt(R2(4)),2))]};
x0 = ax2(1).Position(1);
y0 = ax2(1).Position(2)+ax2(1).Position(4)-0.1;
w0 = 0.2;
d0 = 0.1;
t(1) = annotation('textbox',[x0 y0 w0 d0],'String',strs{1},       ... 
                      'FitBoxToText','on','EdgeColor','none','Color','k', ...
                      'FontSize',14); 
x0 = ax2(2).Position(1);
y0 = ax2(2).Position(2)+ax2(2).Position(4)-0.1;
t(2) = annotation('textbox',[x0 y0 w0 d0],'String',strs{2},       ... 
                      'FitBoxToText','on','EdgeColor','none','Color','k', ...
                      'FontSize',14); 

ro = load('data/runoff_uncertainty.mat');
plot(yr1:yr2-1,ro.annual_grun,'-','Color',cols(3,:),'LineWidth',2); 
plot(yr1+1:yr2-1,ro.annual_lora,'-','Color',cols(4,:),'LineWidth',2); 
plot(yr1:yr2-1,ro.annual_grfr,'-','Color',cols(5,:),'LineWidth',2); 
models = {'GRUN:        ','LORA:         ','GRFR:         '};

strs{3} = ['Calibration: Trend = ' num2str(sens(3)) '0 [%/yr], p-value = ' num2str(pvalue(3))];
t(3) = annotation('textbox',[0.50 0.225 0.3 0.1],'String',strs{3},       ... 
                'FitBoxToText','on','EdgeColor','none','Color',cols(2,:), ...
                'FontSize',13,'FontWeight','bold');
for i = 1 : 3
    sens_slope = num2str(num2str(ro.slope(i+1)));
    pv         = num2str(ro.pvalue(i+1));
    if ro.pvalue(i+1) < 0.001
        pv = ' < 0.05';
    else
        pv = [' = ' pv(1:4)];
    end
    if ro.slope(i+1) < 0
        sens_slope = sens_slope(1:5);
    else
        sens_slope = [' ' sens_slope(1:4)];
    end
    strs{i+3} = [models{i} ' Trend = ' sens_slope ' [%/yr], p-value ' pv];
    t(i+3) = annotation('textbox',[0.50 0.225 - i*0.022 0.3 0.1],'String',strs{i+3},       ... 
                    'FitBoxToText','on','EdgeColor','none','Color',cols(i+2,:),     ...
                    'FontSize',13,'FontWeight','bold'); 
end

add_title(ax2(1), '(a). Calibrated BFI');
add_title(ax2(2), '(b). Calibrated SWF');
add_title(ax2(3), '(c). Annual streamflow');
% All together
figure; set(gcf,'Position',[10 10 950 800]);
w = 0.06;
d = 0.08;

ax(1) = subplot_tight(3,3,1,[w d]);
patch(xv,yv,BFI,'LineStyle','none'); hold on; caxis([0 1]); axis equal; colormap(gca,cmap1);
plot(Sh.X,Sh.Y,'k-','LineWidth',2);
title('(a). GSCD BFI3','FontSize',16,'FontWeight','bold');
ax(2) = subplot_tight(3,3,2,[w d]);
patch(xv,yv,BFI_def,'LineStyle','none'); hold on; caxis([0 1]); axis equal; colormap(gca,cmap1);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]);
title('(b). Default BFI','FontSize',16,'FontWeight','bold');
ax(3) = subplot_tight(3,3,3,[w d]);
patch(xv,yv,BFI_cal,'LineStyle','none'); hold on; caxis([0 1]); 
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]);
title('(c). Calibrated BFI','FontSize',16,'FontWeight','bold');
colormap(gca,cmap1); axis equal;c1 = colorbar('west');

ax(4) = subplot_tight(3,3,4,[w d]);
patch(xv,yv,glad1,'LineStyle','none'); hold on; caxis([0 0.2]); axis equal; colormap(gca,cmap2);
plot(Sh.X,Sh.Y,'k-','LineWidth',2);
title('(d). GLAD SWF','FontSize',16,'FontWeight','bold');
ax(5) = subplot_tight(3,3,5,[w d]);
patch(xv,yv,inund_def,'LineStyle','none'); hold on; caxis([0 0.2]);axis equal; colormap(gca,cmap2);
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]);
title('(e). Default SWF','FontSize',16,'FontWeight','bold');
ax(6) = subplot_tight(3,3,6,[w d]);
patch(xv,yv,inund_cal,'LineStyle','none'); hold on; caxis([0 0.2]); 
plot(Sh.X,Sh.Y,'k-','LineWidth',2); set(gca,'xtick',[],'ytick',[]);
title('(f). Calibrated SWF','FontSize',16,'FontWeight','bold');
colormap(gca,cmap2); axis equal;c2 = colorbar('west');

c1.Position = [ax(3).Position(1) + ax(3).Position(3) + 0.01 ax(3).Position(2) 0.02 ax(3).Position(4)];
c2.Position = [ax(6).Position(1) + ax(6).Position(3) + 0.01 ax(6).Position(2) 0.02 ax(6).Position(4)];
c1.FontSize = 14;
c2.FontSize = 14;

x0 = ax(4).Position(1);
y0 = 0.05;
d0 = ax(4).Position(2) - 0.05 - y0;
w0 = ax(6).Position(1) + ax(6).Position(3) - x0;
axes('Position',[x0 y0 w0 d0]);

yobs_mo = nanmean(reshape(mu_gsim(yr >= yr1 & yr <= yr2),[12,32]),2)';
ydef_mo = nanmean(reshape(nansum(qruno_def,1)',[12,32]),2)';
ycal_mo = nanmean(reshape(nansum(qruno_cal,1)',[12,32]),2)';

plot(yr1:yr2,yobs,'k-','LineWidth',3); hold on; grid on;
plot(yr1:yr2,ydef,'r-','LineWidth',3);
plot(yr1:yr2,ycal,'b-','LineWidth',3);
xlim([yr1 yr2]);
title('(g). Annual streamflow trend','FontSize',16,'FontWeight','bold');

for i = 1 : 4
    strs{i} = {['PBIAS = ' num2str(PBIAS(i)) ' [%]'], ...
               ['\rho = ' num2str(round(sqrt(R2(i)),2))]};
    if i <= 2
        x0 = ax(i+1).Position(1);
        y0 = ax(i+1).Position(2)+ax(i+1).Position(4)-0.1;
        w0 = 0.2;
        d0 = 0.1;
    else
        x0 = ax(i+2).Position(1);
        y0 = ax(i+2).Position(2)+ax(i+2).Position(4)-0.1;
        w0 = 0.2;
        d0 = 0.1;
    end
    t(i) = annotation('textbox',[x0 y0 w0 d0],'String',strs{i},       ... 
                      'FitBoxToText','on','EdgeColor','none','Color','k',     ...
                      'FontSize',14); 
end
strs{5} = ['Observed:  Sen slope = ' num2str(sens(1)) ' [m^{3}/s/yr], p-value = ' num2str(pvalue(1))];
strs{6} = ['Default:      Sen slope = ' num2str(sens(2)) '.0 [m^{3}/s/yr], p-value = ' num2str(pvalue(2))];
strs{7} = ['Calibrated: Sen slope = ' num2str(sens(3)) ' [m^{3}/s/yr], p-value = ' num2str(pvalue(3))];
t(5) = annotation('textbox',[0.52 0.225 0.3 0.1],'String',strs{5},       ... 
                'FitBoxToText','on','EdgeColor','none','Color','k',     ...
                'FontSize',14); 
t(6) = annotation('textbox',[0.52 0.200 0.3 0.1],'String',strs{6},       ... 
                'FitBoxToText','on','EdgeColor','none','Color','r',     ...
                'FontSize',14);
t(7) = annotation('textbox',[0.52 0.175 0.3 0.1],'String',strs{7},       ... 
                'FitBoxToText','on','EdgeColor','none','Color','b',     ...
                'FontSize',14);

figure; set(gcf,'Position',[10 10 1200 500]);
subplot(2,2,1);
patch(xv,yv,fdrai,'LineStyle','none'); hold on; caxis([0.1 5]); colormap(gca,jet(10)); axis equal; colorbar;
title('f_{drain}','FontSize',15,'FontWeight','bold');
subplot(2,2,2);
patch(xv,yv,fover,'LineStyle','none'); hold on; caxis([0.1 5]); colormap(gca,jet(10)); axis equal; colorbar;
title('f_{over}','FontSize',15,'FontWeight','bold');
subplot(2,2,3);
patch(xv,yv,sigma,'LineStyle','none'); hold on;  colormap(gca,jet(10)); axis equal; colorbar;
title('microtopography','FontSize',15,'FontWeight','bold');
subplot(2,2,4);
patch(xv,yv,fc,'LineStyle','none'); hold on; caxis([0 0.5]); colormap(gca,jet(10)); axis equal; colorbar;
title('f_{c}','FontSize',15,'FontWeight','bold');

figure;
plot(yobs_mo,'k-','LineWidth',2); hold on;
plot(ydef_mo,'r-','LineWidth',2); hold on;
plot(ycal_mo,'b-','LineWidth',2); hold on;
plot(reshape(nansum(qruno_cal,1)',[12,32]),'-');

load('zwt_fan.mat');
tmon = get_time_number(yr1,yr2-1);
figure;set(gcf,'Position',[10 10 1200 400]);
plot(tmon,data_def.zwt(1,1:end-12),'r-','LineWidth',2); hold on; grid on;
plot(tmon,nanmean(zwt_cal(:,1:end-12)),'b-','LineWidth',2); 
plot(tmon,ones(length(tmon),1).*nanmean(zwt_fan),'k--','LineWidth',2);
xlim([tmon(1) tmon(end)]); datetick('x','mmmyy','keeplimits');
set(gca, 'YDir','reverse','FontSize',13);
add_title(gca,'Water table depth [m]');
leg1 = legend('Default simulation','Calibrated simulation','Fan et al. 2013');
leg1.FontSize=14;
leg1.FontWeight = 'bold';

[tmp1,lon_gle,lat_gle,~] = get_GLEAM_ET(Sh,1980,2009);
[tmp2,lon_com,lat_com,~] = get_FLUXCOM_ET(Sh,1980,2009);
tmp1 = E1toE2(tmp1.*86400,1);
tmp2 = E1toE2(tmp2.*86400,1);
et_com = NaN(numc,360);
et_gle = NaN(numc,360);
for i = 1 : numc
    ind1 = find(lon_gle == xc(i) & lat_gle == yc(i));
    ind2 = find(lon_com == xc(i) & lat_com == yc(i));
    if ~isnan(ind1)
        et_gle(i,:) = tmp1(ind1,:);
    else
        disp(i);
    end
    if ~isnan(ind2)
        et_com(i,:) = tmp2(ind2,:);
    else
        disp(i);
    end
end

et_def = data_def.etall(1:numc,:);
figure;
subplot(3,2,1);
patch(xv,yv,nanmean(et_def(:,13:end-12),2),'LineStyle','none'); colorbar;
caxis([30 110]);
subplot(3,2,2);
patch(xv,yv,nanmean(et_cal(:,13:end-12),2),'LineStyle','none'); colorbar;
caxis([30 110]);

subplot(3,2,3);
patch(xv,yv,nanmean(et_com,2),'LineStyle','none'); colorbar;
caxis([30 110]);

subplot(3,2,4);
patch(xv,yv,nanmean(et_gle,2),'LineStyle','none'); colorbar;
caxis([30 110]);

subplot(3,2,5:6);
plot(nanmean(reshape(nanmean(et_def(:,13:end-12),1)',[12,30]),2),'r-','LineWidth',2); hold on; grid on;
plot(nanmean(reshape(nanmean(et_cal(:,13:end-12),1)',[12,30]),2),'b-','LineWidth',2);
plot(nanmean(reshape(nanmean(et_com,1)',[12,30]),2),'k--','LineWidth',2);
plot(nanmean(reshape(nanmean(et_gle,1)',[12,30]),2),'g--','LineWidth',2);

figure;
subplot(2,3,1);
subplot(2,3,2);
subplot(2,3,4:6);
plot(nanmean(et_def(:,13:end-12),1),'r-','LineWidth',2); hold on; grid on;
plot(nanmean(et_cal(:,13:end-12),1),'b-','LineWidth',2);
plot(nanmean(et_com,1),'k--','LineWidth',2);

figure;
plot(nanmean(reshape(nanmean(et_com,1)',[12,30]),2),'k--','LineWidth',3);hold on; grid on;
plot(nanmean(reshape(nanmean(et_gle,1)',[12,30]),2),'g--','LineWidth',3);
plot(nanmean(reshape(nanmean(et_def(:,13:end-12),1)',[12,30]),2),'r-','LineWidth',2); 
plot(nanmean(reshape(nanmean(et_cal(:,13:end-12),1)',[12,30]),2),'b-','LineWidth',2);

xlim([1 12]);
ax = gca; set(gca,'FontSize',13);
ax.XTick = 1 : 12;
ax.XTickLabel = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};
leg = legend('FLUXCOM','GLEAM','Default','Calibration');
leg.FontSize = 14;
leg.FontWeight = 'bold';

ylabel('\lambda E [W/m^{2}]','FontSize',15,'FontWeight','bold');