diff --git a/.gitignore b/.gitignore
index bb67a3e..379015d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,7 +23,9 @@ docs/reference
 *.html
 
 # Figures
-*.pdf
+# *.pdf
 .ipynb_checkpoints
 
 __pycache__
+data
+examples
diff --git a/src/photosynthesis.jl b/src/photosynthesis.jl
index 7f180da..b836e2f 100644
--- a/src/photosynthesis.jl
+++ b/src/photosynthesis.jl
@@ -11,26 +11,22 @@ Forests are hypostomatous. Hence, we don't divide the total resistance by 2
 since transfer is going on only one side of a leaf.
 
 # Arguments
-- `ea`     : [kPa]
-- `gb_w`   : leaf laminar boundary layer conductance to H2O, [s m-1]
-- `Vcmax25`: the maximum rate of carboxylation of Rubisco at 25 deg C, [umol m-2
-  s-1]
+- `ea`      : [kPa]
+- `gb_w`    : leaf laminar boundary layer conductance to H2O, [s m-1]
+- `Vcmax25` : the maximum rate of carboxylation of Rubisco at 25℃, [umol m-2 s-1]
 
-- `cii`   : intercellular CO2 concentration (ppm)
-- `g0_h2o`: the minimum stomatal conductance to H2O,      [umol m-2 s-1]
-- `g1_h2o`: the slope of the stomatal conductance to H2O, [unitless]
+- `cii`     : intercellular CO2 concentration (ppm)
+- `g0_h2o`  : the minimum stomatal conductance to H2O,      [umol m-2 s-1]
+- `g1_h2o`  : the slope of the stomatal conductance to H2O, [unitless]
 
-- `LH_leaf`: latent heat of vaporization of water at the leaf temperature, [W
-  m-2]
+- `LH_leaf` : latent heat of vaporization of water at the leaf temperature, [W m-2]
+- `ca`      : atmospheric co2 concentration (ppm)
 
 # Intermediate variables
-- `gb_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)
+- `rh_leaf`    : relative humidity at leaf surface (0-1)
+- `gs_c_mol`   : stomatal conductance to CO2 (umol m-2 s-1)
+- `gs_w_mol`   : stomatal conductance to h2o (umol 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)
@@ -43,14 +39,12 @@ since transfer is going on only one side of a leaf.
 - `Jₓ`         : the flux of electrons through the thylakoid membrane (umol m-2
   s-1)
 """
-@fastmath function photosynthesis_jl(T_leaf_p::Cdouble, Rsn_leaf::Cdouble, ea::Cdouble,
-  gb_w::Cdouble, Vcmax25::Cdouble,
-  β_soil::Cdouble, g0_w::Cdouble, g1_w::Cdouble,
-  cii::Cdouble,
-  T_leaf::Cdouble, LH_leaf::Cdouble)
-
-  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
+@fastmath function photosynthesis_jl(T_leaf_p::T, Rsn_leaf::T, ea::T,
+  gb_w::T, Vcmax25::T,
+  β_soil::T, g0_w::T, g1_w::T, cii::T,
+  T_leaf::T, LH_leaf::T, ca::T=CO2_air) where {T<:Cdouble}
+  # 
+  PPFD::T = 4.55 * 0.5 * Rsn_leaf # incident photosynthetic photon flux density (PPFD) umol m-2 s-1
   (2PPFD < 1) && (PPFD = 0.0)
   T_leaf_K = T_leaf + 273.13
 
@@ -63,20 +57,20 @@ since transfer is going on only one side of a leaf.
   g0_c = g0_w / 1.6
   g1_c = g1_w / 1.6
 
-  RH_leaf = SFC_VPD(T_leaf_K, LH_leaf, λ, rᵥ, ρₐ)
+  rh_leaf = SFC_VPD(T_leaf_K, LH_leaf, λ, rᵥ, ρₐ)
   tprime25 = T_leaf_K - TK25  # temperature difference
 
-  Kc = TEMP_FUNC(kc25, ekc, tprime25, TK25, T_leaf_K)
-  Ko = TEMP_FUNC(ko25, eko, tprime25, TK25, T_leaf_K)
-  tau = TEMP_FUNC(tau25, ektau, tprime25, TK25, T_leaf_K)
+  Kc = kc25 * fTᵥ(ekc, tprime25, TK25, T_leaf_K)
+  Ko = ko25 * fTᵥ(eko, tprime25, TK25, T_leaf_K)
+  tau = tau25 * fTᵥ(ektau, tprime25, TK25, T_leaf_K)
+  Γ = 0.5 * o2 / tau * 1000.0  # [mmol mol-1] -> umol mol-1
 
   K = Kc * (1.0 + o2 / Ko)
-  Γ = 0.5 * o2 / tau * 1000.0  # umol mol-1
   Rd25 = Vcmax25 * 0.004657    # leaf dark respiration (umol m-2 s-1)
 
   # Bin Chen: Reduce respiration by 40% in light according to Amthor
   (2PPFD > 10) && (Rd25 *= 0.4)
-  Rd = TEMP_FUNC(Rd25, erd, tprime25, TK25, T_leaf_K)
+  Rd = Rd25 * fTᵥ(erd, tprime25, TK25, T_leaf_K)
 
   #	jmopt = 29.1 + 1.64*Vcmax25; Chen 1999, Eq. 7
   jmopt = 2.39 * Vcmax25 - 14.2
@@ -84,34 +78,31 @@ since transfer is going on only one side of a leaf.
   Jmax = TBOLTZ(jmopt, ejm, toptjm, T_leaf_K)    # Apply temperature correction to JMAX
   Vcmax = TBOLTZ(Vcmax25, evc, toptvc, T_leaf_K)  # Apply temperature correction to vcmax
 
-  #  * Test for the minimum of Wc and Wj.  Both have the form: `W = (a ci - a d)/(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)
+  # Farquhar and von Cammerer (1981)
   # /*if (jmax > 0) Jₓ = qalpha * iphoton / sqrt(1. +(qalpha2 * iphoton * iphoton / (jmax * jmax)));
   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 * Γ)
-  wc = Vcmax * (cii - Γ) / (cii + K)
+  Wj = Jₓ * (cii - Γ) / (4.0 * cii + 8.0 * Γ)
+  Wc = Vcmax * (cii - Γ) / (cii + K)
 
-  if (wj < wc)
-    # A_guess = wj      # Harley and Farquhar type model for Wj
-    a = Jₓ
+  # Both have the form: `Ag = A + Rd = (a ci - a d)/(e ci + b)`
+  if (Wj < Wc)
+    a = Jₓ    # J limited
     b = 8.0 * Γ
     e = 4.0
   else
-    # A_guess = wc
-    a = Vcmax # x1
+    a = Vcmax # VCmax limited
     b = K
     e = 1.0
   end
   d = Γ
   An = 0.0
 
-  if !(wj <= Rd || wc <= Rd)
-    # g_s =  g0 + g1 * RH_leaf * β_soil * A_g
-    # g_s =  g0 + _c * A_g, (c = g1 * RH_leaf * β_soil)
-    c = g1_c * RH_leaf * β_soil
+  if !(Wj <= Rd || Wc <= Rd)
+    # g_s =  g0 + g1 * rh_leaf * β_soil * A_g
+    # g_s =  g0 + _c * A_g, (c = g1 * rh_leaf * β_soil)
+    c = g1_c * rh_leaf * β_soil
     α = 1.0 + (g0_c / gb_c_mol) - c
     β = ca * (gb_c_mol * c - 2.0 * g0_c - gb_c_mol)
     γ = ca * ca * gb_c_mol * g0_c
@@ -130,7 +121,7 @@ since transfer is going on only one side of a leaf.
   An = max(0.0, An)
   cs = ca - An / gb_c_mol
 
-  gs_w_mol = (β_soil * g1_w * RH_leaf * An / cs) + g0_w  # mol m-2 s-1
+  gs_w_mol = (β_soil * g1_w * rh_leaf * An / cs) + g0_w  # mol m-2 s-1
   gs_c_mol = gs_w_mol / 1.6
 
   ci = cs - An / gs_c_mol
@@ -138,9 +129,12 @@ since transfer is going on only one side of a leaf.
   return gs_w, An, ci
 end
 
+@fastmath function fTᵥ(eact::T, tprime::T, tref::T, t_lk::T)::T where {T<:Real}
+  exp(tprime * eact / (tref * rugc * t_lk))
+end
 
 """
-If `wj` or `wc` are less than Rd then A would probably be less than 0. This
+If `Wj` or `Wc` are less than Rd then A would probably be less than 0. This
 would yield a negative stomatal conductance. 
 
 In this case, assume `gs` equals the cuticular value `g0`. This assumptions
@@ -181,7 +175,7 @@ Rank roots #1, #2 and #3 according to the minimum, intermediate and maximum valu
 ```math
 c_s = c_a - A * 1/ g_b
 c_i = c_a - A * (1/g_b + 1/g_s)
-g_s = g_0 + g_1 RH A / c_s
+g_s = g_0 + g_1 rh A / c_s
 Ag = a(c_i - Gamma) / (e c_i + b})
 A = Ag - Rd
 ```
diff --git a/src/photosynthesis_helper.jl b/src/photosynthesis_helper.jl
index cf3c59a..27101ee 100644
--- a/src/photosynthesis_helper.jl
+++ b/src/photosynthesis_helper.jl
@@ -107,15 +107,11 @@ end
   return rh
 end
 
-@fastmath function TEMP_FUNC(rate::T, eact::T, tprime::T, tref::T, t_lk::T)::T where {T<:Real}
-  rate * exp(tprime * eact / (tref * rugc * t_lk))
-end
-
 function LAMBDA(TK::T)::T where {T<:Real}
   y = 3149000.0 - 2370.0 * TK # J kg-1
   # add heat of fusion for melting ice
   if TK < 273.0
-    y += 333.0
+    y += 333.0 # TODO: unit error, `y += 333_000.0`
   end
   return y
 end
@@ -140,6 +136,5 @@ end
   prodt::T = rugc * topt * tl
   numm::T = hkin * exp(eakin * dtlopt / prodt)
   denom::T = hkin - eakin * (1.0 - exp(hkin * dtlopt / prodt))
-
   return rate * numm / denom
 end