diff --git a/Draft_SPARSE_USERs_Manual.pdf b/Draft_SPARSE_USERs_Manual.pdf
index 33755c2eb956ebb214a1e53b0772a7b3af6cc5c3..494c9ba0308775a54dfca0d0c5e1e8577d500e10 100644
Binary files a/Draft_SPARSE_USERs_Manual.pdf and b/Draft_SPARSE_USERs_Manual.pdf differ
diff --git a/SPARSE.m b/SPARSE.m
index d16eb20486377e54bd654b378be43be6d6239a47..9ae187c10100f671eaf122ff17e464112f507ce0 100644
--- a/SPARSE.m
+++ b/SPARSE.m
@@ -1,4 +1,4 @@
-function [trad,ts,tv,t0,rns,rnv,g,hs,hv,les,lev]=SPARSE(input1,input2,typerun,typeohm,ta,rh,rg,ua,glai,lai,zf,za,albe,rvvmin,xg,sigmoy)
+function [tsurf,ts,tv,t0,rns,rnv,g,hs,hv,les,lev]=SPARSE(input1,input2,typerun,typeohm,albe,emis,ta,rh,rg,ua,glai,lai,zf,za,rvvmin,xg,sigmoy)
% function that computes actual or bounding transpiration and soil evaporation with
% the SPARSE framework, i.e. Shuttleworht-Wallace dual source model with
@@ -12,7 +12,6 @@ function [trad,ts,tv,t0,rns,rnv,g,hs,hv,les,lev]=SPARSE(input1,input2,typerun,ty
if typerun==1 % retrieval:
% imposed observed radiative temperature and view zenith angle
- tsobs=input1;
vza=input2;
betas=1;
betav=1;
@@ -29,6 +28,7 @@ end
n=size(glai,1);
m=size(glai,2);
+% emis = surface emissivity [-]
% ta = air temperature [K]
% rh = air relative humidity [%]
% rg = incoming solar radiation [W/m2]
@@ -43,7 +43,8 @@ m=size(glai,2);
% sigmoy = Beer Lambert law attenuation coefficient
% vza = view zenith angle [rad]
-% RETRIEVAL MODE: tsobs = OBSERVED RADIATIVE surface temperature [°C]
+% RETRIEVAL MODE: input1 = OBSERVED RADIATIVE surface temperature [°C] % A
+% PRECISER
% PRESCRIBED mode: betas: relative stress factor for the soil [-]
% PRESCRIBED mode: betav: relative stress factor for the vegetation [-]
@@ -269,7 +270,12 @@ A3_1=zeros(n,m);
A3_2=-aras-arav;
A3_3=-bras-brav;
-B3=sigma*ones(n,m).*(tsobs+273.15).^4+cras+crav-ratm*ones(n,m);
+% 3 cases according to the type of data:
+%
+
+Mrad=(1-emis).*ratm+sigma*ones(n,m).*emis.*(input1+273.15).^4;
+
+B3=Mrad+cras+crav-ratm*ones(n,m);
det=A1_1.*A2_2.*A3_3 - A1_1.*A2_3.*A3_2 - A1_2.*A2_1.*A3_3 + A1_2.*A2_3.*A3_1 + A1_3.*A2_1.*A3_2 - A1_3.*A2_2.*A3_1;
@@ -346,7 +352,7 @@ A3_1(les<lesmin)=0;
A3_2(les<lesmin)=-aras(les<lesmin)-arav(les<lesmin);
A3_3(les<lesmin)=-bras(les<lesmin)-brav(les<lesmin);
-B3(les<lesmin)=sigma*(tsobs(les<lesmin)+273.15)^4+cras(les<lesmin)+crav(les<lesmin)-ratm;
+B3(les<lesmin)=Mrad(les<lesmin)+cras(les<lesmin)+crav(les<lesmin)-ratm;
det(les<lesmin)=A1_1(les<lesmin).*A2_2(les<lesmin).*A3_3(les<lesmin)...
- A1_1(les<lesmin).*A2_3(les<lesmin).*A3_2(les<lesmin) - A1_2(les<lesmin).*A2_1(les<lesmin).*A3_3(les<lesmin) + A1_2(les<lesmin).*A2_3(les<lesmin)...
@@ -501,7 +507,7 @@ end
end % end stability loop
- trad=((ratm.*ones(n,m)-cras-crav-X2.*(aras+arav)-X3.*(bras+brav))/sigma).^0.25-273.15;
+ tsurf=((emis.*ratm.*ones(n,m)-cras-crav-X2.*(aras+arav)-X3.*(bras+brav))/(sigma.*emis)).^0.25-273.15;
rns=(crns+arns.*X2+brns.*X3)./(1-xg);
rnv=crnv+arnv.*X2+brnv.*X3;
diff --git a/index-1.html b/index-1.html
new file mode 100644
index 0000000000000000000000000000000000000000..ed755f88cbe94020cb528117c222ec7d2d68e902
--- /dev/null
+++ b/index-1.html
@@ -0,0 +1,32 @@
+<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"/><title>Inconnu(e) </title></head><body>
+<h1 id="sparse-model"><strong>SPARSE Model</strong></h1>
+<hr />
+<p>Soil Plant Atmosphere Remote Sensing Evapotranspiration model.</p>
+<h2 id="license">License:</h2>
+<p>This is free software under the GNU General Public License v3.0.
+See https://www.gnu.org/licenses/gpl-3.0.en.html for details or the LICENSE.txt file</p>
+<h2 id="author">Author:</h2>
+<p>Gilles Boulet
+CESBIO UMR5126, (CNRS,UPS,CNES,IRD)
+13 avenue du Colonel Roche
+31401 Toulouse cedex 9
+http://www.cesbio.ups-tlse.fr/</p>
+<h2 id="contact">Contact:</h2>
+<p>gilles.boulet@ird.fr</p>
+<h2 id="acknowledgments">Acknowledgments:</h2>
+<p>The present study received support from the CNES (Centre National d’Etude Spatiale, France).</p>
+<h2 id="reference-for-the-model">Reference for the model:</h2>
+<p>Boulet, G., Mougenot, B., Lhomme, J.P., Fanise, P., Lili-Chabaane, Z., Olioso, A., Bahir, M., Rivalland, V., Jarlan, L., Merlin, O., Coudert, B., Er-Raki, S., Lagouarde, J.P., 2015. The SPARSE model for the prediction of water stress and evapotranspiration components from thermal infra-red data and its evaluation over irrigated and rainfed wheat. Hydrol. Earth Syst. Sci., 19(11): 4653-4672.</p>
+<p>Full-Text: http://www.hydrol-earth-syst-sci.net/19/4653/2015/hess-19-4653-2015.pdf</p>
+<h2 id="summary">Summary:</h2>
+<p>The SPARSE model is coded in Matlab language; it computes the energy budget components for two sources (soil and vegetation) from surface (Leaf Area Index, vegetation height) and meterological forcing data.</p>
+<h2 id="installation-requirements">Installation requirements:</h2>
+<ul>
+<li>respect the basic principles of GPL v3 licence.
+ https://www.gnu.org/licenses/gpl-3.0-standalone.html</li>
+<li>contact the author if bugs are found :
+ mailto:gilles.boulet@ird.fr </li>
+</ul>
+<h2 id="example-file">Example file:</h2>
+<p>main_test.m main sample file to run the model coded as funtion SPARSE.m</p>
+</body></html>
\ No newline at end of file
diff --git a/main_test.m b/main_test.m
index 741f00f0794694bedeb3e86fa4788cc17ba20bdf..251cabb54bb8b387d57692546ece617bc28d6371 100644
--- a/main_test.m
+++ b/main_test.m
@@ -48,7 +48,7 @@ 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)
+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]
@@ -62,109 +62,114 @@ betastress=0.01;
% 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);
+[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
-[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);
+[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
-[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);
+[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
-[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);
+[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;
-[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);
+[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
-[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));
+[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);
-hs(les<lesstress)=hsstress(les<lesstress);
-%hv(les<lesstress)=rnv(les<lesstress)+rnsstress(les<lesstress)-gstress(les<lesstress)-hsstress(les<lesstress)-lev(les<lesstress);
+% % 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);
-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);
+% 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);
-% 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);
+% 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)=rnvstress(hv>hvstress)+rns(hv>hvstress)-g(hv>hvstress)-hs(hv>hvstress)-les(hv>hvstress)-levstress(hv>hvstress);
+%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);
-% 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);
+%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
-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);
+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);
-% 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);
+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;
@@ -187,7 +192,7 @@ 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
+% 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