Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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 (can be adusted according to pot and stress baselines)
% bounding / stress
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
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);
% 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
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 (-)')
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);
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]')
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 [-]')
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]')
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]')
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]')