R Markdown
# detach("package:rpmodel", unload = TRUE)
source("R/rpmodel.R")
source("R/subroutines.R")
out_pmodel <- rpmodel(
tc = 20, # temperature, deg C
vpd = 1000, # Pa,
co2 = 400, # ppm,
fapar = 1, # fraction ,
ppfd = 30, # mol/m2/d,
elv = 0, # m.a.s.l.,
kphio = 0.049977, # quantum yield efficiency as calibrated for setup ORG by Stocker et al. 2020 GMD,
beta = 146, # unit cost ratio a/b,
c4 = FALSE,
method_jmaxlim = "wang17",
do_ftemp_kphio = FALSE, # corresponding to setup ORG
do_soilmstress = FALSE, # corresponding to setup ORG
verbose = TRUE
)
## Warning in rpmodel(tc = 20, vpd = 1000, co2 = 400, fapar = 1, ppfd = 30, :
## Atmospheric pressure (patm) not provided. Calculating it as a function of
## elevation (elv), assuming standard atmosphere (101325 Pa at sea level).
## re-implemented without function calls
tc <- 20 # temperature, deg C
vpd <- 1000 # Pa,
co2 <- 400 # ppm,
fapar <- 1 # fraction ,
ppfd <- 30 # mol/m2/d,
elv <- 0 # m.a.s.l.,
kphio <- 0.049977 # quantum yield efficiency as calibrated for setup ORG by Stocker et al. 2020 GMD,
beta <- 146 # unit cost ratio a/b,
## Parameters
c_molmass <- 12.0107 # molecular mass of carbon (g)
rd_to_vcmax <- 0.015 # Ratio of Rdark to Vcmax25, number from Atkin et al., 2015 for C3 herbaceous
kPo <- 101325.0 # standard atmosphere, Pa (Allen, 1973)
kTo <- 25.0 # base temperature, deg C (Prentice, unpublished)
## standard atmospheric pressure
patm <- calc_patm(elv)
#---- Photosynthesis parameters depending on temperature, pressure, and CO2. -
## ambient CO2 partial pression (Pa)
ca <- co2_to_ca( co2, patm )
## photorespiratory compensation point - Gamma-star (Pa)
gammastar <- calc_gammastar( tc, patm )
## Michaelis-Menten coef. (Pa)
kmm <- calc_kmm( tc, patm ) ## XXX Todo: replace 'NA' here with 'patm'
## viscosity correction factor = viscosity( temp, press )/viscosity( 25 degC, 1013.25 Pa)
ns <- viscosity_h2o( tc, patm ) # Pa s
ns25 <- viscosity_h2o( kTo, kPo ) # Pa s
ns_star <- ns / ns25 # (unitless)
## The heart of the P-model: calculate ci:ca ratio (chi) and additional terms
xi <- sqrt( (beta * ( kmm + gammastar ) ) / ( 1.6 * ns_star ) )
chi <- gammastar / ca + ( 1.0 - gammastar / ca ) * xi / ( xi + sqrt(vpd) )
## alternative variables
gamma <- gammastar / ca
kappa <- kmm / ca
## use chi for calculating mj
mj <- (chi - gamma) / (chi + 2.0 * gamma)
## mc
mc <- (chi - gamma) / (chi + kappa)
## leaf-internal CO2 partial pressure (Pa)
ci <- chi * ca
kc <- 0.41 # Jmax cost coefficient
## original code
mprime <- sqrt(mj^2 - kc^(2.0/3.0) * (mj^(4.0/3.0)))
## Following eq. 17 in Stocker et al., 2020 GMD - EQUIVALENT TO ABOVE
mprime_alt <- mj * sqrt(1.0 - (kc / mj)^(2.0/3.0))
lue = kphio * mprime * c_molmass
vcmax_unitiabs = kphio * mprime / mc
iabs <- fapar * ppfd
# Gross Primary Productivity
gpp <- iabs * lue # in g C m-2 s-1
# Vcmax per unit ground area is the product of the intrinsic quantum
# efficiency, the absorbed PAR, and 'n'
vcmax <- iabs * vcmax_unitiabs
# Jmax using again A_J = A_C, derive the "Jmax limitation factor"
fact_jmaxlim <- vcmax * (ci + 2.0 * gammastar) / (kphio * iabs * (ci + kmm))
# use definition of Jmax limitation factor (L in Eq. 13) and solve for Jmax.
jmax <- 4.0 * kphio * iabs / sqrt( (1.0 / fact_jmaxlim)^2 - 1.0 )
## Alternatively, Jmax can be calculated from Eq. F10 in Stocker et al., 2020 - SAME RESULT!
kc <- 0.41
jmax_alt <- 4.0 * kphio * iabs * sqrt((mj / kc)^(2/3) - 1.0)
fact_jmaxlim_alt <- 1.0 / sqrt(1 + (4.0 * kphio * iabs / jmax_alt)^2)
## check if both m-prime formulations are equivalent
print(all.equal(mprime, mprime_alt, tol = 0.001))
## [1] TRUE
## check if both Jmax limitation factor formulations are equivalent
print(all.equal(fact_jmaxlim, fact_jmaxlim_alt, tol = 0.001))
## [1] TRUE
## check if both Jmax formulations are equivalent
print(all.equal(jmax_alt, jmax, tol = 0.001))
## [1] TRUE
## rpmodel identical output?
print(all.equal(out_pmodel$jmax, jmax, tol = 0.001))
## [1] TRUE