clear;close all;clc;

re = 6.37122e6;% Earth radius

addpath('/Users/xudo627/Developments/mylib/m/');
addpath('/Users/xudo627/Developments/Setup-E3SM-Mac/matlab-scripts-for-mosart/');
addpath('/Users/xudo627/Developments/inpoly');
addpath('/Users/xudo627/Developments/mylib/USGS-download/')
fname = '../inputdata_run/MOSART_Pantanal_all_template_c220126.nc';
area = ncread(fname,'area');
xc   = ncread(fname,'longxy');
yc   = ncread(fname,'latixy');

[lon_gsim,lat_gsim,Sh,yr,mo,da,mu_gsim0,sd,cv] = get_GSIM_discharge('BR_0002835',1);
total_area = max(areaint(Sh.Y,Sh.X)*4*pi.*(re^2)); 
[ioutlet0, icontributing0] = find_mosart_cell(fname,lon_gsim,lat_gsim,total_area);
mu_gsim0 = mu_gsim0(yr>= 1981 & yr <=2010);
mu_gsim0 = reshape(mu_gsim0(:),[12 length(mu_gsim0)/12]);

load('~/DATA/GLAD/GLAD_Pantanal_8th.mat');
axis equal; xlim([-65 -52]); ylim([-22 -14]); 

figure;
cmap = getPyPlot_cMap('BuPu');
imagesc([lon(1,1) lon(end,end)],[lat(1,1) lat(end,end)],nanmean(swf,3)); hold on;
set(gca,'YDir','normal'); colormap(cmap);
downstreamX = [-58 -55 -55 -58 -58];
downstreamY = [-20 -20 -18 -18 -20];
[gauges,in] = search_GSIM_gauge(Sh.X,Sh.Y,[],[],1);  hold on; grid on;
%plot(downstreamX,downstreamY,'r-','LineWidth',2);

% Read MOSART input files
if exist('Default_Streamflow.mat','file')
    load('Default_Streamflow.mat');
else
    files = dir('../data/default/Pantanal*mosart.h0*.nc');
    default = struct([]);
    default(1).qlnd = NaN(191,360);
    default(1).qocn = NaN(191,360);
    for i = 1 : length(files)
        filename  = fullfile(files(i).folder,files(i).name);
        tmp = ncread(filename,'RIVER_DISCHARGE_OVER_LAND_LIQ');
        default(1).qlnd(:,i) = tmp(1:191);
        tmp = ncread(filename,'RIVER_DISCHARGE_TO_OCEAN_LIQ');
        default(1).qocn(:,i) = tmp(1:191);
    end
    save('Default_Streamflow.mat','default');
end

if exist('Calibrated_Streamflow.mat','file')
    load('Calibrated_Streamflow.mat');
else
    files = dir('../data/calibration/Pantanal*mosart.h0*.nc');
    %files = files(133:end);
    qlnd = NaN(191,360);
    qocn = NaN(191,360);
    qsur = NaN(191,360);
    ical = 614;
    for i = 1 : length(files)
        filename  = fullfile(files(i).folder,files(i).name);
        disp(filename);
        tmp = ncread(filename,'RIVER_DISCHARGE_OVER_LAND_LIQ');
        qlnd(:,i) = tmp((ical - 1)*191+1 : ical*191);
        tmp = ncread(filename,'RIVER_DISCHARGE_TO_OCEAN_LIQ');
        qocn(:,i) = tmp((ical - 1)*191+1 : ical*191);
        tmp = ncread(filename,'QSUR_LIQ');
        qsur(:,i) = tmp((ical - 1)*191+1 : ical*191);
    end
    save('Calibrated_Streamflow.mat','qlnd','qocn','qsur');
end

k = 1;
for i = 1981 : 2010
    for j = 1 : 12
        tsim(k) = datenum(i,j,15,1,1,1);
        yrsim(k) = i;
        k = k + 1;
    end
end
k = 0;
for i = 1 : length(gauges)
    [lons(i),lats(i),S,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge(gauges{i},1);
    aream2(i) = max(areaint(S.Y,S.X)*4*pi.*(re^2)); 
    areakm2(i) = max(areaint(S.Y,S.X)*4*pi.*(re^2))/1e6; 
    
    tgsim = NaN(length(yr),1);
    for ii = 1 : length(yr)
        tgsim(ii) = datenum(yr(ii),mo(ii),15,1,1,1);
    end

    if areakm2(i) > 10*50*50
        k = k+1;
        subplot(5,5,k);
        [ioutlet, icontributing] = find_mosart_cell(fname,lons(i),lats(i),aream2(i));
        areamos = nansum(area([ioutlet; icontributing]));
        imagesc([lon(1,1) lon(end,end)],[lat(1,1) lat(end,end)],nanmean(swf,3)); hold on;
        set(gca,'YDir','normal'); colormap(cmap);
        plot(lons(i),lats(i),'kd','LineWidth',1); hold on;
        plot(xc([ioutlet; icontributing]),yc([ioutlet; icontributing]),'rx');
        plot(S.X,S.Y,'-','LineWidth',2);
        plot(Sh.X,Sh.Y,'k-','LineWidth',2);
        %disp([gauges{i} ', i = ' num2str(i) ', area = ' num2str(areakm2(i)) ' [km^2], mosart area = ' num2str(areamos/1e6)]);
        %disp(['Streamflow is avaialble from ' num2str(min(yr)) ' to ' num2str(max(yr))]);
        title(['i = ' num2str(i)]);
    end
end

f101 = figure(101); set(gcf,'Position',[10 10 800  400]);
f102 = figure(102); set(gcf,'Position',[10 10 800  400]);
f103 = figure(103); set(gcf,'Position',[10 10 800  400]);
f104 = figure(104); set(gcf,'Position',[10 10 800  400]);
[~,~,area8th] = xc2xv(lon,lat,1/8,1/8,1);
k = 0; 
gauge_names = {};
for i = [19 65 78 68 92 93]
    [lons(i),lats(i),S,yr,mo,da,mu_gsim,sd,cv] = get_GSIM_discharge(gauges{i},1);
    mu_gsim = mu_gsim(yr >= 1981 & yr <= 2010);
    aream2(i) = max(areaint(S.Y,S.X)*4*pi.*(re^2)); 
    areakm2(i) = max(areaint(S.Y,S.X)*4*pi.*(re^2))/1e6; 
    
    tgsim = NaN(length(yr),1);
    for ii = 1 : length(yr)
        tgsim(ii) = datenum(yr(ii),mo(ii),15,1,1,1);
    end
    tgsim = tgsim(yr >= 1981 & yr <= 2010);
    yr    = yr(yr >= 1981 & yr <= 2010);
    if areakm2(i) > 10 * 50 * 50
        k = k+1;
        figure(101);
        subplot_tight(2,3,k,[0.06 0.04]);
        [ioutlet, icontributing] = find_mosart_cell(fname,lons(i),lats(i),aream2(i));
        areamos = nansum(area([ioutlet; icontributing]));
        plot(Sh.X,Sh.Y,'k-','LineWidth',2);hold on;
        plot(S.X,S.Y,'-','Color',[0.6350 0.0780 0.1840],'LineWidth',2);
        plot(xc([ioutlet; icontributing]),yc([ioutlet; icontributing]), ...
             '.','Color',[0.3010 0.7450 0.9330],'LineWidth',2);
        plot(lons(i),lats(i),'kd','LineWidth',2); 
        disp([gauges{i} ', area = ' num2str(areakm2(i)) ' [km^2], mosart area = ' num2str(areamos/1e6)]);
        disp(['Streamflow is avaialble from ' num2str(min(yr)) ' to ' num2str(max(yr))]);
        add_title(gca,['Subbasin #' num2str(k)]);
        set(gca,'XTick',[],'YTick',[]);
        ind = find(yrsim' >= min(yr) & yrsim' <= max(yr));
%         subplot(3,6,k+6);
%         ind = find(yrsim' >= min(yr) & yrsim' <= max(yr));
%         plot(tgsim,mu_gsim,'k-','LineWidth',2); hold on; grid on;
%         plot(tsim(ind),qlnd(ioutlet,ind),'b--','LineWidth',2);
        [R2(k),~,NSE(k)] = estimate_evaluation_metric(mu_gsim(:),qlnd(ioutlet,ind)');
%         add_title(gca,['\rho = ' num2str(sqrt(R2(k))) ', NSE = ' num2str(NSE(k))]);
        mu_gsim = reshape(mu_gsim(:),[12 length(mu_gsim)/12]);
        
        figure(102);
        subplot(2,3,k);
        yyaxis left
        plot(nanmean(mu_gsim,2),'-','Color',[0, 0.4470, 0.7410],'LineWidth',2); hold on; grid on;
        if k == 1
            ylabel('Discharge [m^{3}/s]','FontSize',14,'FontWeight','bold');
        end
    
        yyaxis right;
        in = inpoly2([lon(:) lat(:)],[S.X' S.Y']);
        tmp = NaN(12,1);
        for jj = 1 : 12
            tmp2      = swf(:,:,jj)./100.*area8th./1e6;
            tmp(jj,1) = nansum(tmp2(in));
        end
        plot(tmp,'-','Color',[0.8500, 0.3250, 0.0980],'LineWidth',2); 
        xlim([1 12]);
        add_title(gca,['Subbasin #' num2str(k)]);

        if k == 3
            ylabel('Surface water [km^{2}]','FontSize',14,'FontWeight','bold');
        end
        if k == 1
            legend('Streamflow','Surface Water','FontSize',11,'color','none','EdgeColor','none');
        end
        set(gca,'XTick',[1:12],'XTickLabel',{'J','F','M','A','M','J','J','S','O','N','D'});
        set(gca,'XTickLabelRotation',0);


        figure(103);
        subplot(2,3,k);
        plot(nanmean(mu_gsim,2),'k-','LineWidth',2); hold on; grid on;
        tmp = nanmean(reshape(qlnd(ioutlet,ind)',[12 length(ind)/12]),2);
        plot(tmp,'b--','LineWidth',2);hold on; grid on;
        tmp = nanmean(reshape(default.qlnd(ioutlet,ind)',[12 length(ind)/12]),2);
        plot(tmp,'r--','LineWidth',2);
        add_title(gca,['Subbasin #' num2str(k)]);
        set(gca,'XTick',[1:12],'XTickLabel',{'J','F','M','A','M','J','J','S','O','N','D'});
        set(gca,'XTickLabelRotation',0);
        if k == 1
            legend('Observation','Calibration','Default','FontSize',11,'Color','none','EdgeColor','none');
        end

        figure(104);
        subplot(2,3,k);
        yyaxis left;
        tmp = nanmean(qsur([ioutlet; icontributing],ind),1);
        tmp = nanmean(reshape(tmp',[12 length(ind)/12]),2);
        plot(tmp,'-','Color',[0, 0.4470, 0.7410],'LineWidth',2);
        yyaxis right;
        in = inpoly2([lon(:) lat(:)],[S.X' S.Y']);
        tmp = NaN(12,1);
        for jj = 1 : 12
            tmp2      = swf(:,:,jj)./100.*area8th./1e6;
            tmp(jj,1) = nansum(tmp2(in));
        end
        plot(tmp,'-','Color',[0.8500, 0.3250, 0.0980],'LineWidth',2); 
        xlim([1 12]);
    end
end


han=axes(f103,'visible','off'); 
han.YLabel.Visible='on';
ylabel(han,'Discharge [m^{3}/s]','FontSize',14,'FontWeight','bold');
han.Position(1) = 0.1;