% Fig 8
clear
close all

%% load
load('CL_TCCIP_c0_end_TRH_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)];

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)];
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

diurnal_TV_e=nan(24,1); diurnal_TV_c=nan(24,1);
diurnal_T2M_e=nan(24,1); diurnal_T2M_c=nan(24,1);
for i=1:24
diurnal_TV_e(i,1)=mean(TV_e(i:24:end),'omitnan');
diurnal_TV_c(i,1)=mean(TV_c(i:24:end),'omitnan');
diurnal_T2M_e(i,1)=mean(T2M_e(i:24:end),'omitnan');
diurnal_T2M_c(i,1)=mean(T2M_c(i:24:end),'omitnan');
end
load('Fig_4_6_7_8_9_10_11_quantile.mat');
time = 0:1:23;

%% plot
% plot
hfig = figure('units','inches','position',[0 0 15 12]);
% subplot 1
sf1 = subplot(2,2,1);
plot(time,ano_CH2O,'linewidth',2,'MarkerSize',5,'Color',[0 0.45 0.74])
hold on
grid on
upper_C = dCH2O_Q3(8,:);
lower_C = dCH2O_Q1(8,:);
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.5)
xlim([0 23])
ylim([-0.05 0.2])
set(sf1,'xTick',(0:2:23),'fontsize',20);
set(sf1,'yTick',-0.1:0.1:0.2,'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(mm)','fontsize',20);
title({'(a) Diurnal Cycle of Canopy Water';'Anomaly'},'fontsize',20,'Interpreter','none');
% subplot 2
sf2 = subplot(2,2,2);
plot(0:1:23,ano_Cevap,'linewidth',2,'MarkerSize',5,'Color','b')
hold on
grid on
plot(0:1:23,ano_Ctrans,'linewidth',2,'MarkerSize',5,'Color',[0 0.5 0])
upper_Ce = dCevap_Q3(8,:);
lower_Ce = dCevap_Q1(8,:);
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(8,:);
lower_Ct = dCtrans_Q1(8,:);
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.5)
xlim([0 23])
ylim([-20 40])
set(sf2,'xTick',(0:2:23),'fontsize',20);
set(sf2,'yTick',[-20:20:40],'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(W/m^2)','fontsize',20);
legend({'Canopy E','Transpir.'},'location','NorthEast','fontsize',15)
title({'(b) Diurnal Cycle of Canopy Evapotranspiration';'Anomaly'},'fontsize',20,'Interpreter','none');
% subplot 3
sf3 = subplot(2,2,3);
plot(0:1:23,ano_Q2M,'--','linewidth',2,'MarkerSize',5,'Color',[0 0 1])
hold on
grid on
plot(0:1:23,ano_QsTV,'linewidth',2,'MarkerSize',5,'Color',[0 0.7 0])
xlim([0 23])
ylim([0 5])
set(sf3,'xTick',(0:2:23),'fontsize',20);
set(sf3,'yTick',(0:1:5),'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(g/kg)','fontsize',20);
legend({'Q @ 2m','Qs of vegetation T'},'location','NorthWest','fontsize',15)
title({'(c) Diurnal Cycle of (Saturated) Specific Humidity';'Anomaly'},'fontsize',20,'Interpreter','none');
% subplot 4
sf4 = subplot(2,2,4);
plot(0:1:23,diurnal_T2M_e-diurnal_T2M_c,'linewidth',2,'color','b')
hold on
grid on
plot(0:1:23,diurnal_TV_e-diurnal_TV_c,'linewidth',2,'color',[0 0.7 0])
xlim([0 23])
ylim([1.5 3.5])
set(sf4,'xTick',(0:2:23),'fontsize',20);
set(sf4,'yTick',(1.5:0.5:3.5),'fontsize',20);
xlabel('hour','fontsize',20);
ylabel('(^oC)','fontsize',20);
legend({'T @ 2m','vegetation T'},'location','NorthWest','fontsize',15)
title({'(d) Diurnal Cycle of Temperature';'Anomaly'},'fontsize',20,'Interpreter','none');
set(gcf,'PaperUnits','inches');
set(gcf,'PaperPosition',[0 0 16 12])
set(gcf,'PaperSize',[16 12])
img = gcf; 
img.PaperUnits = 'inches';  
img.PaperPosition = [0 0 16 12]; 
img.Units = 'inches'; 
img.PaperSize=[16 12]; 
print(img, '-dpdf', '-r600', './Fig_8.pdf')

