FvCB with gm in Ikawa (2019, PPS)

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