Skip to content
Snippets Groups Projects
Commit c636018b authored by gilles.boulet_ird.fr's avatar gilles.boulet_ird.fr
Browse files

Transfer to new github server, few changes to bounding conditions in main-test.m

parent 7c0ed8f3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Draft_SPARSE_USERs_Manual.pdf
+ 0
− 0
View file @ c636018b
No preview for this file type
SPARSE.m
+ 12
− 6
View file @ c636018b
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;
......
index-1.html 0 → 100644
+ 32
− 0
View file @ c636018b
<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
main_test.m
+ 75
− 70
View file @ c636018b
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment