%% FvCB or FvCBL models
%% used in Ikawa et al. (2019) Plant Production Sciecne
%% updated on 5/2/2018 - Hiroki Ikawa
% clear
%% Inputs & example
% Rd25=0.2; % heterotrophic respiration at 25 (degreeC)
% Q=1500*0.85; % absorbed photosynthetic photon flux density(umol m-2 s-1)
% ca=590; % CO2 concentration outside of leaf boundary layer(umol mol-1)
% rh=60; % relative humidity (%)
% tl=25; % leaf temperature (degree C)
% vpdl = 1; % vapor pressure deficit with respect to leaf temperature(kPa)
% * vpdl can be calculated from tl and rh
% ta=25; % air temperature (degree C)
% fb=0.47; % fraction of quanta reaching PSII
% gbw=2; % boundary conductance for H2O (molm-2s-1)
% gsmin=0.006; % minimum gsw 0.02 molm-2s-1 Shimono 2010 Table 2 0.015 for ambient and 0.05 for FACE
% vcmax25=100; % apparent maximum carboxylation rate at 25C
% cbwb=15; % coefficient of BWB or BWBL model
% opt_bl=2; % 1 for BBW and 2 for BBWL
% rjv=1.7; % ratio of Jmax25 to Vcmax25
% gamma=43.8; % CO2 compensation point based on Cs (umol mol-1)
% press=1013; % atmospheric pressure (hPa)
%% Outputs
% A: net leaf photosynthesis (umolm-2 s-1)
% Rd: mitochondrial respiration (umolm-2 s-1)
% E: leaf transpiration (mmolm-2 s-1)
% gsw: stomatal conductance for H2O (mol m-2s-1)
% gm: mesophyll conductance for CO2 (mol m-2s-1)
% cc: CO2 in chloroplast (umol mol-1)
% gstar: CO2 compensation point based on Ci (umol mol-1)
% count: # of iteration
% Av: carboxylation limited A (umolm-2 s-1)
% Aj: electron transport limited A (umolm-2 s-1)
function [A,Rd,E,gsw,gm,cc,count,Av,Aj] = farquhar_gm(Rd25,Q,ca,rh,vpdl,tl,ta,press,fb,cbwb,gbw,vcmax25,gsmin,rjv,opt_bl,gamma,gm25,varargin)
if ~isempty(varargin)
A=varargin{1}; % count can be reduced if a close number to A is known.
else
A=0.1; % any number
end
%% initial values (any numbers)
Ao=1; %umolm-2s-1 initial unreal value
E=1; %molm-2s-1
count = 0;
%% Model parameter
EPS=0.00001;
opt_tern=1; % ternary function
c_max=999; % maximum # of iteration
%% Sharkey (2007) T parameterization based on Bernacchi (2002)
c(1)=18.7145; dH(1)=46.39; % Rd
c(2)=26.355; dH(2)=65.33; % Vcmax
c(3)=17.71; dH(3)=43.9; % Jmax
c(4)=35.9774; dH(4)=80.99; %kc
c(5)=12.3772; dH(5)=23.72; %ko
c(6)=11.187; dH(6)=24.46; % gstar
%% Scafaro (2011) temoerature function for rice
c(7)=18.29; dH(7)=45.22; % temp function for gm
b=0.76+0.018.*tl-3.7.*10^(-4).*tl.^2;
phips2=0.352+0.022.*tl-3.4.*10^(-4).*tl.^2; % quantum yield (umol/mol)
co=210000; % oxygen umol mol-1
%% Internal calculation for model input
tlk=tl+273.15; %K
rh=rh/100;
pressh=press.*100; % pressure (Pa)
R=8.314;
Rd=Rd25*exp(c(1)-dH(1)*1000/R/tlk); % mitochondrial respiration in umolm-2s-1
gbc=gbw/(1.6^(2/3)); % boundary conductance for CO2
esta=100*6.1078*exp(17.2694*ta/(ta+237.3)); % saturated vapour pressure at air temperature (Pa)
estl=100*6.1078*exp(17.2694*tl/(tl+237.3)); % saturated vapour pressure at leaf temperature (Pa)
ea=rh*esta; % vapour pressure in air (Pa)
wi=estl/pressh; % mole fraction of water vapour in leaf (mol mol-1)
wa=ea/pressh; % mole fraction of water vapour in air (mol mol-1)
kc=exp(c(4)-dH(4)*1000/R/tlk)/pressh*10^6; % umol mol-1
ko=exp(c(5)-dH(5)*1000/R/tlk)*1000/pressh*10^6;
gstar=exp(c(6)-dH(6)*1000/R/tlk)/pressh*10^6;
gmtfun=exp(c(7)-dH(7)*1000/R/tlk);
gm=gm25*gmtfun;
vcmax=vcmax25*exp(c(2)-dH(2)*1000/R/tlk);
jmax25=vcmax25*rjv;
jmax=jmax25*exp(c(3)-dH(3)*1000/R/tlk);
j=((jmax+fb*phips2*Q)-((jmax+fb*phips2*Q)^2-4*b*(jmax*fb*phips2*Q))^(0.5))/2/b;
if opt_bl==1
hfunction=rh;
elseif opt_bl==2
hfunction=1/(1+vpdl/1);
end
while abs(Ao-A)>EPS
count = count + 1;
if opt_tern==1
cs=ca-A/gbc-ca*E/gbc; %von Caemmerer81 (B10)
else
cs=ca-A/gbc;
end
if cs<gstar+A/gm+A/gbc
cs=gstar+A/gm+A/gbc;
end
gsw=cbwb*A*hfunction/(cs-gamma)+gsmin;
if gsw<gsmin
gsw=gsmin;
end
if opt_bl==2
gsw=gsw*1.6;
end
if opt_tern==1
E=gbw*gsw/(gbw+gsw)*(wi-wa)/(1-(wi+wa)/2); % mol m-2s-1 von Caemmerer81 (B13)
else
E=gbw*gsw/(gbw+gsw)*(wi-wa);
end
gsc=gsw/1.6; % molar diffusivity of CO2/H2O is 1/1.6
rtc=1/gbc+1/gsc;
gtc=1/rtc;
if opt_tern==1
ci=((gtc-E/2)*ca-A)/(gtc+E/2); % intercellular CO2 (umol mol-1) von Caemmerer81 (B18)
else
ci=cs-A/gsc;
end
if ci<gstar+A/gm
ci=gstar+A/gm;
end
%cc_count(count+1)=cc;
% Vj=j*cc/(4*cc+8*gstar);
% Vc=vcmax*cc/(cc+kc*(1+oi/ko));
Ao = A;
% A=min([Vj,Vc])*(1-gstar/cc)-Rd;
% A=0.1*A+(1-0.9)*Ao;
%% Ethier et al. (2006)
prmv_a=-1/gm;
prmv_b=(vcmax-Rd)/gm+ci+kc*(1+co/ko);
prmv_c=Rd*(ci+kc*(1+co/ko))-vcmax*(ci-gstar);
Av=(-prmv_b+sqrt(prmv_b^2-4*prmv_a*prmv_c))/prmv_a/2;
prmj_a=-1/gm;
prmj_b=(j/4-Rd)/gm+ci+2*gstar;
prmj_c=Rd*(ci+2*gstar)-j/4*(ci-gstar);
Aj=(-prmj_b+sqrt(prmj_b^2-4*prmj_a*prmj_c))/prmj_a/2;
A=min([Av,Aj]);
cc=ci-A/gm;
if cc<gstar
cc=gstar;
end
if count> c_max
break
end
end
end