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',2);
[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,2);

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);
load(['data/' tag '.mat']);
[xtrain,nr_uq] = get_parameters(ex,numc,ntot);

etann = NaN(ntot,32);
qann  = NaN(ntot,32);
qslope = NaN(ntot,1);
rainann = nanmean(reshape(nanmean(rain,1)',[12,32]),1).*86400.*365;
for i = 1 : ntot    
    etann(i,:) = E1toE2(nanmean(reshape(et(i,:)',[12 32]),1),2).*365;
    qann(i,:)= nanmean(reshape(qrunoff(i,:)',[12 32]),1).*86400.*365;
    qslope(i,1) =  sens_slope(qann(i,:)')./nanmean(qrunoff(i,:)');
    e2rlope(i,1) =  sens_slope(etann(i,:)'./rainann');
    a(i,1) = sens_slope(etann(i,:)');
end
tmp = sens_slope(rainann');
figure;
plot(xtrain(:,1), a,'bo'); grid on;
set(gca,'FontSize',13);
ylabel('ET trend [mm/yr]','FontSize',15,'FontWeight','bold');
xlabel('f_{drain}','FontSize',15,'FontWeight','bold');


figure;
plot(xtrain(:,1),nanmean(etann,2)./nanmean(rainann),'bo');

figure;
tmp = repmat(xtrain(:,1),1,32);
scatter(rainann(:),etann(:)./rainann(:),72,tmp(:),'filled'); colormap("jet"); colorbar;
figure;
plot(rainann,'k--','LineWidth',2); hold on;
plot(etann(5,:)','b--','LineWidth',2)
plot(etann(2,:)','r--','LineWidth',2)