Commit 1b6bc4a5b0bb792a6c2f9a085869a7e20ce9b93f

Authored by Gilles Boulet
1 parent 94fdac3a
Exists in master

Add surface temp as input or output instead of the sole radiative temp

Showing 2 changed files with 25 additions and 19 deletions   Show diff stats
SPARSE.m
1   -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)
  1 +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)
2 2  
3 3 % function that computes actual or bounding transpiration and soil evaporation with
4 4 % 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
12 12  
13 13 if typerun==1 % retrieval:
14 14 % imposed observed radiative temperature and view zenith angle
15   - tsobs=input1;
16 15 vza=input2;
17 16 betas=1;
18 17 betav=1;
... ... @@ -29,6 +28,7 @@ end
29 28 n=size(glai,1);
30 29 m=size(glai,2);
31 30  
  31 +% emis = surface emissivity [-]
32 32 % ta = air temperature [K]
33 33 % rh = air relative humidity [%]
34 34 % rg = incoming solar radiation [W/m2]
... ... @@ -43,7 +43,8 @@ m=size(glai,2);
43 43 % sigmoy = Beer Lambert law attenuation coefficient
44 44 % vza = view zenith angle [rad]
45 45  
46   -% RETRIEVAL MODE: tsobs = OBSERVED RADIATIVE surface temperature [°C]
  46 +% RETRIEVAL MODE: input1 = OBSERVED RADIATIVE surface temperature [°C] % A
  47 +% PRECISER
47 48  
48 49 % PRESCRIBED mode: betas: relative stress factor for the soil [-]
49 50 % PRESCRIBED mode: betav: relative stress factor for the vegetation [-]
... ... @@ -269,7 +270,12 @@ A3_1=zeros(n,m);
269 270 A3_2=-aras-arav;
270 271 A3_3=-bras-brav;
271 272  
272   -B3=sigma*ones(n,m).*(tsobs+273.15).^4+cras+crav-ratm*ones(n,m);
  273 +% 3 cases according to the type of data:
  274 +%
  275 +
  276 +Mrad=(1-emis).*ratm+sigma*ones(n,m).*emis.*(input1+273.15).^4;
  277 +
  278 +B3=Mrad+cras+crav-ratm*ones(n,m);
273 279  
274 280 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;
275 281  
... ... @@ -346,7 +352,7 @@ A3_1(les<lesmin)=0;
346 352 A3_2(les<lesmin)=-aras(les<lesmin)-arav(les<lesmin);
347 353 A3_3(les<lesmin)=-bras(les<lesmin)-brav(les<lesmin);
348 354  
349   -B3(les<lesmin)=sigma*(tsobs(les<lesmin)+273.15).^4+cras(les<lesmin)+crav(les<lesmin)-ratm;
  355 +B3(les<lesmin)=Mrad(les<lesmin)+cras(les<lesmin)+crav(les<lesmin)-ratm;
350 356  
351 357 det(les<lesmin)=A1_1(les<lesmin).*A2_2(les<lesmin).*A3_3(les<lesmin)...
352 358 - 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
501 507  
502 508 end % end stability loop
503 509  
504   - trad=((ratm.*ones(n,m)-cras-crav-X2.*(aras+arav)-X3.*(bras+brav))/sigma).^0.25-273.15;
  510 + tsurf=((emis.*ratm.*ones(n,m)-cras-crav-X2.*(aras+arav)-X3.*(bras+brav))/(sigma.*emis)).^0.25-273.15;
505 511  
506 512 rns=(crns+arns.*X2+brns.*X3)./(1-xg);
507 513 rnv=crnv+arnv.*X2+brnv.*X3;
... ...
main_test.m
... ... @@ -48,7 +48,7 @@ rstmin=100; % minimum stomatal resistance (s/m)
48 48 albe=0.25; % surface albedo (BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
49 49 xg=0.4; % soil heat flux to soil net radiation ratio at midday (BEWARE: should vary diurnally / code to be modified in that case)
50 50 za=2.32; % reference height for the climate forcing [m]
51   -emis=0.98; % surface emissivity (BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
  51 +emis=0.99; % surface emissivity (BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
52 52 patm=990; % atmospheric pressure [hPa](BEWARE: should vary diurnally and from day-to-day / code to be modified in that case)
53 53 sigmoy=0.45; % attenuation coef. of incoming radiation within the canopy (0.9 used in the layer case for longwave)
54 54 sigma=5.67e-8; % Stafan-Boltzman constant [W m-2 K-4]
... ... @@ -62,35 +62,35 @@ betastress=0.01;
62 62 % run potential conditions
63 63  
64 64 % series
65   -[tradpot(i),tspot(i),tvpot(i),t0pot(i),rnspot(i),rnvpot(i),gpot(i),hspot(i),hvpot(i),lespot(i),levpot(i)]=...
66   - SPARSE(betapot,betapot,2,1,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
  65 +[tsurfpot(i),tspot(i),tvpot(i),t0pot(i),rnspot(i),rnvpot(i),gpot(i),hspot(i),hvpot(i),lespot(i),levpot(i)]=...
  66 + 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);
67 67  
68 68 % parallel
69 69  
70   -[tradpot2(i),tspot2(i),tvpot2(i),t0pot2(i),rnspot2(i),rnvpot2(i),gpot2(i),hspot2(i),hvpot2(i),lespot2(i),levpot2(i)]=...
71   - SPARSE(betapot,betapot,2,2,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
  70 +[tsurfpot2(i),tspot2(i),tvpot2(i),t0pot2(i),rnspot2(i),rnvpot2(i),gpot2(i),hspot2(i),hvpot2(i),lespot2(i),levpot2(i)]=...
  71 + 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);
72 72  
73 73 % run fully stressed conditions
74 74  
75 75 % series
76   -[tradstress(i),tsstress(i),tvstress(i),t0stress(i),rnsstress(i),rnvstress(i),gstress(i),hsstress(i),hvstress(i),lesstress(i),levstress(i)]=...
77   - SPARSE(betastress,betastress,2,1,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
  76 +[tsurfstress(i),tsstress(i),tvstress(i),t0stress(i),rnsstress(i),rnvstress(i),gstress(i),hsstress(i),hvstress(i),lesstress(i),levstress(i)]=...
  77 + 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);
78 78  
79 79 %parallel
80 80  
81   -[tradstress2(i),tsstress2(i),tvstress2(i),t0stress2(i),rnsstress2(i),rnvstress2(i),gstress2(i),hsstress2(i),hvstress2(i),lesstress2(i),levstress2(i)]...
82   - =SPARSE(betastress,betastress,2,2,ta(i),rh(i),rg(i),ua(i),glai(i),lai(i),zf(i),za,albe,rstmin,xg,sigmoy);
  81 +[tsurfstress2(i),tsstress2(i),tvstress2(i),t0stress2(i),rnsstress2(i),rnvstress2(i),gstress2(i),hsstress2(i),hvstress2(i),lesstress2(i),levstress2(i)]...
  82 + =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);
83 83  
84 84 % run layer/series version
85 85 %tsobs2(i)=((sigma*(tsobs(i)+273.15)^4-(1-emis)*ratmobs(i))/(emis*sigma)).^0.25-273.15;
86 86  
87   -[trad(i),ts(i),tv(i),toto(i),rns(i),rnv(i),g(i),hs(i),hv(i),les(i),lev(i)]=...
88   - 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);
  87 +[tsurf(i),ts(i),tv(i),toto(i),rns(i),rnv(i),g(i),hs(i),hv(i),les(i),lev(i)]=...
  88 + 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);
89 89  
90 90 % run patch/parallel version
91 91  
92   -[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);
93   -%[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));
  92 +[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);
  93 +%[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));
94 94  
95 95 end
96 96  
... ...