From 8d2f9b41e48591bcad02adbb07b9cc773c960dbc Mon Sep 17 00:00:00 2001 From: Alejandro Date: Wed, 17 Jun 2026 12:45:26 +0200 Subject: [PATCH 1/3] Add property-based tests for C3 and C4 photosynthesis models Introduces test_C3_properties.jl and test_C4_properties.jl, which test mathematical invariants of the FvCB model implementations rather than fixed numeric outputs. Each file covers eight (C3) or seven (C4) properties: monotonic response of assimilation to Ca, PAR, and temperature (peaked, with interior maximum in 10-40 C range); monotonic increase of stomatal conductance with RH; dark respiration (An < 0 at PAR = 0); CO2 compensation point (sign change at physiologically appropriate Ca); minimum stomatal conductance (gs >= gso across a grid of PAR x Ca conditions); and O2 inhibition of assimilation via photorespiration (C3 only). Co-Authored-By: Claude Sonnet 4.6 --- test/Photosynthesis/test_C3_properties.jl | 56 +++++++++++++++++++++++ test/Photosynthesis/test_C4_properties.jl | 50 ++++++++++++++++++++ test/runtests.jl | 2 + 3 files changed, 108 insertions(+) create mode 100644 test/Photosynthesis/test_C3_properties.jl create mode 100644 test/Photosynthesis/test_C4_properties.jl diff --git a/test/Photosynthesis/test_C3_properties.jl b/test/Photosynthesis/test_C3_properties.jl new file mode 100644 index 0000000..9e780cf --- /dev/null +++ b/test/Photosynthesis/test_C3_properties.jl @@ -0,0 +1,56 @@ + +using Test +import Ecophys +PH = Ecophys.Photosynthesis + +@testset "C3 photosynthesis properties" begin + + @testset "Monotonic A with Ca" begin + c3 = PH.C3() + A_vals = [PH.photosynthesis(c3, Ca = ca).A for ca in 50.0:50.0:2000.0] + @test all(diff(A_vals) .>= 0) + end + + @testset "Monotonic A with PAR" begin + c3 = PH.C3() + A_vals = [PH.photosynthesis(c3, PAR = par).A for par in 0.0:50.0:2000.0] + @test all(diff(A_vals) .>= 0) + end + + @testset "gs increases with RH" begin + c3 = PH.C3() + gs_vals = [PH.photosynthesis(c3, RH = rh).gs for rh in 0.10:0.05:0.95] + @test all(diff(gs_vals) .>= 0) + end + + @testset "Peaked temperature response (10-40 C)" begin + c3 = PH.C3() + A_vals = [PH.photosynthesis(c3, Tleaf = 273.15 + t).A for t in 10.0:1.0:40.0] + peak_idx = argmax(A_vals) + @test peak_idx > 1 && peak_idx < length(A_vals) + end + + @testset "Dark respiration at PAR = 0" begin + @test PH.photosynthesis(PH.C3(), PAR = 0.0).A < 0 + end + + @testset "CO2 compensation point" begin + c3 = PH.C3() + @test PH.photosynthesis(c3, Ca = 20.0).A < 0 + @test PH.photosynthesis(c3, Ca = 400.0).A > 0 + end + + @testset "gs >= gso" begin + c3 = PH.C3() + results = [PH.photosynthesis(c3, PAR = par, Ca = ca) + for par in 100.0:100.0:2000.0, ca in 200.0:100.0:2000.0] + @test all(r.gs >= c3.gso for r in results) + end + + @testset "O2 inhibition" begin + c3 = PH.C3() + A_vals = [PH.photosynthesis(c3, O2 = o2).A for o2 in 10e3:10e3:210e3] + @test all(diff(A_vals) .<= 0) + end + +end diff --git a/test/Photosynthesis/test_C4_properties.jl b/test/Photosynthesis/test_C4_properties.jl new file mode 100644 index 0000000..e5a3990 --- /dev/null +++ b/test/Photosynthesis/test_C4_properties.jl @@ -0,0 +1,50 @@ + +using Test +import Ecophys +PH = Ecophys.Photosynthesis + +@testset "C4 photosynthesis properties" begin + + @testset "Monotonic A with Ca" begin + c4 = PH.C4() + A_vals = [PH.photosynthesis(c4, Ca = ca).A for ca in 10.0:50.0:2000.0] + @test all(diff(A_vals) .>= 0) + end + + @testset "Monotonic A with PAR" begin + c4 = PH.C4() + A_vals = [PH.photosynthesis(c4, PAR = par).A for par in 0.0:50.0:2000.0] + @test all(diff(A_vals) .>= 0) + end + + @testset "gs increases with RH" begin + c4 = PH.C4() + gs_vals = [PH.photosynthesis(c4, RH = rh).gs for rh in 0.10:0.05:0.95] + @test all(diff(gs_vals) .>= 0) + end + + @testset "Peaked temperature response (10-40 C)" begin + c4 = PH.C4() + A_vals = [PH.photosynthesis(c4, Tleaf = 273.15 + t).A for t in 10.0:1.0:40.0] + peak_idx = argmax(A_vals) + @test peak_idx > 1 && peak_idx < length(A_vals) + end + + @testset "Dark respiration at PAR = 0" begin + @test PH.photosynthesis(PH.C4(), PAR = 0.0).A < 0 + end + + @testset "CO2 compensation point" begin + c4 = PH.C4() + @test PH.photosynthesis(c4, Ca = 1.0).A < 0 + @test PH.photosynthesis(c4, Ca = 400.0).A > 0 + end + + @testset "gs >= gso" begin + c4 = PH.C4() + results = [PH.photosynthesis(c4, PAR = par, Ca = ca) + for par in 100.0:100.0:2000.0, ca in 200.0:100.0:2000.0] + @test all(r.gs >= c4.gso for r in results) + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index dedc1d3..04d56d9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,8 @@ import Aqua include("Photosynthesis/test_components.jl") include("Photosynthesis/test_C3.jl") include("Photosynthesis/test_C4.jl") + include("Photosynthesis/test_C3_properties.jl") + include("Photosynthesis/test_C4_properties.jl") include("Photosynthesis/test_gb.jl") include("Photosynthesis/test_energybalance.jl") end From 83f53716b12afeb54c3b33f6cb45c7b394e23680 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Wed, 17 Jun 2026 12:46:36 +0200 Subject: [PATCH 2/3] Change from calculating gs to calculating Ci for both C3 and C4 --- Project.toml | 2 +- src/Photosynthesis/Components.jl | 6 ++-- src/Photosynthesis/FvCB/C3.jl | 49 ++++++++++++++++---------------- src/Photosynthesis/FvCB/C4.jl | 27 +++++++++--------- test/Photosynthesis/test_C3.jl | 2 +- test/Photosynthesis/test_C4.jl | 2 +- 6 files changed, 45 insertions(+), 43 deletions(-) diff --git a/Project.toml b/Project.toml index ded9049..1dcab03 100644 --- a/Project.toml +++ b/Project.toml @@ -14,4 +14,4 @@ Roots = "2.0.19" SpecialFunctions = "2.3.1" StaticArrays = "v1.9.7" Unitful = "1.17.0" -julia = "1" +julia = "1.12" diff --git a/src/Photosynthesis/Components.jl b/src/Photosynthesis/Components.jl index 1ec9748..7aa3e3c 100644 --- a/src/Photosynthesis/Components.jl +++ b/src/Photosynthesis/Components.jl @@ -42,8 +42,8 @@ Tb(::Type{T}) where {T <: Real} = T(35.86) Tb(::Type{T}) where {T <: Quantity} = 35.86K # es(T = 35.86 K) in Murray's model -es0(::Type{T}) where {T <: Real} = T(610.78) -es0(::Type{T}) where {T <: Quantity} = 610.78Pa +es0(::Type{T}) where {T <: Real} = T(0.61078) # kPa +es0(::Type{T}) where {T <: Quantity} = 0.61078kPa # 0 umol/m2/s zeroflux(::Type{T}) where {T <: Real} = T(0.0) # μmol/m^2/s @@ -126,7 +126,7 @@ function ThermalDiffusivity(Tavg::T) where {T <: Quantity} (0.0001Tavg * Tavg / K^2 + 0.0561 * Tavg / K - 6.1952) * 1e-6m^2 / s # m^2/s end -# Molar volume +# Molar volume MolarVolume(Tavg, P) = GasConstant(typeof(Tavg)) * Tavg / P # Air density diff --git a/src/Photosynthesis/FvCB/C3.jl b/src/Photosynthesis/FvCB/C3.jl index 85ced00..6c8c3ed 100644 --- a/src/Photosynthesis/FvCB/C3.jl +++ b/src/Photosynthesis/FvCB/C3.jl @@ -7,14 +7,14 @@ abstract type FvCB <: Ags end abstract type C3Type <: FvCB end """ - C3(Sco25 = 2800.0, E_Sco = -24.46e3, Kmc25 = 270.0, E_Kmc = 80.99e3, - Kmo25 = 165.0e3, E_Kmo = 23.72e3, Vcmax25 = 120.0, E_Vcmax = 65.33e3, - simpleJ = false, k2ll = 0.35, theta = 0.7, Phi2 = 0.82, sigma2 = 0.5, - beta = 0.85, fcyc = 0.1, - fpseudo = 0.05, Jmax25 = 230.0, E_Jmax = 30.0e3, D_Jmax = 200.0e3, + C3(Sco25 = 2800.0, E_Sco = -24.46e3, Kmc25 = 270.0, E_Kmc = 80.99e3, + Kmo25 = 165.0e3, E_Kmo = 23.72e3, Vcmax25 = 120.0, E_Vcmax = 65.33e3, + simpleJ = false, k2ll = 0.35, theta = 0.7, Phi2 = 0.82, sigma2 = 0.5, + beta = 0.85, fcyc = 0.1, + fpseudo = 0.05, Jmax25 = 230.0, E_Jmax = 30.0e3, D_Jmax = 200.0e3, S_Jmax = 650.0, TPU25 = 12.0, E_TPU = 53.1e3, D_TPU = 20.18e3, - S_TPU = 650.0, Rd25 = 1.2, E_Rd = 46.39e3, gm25 = 0.4, E_gm = 49.6e3, - D_gm = 437.4e3, S_gm = 1400.0, gso = 0.01, a1 = 0.85, b1 = 0.14e-3) + S_TPU = 650.0, Rd25 = 1.2, E_Rd = 46.39e3, gm25 = 0.4, E_gm = 49.6e3, + D_gm = 437.4e3, S_gm = 1400.0, gso = 0.01, a1 = 0.85, b1 = 0.14) Data structure to store all the parameters for the C3 photosynthesis model. @@ -51,7 +51,7 @@ Data structure to store all the parameters for the C3 photosynthesis model. - `S_gm`: Entropy term for gm (K) - `gso`: Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s/Pa) - `a1`: Empirical parameter in gs formula -- `b1`: Empirical parameter in gs formula +- `b1`: Empirical parameter in gs formula (1/kPa) """ Base.@kwdef mutable struct C3{T <: Real} <: C3Type # Rubisco CO2/O2 specificity @@ -95,17 +95,17 @@ Base.@kwdef mutable struct C3{T <: Real} <: C3Type # Stomatal conductance gso::T = 0.01 # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.85 # Empirical parameter in gs formula - b1::T = 0.14e-3 # Empirical parameter in gs formula (1/kPa) + b1::T = 0.14 # Empirical parameter in gs formula (1/kPa) end """ C3Q(Sco25 = 2800.0, E_Sco = -24.46e3J/mol, Kmc25 = 270.0μmol/mol, E_Kmc = 80.99e3J/mol, Kmo25 = 165.0e3μmol/mol, E_Kmo = 23.72e3J/mol, Vcmax25 = 120.0μmol/m^2/s, E_Vcmax = 65.33e3J/mol, - simpleJ = false, k2ll = 0.35, theta = 0.7, Phi2 = 0.82, sigma2 = 0.5, beta = 0.85, fcyc = 0.1, fpseudo = 0.05, - Jmax25 = 230.0μmol/m^2/s, E_Jmax = 30.0e3J/mol, D_Jmax = 200.0e3J/mol, S_Jmax = 650.0J/mol/K, - TPU25 = 12.0μmol/m^2/s, E_TPU = 53.1e3J/mol, D_TPU = 201.8e3J/mol, S_TPU = 650.0K, - Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gm25 = 0.4mol/m^2/s, E_gm = 49.6e3J/mol, - D_gm = 437.4e3J/mol, S_gm = 1400.0K, gso = 0.01mol/m^2/s, a1 = 0.85, b1 = 0.14e-3/Pa) + simpleJ = false, k2ll = 0.35, theta = 0.7, Phi2 = 0.82, sigma2 = 0.5, beta = 0.85, fcyc = 0.1, fpseudo = 0.05, + Jmax25 = 230.0μmol/m^2/s, E_Jmax = 30.0e3J/mol, D_Jmax = 200.0e3J/mol, S_Jmax = 650.0J/mol/K, + TPU25 = 12.0μmol/m^2/s, E_TPU = 53.1e3J/mol, D_TPU = 201.8e3J/mol, S_TPU = 650.0K, + Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gm25 = 0.4mol/m^2/s, E_gm = 49.6e3J/mol, + D_gm = 437.4e3J/mol, S_gm = 1400.0K, gso = 0.01mol/m^2/s, a1 = 0.85, b1 = 0.14/kPa) Data structure to store all the parameters for the C3 photosynthesis model using `Quantity` objects from Unitful.jl. @@ -187,7 +187,7 @@ Base.@kwdef mutable struct C3Q{T <: Real} <: C3Type # Stomatal conductance gso::Quantity{T, dimension(mol / m^2 / s)} = 0.01mol / m^2 / s # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.85 # Empirical parameter in gs formula - b1::Quantity{T, dimension(1 / kPa)} = 0.14e-3 / Pa # Empirical parameter in gs formula (1/kPa) + b1::Quantity{T, dimension(1 / kPa)} = 0.14 / kPa # Empirical parameter in gs formula (1/kPa) end """ @@ -197,10 +197,10 @@ end photosynthesis(par::C4Q, PAR = 1000.0μmol/m^2/s, RH = 0.75, Tleaf = 298.0K, Ca = 400.0μmol/mol, O2 = 210e3μmol/mol, gb = 0.5mol/m^2/s, net = true) Calculate net or gross CO2 assimilation (umol/m2/s) -and stomatal condutance to fluxes of CO2 (mol/m2/s) as a function of +and stomatal condutance to fluxes of CO2 (mol/m2/s) as a function of photosynthetically active radiation (PAR, umol/m2/s), relative humidity (RH), leaf temperature (Tleaf, -K), air CO2 partial pressure (Ca, μmol/mol), oxygen (O2, μmol/mol) and boundary layer +K), air CO2 partial pressure (Ca, μmol/mol), oxygen (O2, μmol/mol) and boundary layer conductance to CO2 (gb, mol/m2/s). Environmental inputs must be scalar. The argument `net` indicates whether the net or gross CO2 assimilation should be returned. """ @@ -259,7 +259,8 @@ function photosynthesis(p::C3Type, PAR, RH, Tleaf, Ca, O2, gb, net) An = min(Ac, min(Aj, Ap)) # μmol/m2/s # Stomatal conductance - gsc = solvegs(p.gso, An, Ca, Ci_star, Rd, fvpd, gb) # mol/m2/s + Ci = CalcCi(p.gso, An, Ca, Ci_star, Rd, fvpd) + gsc = p.gso + ((An + Rd)/(Ci - Ci_star))*fvpd # Choose the right output A = net ? An : An + Rd @@ -282,10 +283,10 @@ function solveAC3(gm, gb, gso, fvpd, x2, x1, gamma_star, Rd, Ca) A = -2 * sqrt(Q) * cos(psi / 3) - p / 3 end -# Calculate gs once A is known -function solvegs(gso, A, Ca, Ci_star, Rd, fvpd, gb) - a = Ca - A / gb - Ci_star - b = -A - Ca * gso + gso * Ci_star - (A + Rd) * fvpd - c = A * gso - A = (-b - sqrt(b * b - 4 * a * c)) / (2 * a) +# Calculation of internal CO2 concentration +function CalcCi(gs0, An, Ca, Ci_star, Rd, fvpd) + a = gs0 + b = An - gs0*Ca - gs0*Ci_star + (An + Rd)*fvpd + c = -An*Ci_star + gs0*Ca*Ci_star - (An + Rd)*Ca*fvpd + Ci = (-b + sqrt(b^2 - 4*a*c))/(2*a) end diff --git a/src/Photosynthesis/FvCB/C4.jl b/src/Photosynthesis/FvCB/C4.jl index 788c500..e0999ef 100644 --- a/src/Photosynthesis/FvCB/C4.jl +++ b/src/Photosynthesis/FvCB/C4.jl @@ -7,11 +7,11 @@ abstract type C4Type <: FvCB end # Data structure to store all the C4 parameters without units """ - C4(Sco25 = 2590.0, E_Sco = -24.46e3, Kmc25 = 650.0, E_Kmc = 79.43e3, Kmo25 = 450e3, - E_Kmo = 36380.0, Vcmax25 = 120.0, E_Vcmax = 65.33, theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, - beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, Jmax25 = 230.0, E_Jmax = 48e3, D_Jmax = 200e3, - S_Jmax = 630.0, x = 0.4, alpha = 0.1, kp25 = 0.7, E_kp = 46.39e3, gbs = 0.003, Rd25 = 1.2, - E_Rd = 46.39e3, gso = 0.01, a1 = 0.9, b1 = 0.15e-3) + C4(Sco25 = 2590.0, E_Sco = -24.46e3, Kmc25 = 650.0, E_Kmc = 79.43e3, Kmo25 = 450e3, + E_Kmo = 36380.0, Vcmax25 = 120.0, E_Vcmax = 65.33, theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, + beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, Jmax25 = 230.0, E_Jmax = 48e3, D_Jmax = 200e3, + S_Jmax = 630.0, x = 0.4, alpha = 0.1, kp25 = 0.7, E_kp = 46.39e3, gbs = 0.003, Rd25 = 1.2, + E_Rd = 46.39e3, gso = 0.01, a1 = 0.9, b1 = 0.15) Data structure to store all the parameters for the C3 photosynthesis model. @@ -81,16 +81,16 @@ Base.@kwdef mutable struct C4{T <: Real} <: C4Type # Stomatal conductance gso::T = 0.01 # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.9 # Empirical parameter in gs formula - b1::T = 0.15e-3 # Empirical parameter in gs formula (1/kPa) + b1::T = 0.15 # Empirical parameter in gs formula (1/kPa) end """ C4(Sco25 = 2590.0, E_Sco = -24.46e3J/mol, Kmc25 = 650.0μmol/mol, E_Kmc = 79.43e3J/mol, Kmo25 = 450e3μmol/mol, E_Kmo = 36380.0J/mol, Vcmax25 = 120.0μmol/m^2/s, E_Vcmax = 65.33J/mol, - theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, - Jmax25 = 230.0μmol/m^2/s, E_Jmax = 48e3J/mol, D_Jmax = 200e3J/mol, S_Jmax = 630.0J/mol/K, - x = 0.4, alpha = 0.1, kp25 = 0.7mol/m^2/s, E_kp = 46.39e3J/mol, gbs = 0.003mol/m^2/s, - Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gso = 0.01mol/m^2/s, a1 = 0.9, b1 = 0.15e-3/Pa) + theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, + Jmax25 = 230.0μmol/m^2/s, E_Jmax = 48e3J/mol, D_Jmax = 200e3J/mol, S_Jmax = 630.0J/mol/K, + x = 0.4, alpha = 0.1, kp25 = 0.7mol/m^2/s, E_kp = 46.39e3J/mol, gbs = 0.003mol/m^2/s, + Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gso = 0.01mol/m^2/s, a1 = 0.9, b1 = 0.15/kPa) Data structure to store all the parameters for the C4 photosynthesis model using `Quantity` objects from Unitful.jl. @@ -161,7 +161,7 @@ Base.@kwdef mutable struct C4Q{T <: Real} <: C4Type # Stomatal conductance gso::Quantity{T, dimension(mol / m^2 / s)} = 0.01mol / m^2 / s # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.9 # Empirical parameter in gs formula - b1::Quantity{T, dimension(1 / kPa)} = 0.15e-3 / Pa # Empirical parameter in gs formula (1/kPa) + b1::Quantity{T, dimension(1 / kPa)} = 0.15 / kPa # Empirical parameter in gs formula (1/kPa) end function photosynthesis(p::C4; @@ -306,8 +306,9 @@ function photosynthesis(p::C4Type, PAR, RH, Tleaf, Ca, O2, gb, net) Aj = min(Aj1, Aj2) An = min(Ac, Aj) - # Stomatal conductance - gsc = solvegs(p.gso, An, Ca, Cs_star, Rd, fvpd, gb) # mol/m2/s + # Stomatal conductance (using Cs_star instead of Ci_star) + Ci = CalcCi(p.gso, An, Ca, Cs_star, Rd, fvpd) + gsc = p.gso + ((An + Rd)/(Ci - Cs_star))*fvpd # Choose the right output A = net ? An : An + Rd diff --git a/test/Photosynthesis/test_C3.jl b/test/Photosynthesis/test_C3.jl index c0710f1..d3a042c 100644 --- a/test/Photosynthesis/test_C3.jl +++ b/test/Photosynthesis/test_C3.jl @@ -23,7 +23,7 @@ let O2 = O2_f, gb = gb_f) @test abs(A_f - 26.15913) < 1e-4 - @test abs(gs_f - 0.0022) < 1e-4 + @test abs(gs_f - 0.31806) < 1e-4 Ag_f, gs_f = PH.photosynthesis(c3, PAR = PAR_f, RH = RH_f, diff --git a/test/Photosynthesis/test_C4.jl b/test/Photosynthesis/test_C4.jl index d4c8d15..77105a1 100644 --- a/test/Photosynthesis/test_C4.jl +++ b/test/Photosynthesis/test_C4.jl @@ -23,7 +23,7 @@ let O2 = O2_f, gb = gb_f) @test abs(A_f - 30.5477) < 0.0001 - @test abs(gs_f - 0.0018) < 1e-4 + @test abs(gs_f - 0.4114188) < 1e-4 Ag_f, gs_f = PH.photosynthesis(c4, PAR = PAR_f, RH = RH_f, From 573830f4539a5f7a17cc5023baeab67f5e9be883 Mon Sep 17 00:00:00 2001 From: Alejandro Date: Wed, 17 Jun 2026 13:29:55 +0200 Subject: [PATCH 3/3] Undo changes to b1: we need in Pa because of the energy balance. --- src/Photosynthesis/Components.jl | 4 ++-- .../EnergyBalance/EnergyBalance.jl | 18 ++++++++--------- src/Photosynthesis/FvCB/C3.jl | 12 +++++------ src/Photosynthesis/FvCB/C4.jl | 12 +++++------ test/Photosynthesis/test_C3.jl | 2 +- test/Photosynthesis/test_components.jl | 2 +- test/Photosynthesis/test_energybalance.jl | 20 +++++++++---------- 7 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/Photosynthesis/Components.jl b/src/Photosynthesis/Components.jl index 7aa3e3c..3a1ee2f 100644 --- a/src/Photosynthesis/Components.jl +++ b/src/Photosynthesis/Components.jl @@ -42,8 +42,8 @@ Tb(::Type{T}) where {T <: Real} = T(35.86) Tb(::Type{T}) where {T <: Quantity} = 35.86K # es(T = 35.86 K) in Murray's model -es0(::Type{T}) where {T <: Real} = T(0.61078) # kPa -es0(::Type{T}) where {T <: Quantity} = 0.61078kPa +es0(::Type{T}) where {T <: Real} = T(610.78) # Pa +es0(::Type{T}) where {T <: Quantity} = 610.78Pa # 0 umol/m2/s zeroflux(::Type{T}) where {T <: Real} = T(0.0) # μmol/m^2/s diff --git a/src/Photosynthesis/EnergyBalance/EnergyBalance.jl b/src/Photosynthesis/EnergyBalance/EnergyBalance.jl index 16ddbb3..b277a14 100644 --- a/src/Photosynthesis/EnergyBalance/EnergyBalance.jl +++ b/src/Photosynthesis/EnergyBalance/EnergyBalance.jl @@ -3,7 +3,7 @@ # energybalance # solve_energy_balance -import Roots: find_zero, Order2 +import Roots: find_zero, A42 abstract type Optical end @@ -40,7 +40,7 @@ Calculate the energy balance of a leaf. - `RH`: Relative humidity - `Tair`: Air temperature (K) - `Ca`: Atmospheric CO2 concentration (μmol/mol) -- `P`: Air pressure (kPa) +- `P`: Air pressure (Pa) - `O2`: Atmospheric O2 concentration (μmol/mol) # Details @@ -94,12 +94,12 @@ end opt = SimpleOptical(), PAR = 1000.0μmol/m^2/s, NIR = 250.0W/m^2, ws = 1.0m/s, RH = 0.75, Tair = 298.0K, Ca = 400.0μmol/mol, P = 101.0kPa, - O2 = 210.0mmol/mol, order = Order2(), xatol = 0.01, + O2 = 210.0mmol/mol, order = A42(), xatol = 0.01, maxfnevals = 100, net = true) solve_energy_balance(Ags::Union{C3, C4}; gb = simplegb(), opt = SimpleOptical(), PAR = 1000.0, NIR = 250.0, ws = 1.0, RH = 0.75, Tair = 298.0, Ca = 400.0, - P = 101.0e3, O2 = 210.0e3, order = Order2(), xatol = 0.01, + P = 101.0e3, O2 = 210.0e3, order = A42(), xatol = 0.01, maxfnevals = 100, net = true) Solve the leaf energy balance coupled to photosynthesis and transpiration. @@ -139,7 +139,7 @@ function solve_energy_balance(Ags::Union{C3Q, C4Q}; gb = simplegbQ(), opt = SimpleOptical(), PAR = 1000.0μmol / m^2 / s, NIR = 250.0W / m^2, ws = 1.0m / s, RH = 0.75, Tair = 298.0K, Ca = 400.0μmol / mol, P = 101.0kPa, - O2 = 210.0mmol / mol, order = Order2(), xatol = 0.01, + O2 = 210.0mmol / mol, order = A42(), xatol = 0.01, maxfnevals = 100, net = true) solve_energy_balance(Ags, gb, opt, PAR, NIR, ws, RH, Tair, Ca, P, O2, order, xatol, maxfnevals, net) @@ -147,14 +147,14 @@ end function solve_energy_balance(Ags::Union{C3, C4}; gb = simplegb(), opt = SimpleOptical(), PAR = 1000.0, NIR = 250.0, ws = 1.0, RH = 0.75, Tair = 298.0, Ca = 400.0, - P = 101.0e3, O2 = 210.0e3, order = Order2(), xatol = 0.01, + P = 101.0e3, O2 = 210.0e3, order = A42(), xatol = 0.01, maxfnevals = 100, net = true) solve_energy_balance(Ags, gb, opt, PAR, NIR, ws, RH, Tair, Ca, P, O2, order, xatol, maxfnevals, net) end function solve_energy_balance(pAgs, pgb, pEb, PAR, NIR, ws, RH, Tair, Ca, P, O2, - order = Order2(), xatol = 0.01, maxfnevals = 100, net = true) + order = A42(), xatol = 0.01, maxfnevals = 100, net = true) # Find the temperature Tleaf = find_zero(x -> energybalance(x, pgb, @@ -168,7 +168,7 @@ function solve_energy_balance(pAgs, pgb, pEb, PAR, NIR, ws, RH, Tair, Ca, P, O2, Ca, P, O2), - (Tair - 10, Tair + 10), order, xatol = xatol, maxfnevals = maxfnevals) + (Tair - 10, Tair + 10), order, xatol = xatol, maxfnevals = maxfnevals, verbose = true) # Boundary layer conductances gbh, gbw, gbc = gb(pgb, ws, Tleaf, Tair, P) # A and gsc @@ -192,7 +192,7 @@ function solve_energy_balance(pAgs, pgb, pEb, PAR, NIR, ws, RH, Tair, Ca, P, O2, end function solve_energy_balance(pAgs, pgb, pEb, PAR::Quantity, NIR::Quantity, ws::Quantity, RH, Tair::Quantity, Ca::Quantity, - P::Quantity, O2::Quantity, order = Order2(), + P::Quantity, O2::Quantity, order = A42(), xatol = 0.01, maxfnevals = 100, net = true) Tau = Tair / 1.0K Tleaf = find_zero(x -> energybalance(x * K, diff --git a/src/Photosynthesis/FvCB/C3.jl b/src/Photosynthesis/FvCB/C3.jl index 6c8c3ed..dd78d37 100644 --- a/src/Photosynthesis/FvCB/C3.jl +++ b/src/Photosynthesis/FvCB/C3.jl @@ -14,7 +14,7 @@ abstract type C3Type <: FvCB end fpseudo = 0.05, Jmax25 = 230.0, E_Jmax = 30.0e3, D_Jmax = 200.0e3, S_Jmax = 650.0, TPU25 = 12.0, E_TPU = 53.1e3, D_TPU = 20.18e3, S_TPU = 650.0, Rd25 = 1.2, E_Rd = 46.39e3, gm25 = 0.4, E_gm = 49.6e3, - D_gm = 437.4e3, S_gm = 1400.0, gso = 0.01, a1 = 0.85, b1 = 0.14) + D_gm = 437.4e3, S_gm = 1400.0, gso = 0.01, a1 = 0.85, b1 = 0.14e-3) Data structure to store all the parameters for the C3 photosynthesis model. @@ -51,7 +51,7 @@ Data structure to store all the parameters for the C3 photosynthesis model. - `S_gm`: Entropy term for gm (K) - `gso`: Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s/Pa) - `a1`: Empirical parameter in gs formula -- `b1`: Empirical parameter in gs formula (1/kPa) +- `b1`: Empirical parameter in gs formula (1/Pa) """ Base.@kwdef mutable struct C3{T <: Real} <: C3Type # Rubisco CO2/O2 specificity @@ -95,7 +95,7 @@ Base.@kwdef mutable struct C3{T <: Real} <: C3Type # Stomatal conductance gso::T = 0.01 # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.85 # Empirical parameter in gs formula - b1::T = 0.14 # Empirical parameter in gs formula (1/kPa) + b1::T = 0.14e-3 # Empirical parameter in gs formula (1/Pa) end """ @@ -105,7 +105,7 @@ end Jmax25 = 230.0μmol/m^2/s, E_Jmax = 30.0e3J/mol, D_Jmax = 200.0e3J/mol, S_Jmax = 650.0J/mol/K, TPU25 = 12.0μmol/m^2/s, E_TPU = 53.1e3J/mol, D_TPU = 201.8e3J/mol, S_TPU = 650.0K, Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gm25 = 0.4mol/m^2/s, E_gm = 49.6e3J/mol, - D_gm = 437.4e3J/mol, S_gm = 1400.0K, gso = 0.01mol/m^2/s, a1 = 0.85, b1 = 0.14/kPa) + D_gm = 437.4e3J/mol, S_gm = 1400.0K, gso = 0.01mol/m^2/s, a1 = 0.85, b1 = 0.14e-3/Pa) Data structure to store all the parameters for the C3 photosynthesis model using `Quantity` objects from Unitful.jl. @@ -143,7 +143,7 @@ Data structure to store all the parameters for the C3 photosynthesis model using - `S_gm`: Entropy term for gm (J/K/mol) - `gso`: Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) - `a1`: Empirical parameter in gs formula -- `b1`: Empirical parameter in gs formula (1/kPa) +- `b1`: Empirical parameter in gs formula (1/Pa) """ Base.@kwdef mutable struct C3Q{T <: Real} <: C3Type # Rubisco CO2/O2 specificity @@ -187,7 +187,7 @@ Base.@kwdef mutable struct C3Q{T <: Real} <: C3Type # Stomatal conductance gso::Quantity{T, dimension(mol / m^2 / s)} = 0.01mol / m^2 / s # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.85 # Empirical parameter in gs formula - b1::Quantity{T, dimension(1 / kPa)} = 0.14 / kPa # Empirical parameter in gs formula (1/kPa) + b1::Quantity{T, dimension(1 / Pa)} = 0.14e-3 / Pa # Empirical parameter in gs formula (1/Pa) end """ diff --git a/src/Photosynthesis/FvCB/C4.jl b/src/Photosynthesis/FvCB/C4.jl index e0999ef..518d247 100644 --- a/src/Photosynthesis/FvCB/C4.jl +++ b/src/Photosynthesis/FvCB/C4.jl @@ -11,7 +11,7 @@ abstract type C4Type <: FvCB end E_Kmo = 36380.0, Vcmax25 = 120.0, E_Vcmax = 65.33, theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, Jmax25 = 230.0, E_Jmax = 48e3, D_Jmax = 200e3, S_Jmax = 630.0, x = 0.4, alpha = 0.1, kp25 = 0.7, E_kp = 46.39e3, gbs = 0.003, Rd25 = 1.2, - E_Rd = 46.39e3, gso = 0.01, a1 = 0.9, b1 = 0.15) + E_Rd = 46.39e3, gso = 0.01, a1 = 0.9, b1 = 0.15e-3) Data structure to store all the parameters for the C3 photosynthesis model. @@ -44,7 +44,7 @@ Data structure to store all the parameters for the C3 photosynthesis model. - `E_Rd`: Activation energy of Rd (J/mol) - `gso`: Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) - `a1`: Empirical parameter in gs formula -- `b1`: Empirical parameter in gs formula (1/kPa) +- `b1`: Empirical parameter in gs formula (1/Pa) """ Base.@kwdef mutable struct C4{T <: Real} <: C4Type # Rubisco @@ -81,7 +81,7 @@ Base.@kwdef mutable struct C4{T <: Real} <: C4Type # Stomatal conductance gso::T = 0.01 # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.9 # Empirical parameter in gs formula - b1::T = 0.15 # Empirical parameter in gs formula (1/kPa) + b1::T = 0.15e-3 # Empirical parameter in gs formula (1/Pa) end """ @@ -90,7 +90,7 @@ end theta = 0.7, Phi2 = 0.83, sigma2 = 0.5, beta = 0.85, fQ = 1.0, fpseudo = 0.1, h = 4.0, Jmax25 = 230.0μmol/m^2/s, E_Jmax = 48e3J/mol, D_Jmax = 200e3J/mol, S_Jmax = 630.0J/mol/K, x = 0.4, alpha = 0.1, kp25 = 0.7mol/m^2/s, E_kp = 46.39e3J/mol, gbs = 0.003mol/m^2/s, - Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gso = 0.01mol/m^2/s, a1 = 0.9, b1 = 0.15/kPa) + Rd25 = 1.2μmol/m^2/s, E_Rd = 46.39e3J/mol, gso = 0.01mol/m^2/s, a1 = 0.9, b1 = 0.15e-3/Pa) Data structure to store all the parameters for the C4 photosynthesis model using `Quantity` objects from Unitful.jl. @@ -124,7 +124,7 @@ Data structure to store all the parameters for the C4 photosynthesis model using - `E_Rd`: Activation energy of Rd (J/mol) - `gso`: Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) - `a1`: Empirical parameter in gs formula -- `b1`: Empirical parameter in gs formula (1/kPa) +- `b1`: Empirical parameter in gs formula (1/Pa) """ Base.@kwdef mutable struct C4Q{T <: Real} <: C4Type # Rubisco @@ -161,7 +161,7 @@ Base.@kwdef mutable struct C4Q{T <: Real} <: C4Type # Stomatal conductance gso::Quantity{T, dimension(mol / m^2 / s)} = 0.01mol / m^2 / s # Minimum stomatal conductance to fluxes of CO2 in darkness (mol/m2/s) a1::T = 0.9 # Empirical parameter in gs formula - b1::Quantity{T, dimension(1 / kPa)} = 0.15 / kPa # Empirical parameter in gs formula (1/kPa) + b1::Quantity{T, dimension(1 / kPa)} = 0.15e-3 / Pa # Empirical parameter in gs formula (1/Pa) end function photosynthesis(p::C4; diff --git a/test/Photosynthesis/test_C3.jl b/test/Photosynthesis/test_C3.jl index d3a042c..d2d0051 100644 --- a/test/Photosynthesis/test_C3.jl +++ b/test/Photosynthesis/test_C3.jl @@ -23,7 +23,7 @@ let O2 = O2_f, gb = gb_f) @test abs(A_f - 26.15913) < 1e-4 - @test abs(gs_f - 0.31806) < 1e-4 + @test abs(gs_f - 0.318063) < 1e-4 Ag_f, gs_f = PH.photosynthesis(c3, PAR = PAR_f, RH = RH_f, diff --git a/test/Photosynthesis/test_components.jl b/test/Photosynthesis/test_components.jl index 6fb002b..9994ca6 100644 --- a/test/Photosynthesis/test_components.jl +++ b/test/Photosynthesis/test_components.jl @@ -39,7 +39,7 @@ let @test PH.peaked(1.0mol, 26900.0J / mol, 2e5J / mol, 650.0J / K / mol, 298.15K) == 1.0mol # Saturated vapour pressure - @test abs(PH.es(298.15) - 3167.69) < 0.01 + @test abs(PH.es(298.15) - 3167.688063635486) < 0.01 @test abs(PH.es(298.15K) - 3.16kPa) < 0.01kPa @test PH.es(35.86) == 0.0 @test PH.es(35.86K) == 0.0kPa diff --git a/test/Photosynthesis/test_energybalance.jl b/test/Photosynthesis/test_energybalance.jl index ac4d1b3..b9d8bf9 100644 --- a/test/Photosynthesis/test_energybalance.jl +++ b/test/Photosynthesis/test_energybalance.jl @@ -27,7 +27,7 @@ O2 = 210e3 # μmol/mol res = PH.energybalance(Tl, gb, c3, opt, PAR, NIR, ws, RH, Tair, Ca, P, O2) # W/m^2 Tl = PH.solve_energy_balance(c3, gb = gb, opt = opt, PAR = PAR, NIR = NIR, ws = ws, RH = RH, Tair = Tair, Ca = Ca, P = P, O2 = O2, - xatol = 0.01, order = Order2()).Tleaf + xatol = 0.01, order = A42()).Tleaf # With units opt = PH.SimpleOptical() @@ -58,8 +58,8 @@ res_q = PH.energybalance(Tl_q, Tl_q = PH.solve_energy_balance(c3_q, gb = gb_q, opt = opt, PAR = PAR_q, NIR = NIR_q, ws = ws_q, RH = RH_q, Tair = Ta_q, Ca = Ca_q, P = P_q, O2 = O2_q, xatol = 0.01, - order = Order2()).Tleaf -@test Tl_q ≈ Tl * K + order = A42()).Tleaf +@test abs(Tl_q - Tl * K) < 0.02K # With effect of angles (horizontal) gb_q = PH.gbAngleQ() @@ -80,8 +80,8 @@ res_q = PH.energybalance(Tl_q, Tl_q = PH.solve_energy_balance(c3_q, gb = gb_q, opt = opt, PAR = PAR_q, NIR = NIR_q, ws = ws_q, RH = RH_q, Tair = Ta_q, Ca = Ca_q, P = P_q, O2 = O2_q, xatol = 0.01, - order = Order2()).Tleaf -@test Tl_q ≈ Tl * K + order = A42()).Tleaf +@test abs(Tl_q - Tl * K) < 0.02K ############################################################################################################# ################################################# C4 leaves ################################################# @@ -93,7 +93,7 @@ Tl = 300.0 res = PH.energybalance(Tl, gb, c4, opt, PAR, NIR, ws, RH, Tair, Ca, P, O2) # W/m^2 Tl = PH.solve_energy_balance(c4, gb = gb, opt = opt, PAR = PAR, NIR = NIR, ws = ws, RH = RH, Tair = Tair, Ca = Ca, P = P, O2 = O2, - xatol = 0.01, order = Order2()).Tleaf + xatol = 0.01, order = A42()).Tleaf # With units c4_q = PH.C4Q() @@ -115,8 +115,8 @@ res_q = PH.energybalance(Tl_q, Tl_q = PH.solve_energy_balance(c4_q, gb = gb_q, opt = opt, PAR = PAR_q, NIR = NIR_q, ws = ws_q, RH = RH_q, Tair = Ta_q, Ca = Ca_q, P = P_q, O2 = O2_q, xatol = 0.01, - order = Order2()).Tleaf -@test Tl_q ≈ Tl * K + order = A42()).Tleaf +@test abs(Tl_q - Tl * K) < 0.02K # With effect of angles (horizontal) gb_q = PH.gbAngleQ() @@ -137,7 +137,7 @@ res_q = PH.energybalance(Tl_q, Tl_q = PH.solve_energy_balance(c4_q, gb = gb_q, opt = opt, PAR = PAR_q, NIR = NIR_q, ws = ws_q, RH = RH_q, Tair = Ta_q, Ca = Ca_q, P = P_q, O2 = O2_q, xatol = 0.01, - order = Order2()).Tleaf -@test Tl_q ≈ Tl * K + order = A42()).Tleaf +@test abs(Tl_q - Tl * K) < 0.02K end