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
-