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

xv = ncread('domain.lnd.r05_oEC60to30v3.190418_GSWP3.nc','xv');
yv = ncread('domain.lnd.r05_oEC60to30v3.190418_GSWP3.nc','yv');
xc = ncread('domain.lnd.r05_oEC60to30v3.190418_GSWP3.nc','xc');
yc = ncread('domain.lnd.r05_oEC60to30v3.190418_GSWP3.nc','yc');
areaTotal = ncread('MOSART_2way_LLR3_c210112_GSWP3.nc','areaTotal');

major_basins_ind = [248;233;188;191;49;78;90;124;114;93;97;99;173;121;14];
load('/Users/xudo627/projects/Global_WM_inund/noinund/global_ind.mat');
load('/Users/xudo627/projects/Global_WM_inund/noinund/major_basins_all.mat');
gsim_nos = {global_ind(major_basins_ind).gsim_no};
gsim_nos{14} = global_ind(121).gsim_no;
area_gsim = [global_ind(major_basins_ind).area_gsim];

scheme1 = 'GSWP3_1way_LLR3_7e85a10';
scheme2 = 'GSWP3_2way_LLR3_7e85a10';
qover2way =load(['DATA/' scheme2 '/QOVER' scheme2 '.mat']);
qover1way =load(['DATA/' scheme1 '/QOVER' scheme1 '.mat']);
qdrai2way =load(['DATA/' scheme2 '/QDRAI' scheme2 '.mat']);
qdrai1way =load(['DATA/' scheme1 '/QDRAI' scheme1 '.mat']);
rdl2way  = load(['DATA/' scheme2 '/RDL' scheme2 '.mat']);
rdl1way  = load(['DATA/' scheme1 '/RDL' scheme1 '.mat']);

qo2mu = nanmean(qover2way.QOVER.QOVER(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,3);
qo1mu = nanmean(qover1way.QOVER.QOVER(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,3);
qo2std = nanstd(qover2way.QOVER.QOVER(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,[],3);
qo1std = nanstd(qover1way.QOVER.QOVER(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,[],3);

qd2mu = nanmean(qdrai2way.QDRAI.QDRAI(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,3);
qd1mu = nanmean(qdrai1way.QDRAI.QDRAI(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,3);
qd2std = nanstd(qdrai2way.QDRAI.QDRAI(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,[],3);
qd1std = nanstd(qdrai1way.QDRAI.QDRAI(:,:,(1993-1971)*12+1 : (2004+1-1971)*12).*86400.*30,[],3);

rdl2mu = nanmean(rdl2way.RDL.RIVER_DISCHARGE_OVER_LAND_LIQ(:,:,(1993-1971)*12+1 : (2004+1-1971)*12),3);
rdl1mu = nanmean(rdl1way.RDL.RIVER_DISCHARGE_OVER_LAND_LIQ(:,:,(1993-1971)*12+1 : (2004+1-1971)*12),3);
rdl2std = nanstd(rdl2way.RDL.RIVER_DISCHARGE_OVER_LAND_LIQ(:,:,(1993-1971)*12+1 : (2004+1-1971)*12),[],3);
rdl1std = nanstd(rdl1way.RDL.RIVER_DISCHARGE_OVER_LAND_LIQ(:,:,(1993-1971)*12+1 : (2004+1-1971)*12),[],3);

i = 12; % 11 Lena
[lon_gsim,lat_gsim,S,yr_gsim,mo_gsim,da_gsim,mu_gsim,sd_gsim,cv_gsim] = ...
                                         get_GSIM_discharge(gsim_nos{i},2);

fname = 'MOSART_2way_LLR3_c210112.nc';
[ioutlet, icontributing] = find_mosart_cell(fname,lon_gsim,lat_gsim,area_gsim(i));

[row,col] = ind2sub([720,360],[ioutlet; icontributing]);
i1 = find(row <= 360);
i2 = find(row > 360);
row(i1) = 360 + row(i1);
row(i2) = row(i2) - 360;

for i = 1 : length(row)
    iall(i) = sub2ind([720,360],row(i),col(i));
end
xv_basin = NaN(4,length(row));
yv_basin = NaN(4,length(row));
aT_basin = areaTotal(iall);
xc_basin = xc(iall);
yc_basin = yc(iall);
k = boundary(xc_basin',yc_basin',1);
for i = 1 : length(row)
    xv_basin(:,i) = xv(:,row(i),col(i));
    yv_basin(:,i) = yv(:,row(i),col(i));
end

qo2mu_basin  = qo2mu(iall);
qo1mu_basin  = qo1mu(iall);
qo2std_basin = qo2std(iall);
qo1std_basin = qo1std(iall);

qd2mu_basin  = qd2mu(iall);
qd1mu_basin  = qd1mu(iall);
qd2std_basin = qd2std(iall);
qd1std_basin = qd1std(iall);

rdl2mu_basin  = rdl2mu(iall);
rdl1mu_basin  = rdl1mu(iall);
rdl2std_basin = rdl2std(iall);
rdl1std_basin = rdl1std(iall);

figure; set(gcf,'Position',[10 10 1200 400]);
subplot(1,2,1);
patch(xv_basin,yv_basin,qo2mu_basin./qo1mu_basin,'LineStyle','none');hold on;
colormap(blue2red(128,'seismic')); caxis([0.8 1.2]);
plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);
title('Average surface runoff 2way/1way','FontSize',15,'FontWeight','bold');
colorbar;
% subplot(2,2,2);
% patch(xv_basin,yv_basin,qo2std_basin./qo1std_basin,'LineStyle','none');colormap(blue2red(128,'seismic')); hold on;
% caxis([0.8 1.2]);
% plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);
subplot(1,2,2);
patch(xv_basin,yv_basin,qd2mu_basin./qd1mu_basin,'LineStyle','none');hold on;
colormap(blue2red(128,'seismic')); caxis([0.8 1.2]);
plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);
title('Average sub-surface runoff 2way/1way','FontSize',15,'FontWeight','bold');
colorbar;
% subplot(2,2,4);
% patch(xv_basin,yv_basin,qd2std_basin./qd1std_basin,'LineStyle','none');colormap(blue2red(128,'seismic')); hold on;
% caxis([0.8 1.2]);
% plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);


figure;
subplot(2,2,1);
patch(xv_basin,yv_basin,rdl2mu_basin-rdl1mu_basin,'LineStyle','none');hold on;
colormap(blue2red(128,'seismic')); colorbar;caxis([-15 15]);
plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);
subplot(2,2,2);
patch(xv_basin,yv_basin,rdl2std_basin./rdl1std_basin,'LineStyle','none');hold on;
colormap(blue2red(128,'seismic')); colorbar;caxis([0.8 1.2]);
plot(xc_basin(k),yc_basin(k),'k-','LineWidth',2);