% Fig 9
clear
close all

%% load
load('CL_TCCIP_c0_end_TtripleRH_transform_raw.mat');
Cevap_e = [nan(7,1);Cevap(1:end-7)];
Ctrans_e = [nan(7,1);Ctrans(1:end-7)];
CH2O_e = [nan(7,1);CH2O(1:end-7)];
Q2M_e = [nan(7,1);Q2M(1:end-7)*1000];
RH2M_e = [nan(7,1);RH2M(1:end-7)];
T2M_e = [nan(7,1);T2M(1:end-7)];
TV_e = [nan(7,1);TV(1:end-7)];
PBOT_e = [nan(7,1);PBOT(1:end-7)];
photo_e = [nan(7,1);photo(1:end-7)];
RSSUN_e = [nan(7,1);RSSUN(1:end-7)];
RSSHA_e = [nan(7,1);RSSHA(1:end-7)];
clear Swnet Lwnet Swdw Swup U10 WIND QDRIP QINTR RSSUN RSSHA
clear Cevap Ctrans Gevap CH2O LE QBOT Q2M RH2M Prec SH TBOT T2M PBOT photo TV TG

load('CL_TCCIP_ctl_transform_raw.mat');
Cevap_c = [nan(7,1);Cevap(1:end-7)];
Ctrans_c = [nan(7,1);Ctrans(1:end-7)];
CH2O_c = [nan(7,1);CH2O(1:end-7)];
Q2M_c = [nan(7,1);Q2M(1:end-7)*1000];
RH2M_c = [nan(7,1);RH2M(1:end-7)];
T2M_c = [nan(7,1);T2M(1:end-7)];
TV_c = [nan(7,1);TV(1:end-7)];
PBOT_c = [nan(7,1);PBOT(1:end-7)];
photo_c = [nan(7,1);photo(1:end-7)];
RSSUN_c = [nan(7,1);RSSUN(1:end-7)];
RSSHA_c = [nan(7,1);RSSHA(1:end-7)];
clear Swnet Lwnet Swdw Swup U10 WIND QDRIP QINTR RSSUN RSSHA
clear Cevap Ctrans Gevap CH2O LE QBOT Q2M RH2M Prec SH TBOT T2M PBOT photo TV TG

day_data = 24;
% Q2M, RH2M--Qs2M [annual diurnal cycle with std]
Lv = 2.5*10^6; Rv = 461; T0 = 273.15;
Qs2M_c = Q2M_c./RH2M_c*100;
Qs2M_e = Q2M_e./RH2M_e*100;
QsTV_c = 6.11*100*exp(-(Lv./Rv)*((1./TV_c)-1/T0))*0.622./PBOT_c; QsTV_c = QsTV_c*1000; % kg/kg-->g/kg
QsTV_e = 6.11*100*exp(-(Lv./Rv)*((1./TV_e)-1/T0))*0.622./PBOT_e; QsTV_e = QsTV_e*1000; % kg/kg-->g/kg
dQ2M = Q2M_e-Q2M_c; dQs2M = Qs2M_e-Qs2M_c; dQsTV =  QsTV_e-QsTV_c; 
ano_Q2M = nan(24,1); ano_Qs2M = nan(24,1); ano_QsTV = nan(24,1);
for i = 1:day_data
    ano_Q2M(i,1) = mean(dQ2M(i:day_data:end),'omitnan');
    ano_Qs2M(i,1) = mean(dQs2M(i:day_data:end),'omitnan');
    ano_QsTV(i,1) = mean(dQsTV(i:day_data:end),'omitnan');
end
dCH2O = CH2O_e-CH2O_c; dCevap = Cevap_e-Cevap_c; dCtrans = Ctrans_e-Ctrans_c;
ano_CH2O = nan(1,24); ano_Cevap = nan(1,24); ano_Ctrans = nan(1,24);
ano_RSSUN = nan(1,24); ano_RSSHA = nan(1,24);
for i = 1:day_data
    ano_CH2O(i) = mean(dCH2O(i:day_data:end),'omitnan');
    ano_Cevap(i) = mean(dCevap(i:day_data:end),'omitnan');
    ano_Ctrans(i) = mean(dCtrans(i:day_data:end),'omitnan');
end

TotalTime = 0;
Time = 0;
photo_c_mon = nan(31*24*3,12);
photo_e_mon = nan(31*24*3,12);
Ctrans_c_mon = nan(31*24*3,12);
Ctrans_e_mon = nan(31*24*3,12);
TV_c_mon = nan(31*24*3,12);
TV_e_mon = nan(31*24*3,12);
T2M_c_mon = nan(31*24*3,12);
T2M_e_mon = nan(31*24*3,12);
Days = [31 28 31 30 31 30 31 31 30 31 30 31];
for yr = 1:3
    for mn = 1:12
        photo_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = photo_c(Time+1:Time+Days(mn)*24,1);
        photo_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = photo_e(Time+1:Time+Days(mn)*24,1);
        Ctrans_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = Ctrans_c(Time+1:Time+Days(mn)*24,1);
        Ctrans_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = Ctrans_e(Time+1:Time+Days(mn)*24,1);
        TV_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = TV_c(Time+1:Time+Days(mn)*24,1);
        TV_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = TV_e(Time+1:Time+Days(mn)*24,1);
        T2M_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = T2M_c(Time+1:Time+Days(mn)*24,1);
        T2M_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = T2M_e(Time+1:Time+Days(mn)*24,1);
        Time = Time+Days(mn)*24;
    end
    TotalTime = TotalTime+31*24;
end
% turn monthly data to seasonal; calculate anomaly
seas_array = [3 4 5 1 2 12; 6 7 8 9 10 11]; % 1=spring&winter; 2=summer&fall
photo_e_seas = nan(31*24*3*6,2);
photo_c_seas = nan(31*24*3*6,2);
Ctrans_e_seas = nan(31*24*3*6,2);
Ctrans_c_seas = nan(31*24*3*6,2);
TV_e_seas = nan(31*24*3*6,2);
TV_c_seas = nan(31*24*3*6,2);
T2M_e_seas = nan(31*24*3*6,2);
T2M_c_seas = nan(31*24*3*6,2);
for seas = 1:2
    photo_e_seas(:,seas) = [photo_e_mon(:,seas_array(seas,1));photo_e_mon(:,seas_array(seas,2));photo_e_mon(:,seas_array(seas,3));...
                     photo_e_mon(:,seas_array(seas,4));photo_e_mon(:,seas_array(seas,5));photo_e_mon(:,seas_array(seas,6))];
    photo_c_seas(:,seas) = [photo_c_mon(:,seas_array(seas,1));photo_c_mon(:,seas_array(seas,2));photo_c_mon(:,seas_array(seas,3));...
                     photo_c_mon(:,seas_array(seas,4));photo_c_mon(:,seas_array(seas,5));photo_c_mon(:,seas_array(seas,6))];
    Ctrans_e_seas(:,seas) = [Ctrans_e_mon(:,seas_array(seas,1));Ctrans_e_mon(:,seas_array(seas,2));Ctrans_e_mon(:,seas_array(seas,3));...
                     Ctrans_e_mon(:,seas_array(seas,4));Ctrans_e_mon(:,seas_array(seas,5));Ctrans_e_mon(:,seas_array(seas,6))];
    Ctrans_c_seas(:,seas) = [Ctrans_c_mon(:,seas_array(seas,1));Ctrans_c_mon(:,seas_array(seas,2));Ctrans_c_mon(:,seas_array(seas,3));...
                     Ctrans_c_mon(:,seas_array(seas,4));Ctrans_c_mon(:,seas_array(seas,5));Ctrans_c_mon(:,seas_array(seas,6))];
    TV_e_seas(:,seas) = [TV_e_mon(:,seas_array(seas,1));TV_e_mon(:,seas_array(seas,2));TV_e_mon(:,seas_array(seas,3));...
                     TV_e_mon(:,seas_array(seas,4));TV_e_mon(:,seas_array(seas,5));TV_e_mon(:,seas_array(seas,6))];
    TV_c_seas(:,seas) = [TV_c_mon(:,seas_array(seas,1));TV_c_mon(:,seas_array(seas,2));TV_c_mon(:,seas_array(seas,3));...
                     TV_c_mon(:,seas_array(seas,4));TV_c_mon(:,seas_array(seas,5));TV_c_mon(:,seas_array(seas,6))];
    T2M_e_seas(:,seas) = [T2M_e_mon(:,seas_array(seas,1));T2M_e_mon(:,seas_array(seas,2));T2M_e_mon(:,seas_array(seas,3));...
                     T2M_e_mon(:,seas_array(seas,4));T2M_e_mon(:,seas_array(seas,5));T2M_e_mon(:,seas_array(seas,6))];
    T2M_c_seas(:,seas) = [T2M_c_mon(:,seas_array(seas,1));T2M_c_mon(:,seas_array(seas,2));T2M_c_mon(:,seas_array(seas,3));...
                     T2M_c_mon(:,seas_array(seas,4));T2M_c_mon(:,seas_array(seas,5));T2M_c_mon(:,seas_array(seas,6))];
end

for i = 1:24
    ano_photo1(i) = mean((photo_e_seas(i:24:end,1)-photo_c_seas(i:24:end,1)),'omitnan');
    ano_photo2(i) = mean((photo_e_seas(i:24:end,2)-photo_c_seas(i:24:end,2)),'omitnan');
    ano_Ctrans1(i) = mean((Ctrans_e_seas(i:24:end,1)-Ctrans_c_seas(i:24:end,1)),'omitnan');
    ano_Ctrans2(i) = mean((Ctrans_e_seas(i:24:end,2)-Ctrans_c_seas(i:24:end,2)),'omitnan');
    ano_TV1(i) = mean((TV_e_seas(i:24:end,1)-TV_c_seas(i:24:end,1)),'omitnan');
    ano_TV2(i) = mean((TV_e_seas(i:24:end,2)-TV_c_seas(i:24:end,2)),'omitnan');
    ano_T2M1(i) = mean((T2M_e_seas(i:24:end,1)-T2M_c_seas(i:24:end,1)),'omitnan');
    ano_T2M2(i) = mean((T2M_e_seas(i:24:end,2)-T2M_c_seas(i:24:end,2)),'omitnan');
end
TotalTime = 0;
Time = 0;
RSSUN_c_mon = nan(31*24*3,12);
RSSUN_e_mon = nan(31*24*3,12);
RSSHA_c_mon = nan(31*24*3,12);
RSSHA_e_mon = nan(31*24*3,12);
Days = [31 28 31 30 31 30 31 31 30 31 30 31];
for yr = 1:3
    for mn = 1:12
        RSSUN_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = RSSUN_c(Time+1:Time+Days(mn)*24,1);
        RSSUN_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = RSSUN_e(Time+1:Time+Days(mn)*24,1);
        RSSHA_c_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = RSSHA_c(Time+1:Time+Days(mn)*24,1);
        RSSHA_e_mon(TotalTime+1:TotalTime+Days(mn)*24,mn) = RSSHA_e(Time+1:Time+Days(mn)*24,1);
        Time = Time+Days(mn)*24;
    end
    TotalTime = TotalTime+31*24;
end

% turn monthly data to seasonal; calculate anomaly

seas_array = [3 4 5 1 2 12; 6 7 8 9 10 11]; % 1=spring&winter; 2=summer&fall
RSSUN_e_seas = nan(31*24*3*6,2);
RSSUN_c_seas = nan(31*24*3*6,2);
RSSHA_e_seas = nan(31*24*3*6,2);
RSSHA_c_seas = nan(31*24*3*6,2);
for seas = 1:2
    RSSUN_e_seas(:,seas) = [RSSUN_e_mon(:,seas_array(seas,1));RSSUN_e_mon(:,seas_array(seas,2));RSSUN_e_mon(:,seas_array(seas,3));...
                     RSSUN_e_mon(:,seas_array(seas,4));RSSUN_e_mon(:,seas_array(seas,5));RSSUN_e_mon(:,seas_array(seas,6))];
    RSSUN_c_seas(:,seas) = [RSSUN_c_mon(:,seas_array(seas,1));RSSUN_c_mon(:,seas_array(seas,2));RSSUN_c_mon(:,seas_array(seas,3));...
                     RSSUN_c_mon(:,seas_array(seas,4));RSSUN_c_mon(:,seas_array(seas,5));RSSUN_c_mon(:,seas_array(seas,6))];
    RSSHA_e_seas(:,seas) = [RSSHA_e_mon(:,seas_array(seas,1));RSSHA_e_mon(:,seas_array(seas,2));RSSHA_e_mon(:,seas_array(seas,3));...
                     RSSHA_e_mon(:,seas_array(seas,4));RSSHA_e_mon(:,seas_array(seas,5));RSSHA_e_mon(:,seas_array(seas,6))];
    RSSHA_c_seas(:,seas) = [RSSHA_c_mon(:,seas_array(seas,1));RSSHA_c_mon(:,seas_array(seas,2));RSSHA_c_mon(:,seas_array(seas,3));...
                     RSSHA_c_mon(:,seas_array(seas,4));RSSHA_c_mon(:,seas_array(seas,5));RSSHA_c_mon(:,seas_array(seas,6))];
end
for i = 1:24
    ano_RSSUN1(i) = mean((RSSUN_e_seas(i:24:end,1)-RSSUN_c_seas(i:24:end,1)),'omitnan');
    ano_RSSUN2(i) = mean((RSSUN_e_seas(i:24:end,2)-RSSUN_c_seas(i:24:end,2)),'omitnan');
    ano_RSSHA1(i) = mean((RSSHA_e_seas(i:24:end,1)-RSSHA_c_seas(i:24:end,1)),'omitnan');
    ano_RSSHA2(i) = mean((RSSHA_e_seas(i:24:end,2)-RSSHA_c_seas(i:24:end,2)),'omitnan');
end
load('Fig_4_6_7_8_9_10_11_quantile.mat')
time = 0:1:23;
%% plot
% plot
hfig = figure('units','inches','position',[5 0 10 20]);
% subplot 1
sf1 = subplot(3,2,1);
plot(time,ano_CH2O,'linewidth',1.5,'MarkerSize',5,'Color',[0 0.45 0.74])
hold on
grid on
upper_C = dCH2O_Q3(9,:);
lower_C = dCH2O_Q1(9,:);
p_C = patch([time fliplr(time)],[upper_C,fliplr(lower_C)],[0 0.45 0.74]);
p_C.EdgeColor = 'none';
p_C.FaceAlpha = 0.2; % transparency of patch
plot([0 23],[0 0],'k--','linewidth',1)
xlim([0 23])
ylim([-0.1 0.7])
set(sf1,'xTick',(0:2:23),'fontsize',16);
set(sf1,'yTick',0:0.2:0.6,'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(mm)','fontsize',16);
title({'(a) Diurnal Cycle of Canopy Water';'Anomaly'},'fontsize',14,'Interpreter','none');
% subplot 2
sf2 = subplot(3,2,2);
plot(0:1:23,ano_Cevap,'linewidth',1.5,'MarkerSize',5,'Color','b')
hold on
grid on
plot(0:1:23,ano_Ctrans,'linewidth',1.5,'MarkerSize',5,'Color',[0 0.5 0])
upper_Ce = dCevap_Q3(9,:);
lower_Ce = dCevap_Q1(9,:);
p_Ce = patch([time fliplr(time)],[upper_Ce,fliplr(lower_Ce)],[0 0.45 0.74]);
p_Ce.EdgeColor = 'none';
p_Ce.FaceAlpha = 0.2; % transparency of patch
upper_Ct = dCtrans_Q3(9,:);
lower_Ct = dCtrans_Q1(9,:);
p_Ct = patch([time fliplr(time)],[upper_Ct,fliplr(lower_Ct)],[0 0.5 0]);
p_Ct.EdgeColor = 'none';
p_Ct.FaceAlpha = 0.3; % transparency of patch
plot([0 23],[0 0],'k--','linewidth',1)
xlim([0 23])
ylim([-40 70])
set(sf2,'xTick',(0:2:23),'fontsize',16);
set(sf2,'yTick',[-40:20:60],'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(W/m^2)','fontsize',16);
lh2 = legend({'Canopy E','Transpir.'},'location','NorthEast','fontsize',9);
title({'(b) Diurnal Cycle of Canopy Evapotranspiration';'Anomaly'},'fontsize',14,'Interpreter','none');
% subplot 3
sf3 = subplot(3,2,3);
plot(0:1:23,ano_Q2M,'--','linewidth',1.5,'MarkerSize',5,'Color',[0 0 1])
hold on
grid on
plot(0:1:23,ano_QsTV,'linewidth',1.5,'MarkerSize',5,'Color',[0 0.7 0])
xlim([0 23])
ylim([4 15])
set(sf3,'xTick',(0:2:23),'fontsize',16);
set(sf3,'yTick',(4:2:14),'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(g/kg)','fontsize',16);
lh3=legend({'Q @ 2m','Qs of vegetation T'},'location','NorthWest','fontsize',9,'Orientation','horizontal'); 
title({'(c) Diurnal Cycle of (Saturated) Specific Humidity';'Anomaly'},'fontsize',14,'Interpreter','none');
% subplot 4
sf4 = subplot(3,2,4);
plot(0:1:23,ano_RSSUN1,'color',[0.6350, 0.0780, 0.1840],'linewidth',1.5)
hold on
grid on
plot(0:1:23,ano_RSSUN2,'r','linewidth',1.5)
plot(0:1:23,ano_RSSHA1,'color',[0.5 0.5 0.5],'linewidth',1.5)
plot(0:1:23,ano_RSSHA2,'k','linewidth',1.5)
plot([0 23],[0 0],'k--','linewidth',1)
xlim([0 23])
ylim([-950 100])
set(sf4,'xTick',(0:2:23),'fontsize',16);
set(sf4,'yTick',[-800:200:0 100],'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(s/m)','fontsize',16);
lh4=legend({'SUNLIT (MAM,DJF)','SUNLIT (JJA,SON)','SHADED (MAM,DJF)','SHADED (JJA,SON)'},'location','South','fontsize',6.5);
title({'(d) Diurnal Cycle of Stomatal Resistance';'Anomaly'},'fontsize',14,'Interpreter','none');
% subplot 5
sf5 = subplot(3,2,5);
plot(0:1:23,ano_Ctrans1,'color',[0 0.35 0],'linewidth',1.5)
hold on
grid on
plot(0:1:23,ano_Ctrans2,'color',[0 0.85 0],'linewidth',1.5)
plot([0 23],[0 0],'k--','linewidth',1)
xlim([0 23])
ylim([-30 20])
set(sf5,'xTick',(0:2:23),'fontsize',16);
set(sf5,'yTick',(-30:10:20),'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(W/m^2)','fontsize',16);
legend({'MAM,DJF','JJA,SON'},'location','SouthEast','fontsize',9)
title({'(e) Diurnal Cycle of Canopy Transpiration';'Anomaly'},'fontsize',14,'Interpreter','none');
% subplot 6
sf6 = subplot(3,2,6);
plot(0:1:23,ano_photo1,'color',[0 0.35 0],'linewidth',1.5)
hold on
grid on
plot(0:1:23,ano_photo2,'color',[0 0.85 0],'linewidth',1.5)
plot([0 23],[0 0],'k--','linewidth',1)
xlim([0 23])
ylim([-3.5 2])
set(gca,'xTick',(0:2:23),'fontsize',16);
set(gca,'yTick',(-3:1:2),'fontsize',16);
xlabel('hour','fontsize',16);
ylabel('(\mumol/m^2s)','fontsize',16);
lh6=legend({'MAM,DJF','JJA,SON'},'location','SouthEast','fontsize',9,'Orientation','horizontal');
title({'(f) Diurnal Cycle of Photosynthesis';'Anomaly'},'fontsize',14,'Interpreter','none');
set(gcf,'PaperUnits','inches');
set(gcf,'PaperPosition',[0 0 18 18])
set(gcf,'PaperSize',[18 18])
img = gcf; 
img.PaperUnits = 'inches';  
img.PaperPosition = [0 0 18 18]; 
img.Units = 'inches'; 
img.PaperSize=[16 20];  
print(img,'-dpdf', '-r600', './Fig_9.pdf')
