% Need to load inputs... 
% Load_Inputs_FINAL;
months = {'J','F','M','A','M','J','J','A','S','O','N','D'};
i = 1;
for y = 1:3
    for m = 1:12
        month_dt(i) = datetime(2013+y,0+m,1);
        i = i+1;
    end
end
month_dt(37) = datetime(2017,1,1);
clearvars i m y;

figure('Units','pixels','Position',[1000 50 950 950]);
plot1 = axes('Units','pixels','Position',[150 600 700 340]);
hold on; box off;
for yeari = 2013:2017
    x1 = datenum(datetime(yeari,06,01));
    x2 = datenum(datetime(yeari,08,31));
    shadeLL([x1 x2],[250 250]',[-250 -250]',[0.7 0.7 0.7],1);
    x1 = datenum(datetime(yeari,12,01));
    x2 = datenum(datetime(yeari+1,02,30));
    shadeLL([x1 x2],[250 250]',[-250 -250]',[0.5 0.5 0.5],1);
end
ERp = bar(plot1,datenum(Monthly_date),ER_solo_Monthly,'r');
%bar(datenum(Monthly_date),ET_Monthly,'b');
GPPp = bar(plot1,datenum(Monthly_date),GPP_solo_Monthly,'g');
NEEp = bar(plot1,datenum(Monthly_date),NEE_solo_Monthly,'k');
xlim(plot1,[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
plot1.XTick = datenum(month_dt);
datetick(plot1,'x','m','keeplimits','keepticks');
plot(plot1,[datenum(2015,01,01) datenum(2015,01,01)],[-250 250],'k');
plot(plot1,[datenum(2016,01,01) datenum(2016,01,01)],[-250 250],'k');
text(datenum(2014,06,01),225,'2014','FontWeight','bold');
text(datenum(2015,06,01),225,'2015','FontWeight','bold');
text(datenum(2016,06,01),225,'2016','FontWeight','bold');
text(datenum(2014,01,15),210,'(a)');
for j = 1:3
    for i = 1:12;
    text(datenum(2013 + j,i,01)-7,ER_solo_Monthly(i+(j-1)*12)+10,months(i));
    end
end
ylabel(plot1,'C fluxes (g C m^-^2 month^-^1)');
legend([ERp GPPp NEEp],'ER','GPP','NEE');
set(plot1,'xticklabel',[]);
set(plot1,'ytick',-200:50:200);

plot2 = axes('Units','pixels','Position',[150 420 700 170]); hold on; box off;
for yeari = 2013:2017
    x1 = datenum(datetime(yeari,06,01));
    x2 = datenum(datetime(yeari,08,31));
    shadeLL([x1 x2],[800 800]',[0 0]',[0.7 0.7 0.7],1);
    x1 = datenum(datetime(yeari,12,01));
    x2 = datenum(datetime(yeari+1,02,30));
    shadeLL([x1 x2],[800 800]',[0 0]',[0.5 0.5 0.5],1);
end
[ax2,p1,p2] = plotyy(datenum(Monthly_date(1:36)),Precip_Monthly(1:36),DateTime_CUP,SWC_s,'bar','plot');
%p1.FaceColor = 'k';
p2.LineWidth = 1; p2.Color = 'b';
set(ax2(2),'YColor','b');
ylim(ax2(1),[0 800]); % left axis
ylim(ax2(2),[-20 45]); % right axis
xlim(ax2(1),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
xlim(ax2(2),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
set(ax2(1),'XTick',datenum(month_dt));
set(ax2(1),'YTick',50:100:400); set(ax2(2),'YTick',10:5:40);
set(ax2(1),'xticklabel',[]);
plot(ax2(1),[datenum(2015,01,01) datenum(2015,01,01)],[0 800],'k');
plot(ax2(1),[datenum(2016,01,01) datenum(2016,01,01)],[0 800],'k');
text(datenum(2014,01,10),700,'(b)');
text(datenum(2014,04,01),300,'SWC','Color','b');
for j = 1:3
    for i = 1:12;
    text(datenum(2013 + j,i,01)-7,Precip_Monthly(i+(j-1)*12)+30,months(i));
    end
end
ylabel(ax2(1),'Rainfall (mm month^-^1)'); % right y-axis
ylabel(ax2(2),'SWC_0_-_8_c_m (%)');
box off;

plot3 = axes('Units','pixels','Position',[150 240 700 170]);hold on; box off;
for yeari = 2013:2017
    x1 = datenum(datetime(yeari,06,01));
    x2 = datenum(datetime(yeari,08,31));
    shadeLL([x1 x2],[50 50]',[0 0]',[0.7 0.7 0.7],1);
    x1 = datenum(datetime(yeari,12,01));
    x2 = datenum(datetime(yeari+1,02,30));
    shadeLL([x1 x2],[50 50]',[0 0]',[0.5 0.5 0.5],1);
end
[ax2,p1,p2] = plotyy(datenum(Monthly_date(1:36)),Tammax_a,datenum(Monthly_date(1:36)),VPDmmax_a,'plot','plot');
p1.LineWidth = 2; p1.Color = 'k';
p2.LineWidth = 2; p2.Color = 'r';
set(ax2(1),'YColor','k');
set(ax2(2),'YColor','r');
ylim(ax2(1),[15 35]); % left axis
ylim(ax2(2),[0.5 3.5]); % right axis
xlim(ax2(1),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
xlim(ax2(2),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
set(ax2(1),'XTick',datenum(month_dt)); set(ax2(1),'xticklabel',[]);
set(ax2(1),'YTick',20:5:30); set(ax2(2),'YTick',1:0.5:3);
plot(ax2(1),[datenum(2015,01,01) datenum(2015,01,01)],[0 500],'k');
plot(ax2(1),[datenum(2016,01,01) datenum(2016,01,01)],[0 500],'k');
text(0.02,0.92,'(c)','Units', 'Normalized', 'VerticalAlignment', 'Top');
text(0.09,0.6,'Ta','Units', 'Normalized', 'VerticalAlignment', 'Top','Color','k');
text(0.06,0.35,'D','Units', 'Normalized', 'VerticalAlignment', 'Top','Color','r');
% for j = 1:3
%     for i = 1:12;
%     text(datenum(2013 + j,i,01)-7,Precip_Monthly(i+(j-1)*12)+20,months(i));
%     end
% end
ylabel(ax2(1),'Ta_m_a_x (°C)'); % right y-axis
ylabel(ax2(2),'D_m_a_x (kPa)');
box off;

plot4 = axes('Units','pixels','Position',[150 60 700 170]); hold on; box off;
for yeari = 2013:2017
    x1 = datenum(datetime(yeari,06,01));
    x2 = datenum(datetime(yeari,08,31));
    shadeLL([x1 x2],[2 2]',[0 0]',[0.7 0.7 0.7],1);
    x1 = datenum(datetime(yeari,12,01));
    x2 = datenum(datetime(yeari+1,02,30));
    shadeLL([x1 x2],[2 2]',[0 0]',[0.5 0.5 0.5],1);
end
text(0.02,0.92,'(d)','Units', 'Normalized', 'VerticalAlignment', 'Top');
text(0.09,0.87,'NDVI','Units', 'Normalized', 'VerticalAlignment', 'Top','Color','r');
text(datenum(2015,2,1),0.85,'LAI','Color','k');
plot(plot4,[datenum(2015,01,01) datenum(2015,01,01)],[0 2],'k');
plot(plot4,[datenum(2016,01,01) datenum(2016,01,01)],[0 2],'k');
LAIdatenum = datenum(DateTime_LAI);
%scatter(LAIdatenum,LAI,10,[0.6 0.6 0.6]); 
[xData, yData] = prepareCurveData( LAIdatenum, LAI );
% Set up fittype and options.
ft = fittype( 'smoothingspline' );
excludedPoints = excludedata( xData, yData, 'Indices', [8 129 131 132] );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 1e-03;
opts.Exclude = excludedPoints;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
beg_LAI = datenum(DateTime_LAI(1)); end_LAI = datenum(DateTime_LAI(length(DateTime_LAI)));
ax7 = gca; ax7.XLim = [beg_LAI end_LAI];
LAI_spline = plot(fitresult); set(LAI_spline,'Color','k'); set(LAI_spline,'LineWidth',2);
legend('hide');
legend off;
ylabel('LAI (m^-^2 m^-^2)'); %xlabel ('2014                       2015                      2016');
ax7 = gca; ax7.XTick = datenum(month_dt);
datetick('x','m','keeplimits','keepticks');

NDVIdatenum = datenum(DateTime_NDVI);
[xData, yData] = prepareCurveData( NDVIdatenum, NDVI );
% Set up fittype and options.
ft = fittype( 'smoothingspline' );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 1e-04;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
s_NDVI = feval(fitresult,NDVIdatenum);
[ax8,p1,p2] = plotyy(LAIdatenum,LAI,datenum(DateTime_NDVI),s_NDVI);
p1.LineWidth = 2; p1.Color = 'none';
p2.LineWidth = 2; p2.Color = 'r';
ylim(ax8(1),[0.65 1.1]); % left axis
ylim(ax8(2),[0.55 0.8]); % right axis
xlim(ax8(1),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
xlim(ax8(2),[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
ylabel(ax8(1),'LAI (m^2 m^-^2)') % left y-axis
ylabel(ax8(2),'NDVI') % right y-axis
set(ax8(2),'YColor', 'r');
set(ax8(1),'YColor','k');
set(ax8(1),'YTick',0.7:0.1:1); set(ax8(2),'YTick',0.6:0.05:0.75);
xlabel('Month');
box off;

% plot5 = axes('Units','pixels','Position',[50 60 2 170]); hold on;
% ylabel('Litterfall');
% delete(plot5);
 
plot6 = axes('Units','pixels','Position',[80 240 700 170]); hold on;
ylabel('PAR_m_a_x');
ylim(plot6,[800 2400]);
plot6.YTick = 1000:400:2200;
plot6.Color = 'none';
plot6.YColor = 'b';
plot6.XColor = 'none';

%delete(plot9);
plot9 = axes('Units','pixels','Position',[150 240 700 170]); hold on;
plot(plot9,datenum(Monthly_date(1:36)),PARmmax_a,'b','LineWidth',2);
xlim(plot9,[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
plot9.XTick = datenum(month_dt);
plot9.Color = 'none';
ylim(plot9,[800 2400]);
plot9.YTick = 1000:400:2200;
set(plot9,'yticklabel',[]);
set(plot9,'xticklabel',[]);
text(datenum(2014,3,15),1000,'PAR','Color','b');

%delete(plot7);
plot7 = axes('Units','pixels','Position',[80 60 700 170]); hold on;
ylabel('LP (m^2 m^-^2 month^-^1)');
ylim(plot7,[0 0.6]);
plot7.YTick = 0.1:0.1:0.5;
plot7.Color = 'none';
plot7.YColor = 'b';
plot7.XColor = 'none';

% delete(plot8);
plot8 = axes('Units','pixels','Position',[150 60 700 170]); hold on;
plot(plot8,datenum(Date_mean),Mist_F_mean + Euc_F_mean,'b','LineWidth',2);
xlim(plot8,[datenum(datetime(2013,12,15)) datenum(datetime(2016,12,15))]);
plot8.XTick = datenum(month_dt);
plot8.Color = 'none';
ylim(plot8,[0 0.6]);
plot8.YTick = 0.1:0.1:0.5;
set(plot8,'yticklabel',[]);
set(plot8,'xticklabel',[]);
plot8.YColor = 'none';
text(datenum(2015,8,12),0.15,'LP','Color','b');