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..3a1ee2f 100644
--- a/src/Photosynthesis/Components.jl
+++ b/src/Photosynthesis/Components.jl
@@ -42,7 +42,7 @@ 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 <: Real} = T(610.78) # Pa
 es0(::Type{T}) where {T <: Quantity} = 610.78Pa
 
 # 0 umol/m2/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/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 85ced00..dd78d37 100644
--- a/src/Photosynthesis/FvCB/C3.jl
+++ b/src/Photosynthesis/FvCB/C3.jl
@@ -7,13 +7,13 @@ 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, 
+        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)
 
 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/Pa)
 """
 Base.@kwdef mutable struct C3{T <: Real} <: C3Type
     # Rubisco CO2/O2 specificity
@@ -95,16 +95,16 @@ 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.14e-3 # Empirical parameter in gs formula (1/Pa)
 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, 
+         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)
 
 Data structure to store all the parameters for the C3 photosynthesis model using
@@ -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.14e-3 / Pa # 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
 
 """
@@ -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..518d247 100644
--- a/src/Photosynthesis/FvCB/C4.jl
+++ b/src/Photosynthesis/FvCB/C4.jl
@@ -7,10 +7,10 @@ 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, 
+    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)
 
 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,15 +81,15 @@ 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.15e-3 # Empirical parameter in gs formula (1/Pa)
 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, 
+       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)
 
 Data structure to store all the parameters for the C4 photosynthesis model using
@@ -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.15e-3 / Pa # 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;
@@ -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..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.0022) < 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_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.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,
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/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
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