Skip to content

Commit

Permalink
second
Browse files Browse the repository at this point in the history
  • Loading branch information
@kongdd
kongdd committed Oct 12, 2024
1 parent a4a759e commit dc8584b
Show file tree
Hide file tree
Showing 3 changed files with 70 additions and 147 deletions.
209 changes: 69 additions & 140 deletions src/photosynthesis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 * Γ)
Expand All @@ -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`
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down
4 changes: 1 addition & 3 deletions test/test-beps_main.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
using Test
using BEPS
using BEPS, DataFrames, Test
using BEPS: path_proj
using DataFrames

function nanmax(x)
x = collect(x)
Expand Down
4 changes: 0 additions & 4 deletions test/test-readparam.jl
View file

This file was deleted.

0 comments on commit dc8584b

Please sign in to comment.