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.

NOTE: some new formatting it seems using JuliaFormatter v1

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

Author: AlexisRenchon

docs/src/APIs/canopy/soil_moisture_stress.md

@@ -0,0 +1,26 @@
+# Soil moisture stress
+
+```@meta
+CurrentModule = ClimaLand.Canopy
+```
+
+## Models and Parameters
+
+```@docs
+ClimaLand.Canopy.TuzetMoistureStressParameters
+ClimaLand.Canopy.TuzetMoistureStressModel
+ClimaLand.Canopy.TuzetMoistureStressModel{FT}()
+ClimaLand.Canopy.PiecewiseMoistureStressParameters
+ClimaLand.Canopy.PiecewiseMoistureStressParameters()
+ClimaLand.Canopy.PiecewiseMoistureStressModel
+ClimaLand.Canopy.PiecewiseMoistureStressModel{FT}()
+ClimaLand.Canopy.PiecewiseMoistureStressParametersFromHydrology()
+ClimaLand.Canopy.NoMoistureStressModel
+```
+
+## Methods
+
+```@docs
+ClimaLand.Canopy.compute_piecewise_moisture_stress
+ClimaLand.Canopy.update_soil_moisture_stress!
+```

docs/src/tutorials/calibration/obs_site_level_calibration.jl

@@ -238,7 +238,7 @@ function G(Vcmax25, g1)
                 simulation.start_date,
                 simulation.start_date + Day(20),
                 stop_date,
-            )
+            ),
         )
     return observation
 end;

docs/src/tutorials/calibration/perfect_model_site_level_calibration.jl

@@ -206,7 +206,7 @@ function G(Vcmax25)
                 lhf,
                 simulation.start_date,
                 simulation.start_date + Day(20),
-            )
+            ),
         )
     return observation
 end

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,7 +278,8 @@ simulation = LandSimulation(
     set_ic!,
     updateat,
     diagnostics = diags,
-);
+)
+
 @time solve!(simulation)
 
 comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset)

experiments/integrated/performance/conservation/ozark_conservation.jl

@@ -86,7 +86,7 @@ for float_type in (Float32, Float64)
             [parent(sv.saveval[k].drivers.T)[1] for k in 1:length(sv.t)]
         @assert mean(
             abs.(
-                cos.(radiation.θs.(sv.t, radiation.start_date)) .- cache_cosθs
+                cos.(radiation.θs.(sv.t, radiation.start_date)) .- cache_cosθs,
             ),
         ) < eps(FT)
         T_mutable = Vector{FT}(undef, 1)
@@ -201,7 +201,7 @@ for float_type in (Float32, Float64)
             eps(FT) .+
             abs.(
                 (soil_mass_change_actual - soil_mass_change_exp) ./
-                soil_mass_change_exp
+                soil_mass_change_exp,
             ),
             label = "Soil Water Balance",
         )
@@ -211,7 +211,7 @@ for float_type in (Float32, Float64)
             eps(FT) .+
             abs.(
                 (canopy_mass_change_actual - canopy_mass_change_exp) ./
-                canopy_mass_change_exp
+                canopy_mass_change_exp,
             ),
             label = "Canopy Water Balance",
         )
@@ -272,7 +272,7 @@ for float_type in (Float32, Float64)
             eps(FT) .+
             abs.(
                 (soil_energy_change_actual - soil_energy_change_exp) ./
-                soil_energy_change_exp
+                soil_energy_change_exp,
             ),
             label = "Soil Energy Balance",
         )

ext/land_sim_vis/leaderboard/leaderboard.jl

@@ -205,7 +205,7 @@ function compute_monthly_leaderboard(
         months =
             Dates.month.(
                 Dates.DateTime(sim_var.attributes["start_date"]) .+
-                Dates.Second.(times)
+                Dates.Second.(times),
             )
         months_split, sim_vec_split, rmse_vec_split, bias_vec_split =
             partition_by_val(12, months, sim_vec, rmse_vec, bias_vec)

ext/neural_snow/DataTools.jl

@@ -761,7 +761,7 @@ function serreze_qc(input::DataFrame, id::Int, state::AbstractString)
     data[!, :tmax] = maxmin[!, :tmax]
     flags =
         ismissing.(
-            data[1:(end - 1), [:SWE, :precip, :air_temp_avg, :tmin, :tmax]]
+            data[1:(end - 1), [:SWE, :precip, :air_temp_avg, :tmin, :tmax]],
         )
     flags[!, :date] = data[1:(end - 1), :date]
 

src/diagnostics/default_diagnostics.jl

@@ -112,7 +112,7 @@ function default_diagnostics(
         diagnostics = get_short_diagnostics(model)
     else
         @assert typeof(output_vars) <: Vector{String}
-        @assert all([var ∈ possible_diags for var in output_vars])
+        @assert all([var in possible_diags for var in output_vars])
         diagnostics = output_vars
     end
 
@@ -178,7 +178,7 @@ function default_diagnostics(
         diagnostics = possible_diags
     else
         @assert typeof(output_vars) <: Vector{String}
-        @assert all([var ∈ possible_diags for var in output_vars])
+        @assert all([var in possible_diags for var in output_vars])
         diagnostics = output_vars
     end
 

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...,