close all;
% read climate forcing etc
meteo=csvread('meteoBS.csv');
n=size(meteo,1);
day(1:n)=meteo(1:n,1); % day of year
hour(1:n)=meteo(1:n,2)+1; % local hour
minute(1:n)=meteo(1:n,3); % minute
doy(1:n)=meteo(1:n,4); % decimal day of year
rg(1:n)=meteo(1:n,5); % global radiation[W/m2]
ta(1:n)=meteo(1:n,6)+273.15; % air temperature at reference level [�C > �K]
rh(1:n)=meteo(1:n,7); % relative humidity at reference level [%]
ua(1:n)=max(meteo(1:n,8),0.1); % wind speed at reference level [m/s]
tsobs(1:n)=meteo(1:n,9); % observed radiative surface temperature [�C]
leobs(1:n)=meteo(1:n,10); % observed latent heat flux (residual correction) [W/m2]
rnobs(1:n)=meteo(1:n,13); % observed net radiation [W/m2]
gobs(1:n)=meteo(1:n,14); % observed soil heat flux [W/m2]
hobs(1:n)=meteo(1:n,15); % observed sensible heat flux[W/m2]
ratmobs(1:n)=meteo(1:n,16); % atmospheric radiation [W/m2]
vza=0; % view zenith angle [rad]
% read vegetation height
zf0=csvread('H_ble.csv');
zfdate=zf0(:,1); % dates de mesure de la hauteur
zfval=zf0(:,2); % valeur de la hauteur � ces dates [m]
zf=interp1(zfdate,zfval,doy,'pchip','extrap');
% read total and green LAI
LAI0=csvread('LAI_ble.csv');
LAIdate=LAI0(:,1); % dates de mesure de l'indice foliaire
gLAIval=LAI0(:,2); % valeur de l'indice foliaire � ces dates
LAIval=LAI0(:,3); % valeur de l'indice foliaire total � ces dates
glai=interp1(LAIdate,gLAIval,doy,'pchip','extrap');
lai=interp1(LAIdate,LAIval,doy,'pchip','extrap');
% read parameters
rstmin=100; % minimum stomatal resistance (s/m)
albe=0.25; % surface albedo (BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
xg=0.4; % soil heat flux to soil net radiation ratio at midday (BEWARE: should vary diurnally / code to be modified in that case)
za=2.32; % reference height for the climate forcing [m]
emis=0.98; % surface emissivity (BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
patm=990; % atmospheric pressure [hPa](BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
sigmoy=0.45; % attenuation coef. of incoming radiation within the canopy (0.9 used in the layer case for longwave)
sigma=5.67e-8; % Stafan-Boltzman constant [W m-2 K-4]
% Efficiencies for fully stressed (0) and potential (1) conditions
betapot=1;
betastress=0.01;
for i=1:n % time loop
% run potential conditions
% series
[tradpot(i),tspot(i),tvpot(i),t0pot(i),rnspot(i),rnvpot(i),gpot(i),hspot(i),hvpot(i),lespot(i),levpot(i)]=...
SPARSE(betapot,betapot,2,1,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
% parallel
[tradpot2(i),tspot2(i),tvpot2(i),t0pot2(i),rnspot2(i),rnvpot2(i),gpot2(i),hspot2(i),hvpot2(i),lespot2(i),levpot2(i)]=...
SPARSE(betapot,betapot,2,2,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
% run fully stressed conditions
% series
[tradstress(i),tsstress(i),tvstress(i),t0stress(i),rnsstress(i),rnvstress(i),gstress(i),hsstress(i),hvstress(i),lesstress(i),levstress(i)]=...
SPARSE(betastress,betastress,2,1,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
%parallel
[tradstress2(i),tsstress2(i),tvstress2(i),t0stress2(i),rnsstress2(i),rnvstress2(i),gstress2(i),hsstress2(i),hvstress2(i),lesstress2(i),levstress2(i)]...
=SPARSE(betastress,betastress,2,2,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
% run layer/series version
%tsobs2(i)=((sigma*(tsobs(i)+273.15)^4-(1-emis)*ratmobs(i))/(emis*sigma)).^0.25-273.15;
[trad(i),ts(i),tv(i),toto(i),rns(i),rnv(i),g(i),hs(i),hv(i),les(i),lev(i)]=...
SPARSE(tsobs(i),vza,1,1,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
% run patch/parallel version
[trad2(i),ts2(i),tv2(i),t02(i),rns2(i),rnv2(i),g2(i),hs2(i),hv2(i),les2(i),lev2(i)]=SPARSE(tsobs(i),vza,1,2,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
%[trad2(i),ts2(i),tv2(i),rns2(i),rnv2(i),g2(i),hs2(i),hv2(i),les2(i),lev2(i)]=TSEBKustas2000(ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,xg,emis,sigmoy,tsobs(i));
end
% bounding / stress
rns(les<lesstress)=rnsstress(les<lesstress);
g(les<lesstress)=gstress(les<lesstress);
hs(les<lesstress)=hsstress(les<lesstress);
%hv(les<lesstress)=rnv(les<lesstress)+rnsstress(les<lesstress)-gstress(les<lesstress)-hsstress(les<lesstress)-lev(les<lesstress);
les(les<lesstress)=lesstress(les<lesstress);
les(hs>hsstress)=lesstress(hs>hsstress);
rns(hs>hsstress)=rnsstress(hs>hsstress);
g(hs>hsstress)=gstress(hs>hsstress);
%hv(hs>hsstress)=rnv(hs>hsstress)+rnsstress(hs>hsstress)-gstress(hs>hsstress)-hsstress(hs>hsstress)-lev(hs>hsstress);
hs(hs>hsstress)=hsstress(hs>hsstress);
% rns2(les2<lesstress2)=rnsstress2(les2<lesstress2);
% g2(les2<lesstress2)=gstress2(les2<lesstress2);
% hs2(les2<lesstress2)=hsstress2(les2<lesstress2);
% %hv2(les2<lesstress2)=rnv2(les2<lesstress2)+rnsstress2(les2<lesstress2)-gstress2(les2<lesstress2)-hsstress2(les2<lesstress2)-lev2(les2<lesstress2);
% les2(les2<lesstress2)=lesstress2(les2<lesstress2);
%
% les2(hs2>hsstress2)=lesstress2(hs2>hsstress2);
% rns2(hs2>hsstress2)=rnsstress2(hs2>hsstress2);
% g2(hs2>hsstress2)=gstress2(hs2>hsstress2);
% %hv2(hs2>hsstress2)=rnv2(hs2>hsstress2)+rnsstress2(hs2>hsstress2)-gstress2(hs2>hsstress2)-hsstress2(hs2>hsstress2)-lev2(hs2>hsstress2);
% hs2(hs2>hsstress2)=hsstress2(hs2>hsstress2);
rnv(lev<levstress)=rnvstress(lev<levstress);
hv(lev<levstress)=hvstress(lev<levstress);
%hv(lev<levstress)=rnvstress(lev<levstress)+rns(lev<levstress)-g(lev<levstress)-hs(lev<levstress)-les(lev<levstress)-levstress(lev<levstress);
lev(lev<levstress)=levstress(lev<levstress);
rnv(hv>hvstress)=rnvstress(hv>hvstress);
lev(hv>hvstress)=levstress(hv>hvstress);
%hv(hv>hvstress)=rnvstress(hv>hvstress)+rns(hv>hvstress)-g(hv>hvstress)-hs(hv>hvstress)-les(hv>hvstress)-levstress(hv>hvstress);
hv(hv>hvstress)=hvstress(hv>hvstress);
% rnv2(lev2<levstress2)=rnvstress2(lev2<levstress2);
% hv2(lev2<levstress2)=hvstress2(lev2<levstress2);
% %hv2(lev2<levstress2)=rnvstress2(lev2<levstress2)+rns2(lev2<levstress2)-g2(lev2<levstress2)-hs2(lev2<levstress2)-les2(lev2<levstress2)-levstress2(lev2<levstress2);
% lev2(lev2<levstress2)=levstress2(lev2<levstress2);
%
% rnv2(hv2>hvstress2)=rnvstress2(hv2>hvstress2);
% lev2(hv2>hvstress2)=levstress2(hv2>hvstress2);
% %hv2(hv2>hvstress2)=rnvstress2(hv2>hvstress2)+rns2(hv2>hvstress2)-g2(hv2>hvstress2)-hs2(hv2>hvstress2)-les2(hv2>hvstress2)-levstress2(hv2>hvstress2);
% hv2(hv2>hvstress2)=hvstress2(hv2>hvstress2);
% bounding / pot
rns(les>lespot)=rnspot(les>lespot);
g(les>lespot)=gpot(les>lespot);
hs(les>lespot)=hspot(les>lespot);
%hv(les>lespot)=rnv(les>lespot)+rnspot(les>lespot)-gpot(les>lespot)-hspot(les>lespot)-lev(les>lespot);
les(les>lespot)=lespot(les>lespot);
% rns2(les2>lespot2)=rnspot2(les2>lespot2);
% g2(les2>lespot2)=gpot2(les2>lespot2);
% hs2(les2>lespot2)=hspot2(les2>lespot2);
% %hv2(les2>lespot2)=rnv2(les2>lespot2)+rnspot2(les2>lespot2)-gpot2(les2>lespot2)-hspot2(les2>lespot2)-lev2(les2>lespot2);
% les2(les2>lespot2)=lespot2(les2>lespot2);
rnv(lev>levpot)=rnvpot(lev>levpot);
hv(lev>levpot)=hvpot(lev>levpot);
%hv(lev>levpot)=rnvpot(lev>levpot)+rns(lev>levpot)-g(lev>levpot)-hs(lev>levpot)-les(lev>levpot)-levpot(lev>levpot);
lev(lev>levpot)=levpot(lev>levpot);
% rnv2(lev2>levpot2)=rnvpot2(lev2>levpot2);
% hv2(lev2>levpot2)=hvpot2(lev2>levpot2);
% %hv2(lev2>levpot2)=rnvpot2(lev2>levpot2)+rns2(lev2>levpot2)-g2(lev2>levpot2)-hs2(lev2>levpot2)-les2(lev2>levpot2)-levpot2(lev2>levpot2);
% lev2(lev2>levpot2)=levpot2(lev2>levpot2);
% computes total fluxes
rn=rns+rnv;
le=les+lev;
h=hs+hv;
rn2=rns2+rnv2;
le2=les2+lev2;
h2=hs2+hv2;
lepot=lespot+levpot;
lepot2=lespot2+levpot2;
hstress=hsstress+hvstress;
hstress2=hsstress2+hvstress2;
hpot=hspot+hvpot;
hpot2=hspot2+hvpot2;
% extract values at midday
a=doy-floor(doy);
time=find(a>0.485-(1/24) & a<0.515-(1/24));
time2=find(a>0.485-(1/24) & a<0.515-(1/24) & leobs >-100 & hobs>-100 & rnobs >-100 & gobs> - 100);
time3=time2;
RMSE_LE = sqrt(mean((leobs(time2) - le(time2)).^2)) % RMSE
RMSE_LE2 = sqrt(mean((leobs(time2) - le2(time2)).^2)) % RMSE 2
RMSE_H = sqrt(mean((hobs(time2) - h(time2)).^2)) % RMSE
RMSE_H2 = sqrt(mean((hobs(time2) - h2(time2)).^2)) % RMSE 2
RMSE_RN = sqrt(mean((rnobs(time2) - rn(time2)).^2)) % RMSE
RMSE_RN2 = sqrt(mean((rnobs(time2) - rn2(time2)).^2)) % RMSE 2
RMSE_G = sqrt(mean((gobs(time2) - g(time2)).^2)) % RMSE
RMSE_G2 = sqrt(mean((gobs(time2) - g2(time2)).^2)) % RMSE 2
figure(1)
plot(doy(time),ts(time),'bo',doy(time),ts2(time),'bs',doy(time),tspot(time),'go',doy(time),tsstress(time),'ro',doy(time),tspot2(time),'gs',doy(time),tsstress2(time),'rs',doy(time),ta(time)-273.15,':k')
legend('simul� s�rie','simul� parall�le','potentiel s�rie','stress s�rie','potentiel parall�le','stress parall�le','air')
title('temp�rature sol')
figure(2)
plot(doy(time),tv(time),'bo',doy(time),tv2(time),'bs',doy(time),tvpot(time),'go',doy(time),tvstress(time),'ro',doy(time),tvpot2(time),'gs',doy(time),tvstress2(time),'rs',doy(time),ta(time)-273.15,':k')
legend('simul� s�rie','simul� parall�le','potentiel s�rie','stress s�rie','potentiel parall�le','stress parall�le','air')
title('temp�rature v�g�tation')
figure(3)
plot(doy(time),trad(time),'bo',doy(time),trad2(time),'bs',doy(time),tradpot(time),'go',doy(time),tradstress(time),'ro',doy(time),tradpot2(time),'gs',doy(time),tradstress2(time),'rs',doy(time),ta(time)-273.15,':k',doy(time),tsobs(time),'ok')
legend('simul� s�rie','simul� parall�le','potentiel s�rie','stress s�rie','potentiel parall�le','stress parall�le','air')
title('temp�rature radiative')
figure(4)
plot(leobs(time2),le(time2),'bo',leobs(time2),le2(time2),'rs')
hold on
lsline
axis([-100 500 -100 500]);
legend('Series model','Parallel model')
xlabel('Observed Latent Heat Flux at midday [W/m2]')
ylabel('Retrieved Latent Heat Flux at midday [W/m2]')
% figure(12)
% %plot(doy(time),lev(time)./levpot(time),'ko',doy(time),lev2(time)./levpot(time),'bo',...
% % doy(time),les(time)./lespot(time),'k*',doy(time),les2(time)./lespot(time),'b*')
% % plot(doy(time),lev(time)./levpot(time),'go',doy(time),lev(time)./levstress(time),'gs',...
% % doy(time),les(time)./lespot(time),'ro',doy(time),les(time)./lesstress(time),'rs')
% plot(doy(time),Hv(time)./Hvpot(time),'go',doy(time),Hv(time)./Hvstress(time),'gs',...
% doy(time),Hs(time)./Hspot(time),'ro',doy(time),Hs(time)./Hsstress(time),'rs')
% axis([0 180 -1 2]);
figure(5)
plot(doy(time2),leobs(time2)./lepot(time2),'rs',doy(time2),le(time2)./lepot(time2),'bo',doy(time2),le2(time2)./lepot2(time2),'g*',doy(time2),lai(time2)/2,'g-')
%plot(doy(time2),leobs(time2).*2./(lepot(time2)+lepot2(time2)),'rs',doy(time2),lereal(time2)./lepot(time2),'ro',doy(time2),lereal2(time2)./lepot2(time2),'r*',doy(time2),le(time2)./lepot(time2),'bo',doy(time2),le2(time2)./lepot2(time2),'g*')
%axis([0 180 -0.2 1.2]);
legend('observed','series','parallel','LAI/2')
xlabel('Date [DOY]')
ylabel('Evaporative Efficiency LE/LEpot [-], LAI/2 [-]')
figure(6)
%plot(1-leobs(time3)./lepot(time3),1-le(time3)./lepot(time3),'bo',1-leobs(time3)./lepot2(time3),1-le2(time3)./lepot2(time3),'rs')
scatter(1-leobs(time3)./lepot(time3),1-le(time3)./lepot(time3),lepot(time3)/2,'bo')
hold on
scatter(1-leobs(time3)./lepot2(time3),1-le2(time3)./lepot2(time3),lepot2(time3)/2,'rs')
%lsline
axis([-0.2 1.2 -0.2 1.2]);
legend('Series model','Parallel model')
xlabel('Observed surface stress at midday (-)')
ylabel('Simulated surface stress at midday (-)')
figure(7)
plot(doy(time2),leobs(time2),'rs',doy(time2),le(time2),'bo',doy(time2),le2(time2),'g*')
axis([0 180 0 500]);
legend('observ�','s�rie','parall�le')
xlabel('Date [DOY]')
ylabel('LE [W/m2]')
save('bounded.mat','doy','time2','le','le2')
weight=0.2*ones(1,5);
t=doy(time3);
es=les(time3);
es2=les2(time3);
esp=lespot(time3);
esp2=lespot2(time3);
esf=filter(weight,1,es);
es2f=filter(weight,1,es2);
espf=filter(weight,1,esp);
esp2f=filter(weight,1,esp2);
figure(8)
%plot(doy(time2),les(time2),'bo',doy(time2),les2(time2),'g*',dateobs+0.5,pluie*10,'k-',doy(time2),lesobs(time2),'r:',doy(time2),leobs(time2),'rs')
plot(t,esf,'bo:',t,es2f,'rs:',t,espf,'g:',t,esp2f,'g-')
axis([0 180 0 250]);
legend('Series model','Parallel model','Potential serie','potential parallel')
xlabel('Date [DOY]')
ylabel('Soil Evaporation, 5 days running average [W/m2]')
figure(12)
plot(t,esf./espf,'bo:',t,es2f./esp2f,'rs:')
axis([0 180 -0.1 1.1]);
legend('Series model','Parallel model')
xlabel('Date [DOY]')
ylabel('Soil Evaporation Efficiency, 5 days running average [-]')
% figure(9)
% plot(doy(time2),lev2(time2),'bo',doy(time2),lev(time2),'b-',doy(time2),levpot(time2),'g:')
% legend('Old Srface Energy Balance Model TSEB','New Surface Energy Balance Model TSEB2','potentiel')
% xlabel('Date [DOY]')
% ylabel('LE [W/m2]')
figure(9)
plot(doy(time2),(h(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'b-',doy(time2),(hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'bo',...
doy(time2),1-le(time2)./lepot(time2),'r-',doy(time2),1-leobs(time2)./lepot(time2),'rs')
%axis([0 180 -0.2 1.2]);
legend('Stress H simul�','Stress H observ�','Stress LE simul�','Stress LE observ�')
xlabel('Date [DOY]')
ylabel('Stress [-]')
figure(10)
plot(doy(time2),hv2(time2),'b-',doy(time2),hvpot2(time2),'r-',doy(time2),hvstress2(time2),'g-')
axis([0 180 -200 500]);
figure(11)
plot(doy(time2),hv(time2),'b-',doy(time2),hvpot(time2),'r-',doy(time2),hvstress(time2),'g-')
axis([0 180 -200 500]);
figure(13)
plot((hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),(h(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'bo',...
(hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),(h2(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'rs')
hold on
lsline
axis([-0.2 1.2 -0.2 1.2]);
legend('Series model','Parallel model')
xlabel('Observed surface stress at midday (-)')
ylabel('Simulated surface stress at midday (-)')
% figure(13)
% plot(doy(time2),leobs(time2),'rs',doy(time2),lev(time2),'bo',doy(time2),lev2(time2),'g*',dateobs+0.5,pluie*10,'k:')
% axis([0 180 -200 500]);
% legend('observ� total','v�g�tation s�rie','v�g�tation parall�le','pluie')
% xlabel('Date [DOY]')
% ylabel('LE [W/m2], pluie (*0.1 mm)')
% figure(14)
% plot(lesobs(time2),les(time2),'bo',lesobs(time2),les2(time2),'g*')
% hold on
% lsline
% axis([-100 400 -100 400]);
% title('�vaporation du sol')
% legend('s�rie','parall�le')
% xlabel('Flux observ� � midi (W/m2)')
% ylabel('Flux simul� � midi (W/m2)')
figure(15)
plot(hobs(time2),h(time2),'bo',hobs(time2),h2(time2),'g*')
hold on
lsline
axis([-100 500 -100 500]);
legend('series','parallel')
xlabel('Observed Sensible Heat Flux at Midday [W/m2]')
ylabel('Simulated Sensible Heat Flux at Midday [W/m2]')
figure(16)
plot(rnobs(time2),rn(time2),'bo',rnobs(time2),rn2(time2),'g*')
hold on
lsline
axis([-100 700 -100 700]);
legend('series','parallel')
xlabel('Observed Net Radiation at Midday [W/m2]')
ylabel('Simulated Net Radition at Midday [W/m2]')
figure(18)
plot(gobs(time2),g(time2),'bo',gobs(time2),g2(time2),'g*')
hold on
lsline
axis([-50 300 -50 300]);
legend('series','parallel')
xlabel('Observed Soil Heat Flux at Midday [W/m2]')
ylabel('Observed Soil Heat Flux at Midday [W/m2]')
figure(21)
plot(doy(time3),1-leobs(time3)./lepot(time3),'rs',doy(time3),(tsobs(time3)-tradpot(time3))./(tradstress(time3)-tradpot(time3)),'b-',...
doy(time3),(tsobs(time3)-tradpot2(time3))./(tradstress2(time3)-tradpot2(time3)),'g-')
axis([0 180 -0.2 1.2]);
legend('Stress observ�','Stress tso-tsp/ts0-tsp serie','Stress tso-tsp/ts0-tsp parallel')
xlabel('Date [DOY]')
ylabel('Stress [-]')
figure(22)
plot(1-leobs(time2)./lepot(time2),(tsobs(time2)-tradpot(time2))./(tradstress(time2)-tradpot(time2)),'bo',...
1-leobs(time2)./lepot(time2),1-le(time2)./lepot(time2),'go')
hold on
lsline
axis([-0.2 1.2 -0.2 1.2]);
legend('Stress tso-tsp/ts0-tsp serie','Stress 1-LE/LEpot')
xlabel('Stress observ� [-]')
ylabel('Stress simul� [-]')
figure(17)
plot(1-leobs(time2)./lepot(time2),(hobs(time2)-hpot(time2))./(hstress(time2)-hpot(time2)),'rs')
hold on
lsline
axis([-0.2 1.2 -0.2 1.2]);
legend('Stress LE','Stress H')
xlabel('Stress [-]')
ylabel('Stress [-]')