use new LandModel constructors in experiments

Author: juliasloan25

experiments/benchmarks/snowy_land.jl

@@ -77,10 +77,10 @@ outdir = "snowy_land_benchmark_$(device_suffix)"
 function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     earth_param_set = LP.LandParameters(FT)
     domain = ClimaLand.Domains.global_domain(FT; nelements = nelements)
+    surface_domain = ClimaLand.Domains.obtain_surface_domain(domain)
     surface_space = domain.space.surface
-    subsurface_space = domain.space.subsurface
-
     start_date = DateTime(2008)
+
     # Forcing data
     era5_ncdata_path = ClimaLand.Artifacts.era5_land_forcing_data2008_path(;
         context,
@@ -94,216 +94,69 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         FT;
         time_interpolation_method = time_interpolation_method,
     )
+    forcing = (; atmos, radiation)
 
-    (; ν_ss_om, ν_ss_quartz, ν_ss_gravel) =
-        ClimaLand.Soil.soil_composition_parameters(subsurface_space, FT)
-    (; ν, hydrology_cm, K_sat, θ_r) =
-        ClimaLand.Soil.soil_vangenuchten_parameters(subsurface_space, FT)
-    soil_albedo = Soil.CLMTwoBandSoilAlbedo{FT}(;
-        ClimaLand.Soil.clm_soil_albedo_parameters(surface_space)...,
-    )
-    S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3)
-    soil_params = Soil.EnergyHydrologyParameters(
-        FT;
-        ν,
-        ν_ss_om,
-        ν_ss_quartz,
-        ν_ss_gravel,
-        hydrology_cm,
-        K_sat,
-        S_s,
-        θ_r,
-        albedo = soil_albedo,
-    )
-    f_over = FT(3.28) # 1/m
-    R_sb = FT(1.484e-4 / 1000) # m/s
-    runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
-        f_over = f_over,
-        f_max = ClimaLand.Soil.topmodel_fmax(surface_space, FT),
-        R_sb = R_sb,
-    )
-
-    # Spatially varying canopy parameters from CLM
-    g1 = ClimaLand.Canopy.clm_medlyn_g1(surface_space)
-    rooting_depth = ClimaLand.Canopy.clm_rooting_depth(surface_space)
-    (; is_c3, Vcmax25) =
-        ClimaLand.Canopy.clm_photosynthesis_parameters(surface_space)
-    (; Ω, G_Function, α_PAR_leaf, τ_PAR_leaf, α_NIR_leaf, τ_NIR_leaf) =
-        ClimaLand.Canopy.clm_canopy_radiation_parameters(surface_space)
-
-    # Energy Balance model
-    ac_canopy = FT(2.5e3)
-    # Plant Hydraulics and general plant parameters
-    SAI = FT(0.0) # m2/m2
-    f_root_to_shoot = FT(3.5)
-    RAI = FT(1.0)
-    K_sat_plant = FT(5e-9) # m/s # seems much too small?
-    ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa
-    Weibull_param = FT(4) # unitless, Holtzman's original c param value
-    a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
-    conductivity_model =
-        Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
-    retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
-    plant_ν = FT(1.44e-4)
-    plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    n_stem = 0
-    n_leaf = 1
-    h_stem = FT(0.0)
-    h_leaf = FT(1.0)
-    zmax = FT(0.0)
-    h_canopy = h_stem + h_leaf
-    compartment_midpoints =
-        n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2]
-    compartment_surfaces =
-        n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf]
-
-    z0_m = FT(0.13) * h_canopy
-    z0_b = FT(0.1) * z0_m
-
-    soil_args = (domain = domain, parameters = soil_params)
-    soil_model_type = Soil.EnergyHydrology{FT}
-
-    # Soil microbes model
-    soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT}
-    soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters(FT)
-    Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5))
-
-    soilco2_args = (; domain = domain, parameters = soilco2_ps)
-
-    # Now we set up the canopy model, which we set up by component:
-    # Component Types
-    canopy_component_types = (;
-        autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT},
-        radiative_transfer = Canopy.TwoStreamModel{FT},
-        photosynthesis = Canopy.FarquharModel{FT},
-        conductance = Canopy.MedlynConductanceModel{FT},
-        hydraulics = Canopy.PlantHydraulicsModel{FT},
-        energy = Canopy.BigLeafEnergyModel{FT},
-    )
-    # Individual Component arguments
-    # Set up autotrophic respiration
-    autotrophic_respiration_args =
-        (; parameters = Canopy.AutotrophicRespirationParameters(FT))
-    # Set up radiative transfer
-    radiative_transfer_args = (;
-        parameters = Canopy.TwoStreamParameters(
-            FT;
-            Ω,
-            α_PAR_leaf,
-            τ_PAR_leaf,
-            α_NIR_leaf,
-            τ_NIR_leaf,
-            G_Function,
-        )
-    )
-    # Set up conductance
-    conductance_args =
-        (; parameters = Canopy.MedlynConductanceParameters(FT; g1))
-    # Set up photosynthesis
-    photosynthesis_args =
-        (; parameters = Canopy.FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
-    # Set up plant hydraulics
+    # Read in LAI from MODIS data
     modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(;
         context = nothing,
         start_date = start_date + Second(t0),
         end_date = start_date + Second(t0) + Second(tf),
     )
-    LAIfunction = ClimaLand.prescribed_lai_modis(
+    LAI = ClimaLand.prescribed_lai_modis(
         modis_lai_ncdata_path,
         surface_space,
         start_date;
         time_interpolation_method = time_interpolation_method,
     )
-    ai_parameterization =
-        Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
-
-    plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
-        ai_parameterization = ai_parameterization,
-        ν = plant_ν,
-        S_s = plant_S_s,
-        rooting_depth = rooting_depth,
-        conductivity_model = conductivity_model,
-        retention_model = retention_model,
-    )
-    plant_hydraulics_args = (
-        parameters = plant_hydraulics_ps,
-        n_stem = n_stem,
-        n_leaf = n_leaf,
-        compartment_midpoints = compartment_midpoints,
-        compartment_surfaces = compartment_surfaces,
-    )
-
-    energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),)
-
-    # Canopy component args
-    canopy_component_args = (;
-        autotrophic_respiration = autotrophic_respiration_args,
-        radiative_transfer = radiative_transfer_args,
-        photosynthesis = photosynthesis_args,
-        conductance = conductance_args,
-        hydraulics = plant_hydraulics_args,
-        energy = energy_args,
-    )
 
-    # Other info needed
-    shared_params = Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}(
-        z0_m,
-        z0_b,
-        earth_param_set,
-    )
+    # Overwrite some defaults for the canopy model
+    # Energy model
+    ac_canopy = FT(2.5e3)
+    energy = Canopy.BigLeafEnergyModel{FT}(; ac_canopy)
 
-    canopy_model_args = (;
-        parameters = shared_params,
-        domain = ClimaLand.obtain_surface_domain(domain),
-    )
+    # Roughness lengths
+    hydraulics = Canopy.PlantHydraulicsModel{FT}(surface_domain, LAI)
+    h_canopy = hydraulics.compartment_surfaces[end]
+    z0_m = FT(0.13) * h_canopy
+    z0_b = FT(0.1) * z0_m
 
-    # Snow model
-    snow_parameters = SnowParameters{FT}(Δt; earth_param_set = earth_param_set)
-    snow_args = (;
-        parameters = snow_parameters,
-        domain = ClimaLand.obtain_surface_domain(domain),
+    ground = ClimaLand.PrognosticGroundConditions{FT}()
+    canopy_forcing = (; atmos, radiation, ground)
+
+    canopy = ClimaLand.Canopy.CanopyModel{FT}(
+        surface_domain,
+        canopy_forcing,
+        LAI,
+        earth_param_set;
+        prognostic_land_components = (:canopy, :snow, :soil, :soilco2),
+        energy,
+        hydraulics,
+        z_0m = z0_m,
+        z_0b = z0_b,
     )
-    snow_model_type = Snow.SnowModel
 
-    land_input = (
-        atmos = atmos,
-        radiation = radiation,
-        runoff = runoff_model,
-        soil_organic_carbon = Csom,
-    )
-    land = LandModel{FT}(;
-        soilco2_type = soilco2_type,
-        soilco2_args = soilco2_args,
-        land_args = land_input,
-        soil_model_type = soil_model_type,
-        soil_args = soil_args,
-        canopy_component_types = canopy_component_types,
-        canopy_component_args = canopy_component_args,
-        canopy_model_args = canopy_model_args,
-        snow_args = snow_args,
-        snow_model_type = snow_model_type,
-    )
+    # Construct land model with all default components
+    land = LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; canopy)
 
     Y, p, cds = initialize(land)
 
+    # Set initial conditions
+    (; θ_r, ν, ρc_ds) = land.soil.parameters
     @. Y.soil.ϑ_l = θ_r + (ν - θ_r) / 2
     Y.soil.θ_i .= 0
     T = FT(276.85)
     ρc_s =
         Soil.volumetric_heat_capacity.(
             Y.soil.ϑ_l,
             Y.soil.θ_i,
-            soil_params.ρc_ds,
-            soil_params.earth_param_set,
+            ρc_ds,
+            earth_param_set,
         )
     Y.soil.ρe_int .=
-        Soil.volumetric_internal_energy.(
-            Y.soil.θ_i,
-            ρc_s,
-            T,
-            soil_params.earth_param_set,
-        )
+        Soil.volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T, earth_param_set)
     Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air
+
+    plant_ν = land.canopy.hydraulics.parameters.ν
     Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν
     evaluate!(Y.canopy.energy.T, atmos.T, t0)