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
% full linearization around air temperature

%License:
%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

% inputs

if typerun==1 % retrieval: 
    % imposed observed radiative temperature and view zenith angle
    vza=input2;
    betas=1;
    betav=1;
else % precribed:
    % imposed relative stress "beta" for the soil "s" and the vegetation "v"
    betas=input1;
    betav=input2;
    vza=0;
end
    
% n: row number
% m: line number

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]
% ua = wind speed [m/s]
% glai = green Leaf Are Index [m2/m2]
% lai = total Leaf Are Index [m2/m2]
% zf = vegetation height [m]
% za = wind speed meas. height [m]
% albe = OBSERVED albedo [-]
% rvvmin = mimimum stomatal resistance [s/m]
% xg = soil to soil net radiation ratio
% sigmoy = Beer Lambert law attenuation coefficient
% vza = view zenith angle [rad]

% 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 [-]

% outputs

% trad = radiative surface temperature (should be equal to tsobs in retrieval mode)
% ts = soil temperature in precribed conditions [K then �C]
% tv = vegetation temperature [K then �C]
% t0 = aerodynamic temperature [K then �C]
% rns = soil net radiation [W/m2]
% rnv = vegetation net radiation [W/m2]
% g = soil heat flux [W/m2]
% hs = soil sensible heat flux [W/m2]
% hv = vegetation sensible heat flux [W/m2]
% les = soil evaporation [W/m2]
% lev = transpiration [W/m2]

% constants

    rcp=1170; % product of air density and mass heat constant of the air
    gamma=0.66; % psychrometric constant
    sigma=5.66e-8; % Setfan Boltzmann constant
	alfo=5e-3; % ?
    xn=2.5;  
    zoms=5e-3; % bare soil roughness length
    wl=0.05; % mean leaf width [m] 
    albv=0.2; % vegetation albedo
    emisv=0.98; % vegetation emissivity
    emiss=0.96; % soil emissivity
        
% initialisation   
            
    X0=5*ones(n,m); % aerodynamic temperature initialization
    error=10; % criteria for convergence of aerodynamic temperature (stability loop)
    k=0; % index for stability loop (stop  100)
    
    lesmin=0;
               
%	computes partial pressure from rh

	ea=0.01*rh*6.1878.*exp(17.269*(ta-273)./(ta-35.86));
    
%	fraction cover

    fracover=(1-exp(-sigmoy*lai./cos(vza))); % fraction cover
    gfracover=1-exp(-sigmoy*glai./cos(vza)); % green fraction cover
    fracoverl=(1-exp(-0.9*lai)); % fraction cover longwave   

%   patch model: green LAI transformed into a clump green LAI:

if typeohm>1 
    glai=glai./gfracover;
end

%	roughness length zom and displacement height d (several formulations)

	d = 0.65*zf;  
	zom = max(0.13*zf,zoms);
    %d= zf*(1-2*(1-exp(-lai/2))/lai); % shaw and pereira ?
    %zom= zf*exp(-lai/2)*(1-exp(-lai/2));
    
%	soil-air aerodynamic conductance

	xkzf=0.4*0.4*ua.*(zf-d)./log((za-d)./zom);
	ras=zf.*exp(xn).*(exp((-xn*zoms)./zf)-exp((-xn*(d+zom))./zf))./(xn*xkzf);

%	leaf-air aerodynamic conductance

	uzf=ua.*log((zf-d)./zom)./log((za-d)./zom);
	%rav=xn*sqrt(wl/uzf)/(4*alfo*lai*(1-exp(-xn/2)));
    
    rav=xn*sqrt(wl./uzf)./(4*alfo*glai*(1-exp(-xn/2)));
        
%	emissivity of the air in clear sky conditions

	emisa = 1.24*(ea./(ta)).^(1/7);
    
%	Radiation components (longwave and total net radiation)

	ratm=emisa.*sigma*ta.^4;
    
    if typeohm==1 % layer approach
   
    [albs,ans,bns,cns,anv,bnv,cnv,cnas,cnav]=calcrn_vect(rg,ratm,emiss,emisv,albe,albv,fracover,fracoverl);
      
	arns=(ans.*sigma*4*ta.^3)*(1-xg);
    aras=ans.*sigma*4*ta.^3;
    brns=(bns.*sigma*4*ta.^3)*(1-xg);
    bras=bns.*sigma*4*ta.^3;
    arnv=anv.*sigma*4*ta.^3;
    arav=arnv;
    brnv=bnv.*sigma*4*ta.^3;
    brav=brnv;
    crns=((ans+bns).*sigma*ta.^4+cns)*(1-xg);
    crnv=(anv+bnv).*sigma*ta.^4+cnv;
    cras=(ans+bns).*sigma*ta.^4+cnas;
    crav=(anv+bnv).*sigma*ta.^4+cnav;
    
    else % patch approach
    
    albs=(albe-fracover*albv)./(1-fracover);
    arns=-(1-fracover).*emiss*(sigma*4*ta.^3)*(1-xg);
    brns=zeros(n,m);
    bras=zeros(n,m);
    arnv=zeros(n,m);
    arav=zeros(n,m);
    brnv=-fracover.*emisv*sigma*4*ta.^3;
    crns=(1-fracover).*((1-albs).*rg+emiss*(ratm-sigma*ta.^4))*(1-xg);
    crnv=fracover.*((1-albv)*rg+emisv*(ratm-sigma*ta.^4));
    aras=-(1-fracover).*emiss*sigma*4*ta.^3;
    brav=brnv;
    cras=(1-fracover).*emiss*(ratm-sigma*ta.^4);
    crav=fracover.*emisv*(ratm-sigma*ta.^4);

    end
    
%	various terms appearing in the linearization

	da=6.1878*exp(17.269*(ta-273)./(ta-35.86))-ea; % vapour pressure deficit VPD
	delta=((da+ea).*4097.9337).*(ta-35.86).^-2; % slope of the saturation curve
    rcpg=rcp/gamma;
	rcpgd=rcpg*delta;
    
%	leaf-air stomatal conductance

    f=(0.0055*max(rg,10)*2)./glai; 
	frg=(1+f)./(f+rvvmin/5000); % stress factor / solar radiation 
    fea=1+0.04*da; % stress factor / VPD
    
%	aerodynamic resistance in neutral conditions
		
   	alg = log((za-d)./zom);
  	ra0 = alg.*alg/(.4*.4*ua);
    
    X0old=X0;

while error>1e-2 && k<100

        k=k+1;
                
%	air-air aerodynamic conductance

%	Richardson number:

   	ri = 5*(za-d).*9.81.*X0./(ua.*ua.*ta);
        if rg<100
        ri=0;
        end
        
    p(X0>0) = 0.75; % unstable conditions
    p(X0<0) = 2; % stable conditions

    ra = ra0./(1+ri).^p;
    ga=1./ra;
        
% layer approach
    
if typeohm==1

% aggregated conductances for the layer network

gav=1./rav;
gvv=betav./(rvvmin*fea*frg./glai+rav);
gas=1./ras;
gss=betas./ras;
g3a=ga+gas+gav;
g3=ga+gss+gvv;

else
    
% and for the patch network

ga2s=(1-fracover)./(ras+ra);
ga2v=fracover./(rav+ra);
g3s=betas.*(1-fracover)./(ras+ra);
g3v=betav.*fracover./(rav+rvvmin*fea*frg./glai+ra);

end

    % RETRIEVAL MODE

if typerun==1    

% MODE 1: potential rate transpiration

% var 1: LEs, var 2: Xs=Ts-Ta, var 3: Xv=Tv-Ta

if typeohm==1 % layer
    
A1_1=ones(n,m);
A1_2=-rcp*gas.*gas./g3a-arns+rcp*gas;
A1_3=-rcp*gas.*gav./g3a-brns;

B1=crns;

A2_1=-gvv./(gvv+ga);
A2_2=-rcp*gav.*gas./g3a-arnv;
A2_3=-rcpgd.*gvv.*gvv./(gvv+ga)-rcp*gav.*gav./g3a+rcpgd.*gvv+rcp.*gav-brnv;

B2=crnv-rcpg.*gvv.*ga.*da./(gvv+ga);

else % patch
    
A1_1=ones(n,m);
A1_2=-arns+rcp*ga2s;
A1_3=zeros(n,m);

B1=crns;

A2_1=zeros(n,m);
A2_2=zeros(n,m);
A2_3=rcpgd.*g3v+rcp*ga2v-brnv;

B2=crnv-rcpg.*g3v.*da;
    
end

% same equation for both approaches for the link between radiative and
% source surface temperatures

A3_1=zeros(n,m);
A3_2=-aras-arav;
A3_3=-bras-brav;

% 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;

IA1_1=(A2_2.*A3_3 - A2_3.*A3_2)./det;
IA1_2=-(A1_2.*A3_3 - A1_3.*A3_2)./det;
IA1_3=(A1_2.*A2_3 - A1_3.*A2_2)./det;

IA2_1=-(A2_1.*A3_3 - A2_3.*A3_1)./det;
IA2_2=(A1_1.*A3_3 - A1_3.*A3_1)./det;
IA2_3=-(A1_1.*A2_3 - A1_3.*A2_1)./det;

IA3_1=(A2_1.*A3_2 - A2_2.*A3_1)./det;
IA3_2=-(A1_1.*A3_2 - A1_2.*A3_1)./det;
IA3_3=(A1_1.*A2_2 - A1_2.*A2_1)./det;

X1=IA1_1.*B1+IA1_2.*B2+IA1_3.*B3;
X2=IA2_1.*B1+IA2_2.*B2+IA2_3.*B3;
X3=IA3_1.*B1+IA3_2.*B2+IA3_3.*B3;

% outputs
 
    les=X1;   
    
if typeohm==1 % layer
    
    X0=(gas.*X2+gav.*X3)./g3a;
	d0=(rcpg.*ga.*da-X1-rcpgd.*gvv.*X3)./(rcpg.*(gvv+ga));
	lev=rcpg.*gvv.*(d0+delta.*X3);
    
else % patch
    
    X0=(ga2s.*X2+ga2v.*X3)./ga;
	d0=da-X1./(ga.*rcpg)-g3v.*(da+delta.*X3)./ga;
	lev=rcpg.*g3v.*(da+delta.*X3);
    
end
    
% MODE 2: case LES<lesmin: var 1 lev; var 2 Xs=ts-ta; var 3 Xv=tv-ta

if typeohm==1 % layer

A1_1(les<lesmin)=0;
A1_2(les<lesmin)=-rcp*gas(les<lesmin).*gas(les<lesmin)./g3a(les<lesmin)-arns(les<lesmin)+rcp*gas(les<lesmin);
A1_3(les<lesmin)=-rcp*gas(les<lesmin).*gav(les<lesmin)./g3a(les<lesmin)-brns(les<lesmin);

B1(les<lesmin)=crns(les<lesmin)-lesmin;

A2_1(les<lesmin)=1;
A2_2(les<lesmin)=-rcp*gav(les<lesmin).*gas(les<lesmin)./g3a(les<lesmin)-arnv(les<lesmin);
A2_3(les<lesmin)=-rcp*gav(les<lesmin).*gav(les<lesmin)./g3a(les<lesmin)+rcp*gav(les<lesmin)-brnv(les<lesmin);

B2(les<lesmin)=crnv(les<lesmin);

else % patch
   
A1_1(les<lesmin)=0;
A1_2(les<lesmin)=rcp*ga2s(les<lesmin)-arns(les<lesmin);
A1_3(les<lesmin)=0;

B1(les<lesmin)=crns(les<lesmin)-lesmin;

A2_1(les<lesmin)=1;
A2_2(les<lesmin)=0;
A2_3(les<lesmin)=rcp*ga2v(les<lesmin)-brnv(les<lesmin);

B2(les<lesmin)=crnv(les<lesmin);
    
end

% same equation for both approaches for the link between radiative and
% source surface temperatures

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)=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)...
    .*A3_1(les<lesmin) + A1_3(les<lesmin).*A2_1(les<lesmin).*A3_2(les<lesmin) - A1_3(les<lesmin).*A2_2(les<lesmin).*A3_1(les<lesmin);

IA1_1(les<lesmin)=(A2_2(les<lesmin).*A3_3(les<lesmin) - A2_3(les<lesmin).*A3_2(les<lesmin))./det(les<lesmin);
IA1_2(les<lesmin)=-(A1_2(les<lesmin).*A3_3(les<lesmin) - A1_3(les<lesmin).*A3_2(les<lesmin))./det(les<lesmin);
IA1_3(les<lesmin)=(A1_2(les<lesmin).*A2_3(les<lesmin) - A1_3(les<lesmin).*A2_2(les<lesmin))./det(les<lesmin);

IA2_1(les<lesmin)=-(A2_1(les<lesmin).*A3_3(les<lesmin) - A2_3(les<lesmin).*A3_1(les<lesmin))./det(les<lesmin);
IA2_2(les<lesmin)=(A1_1(les<lesmin).*A3_3(les<lesmin) - A1_3(les<lesmin).*A3_1(les<lesmin))./det(les<lesmin);
IA2_3(les<lesmin)=-(A1_1(les<lesmin).*A2_3(les<lesmin) - A1_3(les<lesmin).*A2_1(les<lesmin))./det(les<lesmin);

IA3_1(les<lesmin)=(A2_1(les<lesmin).*A3_2(les<lesmin) - A2_2(les<lesmin).*A3_1(les<lesmin))./det(les<lesmin);
IA3_2(les<lesmin)=-(A1_1(les<lesmin).*A3_2(les<lesmin) - A1_2(les<lesmin).*A3_1(les<lesmin))./det(les<lesmin);
IA3_3(les<lesmin)=(A1_1(les<lesmin).*A2_2(les<lesmin) - A1_2(les<lesmin).*A2_1(les<lesmin))./det(les<lesmin);

X1(les<lesmin)=IA1_1(les<lesmin).*B1(les<lesmin)+IA1_2(les<lesmin).*B2(les<lesmin)+IA1_3(les<lesmin).*B3(les<lesmin);
X2(les<lesmin)=IA2_1(les<lesmin).*B1(les<lesmin)+IA2_2(les<lesmin).*B2(les<lesmin)+IA2_3(les<lesmin).*B3(les<lesmin);
X3(les<lesmin)=IA3_1(les<lesmin).*B1(les<lesmin)+IA3_2(les<lesmin).*B2(les<lesmin)+IA3_3(les<lesmin).*B3(les<lesmin);

lev(les<lesmin)=real(X1(les<lesmin));
les(les<lesmin)=lesmin;

% outputs

d0(les<lesmin)=da-(lesmin+X1(les<lesmin))./(rcpg.*ga(les<lesmin));

if typeohm==1 % layer
    
    X0(les<lesmin)=(gas(les<lesmin).*X2(les<lesmin)+gav(les<lesmin).*X3(les<lesmin))./g3a(les<lesmin);
                   
else % patch
    
    X0(les<lesmin)=(ga2s(les<lesmin).*X2(les<lesmin)+ga2v(les<lesmin).*X3(les<lesmin))./ga(les<lesmin);
               
end

% end MODE 2

% MODE 3: LEv<0

if typeohm==1 % layer

		A1_1(lev<0)=-arns(lev<0)+rcp*gas(lev<0)-rcp*gas(lev<0).*gas(lev<0)./g3a(lev<0);
		A1_2(lev<0)=-brns(lev<0)-rcp*gas(lev<0).*gav(lev<0)./g3a(lev<0);

        B1(lev<0)=crns(lev<0)-lesmin(lev<0);

		A2_1(lev<0)=-arnv(lev<0)-rcp*gav(lev<0).*gas(lev<0)./g3a(lev<0);
		A2_2(lev<0)=-brnv(lev<0)+rcp*gav(lev<0)-rcp*gav(lev<0).*gav(lev<0)./g3a(lev<0);
		B2(lev<0)=crnv(lev<0);
   
        det(lev<0)=A1_1(lev<0).*A2_2(lev<0)-A1_2(lev<0).*A2_1(lev<0);
        
        X2(lev<0)=(1./det(lev<0)).*(A2_2(lev<0).*B1(lev<0)-A1_2(lev<0).*B2(lev<0));
        X3(lev<0)=(1./det(lev<0)).*(-A2_1(lev<0).*B1(lev<0)+A1_1(lev<0).*B2(lev<0));
           
        X0(lev<0)=(gas(lev<0)./g3a(lev<0)).*X2(lev<0)+(gav(lev<0)./g3a(lev<0)).*X3(lev<0);
        d0(lev<0)=(ga(lev<0).*da-gss(lev<0).*delta.*X2(lev<0)-gvv(lev<0).*delta.*X3(lev<0))./g3(lev<0);

else % patch
    
    X2(lev<0)=(crns(lev<0)-lesmin)./(-arns(lev<0)+rcp*ga2s(lev<0));
    X3(lev<0)=crnv(lev<0)./(-brnv(lev<0)+rcp*ga2v(lev<0));
    X0(lev<0)=ga2s(lev<0).*X2(lev<0)./ga(lev<0)+ga2v(lev<0).*X3(lev<0)./ga(lev<0);
	d0(lev<0)=da-lesmin./(ga(lev<0).*rcpg);
    
end

les(lev<0)=lesmin;
lev(lev<0)=0;

else
    
    % PRESCRIBED MODE
    
    if typeohm==1 % layer
        
		A1_1=-arns+rcp.*gas+rcpgd.*gss-rcp.*gas.*gas./g3a;
			A1_1=A1_1-rcpgd.*gss.*gss./g3;
		A1_2=-brns-rcp.*gas.*gav./g3a;
			A1_2=A1_2-rcpgd.*gss.*gvv./g3;
		B1=crns-rcpg.*gss.*ga.*da./g3;

		A2_1=-arnv-rcp.*gav.*gas./g3a;
			A2_1=A2_1-rcpgd.*gvv.*gas./g3;
		A2_2=-brnv+rcp.*gav+rcpgd.*gvv-rcp.*gav.*gav./g3a;
			A2_2=A2_2-rcpgd.*gvv.*gvv./g3;
		B2=crnv-rcpg.*gvv.*ga.*da./g3;
   
        det=A1_1.*A2_2-A1_2.*A2_1;
        
        X2=(1./det).*(A2_2.*B1-A1_2.*B2);
        X3=(1./det).*(-A2_1.*B1+A1_1.*B2);
           
        X0=(gas./g3a).*X2+(gav./g3a).*X3;
        d0=(ga.*da-gss.*delta.*X2-gvv.*delta.*X3)./g3;
        
        les=real(rcpg.*gss.*(d0+delta.*X2));
        lev=real(rcpg.*gvv.*(d0+delta.*X3));
        
    else % patch
        
        A1_1=-arns+rcp*ga2s+rcpgd.*g3s;
        A1_2=zeros(n,m);
        B1=crns-rcpg.*g3s.*da;
        
        A2_1=zeros(n,m);
        A2_2=-brnv+rcp*ga2v+rcpgd.*g3v;
        B2=crnv-rcpg.*g3v.*da;
        
        det=A1_1.*A2_2-A1_2.*A2_1;
        
        X2=(1./det).*(A2_2.*B1-A1_2.*B2);
        X3=(1./det).*(-A2_1.*B1+A1_1.*B2);
           
        X0=(ga2s.*X2+ga2v.*X3)./ga;
        %d0=da-(g3s.*(da+delta*X2)+g3v.*(da+delta*X3))./ga;
        
        les=real(rcpg.*g3s.*(da+delta.*X2));
        lev=real(rcpg.*g3v.*(da+delta.*X3));
                
    end
            
end     

% outputs whose formulation is the same for prescribed and retrieval modes
    
    if typeohm==1
        
        hs=real(rcp*gas.*(X2-X0));
        hv=real(rcp*gav.*(X3-X0));
        
    else
        
        hs=real(rcp*ga2s.*X2);
        hv=real(rcp*ga2v.*X3);

    end
        
if k==100
     'STOP k=100';
end
          
    error=max(max(abs(X0-X0old)));    
    X0old=X0;
    
         if rg<20
            error=1e-3;
        end
        
end % end stability loop 
 
    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;

    t0=X0+ta;
        
    ts=X2+ta;
	tv=X3+ta;
       
    g=xg*rns;
                 
    ts=ts-273.15;
    tv=tv-273.15;
    t0=t0-273.15;

end