fixing bug, adding test

src/integrated/soil_canopy_model.jl

@@ -173,12 +173,12 @@ function SoilCanopyModel{FT}(;
     end
 
     co2_prognostic_soil = Soil.Biogeochemistry.PrognosticMet(soil.parameters)
-    soil_co2_drivers = Soil.Biogeochemistry.SoilDrivers(
+    soilco2_drivers = Soil.Biogeochemistry.SoilDrivers(
         co2_prognostic_soil,
         Soil.Biogeochemistry.PrescribedSOC{FT}(soil_organic_carbon),
         atmos,
     )
-    soilco2 = soil_co2_model(; soilco2_args..., drivers = soil_co2_drivers)
+    soilco2 = soilco2_type(; soilco2_args..., drivers = soilco2_drivers)
 
     return SoilCanopyModel{FT, typeof(soilco2), typeof(soil), typeof(canopy)}(
         soilco2,

src/standalone/Soil/soil_hydrology_parameterizations.jl

@@ -69,11 +69,13 @@ end
 """
      approximate_ψ_S_slope(cm::vanGenuchten)
 
-An estimate of the slop of the logψ-logS curve, following
-Lehmann, Assouline, and Or (2008).
+An estimate of the slope of the absolute value of the logψ-logS curve.
+Following Lehmann, Assouline, and Or (2008), we linearize the ψ(S) curve about the inflection point (where d²ψ/dS² = 0, at S = (1+m)^(-m)). 
 """
 function approximate_ψ_S_slope(cm::vanGenuchten)
-    return cm.n
+    m = cm.m
+    n = cm.n
+    return (1 + m) / (n * m * m)
 end
 
 

test/integrated/soil_energy_hydrology_biogeochemistry.jl

@@ -4,7 +4,7 @@ import ClimaComms
 using ClimaCore
 import ClimaParams
 using ClimaLand
-using ClimaLand.Domains: Column
+using ClimaLand.Domains: Column, HybridBox
 using ClimaLand.Soil
 using ClimaLand.Soil.Biogeochemistry
 using Dates
@@ -109,6 +109,8 @@ for FT in (Float32, Float64)
             soil_args = soil_args,
             soilco2_args = soilco2_args,
         )
+        @test model.soilco2.drivers.met.ν == model.soil.parameters.ν
+        @test model.soilco2.drivers.met.ν isa ClimaLand.PrognosticMet
         Y, p, coords = initialize(model)
         @test propertynames(p.drivers) ==
               (:P_liq, :P_snow, :T, :P, :u, :q, :c_co2, :thermal_state)
@@ -173,4 +175,92 @@ for FT in (Float32, Float64)
             z,
         )
     end
+    @testset "PrognosticMet, FT = $FT" begin
+        earth_param_set = LP.LandParameters(FT)
+        zmax = FT(0)
+        zmin = FT(-1)
+        nelems = 10
+        xmin = ymin = zmin
+        xmax = ymax = zmax
+        col = Column(; zlim = (zmin, zmax), nelements = nelems)
+        box = HybridBox(;
+            xlim = (xmin, xmax),
+            ylim = (ymin, ymax),
+            zlim = (zmin, zmax),
+            nelements = (nelems, nelems, nelems),
+            npolynomial = 2,
+        )
+
+        ν = FT(0.556)
+        K_sat = FT(0.0443 / 3600 / 100)
+        S_s = FT(1e-3)
+        vg_n = FT(2.0)
+        vg_α = FT(2.6)
+        hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
+        θ_r = FT(0.1)
+        ν_ss_om = FT(0.0)
+        ν_ss_quartz = FT(1.0)
+        ν_ss_gravel = FT(0.0)
+
+
+        soil_ps_col = Soil.EnergyHydrologyParameters(
+            FT;
+            ν,
+            ν_ss_om,
+            ν_ss_quartz,
+            ν_ss_gravel,
+            K_sat,
+            hydrology_cm,
+            S_s,
+            θ_r,
+        )
+
+        zero_field = ClimaCore.Fields.zeros(box.space.subsurface)
+        vg_fields_to_hcm_field(α::FT, n::FT) where {FT} =
+            ClimaLand.Soil.vanGenuchten{FT}(;
+                @NamedTuple{α::FT, n::FT}((α, n))...,
+            )
+        hydrology_cm_field =
+            vg_fields_to_hcm_field.(vg_α .+ zero_field, vg_n .+ zero_field)
+        ν_field = ν .+ zero_field
+        ν_ss_om_field = ν_ss_om .+ zero_field
+        ν_ss_quartz_field = ν_ss_quartz .+ zero_field
+        ν_ss_gravel_field = ν_ss_gravel .+ zero_field
+        K_sat_field = K_sat .+ zero_field
+        S_s_field = S_s .+ zero_field
+        θ_r_field = θ_r .+ zero_field
+
+        soil_ps_box = Soil.EnergyHydrologyParameters(
+            FT;
+            ν = ν_field,
+            ν_ss_om = ν_ss_om_field,
+            ν_ss_quartz = ν_ss_quartz_field,
+            ν_ss_gravel = ν_ss_gravel_field,
+            K_sat = K_sat_field,
+            hydrology_cm = hydrology_cm_field,
+            S_s = S_s_field,
+            θ_r = θ_r_field,
+        )
+
+        met_col = ClimaLand.PrognosticMet(soil_ps_col)
+        met_box = ClimaLand.PrognosticMet(soil_ps_box)
+        @test met_box.ν == soil_ps_box.ν
+        @test met_col.ν == soil_ps_col.ν
+        @test met_box.θ_a100 == @. Soil.inverse_matric_potential(
+            soil_ps_box.hydrology_cm,
+            -FT(1),
+        ) * (soil_ps_box.ν - soil_ps_box.θ_r) + soil_ps_box.θ_r
+        @test met_col.θ_a100 ==
+              Soil.inverse_matric_potential(soil_ps_col.hydrology_cm, -FT(1)) *
+              (soil_ps_col.ν - soil_ps_col.θ_r) + soil_ps_col.θ_r
+        @test met_box.b ==
+              @. Soil.approximate_ψ_S_slope(soil_ps_box.hydrology_cm)
+        @test met_col.b == Soil.approximate_ψ_S_slope(soil_ps_col.hydrology_cm)
+        vg_m = 1 - 1 / vg_n
+        @test Soil.approximate_ψ_S_slope(soil_ps_col.hydrology_cm) ==
+              (1 + vg_m) / (vg_n * vg_m^2)
+        @test Soil.approximate_ψ_S_slope(
+            BrooksCorey{FT}(; c = FT(1.0), ψb = FT(1.0)),
+        ) == FT(1)
+    end
 end