Add a AbstractSoilMoistureStressModel, 3 methods

This commits adds a portion of Yuchen Li PR #1330,
a new Canopy component: AbstractSoilMoistureStressModel,
with 3 methods: the existing one (Tuzet) which was previously
in PlantHydraulics, a Piecewise method using dry to wet soil moisture,
and no moisture stress.

Co-authored-by: AlexisRenchon <[email protected]>
Co-authored-by: ychnli <[email protected]>

[skip ci]

Author: AlexisRenchon

experiments/integrated/fluxnet/ozark_pmodel.jl

@@ -112,7 +112,8 @@ dt = Float64(450) # 7.5 minutes
     atmos_h,
     start_date,
     earth_param_set,
-    FT,
+    FT;
+    construct_prescribed_soil = false,
 )
 
 
@@ -179,6 +180,17 @@ conductance = PModelConductance{FT}()
 # Set up photosynthesis
 photosynthesis = PModel{FT}()
 
+# Set up soil moisture stress
+soil_moisture_stress_params = PiecewiseMoistureStressParameters(
+                                                                FT;
+                                                                θ_c = FT(0.60),
+                                                                θ_w = FT(0.10),
+                                                                c = FT(1.0),
+                                                                β0 = FT(1.0),
+                                                               )
+soil_moisture_stress =
+PiecewiseMoistureStressModel{FT}(soil_moisture_stress_params)
+
 # Set up plant hydraulics
 # Read in LAI from MODIS data
 surface_space = land_domain.space.surface;
@@ -222,6 +234,7 @@ canopy = Canopy.CanopyModel{FT}(
     radiative_transfer,
     photosynthesis,
     conductance,
+    soil_moisture_stress,
     hydraulics,
     energy,
 )
@@ -265,8 +278,9 @@ simulation = LandSimulation(
     set_ic!,
     updateat,
     diagnostics = diags,
-);
-@time solve!(simulation)
+)
+
+@time sol = solve!(simulation)
 
 comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset)
 savedir = joinpath(

src/diagnostics/land_compute_methods.jl

@@ -115,7 +115,14 @@ function compute_stomatal_conductance!(
 end
 
 function compute_stomatal_conductance!(out, Y, p, t, land_model::CanopyModel)
-    compute_stomatal_conductance!(out, Y, p, t, land_model, canopy.conductance)
+        compute_stomatal_conductance!(
+        out,
+        Y,
+        p,
+        t,
+        land_model,
+        land_model.conductance,
+    )
 end
 
 function compute_stomatal_conductance!(
@@ -273,61 +280,6 @@ function compute_leaf_water_potential!(out, Y, p, t, land_model::CanopyModel)
     end
 end
 
-function compute_moisture_stress_factor!(
-    out,
-    Y,
-    p,
-    t,
-    land_model::Union{SoilCanopyModel, LandModel},
-)
-    canopy = land_model.canopy
-    hydraulics = canopy.hydraulics
-    n_stem = hydraulics.n_stem
-    n_leaf = hydraulics.n_leaf
-    n = n_stem + n_leaf
-
-    earth_param_set = canopy.parameters.earth_param_set
-    grav = LP.grav(earth_param_set)
-    ρ_l = LP.ρ_cloud_liq(earth_param_set)
-    (; sc, pc) = canopy.photosynthesis.parameters
-    ψ = p.canopy.hydraulics.ψ
-    if isnothing(out)
-        out = zeros(land_model.canopy.domain.space.surface) # Allocates
-        fill!(field_values(out), NaN)
-        @. out = moisture_stress(ψ.:($$n) * ρ_l * grav, sc, pc)
-        return out
-    else
-        @. out = moisture_stress(ψ.:($$n) * ρ_l * grav, sc, pc)
-    end
-end
-function compute_moisture_stress_factor!(
-    out,
-    Y,
-    p,
-    t,
-    land_model::Union{CanopyModel},
-)
-    canopy = land_model
-    hydraulics = canopy.hydraulics
-    n_stem = hydraulics.n_stem
-    n_leaf = hydraulics.n_leaf
-    n = n_stem + n_leaf
-
-    earth_param_set = canopy.parameters.earth_param_set
-    grav = LP.grav(earth_param_set)
-    ρ_l = LP.ρ_cloud_liq(earth_param_set)
-    (; sc, pc) = canopy.photosynthesis.parameters
-    ψ = p.canopy.hydraulics.ψ
-    if isnothing(out)
-        out = zeros(land_model.canopy.domain.space.surface) # Allocates
-        fill!(field_values(out), NaN)
-        @. out = moisture_stress(ψ.:($$n) * ρ_l * grav, sc, pc)
-        return out
-    else
-        @. out = moisture_stress(ψ.:($$n) * ρ_l * grav, sc, pc)
-    end
-end
-
 # @diagnostic_compute "flux_per_ground_area" Union{SoilCanopyModel, LandModel} p.canopy.hydraulics.fa # return a Tuple
 @diagnostic_compute "root_flux_per_ground_area" Union{
     SoilCanopyModel,
@@ -340,6 +292,13 @@ end
     CanopyModel,
 } p.canopy.hydraulics.area_index.leaf
 
+# Canopy - Soil moisture stress
+@diagnostic_compute "moisture_stress_factor" Union{
+    SoilCanopyModel,
+    LandModel,
+    CanopyModel,
+} p.canopy.soil_moisture_stress.βm
+
 # Canopy - Hydraulics
 @diagnostic_compute "root_area_index" Union{
     SoilCanopyModel,
@@ -574,8 +533,8 @@ end # Convert from kg C to mol CO2.
 
 ## Other ##
 @diagnostic_compute "sw_albedo" Union{SoilCanopyModel, LandModel} p.α_sfc
-@diagnostic_compute "lw_up" Union{SoilCanopyModel, LandModel} p.LW_u
-@diagnostic_compute "sw_up" Union{SoilCanopyModel, LandModel} p.SW_u
+@diagnostic_compute "lw_up" Union{SoilCanopyModel, LandModel, CanopyModel} p.LW_u
+@diagnostic_compute "sw_up" Union{SoilCanopyModel, LandModel, CanopyModel} p.SW_u
 @diagnostic_compute "surface_runoff" Union{SoilCanopyModel, LandModel} p.soil.R_s
 @diagnostic_compute "subsurface_runoff" Union{SoilCanopyModel, LandModel} p.soil.R_ss
 

src/integrated/land.jl

@@ -184,6 +184,17 @@ function LandModel{FT}(;
         energy_model = PrescribedCanopyTempModel{FT}()
     end
 
+    if :soil_moisture_stress in propertynames(canopy_component_types)
+        soil_moisture_stress_model =
+        canopy_component_types.soil_moisture_stress(
+                                                    canopy_component_args.soil_moisture_stress...,
+                                                   )
+    else
+        @info "No soil moisture stress model provided, using NoMoistureStressModel (βm = 1.0)"
+        soil_moisture_stress_model =
+        ClimaLand.Canopy.NoMoistureStressModel{FT}()
+    end
+
     canopy = Canopy.CanopyModel{FT}(;
         autotrophic_respiration = canopy_component_types.autotrophic_respiration(
             canopy_component_args.autotrophic_respiration...,
@@ -197,6 +208,7 @@ function LandModel{FT}(;
         conductance = canopy_component_types.conductance(
             canopy_component_args.conductance...,
         ),
+        soil_moisture_stress = soil_moisture_stress_model,
         hydraulics = canopy_component_types.hydraulics(;
             transpiration = transpiration,
             canopy_component_args.hydraulics...,