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
39
40
41
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
97
98
99
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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
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 srie','simul parallle','potentiel srie','stress srie','potentiel parallle','stress parallle','air')
title('temprature 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 srie','simul parallle','potentiel srie','stress srie','potentiel parallle','stress parallle','air')
title('temprature vgtation')
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 srie','simul parallle','potentiel srie','stress srie','potentiel parallle','stress parallle','air')
title('temprature 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','srie','parallle')
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','vgtation srie','vgtation parallle','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('srie','parallle')
% 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 [-]')