FvCB in Ikawa (2018, GCB)
%% FvCB or FvCBL models with infinite mesophyll conductance
%% used in Ikawa et al. (2018) Global Change Biology
%% updated on 02/21/2017 - Hiroki Ikawa
% clear
%% Inputs & example
% Rd25=0.2; % heterotrophic respiration at 25 (degreeC)
% Q=1500*0.85; % absorbed PPFD (umol m-2s-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; % initial slope of J to Q (umolCO2m-2s-1/umolEm-2s-1)
% 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)
% ci: interceulluar CO2 (umol mol-1)
% gstar: CO2 compensation point based on Ci (umol mol-1)
% count: # of iteration
function [A,Rd,E,gsw,ci,gstar,count,Vj,Vc] = farquhar(Rd25,Q,ca,rh,vpdl,tl,ta,press,fb,cbbw,gbw,vcmax25,gsmin,rjv,opt_bl,gamma,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.1; %molm-2s-1
count = 0;
%% Model parameter
EPS=0.00001; % setting for iteration
opt_tern=1; % ternary function
c_max=10; % maximum # of iteration
%% Bernacchi (2001) parameterization for apparent Vcmax
c(1)=18.72; dH(1)=46.39; % Rd
c(2)=26.35; dH(2)=65.33; % Vcmax
c(3)=17.57; dH(3)=43.54; % Jmax
c(4)=38.05; dH(4)=79.43; % Kc 404.9 at 25 C
c(5)=20.30; dH(5)=36.38; % Ko 278.4 at 25 C
b=0.76+0.018.*tl-3.7.*10^(-4).*tl.^2; % convexity
phips2=0.352+0.022.*tl-3.4.*10^(-4).*tl.^2;
vcor=0.255; % Vomax/Vcmax Hikosaka (1998)
oi=210000; % oxygen umol mol-1 Hikosaka (1998)
%% Internal calculation for model input
tlk=tl+273.15; % K
rh=rh/100;
pressh=press.*100; % pressure (Pa)
R=8.314; % gas constant
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);
ko=exp(c(5)-dH(5)*1000/R./tlk)*1000;
gstar=kc*oi/2/ko*vcor; % 39 at 25 degC
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
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 (ppm) % umol mol-1 von Caemmerer81 (B18)
else
ci=cs-A/gsc;
end
if ci<gstar
ci=gstar;
end
Vj=j*ci/(4*ci+8*gstar);
Vc=vcmax*ci/(ci+kc*(1+oi/ko));
Ao = A;
A=min([Vj,Vc])*(1-gstar/ci)-Rd;
if count> c_max
break
end
end
end