main_test.m

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
    gLAIval=LAI0(:,2); % Grean LAI at those dates
    LAIval=LAI0(:,3); % Total LAI at those 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.99; % 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
[tsurfpot(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,rstmin,xg,sigmoy);

    % parallel
    
[tsurfpot2(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,rstmin,xg,sigmoy);

% run fully stressed conditions

    % series
[tsurfstress(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,rstmin,xg,sigmoy);

    %parallel
    
[tsurfstress2(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,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;

[tsurf(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,rstmin,xg,sigmoy);

% run patch/parallel version

[tsurf2(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,albe,emis,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,rstmin,xg,sigmoy);
%[tsurf2(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 (can be adusted according to pot and stress baselines)
 
% % bounding / stress 
 
%  rns(les<lesstress)=rnsstress(les<lesstress);
%  g(les<lesstress)=gstress(les<lesstress);
les(les<lesstress)=lesstress(les<lesstress);
hs(les<lesstress)=rns(les<lesstress)-g(les<lesstress)-lesstress(les<lesstress);

% rns(hs>hsstress)=rnsstress(hs>hsstress);
% g(hs>hsstress)=gstress(hs>hsstress);
hs(hs>hsstress)=hsstress(hs>hsstress);
les(hs>hsstress)=rns(hs>hsstress)-g(hs>hsstress)-hsstress(hs>hsstress);

% rnv(lev<levstress)=rnvstress(lev<levstress);
lev(lev<levstress)=levstress(lev<levstress);
hv(lev<levstress)=rnv(lev<levstress)-levstress(lev<levstress);

%rnv(hv>hvstress)=rnvstress(hv>hvstress);
%lev(hv>hvstress)=levstress(hv>hvstress);
hv(hv>hvstress)=hvstress(hv>hvstress);
lev(hv>hvstress)=rnv(hv>hvstress)-hvstress(hv>hvstress);

%***
%rns2(les2<lesstress)=rnsstress(les2<lesstress);
%g2(les2<lesstress)=gstress(les2<lesstress);
les2(les2<lesstress)=lesstress(les2<lesstress);
hs2(les2<lesstress)=rns2(les2<lesstress)-g2(les2<lesstress)-lesstress(les2<lesstress);

%rns2(hs2>hsstress)=rnsstress(hs2>hsstress);
%g2(hs2>hsstress)=gstress(hs2>hsstress);
hs2(hs2>hsstress)=hsstress(hs2>hsstress);
les2(hs2>hsstress)=rns2(hs2>hsstress)-g2(hs2>hsstress)-hsstress(hs2>hsstress);

 %rnv2(lev2<levstress)=rnvstress(lev2<levstress);
lev2(lev2<levstress)=levstress(lev2<levstress);
hv2(lev2<levstress)=rnv2(lev2<levstress)-levstress(lev2<levstress);

%rnv2(hv2>hvstress)=rnvstress(hv2>hvstress);
hv2(hv2>hvstress)=hvstress(hv2>hvstress);
lev2(hv2>hvstress)=rnv2(hv2>hvstress)-hvstress(hv2>hvstress);

% bounding / pot

lespot(lespot<10)=10;
hspot(lespot<10)=rnspot(lespot<10)-lespot(lespot<10);

%rns(les>lespot)=rnspot(les>lespot);
%g(les>lespot)=gpot(les>lespot);
les(les>lespot)=max(lespot(les>lespot),0);
hs(les>lespot)=rns(les>lespot)-g(les>lespot)-max(lespot(les>lespot),0);

%rnv(lev>levpot)=rnvpot(lev>levpot);
lev(lev>levpot)=levpot(lev>levpot);
hv(lev>levpot)=rnv(lev>levpot)-levpot(lev>levpot);

%rns2(les2>lespot)=rnspot(les2>lespot);
%g2(les2>lespot)=gpot(les2>lespot);
hs2(les2>lespot)=rns2(les2>lespot)-g2(les2>lespot)-max(lespot(les2>lespot),0);
les2(les2>lespot)=max(lespot(les2>lespot),0);

%rnv2(lev2>levpot)=rnvpot(lev2>levpot);
hv2(lev2>levpot)=rnv2(lev2>levpot)-levpot(lev2>levpot);
lev2(lev2>levpot)=levpot(lev2>levpot);

hs(hs<hspot)=max(hspot(hs<hspot),-100);
les(hs<hspot)=rns(hs<hspot)-g(hs<hspot)-max(hspot(hs<hspot),-100);

hv(hv<hvpot)=max(hvpot(hv<hvpot),-100);
lev(hv<hvpot)=rnv(hv<hvpot)-g(hv<hvpot)-max(hspot(hv<hvpot),-100);

hs2(hs2<hspot)=max(hspot(hs2<hspot),-100);
les2(hs2<hspot)=rns2(hs2<hspot)-g2(hs2<hspot)-max(hspot(hs2<hspot),-100);

hv2(hv2<hvpot)=max(hvpot(hv2<hvpot),-100);
lev2(hv2<hvpot)=rnv2(hv2<hvpot)-g(hv2<hvpot)-max(hspot(hv2<hvpot),-100);

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

% Rott Mean Square Error on Latent Heat LE, Sensible Heat H, Net radiation
% Rn and Soil Heat Flux G

RMSE_LE = sqrt(mean((leobs(time2) - le(time2)).^2))  % RMSE Layer
RMSE_LE2 = sqrt(mean((leobs(time2) - le2(time2)).^2)) % RMSE Patch

RMSE_H = sqrt(mean((hobs(time2) - h(time2)).^2))  % RMSE Layer 
RMSE_H2 = sqrt(mean((hobs(time2) - h2(time2)).^2)) % RMSE Patch

RMSE_RN = sqrt(mean((rnobs(time2) - rn(time2)).^2))  % RMSE Layer 
RMSE_RN2 = sqrt(mean((rnobs(time2) - rn2(time2)).^2)) % RMSE Patch

RMSE_G = sqrt(mean((gobs(time2) - g(time2)).^2))  % RMSE Layer 
RMSE_G2 = sqrt(mean((gobs(time2) - g2(time2)).^2)) % RMSE Patch

figure(1)
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(2)
scatter(1-leobs(time2)./lepot(time2),1-le(time2)./lepot(time2),'bo')
hold on
scatter(1-leobs(time2)./lepot2(time2),1-le2(time2)./lepot2(time2),'rs')
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(3)
plot(doy(time2),leobs(time2),'rs',doy(time2),le(time2),'bo',doy(time2),le2(time2),'g*')
axis([0 180 0 500]);
legend('observ','srie','parallle')
xlabel('Date [DOY]')
ylabel('LE [W/m2]')

weight=0.2*ones(1,5);

t=doy(time2);
es=les(time2);
es2=les2(time2);
esp=lespot(time2);
esp2=lespot2(time2);

esf=filter(weight,1,es);
es2f=filter(weight,1,es2);
espf=filter(weight,1,esp);
esp2f=filter(weight,1,esp2);

figure(4)
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(5)
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(6)
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(7)
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(8)
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]')