P Model

The R script below corresponds to the Pmodel version used in Wang et al. Nature Plants (2017) paper. The script includes four sections: (1) specifying the values of two parameters in Pmodel; (2) defining the used functions such as calculating air pressure in the unit of Pascal as a function of elevation; (3) reading the input data of temperature, elevation, photosynthetic photon flux density, vapor pressure deficit, fractional absorbed photosynthetic active radiation and atmospheric CO2 concentration. Here using the standard values as example; (4) the Pmodel code for estimating GPP from inputs data
The R software (https://www.r-project.org) is required to run the script.

##===================================
## 1. parameterization
##===================================

# the ratio of cost factor b to a at reference temperature
beta <- 240

# the cost factor of maintaining Jmax
c <- 0.41

##===================================
## 2. define functions
##===================================
{
  # calculate air pressure in Pa
  calc_patm <- function( elv ){
    #-----------------------------------------------------------------------
    # Input:    - elevation, m (elv)
    # Output:   - float, atmospheric pressure at elevation 'elv', Pa (patm)
    # Features: Returns the atmospheric pressure as a function of elevation
    #           and standard atmosphere (1013.25 hPa)
    # Depends:  - connect_sql
    #           - flux_to_grid
    #           - get_data_point
    #           - get_msvidx
    # Ref:      Allen et al. (1998)
    #-----------------------------------------------------------------------
    
    # Define constants:
    kPo <- 101325   # standard atmosphere, Pa (Allen, 1973)
    kTo <- 298.15   # base temperature, K (Prentice, unpublished)
    kL <- 0.0065    # temperature lapse rate, K/m (Allen, 1973)
    kG <- 9.80665   # gravitational acceleration, m/s^2 (Allen, 1973)
    kR <- 8.3143    # universal gas constant, J/mol/K (Allen, 1973)
    kMa <- 0.028963 # molecular weight of dry air, kg/mol (Tsilingiris, 2008)
    
    # Convert elevation to pressure, Pa:
    patm <- kPo*(1.0 - kL*elv/kTo)**(kG*kMa/(kR*kL))
    
    return (patm)
  }
  
  # calculate K (MM coefficient of Rubisco) in Pa
  calc_k <- function(temp, patm) {
    #-----------------------------------------------------------------------
    # Input:    - float, air temperature, deg C (temp)
    #           - float, atmospheric pressure, Pa (patm)
    # Output:   float, Pa (mmk)
    # Features: Returns the temperature & pressure dependent Michaelis-Menten
    #           coefficient, K (Pa).
    # Ref:      Bernacchi et al. (2001), Improved temperature response 
    #           functions for models of Rubisco-limited photosynthesis, 
    #           Plant, Cell and Environment, 24, 253--259.
    #-----------------------------------------------------------------------
    
    # Define constants
    kc25 <- 39.97      # Pa, assuming 25 deg C & 98.716 kPa
    ko25 <- 2.748e4    # Pa, assuming 25 deg C & 98.716 kPa
    dhac <- 79430      # J/mol
    dhao <- 36380      # J/mol
    kR   <- 8.3145     # J/mol/K
    kco  <- 2.09476e5  # ppm, US Standard Atmosphere
    
    vc <- kc25*exp(dhac*(temp - 25.0)/(298.15*kR*(temp + 273.15)))
    vo <- ko25*exp(dhao*(temp - 25.0)/(298.15*kR*(temp + 273.15)))
    k  <- vc*(1 + kco*(1e-6)*patm/vo)
    
    return(k)
    
  }
  
  # calculate Gstar (CO2 compensation point) in Pa
  calc_gstar_gepisat <- function( temp ) {
    #-----------------------------------------------------------------------
    # Input:    float, air temperature, degrees C (tc)
    # Output:   float, gamma-star, Pa (gs)
    # Features: Returns the temperature-dependent photorespiratory 
    #           compensation point, Gamma star (Pascals), based on constants 
    #           derived from Bernacchi et al. (2001) study.
    # Ref:      Bernacchi et al. (2001), Improved temperature response 
    #           functions for models of Rubisco-limited photosynthesis, 
    #           Plant, Cell and Environment, 24, 253--259.
    #-----------------------------------------------------------------------
    
    # Define constants
    gs25 <- 4.220    # Pa, assuming 25 deg C & 98.716 kPa)
    dha  <- 37830    # J/mol
    kR   <- 8.3145   # J/mol/K
    
    gs <- gs25 * exp( dha * ( temp - 25.0 ) / ( 298.15 * kR * ( temp + 273.15 ) ) )
    
    return( gs )
    
  }
  
  # conver CO2 from ppm to Pa
  co2_to_ca <- function( co2, patm ){
    #-----------------------------------------------------------------------
    # Input:    - float, annual atm. CO2, ppm (co2)
    #           - float, monthly atm. pressure, Pa (patm)
    # Output:   - ca in units of Pa
    # Features: Converts ca (ambient CO2) from ppm to Pa.
    #-----------------------------------------------------------------------
    ca   <- ( 1.e-6 ) * co2 * patm         # Pa, atms. CO2
    return( ca )
  }
}

##===================================
## 3. read monthly input data 
## (replace this with your input data)
## daily/weekly/monthly/annual
##===================================

# example of driving forces at the standard conditions
Tc.deg_C  <- 25 # temperature in degreeC
elv.m  <- 0 # elevation in meter
PPFD.mol_m2  <- 1000 # PAR in mol photon m-2 month-1
VPD.kPa  <- 1 # vapour pressure deficit in kPa
fAPAR  <- 1 # factional absored PAR
CO2.ppm <- 400 # ppm

##==========================================================
## 4. running Pmodel at daily/weekly/monthly/annual step
##==========================================================

# instrinsic quantum yield, based on Cozettle et al.1998, unit: g C/ mol photon
phi0 <- 0.085 # mol C/ mol photon

maxQE <- phi0*12 

# light use efficiency
K <- calc_k(Tc.deg_C,calc_patm(elv.m))
Gstar <- calc_gstar_gepisat(Tc.deg_C)
f1 <- exp(-0.0227*(Tc.deg_C-25)) # the viscosity of water relative to its value at 25˚C
pa <- co2_to_ca(CO2.ppm,calc_patm(elv.m))

m <- (pa - Gstar)/(pa + 2*Gstar + 3*Gstar*sqrt(1.6*VPD.kPa*1000*f1/(K + Gstar)/(beta)))
M <- m*sqrt(1-(c/m)^(2/3))

LUE <- M*maxQE

Iabs <- fAPAR*PPFD.mol_m2
GPP <- LUE*Iabs

Contributed by Han Wang

Close