clear;close all;clc;

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

yr1 = 1981;
yr2 = 2010;
tmon = get_time_number(yr1,yr2,'mon');

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

mosart_input = './inputdata_run/MOSART_Pantanal_default_c211116.nc';
domain_input = './inputdata_run/domain_lnd_Pantanal_default.nc';

numc = 3194;

areaTotal = ncread(mosart_input,'areaTotal2');
rslp      = ncread(mosart_input,'rslp');
xv = ncread(domain_input,'xv');
yv = ncread(domain_input,'yv');
xc = ncread(domain_input,'xc');
yc = ncread(domain_input,'yc');

ioutlet = find(areaTotal == max(areaTotal));
cases = fieldnames(data);

colors = distinguishable_colors(5);

figure(1);
ind = find(yr >= yr1 & yr <= yr2);
mu_gsim = mu_gsim(ind);
plot(tmon,mu_gsim,'k-','LineWidth',2); hold on; grid on;
xlim([tmon(1) tmon(end)]);
datetick('x','mmmyy','keeplimits');
figure(2);
plot(1:12,nanmean(reshape(mu_gsim,[12,30]),2),'k-','LineWidth',2);hold on; grid on;

figure(3); set(gcf,'Position',[10 10 1400 500]);
for i = 1 : length(cases)
    if i < 4
        figure(1);
        q = data.(cases{i}).q2(ioutlet,(yr1-1972)*12+1:(yr2-1972+1)*12);
        plot(tmon,q,'-','Color',colors(i,:),'LineWidth',1.5); hold on; grid on;

        figure(2);
        plot(1:12,nanmean(reshape(q',[12,30]),2),'--','Color',colors(i,:),'LineWidth',1.5);
    elseif i == 4
        qall = data.(cases{i}).q2(ioutlet:numc:end,(yr1-1972)*12+1:(yr2-1972+1)*12);
        
        
        for j = 1 : 10
            q = qall(j,:);
            figure(1);
            plot(tmon,q,'--','LineWidth',1.5); hold on; grid on;
            figure(2);
            plot(1:12,nanmean(reshape(q',[12,30]),2),'--','LineWidth',1.5);
        end
        
    end
    figure(3);
    if i < 3
        subplot(1,2,i);
        f = data.(cases{i}).f2(:,(yr1-1972)*12+1:(yr2-1972+1)*12);
        patch(xv,yv,nanmean(f,2),'LineStyle','none'); colorbar; hold on; caxis([0 0.4]);       
    end
%     ind = find(rslp <= 0);
%     plot(xc(ind),yc(ind),'r.');
end

figure;
plot(Sh.X,Sh.Y,'k-','LineWidth',2); hold on;
plot(xc,yc,'b.');
plot(xc(ioutlet),yc(ioutlet),'ro','LineWidth',2);

figure;
ind = find(rslp <= 0);
rslp(ind) = 0.0001;
patch(xv,yv,rslp,'LineStyle','none'); colorbar; hold on;

plot(xc(ind),yc(ind),'r.');
caxis([0 0.001]); colormap(jet);

