% PAR saturated (> 800) ET, GEP and WUE vs VPD, for SWC quantiles (n quantiles)
figure; hold on;
for s = 1:2 % 2 seasons, winter or summer
    if s == 1
        %this_season = find(Precip_2daysago < 1 & Precip_1daysago < 0.5 & Precip_halfdayago < 0.2 & Month_t >= 5 & Month_t <= 8 & daytime == 1 & qc == 0 & qc_h20_flux == 0 & u > 0.2 & AGC_c == 0 & qc_Sc == 0 & NEE_c > -20 & NEE_c < 15 & PAR > 1000); % winter daytime
        this_season = find(Precip_2daysago < 1 & Precip_1daysago < 0.5 & Precip_halfdayago < 0.2 & Month_t >= 6 & Month_t <= 8 & daytime == 1 & qc == 0 & qc_h2o_flux == 0 & u > 0.2 & AGC_c == 0 & qc_Sc == 0 & NEE_c > -20 & NEE_c < 15 & PAR > 800); % winter daytime
    elseif s == 2
        %this_season = find(Precip_2daysago < 1 & Precip_1daysago < 0.5 & Precip_halfdayago < 0.2 & (Month_t >= 10 | Month_t <= 2) & daytime == 1 & qc == 0  & qc_h20_flux == 0 & u > 0.2 & AGC_c == 0 & qc_Sc == 0 & NEE_c > -20 & NEE_c < 15 & PAR > 1000); % summer daytime
        this_season = find(Precip_2daysago < 1 & Precip_1daysago < 0.5 & Precip_halfdayago < 0.2 & (Month_t >= 12 | Month_t <= 2) & daytime == 1 & qc == 0  & qc_h2o_flux == 0 & u > 0.2 & AGC_c == 0 & qc_Sc == 0 & NEE_c > -20 & NEE_c < 15 & PAR > 800); % summer daytime
    end
    GEP_s = -GPP(this_season); 
    ET_s = h2o_flux(this_season);
    WUE_s = GEP_s./(ET_s);
    PAR_s = PAR(this_season);
    SWCs = SWC_s(this_season);
    VPD_s = VPD(this_season);
    gc_s = gss(this_season);
    %VPD_binedges = quantile(VPD_s,0:1/n:1); %* n bins of driver class, by quantile (n bins with same amount of data) 
    SWC_binedges = quantile(SWCs,0:1/3:1); % 3 quantiles of SWC
    c = 0.8; % dark blue/red is wet soil
    for j = 1:3 % SWC
        if s == 1
            colors = [c c 1];
        elseif s == 2
            colors = [1 c c];
        end
        c = c - 0.4;
        this_bin = find(SWCs >= SWC_binedges(j) & SWCs <= SWC_binedges(j+1)); %* iteration change SWC
        subplot(2,2,1); 
        binplot_v3(VPD_s(this_bin),GEP_s(this_bin),4,colors,6);
        subplot(2,2,2); 
        binplot_v3(VPD_s(this_bin),ET_s(this_bin),4,colors,6);
        subplot(2,2,3);
        binplot_v3(VPD_s(this_bin),WUE_s(this_bin),4,colors,6);
        subplot(2,2,4);
        binplot_v3(VPD_s(this_bin),gc_s(this_bin),4,colors,6);
    end
end

subplot(2,2,1); hold on; ylabel('-GPP or A_m_a_x (\mumol m^-^2 s^-^1)'); box off;
ax = gca; ax.XTick = 0:1:5; ax.YTick = 4:4:24; %ax.FontSize = 12;
set(gca,'xticklabel',[]); ylim([4 24]); xlim([0 5]);
binplot_v3(VpdL_gim,Photo_gim,8,'k',6);

subplot(2,2,2); hold on; ylabel('ET or T (mmol m^-^2 s^-^1)'); box off;
ax = gca; ax.XTick = 0:1:5; %ax.YTick = 0.05:0.05:0.2; 
%ax.FontSize = 12;
set(gca,'xticklabel',[]); %ylim([0.05 0.23]);
xlim([0 5]);
binplot_v3(VpdL_gim,T_gim,8,'k',6);

subplot(2,2,3); hold on; ylabel('WUE (\mumol mmol^-^1)'); xlabel('D (kPa)'); box off;
ax = gca; ax.XTick = 0:1:5; % ax.YTick = 1000:2000:7000;
%ax.FontSize = 12;
%ylim([1000 7000]);
xlim([0 5]);
binplot_v3(VpdL_gim,Photo_gim./T_gim,8,'k',6);

subplot(2,2,4); hold on; ylabel('G_s or g_s (mmol m^-^2 s^-^1)'); xlabel('D (kPa)'); box off;
ax = gca; ax.XTick = 0:1:5; ax.YTick = 0:0.1:0.5; ylim([0 0.5]); 
%ax.FontSize = 12; 
xlim([0 5]);
binplot_v3(VpdL_gim,gs_gim,8,'k',6);

labelp = {'(a)','(b)','(c)','(d)'};
for i = 1:4
    subplot(2,2,i); hold on;
    text(0.90,0.98,labelp(i),'Units', 'Normalized', 'VerticalAlignment', 'Top');
    %sub_pos = get(gca,'position'); % get subplot axis position
    %set(gca,'position',sub_pos.*[1.2 1.2 1 1]) % stretch its width and height
end

% subplot(2,2,4); hold on; 
% h = legend('0-33%', '33-67%', '67-100%');
% set(h,'FontSize',6);legend('boxoff');