Minor fixes

- Removed the PModelDrivers struct due to redundancy and broadcasting issues. Now we take in the drivers as separate arguments
- Clarified units of Kc, Ko, and Gammastar
- Removed T, P-dependent density of water. Now we only use constant density
- Fixed accidental deletions during rebase

Author: ychnli

src/standalone/Vegetation/Canopy.jl

@@ -46,6 +46,7 @@ include("./canopy_energy.jl")
 include("./canopy_parameterizations.jl")
 using Dates
 include("./autotrophic_respiration.jl")
+include("./spatially_varying_parameters.jl")
 
 """
     SharedCanopyParameters{FT <: AbstractFloat, PSE}

src/standalone/Vegetation/canopy_parameterizations.jl

@@ -845,7 +845,7 @@ end
           R::FT) where {FT}
 
 Computes the Michaelis-Menten coefficient for CO2 (`Kc`),
-in units of mol/mol,
+in units of mol/mol or Pa (depending on the units of `Kc25`),
 as a function of its value at 25 °C (`Kc25`),
 a constant (`ΔHkc`), a standard temperature (`To`),
 the unversal gas constant (`R`), and the temperature (`T`).
@@ -866,7 +866,7 @@ end
           R::FT) where {FT}
 
 Computes the Michaelis-Menten coefficient for O2 (`Ko`),
-in units of mol/mol,
+in units of mol/mol or Pa (depending on the units of `Ko25`),
 as a function of its value at 25 °C (`Ko25`),
 a constant (`ΔHko`), a standard temperature (`To`),
 the universal gas constant (`R`), and the temperature (`T`).
@@ -1141,59 +1141,6 @@ function intrinsic_quantum_yield(
 end
 
 
-"""
-    density_h20(
-        T::FT, 
-        p::FT
-    ) where {FT}
-
-    Computes the density of water in kg/m^3 given temperature T (K) and pressure p (Pa)
-    according to F.H. Fisher and O.E Dial, Jr. (1975) Equation of state of pure water and 
-    sea water, Tech. Rept., Marine Physical Laboratory, San Diego, CA.
-
-    (can consider removing the higher order terms?) 
-"""
-function density_h2o(T::FT, p::FT) where {FT}
-    # Convert temperature to Celsius
-    T = T - FT(273.15)
-    # Convert pressure to bar (1 bar = 1e5 Pa)
-    pbar = FT(1e-5) * p
-
-    # λ(T): bar cm3/g
-    λ =
-        FT(1788.316) + FT(21.55053) * T - FT(0.4695911) * T^2 +
-        FT(3.096363e-3) * T^3 - FT(7.341182e-6) * T^4
-
-    # p0(T): bar
-    p0 =
-        FT(5918.499) + FT(58.05267) * T - FT(1.1253317) * T^2 +
-        FT(6.6123869e-3) * T^3 - FT(1.4661625e-5) * T^4
-
-    # v_inf(T): cm3/g
-    T_powers = (T .^ (0:9))  # T^0 to T^9
-    coeffs_vinf =
-        FT.([
-            0.6980547,
-            -7.435626e-4,
-            3.704258e-5,
-            -6.315724e-7,
-            9.829576e-9,
-            -1.197269e-10,
-            1.005461e-12,
-            -5.437898e-15,
-            1.69946e-17,
-            -2.295063e-20,
-        ])
-    v_inf = dot(coeffs_vinf, T_powers)
-
-    # Specific volume [cm3/g]
-    v = v_inf + λ / (p0 + pbar)
-
-    # Convert to density [kg/m3]
-    return FT(1e3) / v
-end
-
-
 """
     viscosity_h20(
         T::FT, 
@@ -1206,16 +1153,15 @@ end
 
     Can consider simplifying if this level of precision is not needed
 """
-function viscosity_h2o(T::FT, p::FT, constant_density::Bool) where {FT}
+function viscosity_h2o(T::FT, p::FT) where {FT}
     # Reference constants
     tk_ast = FT(647.096)    # K
     ρ_ast = FT(322.0)      # kg/m^3
     μ_ast = FT(1e-6)       # Pa s
 
     # Get density of water [kg/m^3]
     earth_param_set = LP.LandParameters(FT)
-    ρ0 = LP.ρ_cloud_liq(earth_param_set)
-    constant_density ? ρ = ρ0 : ρ = density_h2o(T, p)
+    ρ = LP.ρ_cloud_liq(earth_param_set)
 
     # Dimensionless variables
     tbar = T / tk_ast
@@ -1278,10 +1224,9 @@ end
 function compute_viscosity_ratio(
     T::FT,
     p::FT,
-    constant_density::Bool,
 ) where {FT}
-    η25 = viscosity_h2o(FT(298.15), FT(101325.0), constant_density)
-    ηstar = viscosity_h2o(T, p, constant_density) / η25
+    η25 = viscosity_h2o(FT(298.15), FT(101325.0))
+    ηstar = viscosity_h2o(T, p) / η25
     return FT(ηstar)
 end
 
@@ -1434,7 +1379,6 @@ end
     is constant. cstar is defined as 4c, a free parameter. Wang etal (2017) derive cstar = 0.412
     at STP and using Vcmax/Jmax = 1.88. 
 """
-# TODO: test if cstar should be made a free parameter? 
 function compute_mj_with_jmax_limitation(mj::FT, cstar::FT) where {FT}
     return FT(mj * sqrt(1 - (cstar / mj)^(FT(2 / 3))))
 end

src/standalone/Vegetation/pmodel.jl

@@ -1,5 +1,7 @@
-export PModelParameters,
-    PModelDrivers, PModelConstants, compute_full_pmodel_outputs, PModel
+export PModelParameters, 
+    PModelConstants, 
+    compute_full_pmodel_outputs, 
+    PModel
 
 """
     PModelParameters{FT<:AbstractFloat}
@@ -29,42 +31,26 @@ Base.@kwdef struct PModelParameters{FT <: AbstractFloat}
     ϕa2::FT
 end
 
-"""
-    PModelDrivers{FT<:AbstractFloat}
-
-The required drivers for P-model (Stocker et al. 2020). Drivers are defined as
-external variables used to compute the optimal photosynthetic capacities. 
-$(DocStringExtensions.FIELDS)
-"""
-Base.@kwdef struct PModelDrivers{FT <: AbstractFloat}
-    "Canopy temperature (K)"
-    T_canopy::FT
-    "Absorbed PAR in moles of photons (mol m^-2 s^-1)"
-    I_abs::FT
-    "Ambient CO2 partial pressure (Pa)"
-    ca::FT
-    "Ambient air pressure (Pa)"
-    P_air::FT
-    "Vapor pressure deficit (Pa)"
-    VPD::FT
-    """Soil moisture stress factor (unitless)"""
-    βm::FT
-end
-
 
 """
     PModelConstants{FT<:AbstractFloat}
 
 The required constants for P-model (Stocker et al. 2020). These are physical
 or biochemical constants that are not expected to change with time or space.
+
+IMPORTANT NOTE: the dimensions of the Michaelis-Menten parameters Kc and Ko and compensation
+point Γstar correspond to gaseous concentrations, so they can either be expressed as a mixing
+ratio (mol mol^-1) or as a partial pressure (Pa) of 1 atm. These are equivalent up to a constant
+factor of P_atm = 101325 Pa. As long as you are consistent with which convention you use,
+the model will work correctly. By default, these constants are set using the Pa convention. 
 $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef struct PModelConstants{FT}
     """Gas constant (J mol^-1 K^-1)"""
     R::FT
-    """Michaelis-Menten parameter for carboxylation at 25°C (μmol mol^-1)"""
+    """Michaelis-Menten parameter for carboxylation at 25°C (mol mol^-1 or Pa)"""
     Kc25::FT
-    """Michaelis-Menten parameter for oxygenation at 25°C (μmol mol^-1)"""
+    """Michaelis-Menten parameter for oxygenation at 25°C (mol mol^-1 or Pa)"""
     Ko25::FT
     """Reference temperature equal to 25˚C (K)"""
     To::FT
@@ -76,7 +62,7 @@ Base.@kwdef struct PModelConstants{FT}
     Drel::FT
     """Effective energy of activation for Γstar (J mol^-1)"""
     ΔHΓstar::FT
-    """Γstar at 25 °C (Pa)"""
+    """Γstar at 25 °C (mol mol^-1 or Pa)"""
     Γstar25::FT
     """Effective energy of activation for Vcmax (J mol^-1)"""
     Ha_Vcmax::FT
@@ -96,7 +82,7 @@ Base.@kwdef struct PModelConstants{FT}
     bS_Jmax::FT
     """Molar mass of carbon (kg mol^-1)"""
     Mc::FT
-    """Intercellular O2 mixing ratio (unitless)"""
+    """Intercellular O2 mixing ratio (mol mol^-1)"""
     oi::FT
     """First order coefficient for temp-dependent Rd (K^-1)"""
     aRd::FT
@@ -107,11 +93,9 @@ Base.@kwdef struct PModelConstants{FT}
 end
 
 Base.eltype(::PModelParameters{FT}) where {FT} = FT
-Base.eltype(::PModelDrivers{FT}) where {FT} = FT
 Base.eltype(::PModelConstants{FT}) where {FT} = FT
 
 Base.broadcastable(m::PModelParameters) = tuple(m)
-Base.broadcastable(m::PModelDrivers) = tuple(m)
 Base.broadcastable(m::PModelConstants) = tuple(m)
 
 """
@@ -178,10 +162,15 @@ ClimaLand.auxiliary_domain_names(::PModel) =
 
 """
     compute_full_pmodel_outputs(
-        parameters::PModelParameters, 
-        drivers::PModelDrivers, 
-        constants::PModelConstants
-    )
+        parameters::PModelParameters{FT},
+        constants::PModelConstants{FT},
+        T_canopy::FT,
+        I_abs::FT,
+        ca::FT,
+        P_air::FT,
+        VPD::FT,
+        βm::FT
+    ) where {FT}
 
 Performs the P-model computations as defined in Stocker et al. (2020) 
 and returns a dictionary of outputs. See https://github.com/geco-bern/rpmodel
@@ -207,15 +196,17 @@ Output name      Description (units)
 """
 function compute_full_pmodel_outputs(
     parameters::PModelParameters{FT},
-    drivers::PModelDrivers{FT},
     constants::PModelConstants{FT},
+    T_canopy::FT,
+    I_abs::FT,
+    ca::FT,
+    P_air::FT,
+    VPD::FT,
+    βm::FT
 ) where {FT}
     # Unpack parameters
     (; cstar, β, ϕc, ϕ0, ϕa0, ϕa1, ϕa2) = parameters
 
-    # Unpack drivers
-    (; T_canopy, I_abs, ca, P_air, VPD, βm) = drivers
-
     # Unpack constants
     (;
         R,
@@ -246,7 +237,7 @@ function compute_full_pmodel_outputs(
     ϕ0 = isnan(ϕ0) ? intrinsic_quantum_yield(T_canopy, ϕc, ϕa0, ϕa1, ϕa2) : ϕ0
 
     Γstar = co2_compensation_p(T_canopy, To, P_air, R, ΔHΓstar, Γstar25)
-    ηstar = compute_viscosity_ratio(T_canopy, P_air, true)
+    ηstar = compute_viscosity_ratio(T_canopy, P_air)
     Kmm = compute_Kmm(T_canopy, P_air, Kc25, Ko25, ΔHkc, ΔHko, To, R, oi)
     χ, ξ, mj, mc = optimal_co2_ratio_c3(Kmm, Γstar, ηstar, ca, VPD, β, Drel)
     ci = χ * ca

test/Artifacts.jl

@@ -54,7 +54,8 @@ end
 """
     pmodel_unittests_path(; context = nothing)
 
-Returns the filepaths for input-output pairs for the P-model. 
+Returns the filepaths for input-output pairs for the P-model. This data is generated from the R
+implementation of the P-model and is used to test the ClimaLand P-model.
     
 Experiment details are in `test/standalone/Vegetation/test_pmodel.jl`.
 """

test/runtests.jl

@@ -83,6 +83,9 @@ end
 @safetestset "Soil integrated water and energy content" begin
     include("standalone/Soil/conservation.jl")
 end
+@safetestset "Soil spatial parameters" begin
+    include("standalone/Soil/spatial_parameters.jl")
+end
 
 # Standalone Surface Water model tests
 @safetestset "Pond module tests" begin
@@ -105,6 +108,9 @@ end
 @safetestset "Canopy integrated water and energy content" begin
     include("standalone/Vegetation/conservation.jl")
 end
+@safetestset "Canopy spatial parameters" begin
+    include("standalone/Vegetation/spatial_parameters.jl")
+end
 @safetestset "P model tests" begin
     include("standalone/Vegetation/test_pmodel.jl")
 end

test/standalone/Vegetation/test_pmodel.jl

@@ -190,14 +190,23 @@ end
                 verbose &&
                     println("Running test case: $testcase_name with FT = $FT")
 
-                # Create constants, drivers, and parameters for the current FT
+                # prepare constants, parameters, and drivers for the current FT
                 constants = PModelConstants(FT)
-                drivers = PModelDrivers(inputs, FT)
                 parameters = PModelParameters(inputs, FT)
 
+                T_canopy = FT(inputs["tc"] + 273.15)  # Convert from Celsius to Kelvin
+                VPD = FT(inputs["vpd"])
+                ca = FT(inputs["co2"]) * FT(1e-6) * FT(101325.0)  # Convert ppm to Pa
+                P_air = FT(get(inputs, "patm", 101325.0))
+                I_abs = FT(inputs["fapar"]) * FT(inputs["ppfd"])
+                βm =
+                    Bool(inputs["do_soilmstress"]) ?
+                    quadratic_soil_moisture_stress(FT(inputs["soilm"])) : FT(1.0)
+
                 # Run the model
                 outputs =
-                    compute_full_pmodel_outputs(parameters, drivers, constants)
+                    compute_full_pmodel_outputs(parameters, constants, 
+                        T_canopy, I_abs, ca, P_air, VPD, βm)
 
                 # Compare each output field
                 for key in keys(ref_outputs_typed)