diff --git a/README.md b/README.md index 4a3be4a..604923f 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,8 @@ # BEPS.jl -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://CUG-hydro.github.io/BEPS.jl/dev) -[![CI](https://github.com/CUG-hydro/BEPS.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/CUG-hydro/BEPS.jl/actions/workflows/CI.yml) -[![Codecov](https://codecov.io/gh/CUG-hydro/BEPS.jl/branch/master/graph/badge.svg)](https://app.codecov.io/gh/CUG-hydro/BEPS.jl/tree/master) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://jl-pkgs.github.io/BEPS.jl/dev) +[![CI](https://github.com/jl-pkgs/BEPS.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/jl-pkgs/BEPS.jl/actions/workflows/CI.yml) +[![Codecov](https://codecov.io/gh/jl-pkgs/BEPS.jl/branch/master/graph/badge.svg)](https://app.codecov.io/gh/jl-pkgs/BEPS.jl/tree/master) Boreal Ecosystem Productivity Simulator in Julia @@ -22,16 +22,16 @@ Boreal Ecosystem Productivity Simulator in Julia - For developers ```bash - git clone https://github.com/CUG-hydro/BEPS.jl + git clone https://github.com/jl-pkgs/BEPS.jl cd BEPS.jl/deps - git clone https://github.com/CUG-hydro/BEPS.c + git clone https://github.com/jl-pkgs/BEPS.c ``` - For users ```bash # In Julia - ] add https://github.com/CUG-hydro/BEPS.jl + ] add https://github.com/jl-pkgs/BEPS.jl ``` ## Usage diff --git a/src/photosynthesis.jl b/src/photosynthesis.jl index 394ceeb..862e717 100644 --- a/src/photosynthesis.jl +++ b/src/photosynthesis.jl @@ -32,138 +32,96 @@ function photosynthesis(Tc_old::Leaf, R::Leaf, Ci_old::Leaf, leleaf::Leaf, end # 这里依赖的变量过多，不易核对 +""" +# Intermediate variables +- `g_lb_c` : leaf laminar boundary layer condunctance to CO2 (mol m-2 s-1) +- `RH_leaf` : relative humidity at leaf surface (0-1) +- `gs_co2_mole`: stomatal conductance to CO2 (mol m-2 s-1) +- `gs_h2o_mole`: stomatal conductance to h2o (mol m-2 s-1) +- `cs` : CO2 concentration at leaf surface (ppm) +- `b_co2` : the intercept term in BWB model (mol CO2 m-2 s-1): b_h2o/1.6 +- `m_co2` : the slope in BWB model: m_h2o/1.6 +- `Γ` : CO2 compensation point (ppm) +- `jmopt` : the maximum potential electron transport rate at 25 deg C (umol m-2 s-1) +- `Jmax` : the maximum potential electron transport rate (umol m-2 s-1) +- `Vcmax` : the maximum velocities of carboxylation of Rubisco (umol m-2 s-1) +- `km_co2` : Michaelis-Menten constant for CO2 (µmol mol-1) +- `km_o2` : Michaelis-Menten constant for O2 (mmol mol-1) +- `tau` : the specifity of Rubisco for CO2 compared with O2 +- `Jₓ` : the flux of electrons through the thylakoid membrane (umol m-2 s-1) +""" @fastmath function photosynthesis_jl(temp_leaf_p::Cdouble, Rsn_leaf::Cdouble, e_air::Cdouble, g_lb_w::Cdouble, vc_opt::Cdouble, - f_soilwater::Cdouble, b_h2o::Cdouble, m_h2o::Cdouble, + β_soil::Cdouble, b_h2o::Cdouble, m_h2o::Cdouble, cii::Cdouble, T_leaf::Cdouble, LH_leaf::Cdouble) - Gs_w::Cdouble = 0.0 - An::Cdouble = 0.0 - ci::Cdouble = 0.0 - - g_lb_c::Cdouble = 0.0 # leaf laminar boundary layer condunctance to CO2 (mol m-2 s-1) - RH_leaf::Cdouble = 0.0 # relative humidity at leaf surface (0-1) - temp_leaf_K::Cdouble = 0.0 # leaf temperature (K) - gs_co2_mole::Cdouble = 0.0 # stomatal conductance to CO2 (mol m-2 s-1) - gs_h2o_mole::Cdouble = 0.0 # stomatal conductance to h2o (mol m-2 s-1) - K::Cdouble = 0.0 # temporary variable - cs::Cdouble = 0.0 # CO2 concentration at leaf surface (ppm) - b_co2::Cdouble = 0.0 # the intercept term in BWB model (mol CO2 m-2 s-1): b_h2o/1.6 - m_co2::Cdouble = 0.0 # the slope in BWB model: m_h2o/1.6 - Γ::Cdouble = 0.0 # CO2 compensation point (ppm) - jmopt::Cdouble = 0.0 # the maximum potential electron transport rate at 25 deg C (umol m-2 s-1) - Jmax::Cdouble = 0.0 # the maximum potential electron transport rate (umol m-2 s-1) - Vcmax::Cdouble = 0.0 # the maximum velocities of carboxylation of Rubisco (umol m-2 s-1) - km_co2::Cdouble = 0.0 # Michaelis-Menten constant for CO2 (µmol mol-1) - km_o2::Cdouble = 0.0 # Michaelis-Menten constant for O2 (mmol mol-1) - tau::Cdouble = 0.0 # the specifity of Rubisco for CO2 compared with O2 - Rd::Cdouble = 0.0 # leaf dark respiration (umol m-2 s-1) - resp_ld25::Cdouble = 0.0 # leaf dark respiration at 25 deg C (umol m-2 s-1) - Jₓ::Cdouble = 0.0 # the flux of electrons through the thylakoid membrane (umol m-2 s-1) - α_ps::Cdouble = 0.0 - β_ps::Cdouble = 0.0 - γ_ps::Cdouble = 0.0 - θ_ps::Cdouble = 0.0 - denom::Cdouble = 0.0 - p_cubic::Cdouble = 0.0 - q_cubic::Cdouble = 0.0 - r_cubic::Cdouble = 0.0 - Qroot::Cdouble = 0.0 - Rroot::Cdouble = 0.0 - - # double root1, root2, root3; - # double root_min = 0, root_max = 0, root_mid = 0; - ψ::Cdouble = 0.0 - j_sucrose::Cdouble = 0.0 # net photosynthesis rate limited by sucrose synthesis (umol m-2 s-1) - # double wc, wj, psguess; # gross photosynthesis rate limited by light (umol m-2 s-1) - # double Aquad, Bquad, Cquad; - # double b_ps, a_ps, e_ps, d_ps; - product::Cdouble = 0.0 - ps_1::Cdouble = 0.0 - delta_1::Cdouble = 0.0 - r3q::Cdouble = 0.0 - tprime25::Cdouble = 0.0 - ca::Cdouble = CO2_air # atmospheric co2 concentration (ppm) PPFD::Cdouble = 4.55 * 0.5 * Rsn_leaf # incident photosynthetic photon flux density (PPFD) umol m-2 s-1 - if (2 * PPFD < 1) - PPFD = 0.0 - end - - temp_leaf_K = T_leaf + 273.13 + (2PPFD < 1) && (PPFD = 0.0) + T_leaf_K = T_leaf + 273.13 fact_latent = LAMBDA(temp_leaf_p) bound_vapor = 1.0 / g_lb_w # air_pres::Cdouble = 101.325 # air pressure (kPa) # g_lb_c = (g_lb_w/1.6)*air_pres/(temp_leaf_K*rugc); # (mol m-2 s-1) - press_bars = 1.013 - pstat273 = 0.022624 / (273.16 * press_bars) + atm = 1.013 # 1 atm, [1013.25 hPa] -> 1.013 bar + pstat273 = 0.022624 / (273.16 * atm) T_Kelvin = T_leaf + 273.13 rhova_g = e_air * 2165 / T_Kelvin # absolute humidity, g m-3 rhova_kg = rhova_g / 1000.0 # absolute humidity, kg m-3 - g_lb_c = 1.0 / (1.0 / g_lb_w * 1.6 * temp_leaf_K * (pstat273)) + g_lb_c = 1.0 / (1.0 / g_lb_w * 1.6 * T_leaf_K * (pstat273)) - m_co2 = m_h2o / 1.6 + m_co2 = m_h2o / 1.6 # b_co2 = b_h2o / 1.6 - RH_leaf = SFC_VPD(temp_leaf_K, LH_leaf, fact_latent, bound_vapor, rhova_kg) - tprime25 = temp_leaf_K - tk_25 # temperature difference + RH_leaf = SFC_VPD(T_leaf_K, LH_leaf, fact_latent, bound_vapor, rhova_kg) + tprime25 = T_leaf_K - tk_25 # temperature difference #/***** Use Arrhenius Eq. to compute KC and km_o2 *****/ - km_co2 = TEMP_FUNC(kc25, ekc, tprime25, tk_25, temp_leaf_K) - km_o2 = TEMP_FUNC(ko25, eko, tprime25, tk_25, temp_leaf_K) - tau = TEMP_FUNC(tau25, ektau, tprime25, tk_25, temp_leaf_K) + km_co2 = TEMP_FUNC(kc25, ekc, tprime25, tk_25, T_leaf_K) + km_o2 = TEMP_FUNC(ko25, eko, tprime25, tk_25, T_leaf_K) + tau = TEMP_FUNC(tau25, ektau, tprime25, tk_25, T_leaf_K) K = km_co2 * (1.0 + o2 / km_o2) Γ = 0.5 * o2 / tau * 1000.0 # umol mol-1 - resp_ld25 = vc_opt * 0.004657 + Rd25 = vc_opt * 0.004657 #leaf dark respiration (umol m-2 s-1) # Bin Chen: check this later. reduce respiration by 40% in light according to Amthor - if (2.0 * PPFD > 10) - resp_ld25 *= 0.4 - end - Rd = TEMP_FUNC(resp_ld25, erd, tprime25, tk_25, temp_leaf_K) + (2PPFD > 10) && (Rd25 *= 0.4) + Rd = TEMP_FUNC(Rd25, erd, tprime25, tk_25, T_leaf_K) # jmopt = 29.1 + 1.64*vc_opt; Chen 1999, Eq. 7 jmopt = 2.39 * vc_opt - 14.2 - Jmax = TBOLTZ(jmopt, ejm, toptjm, temp_leaf_K) # Apply temperature correction to JMAX - Vcmax = TBOLTZ(vc_opt, evc, toptvc, temp_leaf_K) # Apply temperature correction to vcmax + Jmax = TBOLTZ(jmopt, ejm, toptjm, T_leaf_K) # Apply temperature correction to JMAX + Vcmax = TBOLTZ(vc_opt, evc, toptvc, T_leaf_K) # Apply temperature correction to vcmax # /* # * APHOTO = PG - resp_ld, net photosynthesis is the difference # * between gross photosynthesis and dark respiration. Note # * photorespiration is already factored into PG. - # * ********************************************************** - # * + # # * Gs from Ball-Berry is for water vapor. It must be divided # * by the ratio of the molecular diffusivities to be valid for A # */ - α_ps = 1.0 + (b_co2 / g_lb_c) - m_co2 * RH_leaf * f_soilwater - β_ps = ca * (g_lb_c * m_co2 * RH_leaf * f_soilwater - 2.0 * b_co2 - g_lb_c) - γ_ps = ca * ca * g_lb_c * b_co2 - θ_ps = g_lb_c * m_co2 * RH_leaf * f_soilwater - b_co2 + α = 1.0 + (b_co2 / g_lb_c) - m_co2 * RH_leaf * β_soil + β = ca * (g_lb_c * m_co2 * RH_leaf * β_soil - 2.0 * b_co2 - g_lb_c) + γ = ca * ca * g_lb_c * b_co2 + θ = g_lb_c * m_co2 * RH_leaf * β_soil - b_co2 - # /* # * Test for the minimum of Wc and Wj. Both have the form: - # * - # * W = (a ci - ad)/(e ci + b) + # * `W = (a ci - ad)/(e ci + b)` # * # * after the minimum is chosen set a, b, e and d for the cubic solution. - # * # * estimate of J according to Farquhar and von Cammerer (1981) - # */ # /*if (jmax > 0) # Jₓ = qalpha * iphoton / sqrt(1. +(qalpha2 * iphoton * iphoton / (jmax * jmax))); - # else - # Jₓ = 0;*/ - # J photon from Harley - Jₓ = Jmax * PPFD / (PPFD + 2.1 * Jmax) # chen1999, eq.6 + Jₓ = Jmax * PPFD / (PPFD + 2.1 * Jmax) # chen1999, eq.6, J photon from Harley # initial guess of intercellular CO2 concentration to estimate Wc and Wj: wj = Jₓ * (cii - Γ) / (4.0 * cii + 8.0 * Γ) @@ -175,14 +133,13 @@ end a_ps = Jₓ b_ps = 8.0 * Γ e_ps = 4.0 - d_ps = Γ else ps_guess = wc a_ps = Vcmax b_ps = K e_ps = 1.0 - d_ps = Γ end + d_ps = Γ # If `wj` or `wc` are less than resp_ld then A would probably be less than 0. # This would yield a negative stomatal conductance. In this case, assume `gs` @@ -218,10 +175,11 @@ end An = findroot(root1, root2, root3) return An end - An = solve_cubic(α_ps, β_ps, γ_ps, θ_ps, Rd, a_ps, b_ps, d_ps, e_ps, ca) + An = solve_cubic(α, β, γ, θ, Rd, a_ps, b_ps, d_ps, e_ps, ca) isnan(An) && (@goto quad) - # also test for sucrose limitation of photosynthesis, as suggested by Collatz. Js=Vmax/2 + # Sucrose limitation of photosynthesis, as suggested by Collatz. `Js=Vmax/2` + # net photosynthesis rate limited by sucrose synthesis (umol m-2 s-1) j_sucrose = Vcmax / 2.0 - Rd An = min(An, j_sucrose) @@ -236,17 +194,15 @@ end @goto OUTDAT end # if aphoto < 0 set stomatal conductance to cuticle value - # /* - # * a quadratic solution of A is derived if gs=b, but a cubic form occur - # * if gs = ax + b. Use quadratic case when A <=0 - # * - # * Bin Chen: - # * r_tot = 1.0/b_co2 + 1.0/g_lb_c; # total resistance to CO2 (m2 s mol-1) - # * denom = g_lb_c * b_co2; - # * Aquad = r_tot * e_ps; - # * Bquad = (e_ps*resp_ld + a_ps)*r_tot - b_ps - e_ps*ca; - # * Cquad = a_ps*(ca-d_ps) - resp_ld*(e_ps*ca+b_ps); - # */ + # A quadratic solution of A is derived if gs=b, but a cubic form occur + # if gs = ax + b. Use quadratic case when A <=0 + # + # Bin Chen: + # r_tot = 1.0/b_co2 + 1.0/g_lb_c; # total resistance to CO2 (m2 s mol-1) + # denom = g_lb_c * b_co2; + # Aquad = r_tot * e_ps; + # Bquad = (e_ps*resp_ld + a_ps)*r_tot - b_ps - e_ps*ca; + # Cquad = a_ps*(ca-d_ps) - resp_ld*(e_ps*ca+b_ps); # original version @label quad ps_1 = ca * g_lb_c * b_co2 @@ -267,12 +223,11 @@ end An = max(0.0, An) cs = ca - An / g_lb_c - gs_h2o_mole = (f_soilwater * m_h2o * RH_leaf * An / cs) + b_h2o # mol m-2 s-1 + gs_h2o_mole = (β_soil * m_h2o * RH_leaf * An / cs) + b_h2o # mol m-2 s-1 gs_co2_mole = gs_h2o_mole / 1.6 ci = cs - An / gs_co2_mole - Gs_w = gs_h2o_mole * temp_leaf_K * (pstat273) # m s-1 - + Gs_w = gs_h2o_mole * T_leaf_K * (pstat273) # m s-1 return Gs_w, An, ci end @@ -327,55 +282,29 @@ end return rate * numm / denom end - -function findroot(root1::Float64, root2::Float64, root3::Float64)::Float64 - root_min::Float64 = 0.0 - root_mid::Float64 = 0.0 - root_max::Float64 = 0.0 - A::Float64 = 0.0 - - if root1 <= root2 && root1 <= root3 - root_min = root1 - if root2 <= root3 - root_mid = root2 - root_max = root3 - else - root_mid = root3 - root_max = root2 - end +function sort3(a::T, b::T, c::T) where {T<:Real} + if a > b + a, b = b, a end - - if root2 <= root1 && root2 <= root3 - root_min = root2 - if root1 <= root3 - root_mid = root1 - root_max = root3 - else - root_mid = root3 - root_max = root1 - end + if b > c + b, c = c, b end - - if root3 <= root1 && root3 <= root2 - root_min = root3 - if root1 < root2 - root_mid = root1 - root_max = root2 - else - root_mid = root2 - root_max = root1 - end + if a > b + a, b = b, a end + return a, b, c +end +function findroot(root1::Float64, root2::Float64, root3::Float64)::Float64 + A::Float64 = 0.0 + root_min, root_mid, root_max = sort3(root1, root2, root3) # find out where roots plop down relative to the x-y axis if root_min > 0 && root_mid > 0 && root_max > 0 A = root_min end - if root_min < 0 && root_mid < 0 && root_max > 0 A = root_max end - if root_min < 0 && root_mid > 0 && root_max > 0 A = root_mid end diff --git a/test/test-beps_main.jl b/test/test-beps_main.jl index ec7de93..6077f13 100644 --- a/test/test-beps_main.jl +++ b/test/test-beps_main.jl @@ -1,7 +1,5 @@ -using Test -using BEPS +using BEPS, DataFrames, Test using BEPS: path_proj -using DataFrames function nanmax(x) x = collect(x) diff --git a/test/test-readparam.jl b/test/test-readparam.jl deleted file mode 100644 index 098be73..0000000 --- a/test/test-readparam.jl +++ /dev/null @@ -1,4 +0,0 @@ -using Test -using BEPS -# import BEPS: readparam, readcoef -