diff --git a/docs/src/APIs/canopy/PlantHydraulics.md b/docs/src/APIs/canopy/PlantHydraulics.md index 4e7ab468d3..aaeecf693d 100644 --- a/docs/src/APIs/canopy/PlantHydraulics.md +++ b/docs/src/APIs/canopy/PlantHydraulics.md @@ -6,7 +6,7 @@ CurrentModule = ClimaLand.PlantHydraulics ## Models ```@docs -ClimaLand.PlantHydraulics.PlantHydraulicsModel +ClimaLand.PlantHydraulics.BigLeafHydraulicsModel ``` ## Plant Hydraulics Diagnostic Variables diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl index 3309fe6beb..120c63368f 100644 --- a/docs/tutorials/integrated/soil_canopy_tutorial.jl +++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl @@ -68,12 +68,8 @@ f_root_to_shoot = FT(3.5) SAI = FT(0.00242) maxLAI = FT(4.2) plant_ν = FT(2.46e-4) -n_stem = Int64(1) -n_leaf = Int64(1) h_stem = FT(9) h_leaf = FT(9.5) -compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2] -compartment_surfaces = [zmax, h_stem, h_stem + h_leaf] land_domain = Column(; zlim = (zmin, zmax), nelements = nelements); # - We will be using prescribed atmospheric and radiative drivers from the @@ -215,7 +211,7 @@ canopy_component_types = (; radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, ); # Then provide arguments to the canopy radiative transfer, stomatal conductance, @@ -270,13 +266,8 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; 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, -); +plant_hydraulics_args = + (parameters = plant_hydraulics_ps, h_stem = h_stem, h_leaf = h_leaf); # We may now collect all of the canopy component argument tuples into one # arguments tuple for the canopy component models. diff --git a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl index e84449b1d7..58b65af0a7 100644 --- a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl +++ b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl @@ -84,12 +84,8 @@ f_root_to_shoot = FT(3.5) SAI = FT(0.00242) maxLAI = FT(4.2) plant_ν = FT(2.46e-4) # kg/m^2 -n_stem = Int64(1) -n_leaf = Int64(1) h_stem = FT(9) h_leaf = FT(9.5) -compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2] -compartment_surfaces = [zmax, h_stem, h_stem + h_leaf] land_domain = Point(; z_sfc = FT(0.0)) # - We will be using prescribed atmospheric and radiative drivers from the @@ -233,12 +229,10 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; # Define the remaining variables required for the plant hydraulics model. -plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; +plant_hydraulics = PlantHydraulics.BigLeafHydraulicsModel{FT}(; parameters = plant_hydraulics_ps, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_surfaces, - compartment_midpoints = compartment_midpoints, + h_stem = h_stem, + h_leaf = h_leaf, ); # Now, instantiate the canopy model, using the atmospheric and radiative diff --git a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl index 6afcd4670d..1eb17ad89d 100644 --- a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl +++ b/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl @@ -9,8 +9,6 @@ dz_bottom = FT(1.5) dz_top = FT(0.025) # Stem and leaf compartments and their heights: -n_stem = Int64(1) -n_leaf = Int64(1) h_leaf = FT(12) # m h_stem = FT(14) # m diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl index 993e855acc..115144cce8 100644 --- a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl +++ b/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl @@ -9,8 +9,6 @@ dz_bottom = FT(1.5) dz_top = FT(0.025) # Stem and leaf compartments and their heights: -n_stem = Int64(1) -n_leaf = Int64(1) h_stem = FT(9) # m h_leaf = FT(9.5) # m diff --git a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl index e27c0b4827..3ad13a9079 100644 --- a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl +++ b/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl @@ -9,8 +9,6 @@ dz_bottom = FT(1.25) dz_top = FT(0.05) # Stem and leaf compartments and their heights: -n_stem = Int64(1) -n_leaf = Int64(1) h_leaf = FT(6.5) # m h_stem = FT(7.5) # m diff --git a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl index c19746d33c..093b73a16a 100644 --- a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl +++ b/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl @@ -9,8 +9,6 @@ dz_bottom = FT(1.0) dz_top = FT(0.05) # Stem and leaf compartments and their heights: -n_stem = Int64(0) -n_leaf = Int64(1) h_leaf = FT(0.7) # m h_stem = FT(0) # m diff --git a/experiments/integrated/fluxnet/fluxnet_domain.jl b/experiments/integrated/fluxnet/fluxnet_domain.jl index b5d1be8012..0f969ad59c 100644 --- a/experiments/integrated/fluxnet/fluxnet_domain.jl +++ b/experiments/integrated/fluxnet/fluxnet_domain.jl @@ -8,9 +8,6 @@ nelements = 20 zmin = FT(-10) zmax = FT(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] land_domain = Column(; zlim = (zmin, zmax), diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl index 1537c3c77b..baae661d20 100644 --- a/experiments/integrated/fluxnet/ozark_pft.jl +++ b/experiments/integrated/fluxnet/ozark_pft.jl @@ -168,7 +168,7 @@ canopy_component_types = (; radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, energy = Canopy.BigLeafEnergyModel{FT}, ) # Individual Component arguments @@ -209,13 +209,8 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; 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, -) +plant_hydraulics_args = + (parameters = plant_hydraulics_ps, h_stem = h_stem, h_leaf = h_leaf) energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) @@ -276,12 +271,12 @@ Y.soil.ρe_int = Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m] ψ_leaf_0 = FT(-2e5 / 9800) -ψ_comps = n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 +ψ_comps = h_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 S_l_ini = inverse_water_retention_curve.(retention_model, ψ_comps, plant_ν, plant_S_s) -for i in 1:(n_stem + n_leaf) +for i in 1:(h_stem < 0 + 1) Y.canopy.hydraulics.ϑ_l.:($i) .= augmented_liquid_fraction.(plant_ν, S_l_ini[i]) end diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index 321a2f0d46..69cc9a60dd 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -130,7 +130,7 @@ canopy_component_types = (; radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, energy = Canopy.BigLeafEnergyModel{FT}, ) # Individual Component arguments @@ -169,13 +169,8 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; 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, -) +plant_hydraulics_args = + (parameters = plant_hydraulics_ps, h_stem = h_stem, h_leaf = h_leaf) energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) @@ -236,12 +231,12 @@ Y.soil.ρe_int = Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m] ψ_leaf_0 = FT(-2e5 / 9800) -ψ_comps = n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 +ψ_comps = h_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 S_l_ini = inverse_water_retention_curve.(retention_model, ψ_comps, plant_ν, plant_S_s) -for i in 1:(n_stem + n_leaf) +for i in 1:(h_stem > 0 + 1) Y.canopy.hydraulics.ϑ_l.:($i) .= augmented_liquid_fraction.(plant_ν, S_l_ini[i]) end diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl index ebf0e98ae2..df788f1b18 100644 --- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl +++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl @@ -6,8 +6,6 @@ dz_bottom = FT(1.5) dz_top = FT(0.025) # Stem and leaf compartments and their heights: -n_stem = Int64(1) -n_leaf = Int64(1) h_stem = FT(9) # m h_leaf = FT(9.5) # m @@ -112,7 +110,7 @@ canopy_component_types = (; radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, energy = Canopy.BigLeafEnergyModel{FT}, ) # Individual Component arguments @@ -153,13 +151,8 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; 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, -) +plant_hydraulics_args = + (parameters = plant_hydraulics_ps, h_stem = h_stem, h_leaf = h_leaf) # Canopy component args canopy_component_args = (; diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl index bcc5e675fd..214cd40227 100644 --- a/experiments/integrated/performance/profile_allocations.jl +++ b/experiments/integrated/performance/profile_allocations.jl @@ -86,14 +86,10 @@ dz_bottom = FT(2.0) dz_top = FT(0.2) soil_depth = FT(5) z_sfc = FT(0) -n_stem = Int64(1) -n_leaf = Int64(1) h_stem = FT(9) # m h_leaf = FT(9.5) # m h_canopy = h_stem + h_leaf -compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2] -compartment_surfaces = [z_sfc, h_stem, h_canopy] land_domain = ClimaLand.Domains.SphericalShell(; radius = FT(6.3781e6), @@ -258,7 +254,7 @@ canopy_component_types = (; radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, energy = Canopy.BigLeafEnergyModel{FT}, ) # Individual Component arguments @@ -287,13 +283,8 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; 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, -) +plant_hydraulics_args = + (parameters = plant_hydraulics_ps, h_stem = h_stem, h_leaf = h_leaf) # Canopy component args canopy_component_args = (; diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl index bccdbc01ef..455114293f 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_simulation.jl @@ -4,8 +4,6 @@ export harvard_default_args function harvard_default_args(; dz_bottom = FT(1.5), dz_top = FT(0.025), - n_stem = Int64(1), - n_leaf = Int64(1), h_leaf = FT(12), # m h_stem = FT(14), # m t0 = Float64(120 * 3600 * 24), @@ -15,8 +13,6 @@ export harvard_default_args return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, @@ -30,8 +26,6 @@ Named Tuple containing default simulation parameter for the Harvard site. function harvard_default_args(; dz_bottom = FT(1.5), dz_top = FT(0.025), - n_stem = Int64(1), - n_leaf = Int64(1), h_leaf = FT(12), # m h_stem = FT(14), # m t0 = Float64(120 * 3600 * 24), @@ -41,8 +35,6 @@ function harvard_default_args(; return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl index dcbff11279..23d551438b 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-MOz/US-MOz_simulation.jl @@ -4,8 +4,6 @@ export ozark_default_args function ozark_default_args(; dz_bottom = FT(1.5), dz_top = FT(0.025), - n_stem = Int64(1), - n_leaf = Int64(1), h_stem = FT(9), h_leaf = FT(9.5), t0 = Float64(120 * 3600 * 24), @@ -15,8 +13,6 @@ export ozark_default_args return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, @@ -30,8 +26,6 @@ Named Tuple containing default simulation parameter for the Ozark site. function ozark_default_args(; dz_bottom = FT(1.5), dz_top = FT(0.025), - n_stem = Int64(1), - n_leaf = Int64(1), h_stem = FT(9), h_leaf = FT(9.5), t0 = Float64(120 * 3600 * 24), @@ -41,8 +35,6 @@ function ozark_default_args(; return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl index 0d3acf86be..4227a98414 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_simulation.jl @@ -4,8 +4,6 @@ export niwotridge_default_args function niwotridge_default_args(; dz_bottom = FT(1.25), dz_top = FT(0.05), - n_stem = Int64(1), - n_leaf = Int64(1), h_leaf = FT(6.5), # m, h_stem = FT(7.5), # m, t0 = Float64(120 * 3600 * 24), # start mid year to avoid snow, @@ -15,8 +13,6 @@ export niwotridge_default_args return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, @@ -30,8 +26,6 @@ Named Tuple containing default simulation parameter for the Niwot Ridge site. function niwotridge_default_args(; dz_bottom = FT(1.25), dz_top = FT(0.05), - n_stem = Int64(1), - n_leaf = Int64(1), h_leaf = FT(6.5), # m, h_stem = FT(7.5), # m, t0 = Float64(120 * 3600 * 24), # start mid year to avoid snow, @@ -41,8 +35,6 @@ function niwotridge_default_args(; return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl index 66dc5735c3..360ec63b2a 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_simulation.jl @@ -4,8 +4,6 @@ export vairaranch_default_args function vairaranch_default_args(; dz_bottom = FT(1.0), dz_top = FT(0.05), - n_stem = Int64(0), - n_leaf = Int64(1), h_leaf = FT(0.7), # m, h_stem = FT(0), # m, t0 = Float64(21 * 3600 * 24), # start day 21 of the year, @@ -15,8 +13,6 @@ export vairaranch_default_args return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, @@ -30,8 +26,6 @@ Named Tuple containing default simulation parameter for the Vaira Ranch site. function vairaranch_default_args(; dz_bottom = FT(1.0), dz_top = FT(0.05), - n_stem = Int64(0), - n_leaf = Int64(1), h_leaf = FT(0.7), # m, h_stem = FT(0), # m, t0 = Float64(21 * 3600 * 24), # start day 21 of the year, @@ -41,8 +35,6 @@ function vairaranch_default_args(; return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl index 5f80beee14..9fb3738ecd 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl @@ -6,30 +6,16 @@ export make_setup setup site specific domain and time stepping: - dz_bottom, the bottom of the soil domain (m) - dz_top, the top of the soil domain (m) -- n_stem, the number of stems -- n_leaf, the number of leaves - h_stem, the height of each stem (m) - h_leaf, the height of each leaf (m) - t0, the start time of the simulation (s) - dt, the time step (s) - n, the number of time step between saving outputs """ -function make_setup(; - dz_bottom, - dz_top, - n_stem, - n_leaf, - h_stem, - h_leaf, - t0, - dt, - n, -) +function make_setup(; dz_bottom, dz_top, h_stem, h_leaf, t0, dt, n) return ( dz_bottom = dz_bottom, dz_top = dz_top, - n_stem = n_stem, - n_leaf = n_leaf, h_stem = h_stem, h_leaf = h_leaf, t0 = t0, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl index 0c4a36fa11..b49356f4a3 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl @@ -93,7 +93,7 @@ function run_fluxnet( radiative_transfer = Canopy.TwoStreamModel{FT}, photosynthesis = Canopy.FarquharModel{FT}, conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, + hydraulics = Canopy.BigLeafHydraulicsModel{FT}, energy = Canopy.BigLeafEnergyModel{FT}, ) # Individual Component arguments @@ -180,13 +180,7 @@ function run_fluxnet( conductivity_model = params.plant_hydraulics.conductivity_model, retention_model = params.plant_hydraulics.retention_model, ) - plant_hydraulics_args = ( - parameters = plant_hydraulics_ps, - n_stem = setup.n_stem, - n_leaf = setup.n_leaf, - compartment_midpoints = domain.compartment_midpoints, - compartment_surfaces = domain.compartment_surfaces, - ) + plant_hydraulics_args = (parameters = plant_hydraulics_ps, h_stem, h_leaf) energy_args = ( parameters = Canopy.BigLeafEnergyParameters{FT}( @@ -261,7 +255,7 @@ function run_fluxnet( Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m] ψ_leaf_0 = FT(-2e5 / 9800) - ψ_comps = setup.n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 + ψ_comps = setup.h_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 S_l_ini = inverse_water_retention_curve.( @@ -271,7 +265,7 @@ function run_fluxnet( params.plant_hydraulics.S_s, ) - for i in 1:(setup.n_stem + setup.n_leaf) + for i in 1:(setup.h_stem > 0 + 1) Y.canopy.hydraulics.ϑ_l.:($i) .= augmented_liquid_fraction.(drivers.plant_ν, S_l_ini[i]) end diff --git a/lib/ClimaLandSimulations/src/utilities/make_domain.jl b/lib/ClimaLandSimulations/src/utilities/make_domain.jl index 42755bb4f5..45b6a86cf9 100644 --- a/lib/ClimaLandSimulations/src/utilities/make_domain.jl +++ b/lib/ClimaLandSimulations/src/utilities/make_domain.jl @@ -9,8 +9,6 @@ export make_domain setup the soil and canopy domain of the simulation. Returns: - h_canopy, -- compartment_midpoints, -- compartment_surfaces, - land_domain, - canopy_domain, """ @@ -18,13 +16,6 @@ function make_domain(setup, FT; nelements = 20, zmin = FT(-10), zmax = FT(0)) h_canopy = setup.h_stem + setup.h_leaf - compartment_midpoints = - setup.n_stem > 0 ? [setup.h_stem / 2, setup.h_stem + setup.h_leaf / 2] : - [setup.h_leaf / 2] - - compartment_surfaces = - setup.n_stem > 0 ? [zmax, setup.h_stem, h_canopy] : [zmax, setup.h_leaf] - land_domain = Column(; zlim = (zmin, zmax), nelements = nelements, @@ -38,8 +29,6 @@ function make_domain(setup, FT; nelements = 20, zmin = FT(-10), zmax = FT(0)) zmin = zmin, zmax = zmax, h_canopy = h_canopy, - compartment_midpoints = compartment_midpoints, - compartment_surfaces = compartment_surfaces, land_domain = land_domain, canopy_domain = canopy_domain, ) diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl index 754f9184b3..92b13171a3 100644 --- a/src/integrated/soil_canopy_model.jl +++ b/src/integrated/soil_canopy_model.jl @@ -224,13 +224,14 @@ lsm_aux_domain_names(m::SoilCanopyModel) = FT, MM <: Soil.Biogeochemistry.SoilCO2Model{FT}, SM <: Soil.RichardsModel{FT}, - RM <: Canopy.CanopyModel{FT} + RM <: Canopy.CanopyModel{FT, PlantHydraulics.BigLeafHydraulicsModel{FT}} } A method which makes a function; the returned function updates the additional auxiliary variables for the integrated model, as well as updates the boundary auxiliary variables for all component -models. +models. This version is for an integrated canopy model which suports a +multi-layered hydraulics model. This function is called each ode function evaluation, prior to the tendency function evaluation. @@ -241,7 +242,7 @@ function make_update_boundary_fluxes( FT, MM <: Soil.Biogeochemistry.SoilCO2Model{FT}, SM <: Soil.EnergyHydrology{FT}, - RM <: Canopy.CanopyModel{FT}, + RM <: Canopy.CanopyModel{FT, PlantHydraulics.BigLeafHydraulicsModel{FT}}, } update_soil_bf! = make_update_boundary_fluxes(land.soil) update_soilco2_bf! = make_update_boundary_fluxes(land.soilco2) @@ -256,16 +257,20 @@ function make_update_boundary_fluxes( area_index = p.canopy.hydraulics.area_index - above_ground_area_index = getproperty( - area_index, - land.canopy.hydraulics.compartment_labels[1], - ) + labels = land.canopy.hydraulics.h_stem > 0 ? [:stem, :leaf] : [:leaf] + + above_ground_area_index = getproperty(area_index, labels[1]) + + h_stem = land.canopy.hydraulics.h_stem + h_leaf = land.canopy.hydraulics.h_leaf + + midpoint = labels[1] == :stem ? h_stem / 2 : model.h_leaf / 2 @. p.root_extraction = (area_index.root + above_ground_area_index) / 2 * PlantHydraulics.flux( z, - land.canopy.hydraulics.compartment_midpoints[1], + midpoint, p.soil.ψ, p.canopy.hydraulics.ψ.:1, PlantHydraulics.hydraulic_conductivity( @@ -300,6 +305,87 @@ function make_update_boundary_fluxes( return update_boundary_fluxes! end +""" + make_update_boundary_fluxes( + land::SoilCanopyModel{FT, MM, SM, RM}, + ) where { + FT, + MM <: Soil.Biogeochemistry.SoilCO2Model{FT}, + SM <: Soil.RichardsModel{FT}, + RM <: Canopy.CanopyModel{FT, PlantHydraulics.BigLeafHydraulicsModel{FT}} + } + +A method which makes a function; the returned function updates the additional +auxiliary variables for the integrated model, as well as updates the boundary +auxiliary variables for all component models. This version is for an integrated +canopy model which supports only a big leaf hydraulics model. + +This function is called each ode function evaluation, prior to the tendency +function evaluation. +""" +function make_update_boundary_fluxes( + land::SoilCanopyModel{FT, MM, SM, RM}, +) where { + FT, + MM <: Soil.Biogeochemistry.SoilCO2Model{FT}, + SM <: Soil.EnergyHydrology{FT}, + RM <: Canopy.CanopyModel{FT, PlantHydraulics.BigLeafHydraulicsModel{FT}}, +} + update_soil_bf! = make_update_boundary_fluxes(land.soil) + update_soilco2_bf! = make_update_boundary_fluxes(land.soilco2) + update_canopy_bf! = make_update_boundary_fluxes(land.canopy) + function update_boundary_fluxes!(p, Y, t) + # update root extraction + z = + ClimaCore.Fields.coordinate_field( + land.soil.domain.space.subsurface, + ).z + (; h_stem, h_leaf, conductivity_model) = + land.canopy.hydraulics.parameters + (label, midpoint) = + h_stem > 0 ? (:stem, h_stem / 2) : (:leaf, h_leaf / 2) + + area_index = p.canopy.hydraulics.area_index + + above_ground_area_index = getproperty(area_index, label) + + @. p.root_extraction = + PlantHydraulics.flux( + z, + midpoint, + p.soil.ψ, + p.canopy.hydraulics.ψ.:1, + PlantHydraulics.hydraulic_conductivity( + conductivity_model, + p.soil.ψ, + ), + PlantHydraulics.hydraulic_conductivity( + conductivity_model, + p.canopy.hydraulics.ψ.:1, + ), + ) * (land.canopy.hydraulics.parameters.root_distribution(z)) + @. p.root_energy_extraction = + p.root_extraction * ClimaLand.Soil.volumetric_internal_energy_liq( + p.soil.T, + land.soil.parameters, + ) + + # Radiation + lsm_radiant_energy_fluxes!( + p, + land.canopy.radiative_transfer, + land.canopy, + land.soil, + Y, + t, + ) + update_soil_bf!(p, Y, t) + update_canopy_bf!(p, Y, t) + update_soilco2_bf!(p, Y, t) + end + return update_boundary_fluxes! +end + """ lsm_radiant_energy_fluxes!(p, canopy_radiation::Canopy.AbstractRadiationModel{FT}, @@ -456,7 +542,7 @@ end PlantHydraulics.root_water_flux_per_ground_area!( fa::ClimaCore.Fields.Field, s::PrognosticSoil, - model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, + model::Canopy.PlantHydraulics.BigLeafHydraulicsModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, @@ -475,7 +561,7 @@ roots and soil at each soil layer. function PlantHydraulics.root_water_flux_per_ground_area!( fa::ClimaCore.Fields.Field, s::PrognosticSoil, - model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, + model::Canopy.PlantHydraulics.BigLeafHydraulicsModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index 1e4baf9194..9b22bc748e 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/src/standalone/Vegetation/Canopy.jl @@ -87,7 +87,7 @@ running the canopy model in standalone mode), or prognostic (of type PrognosticSoil, for running integrated soil+canopy models) Note that the canopy height is specified as part of the -PlantHydraulicsModel, along with the area indices of the leaves, roots, and +BigLeafHydraulicsModel, along with the area indices of the leaves, roots, and stems. Eventually, when plant biomass becomes a prognostic variable (by integrating with a carbon model), some parameters specified here will be treated differently. @@ -143,7 +143,7 @@ end An outer constructor for the `CanopyModel`. The primary constraints this applies are (1) ensuring that the domain is 1d or 2d (a ``surface" domain of a column, box, or sphere) and (2) ensuring -consistency between the PlantHydraulics model and the general canopy model, +consistency between the BigLeafHydraulics model and the general canopy model, since these are built separately. """ function CanopyModel{FT}(; @@ -362,13 +362,13 @@ end <:Union{BeerLambertModel, TwoStreamModel}, <:FarquharModel, <:MedlynConductanceModel, - <:PlantHydraulicsModel,}, + <:BigLeafHydraulicsModel,}, ) where {FT} Creates the `update_aux!` function for the `CanopyModel`; a specific method for `update_aux!` for the case where the canopy model components are of the type in the parametric type signature: `AutotrophicRespirationModel`, `AbstractRadiationModel`, -`FarquharModel`, `MedlynConductanceModel`, and `PlantHydraulicsModel`. +`FarquharModel`, `MedlynConductanceModel`, and `BigLeafHydraulicsModel`. Please note that the plant hydraulics model has auxiliary variables that are updated in its prognostic `compute_exp_tendency!` function. @@ -386,7 +386,7 @@ function ClimaLand.make_update_aux( <:Union{BeerLambertModel, TwoStreamModel}, <:Union{FarquharModel, OptimalityFarquharModel}, <:MedlynConductanceModel, - <:PlantHydraulicsModel, + <:BigLeafHydraulicsModel, <:AbstractCanopyEnergyModel, }, ) where {FT} @@ -488,9 +488,7 @@ function ClimaLand.make_update_aux( # update plant hydraulics aux hydraulics = canopy.hydraulics - n_stem = hydraulics.n_stem - n_leaf = hydraulics.n_leaf - PlantHydraulics.lai_consistency_check.(n_stem, n_leaf, area_index) + PlantHydraulics.lai_consistency_check.(hydraulics, area_index) (; retention_model, conductivity_model, S_s, ν) = hydraulics.parameters # We can index into a field of Tuple{FT} to extract a field of FT # using the following notation: field.:index @@ -506,8 +504,10 @@ function ClimaLand.make_update_aux( # for broadcasted expressions using the macro @. # field.:($index) .= value # works # @ field.:($$index) = value # works - @inbounds for i in 1:(n_stem + n_leaf - 1) - ip1 = i + 1 + labels = hydraulics.h_stem > 0 ? [:stem, :leaf] : [:leaf] + if hydraulics.h_stem > 0 + i = 1 + ip1 = 2 @. ψ.:($$ip1) = PlantHydraulics.water_retention_curve( retention_model, PlantHydraulics.effective_saturation(ν, ϑ_l.:($$ip1)), @@ -515,16 +515,18 @@ function ClimaLand.make_update_aux( S_s, ) - areai = getproperty(area_index, hydraulics.compartment_labels[i]) - areaip1 = - getproperty(area_index, hydraulics.compartment_labels[ip1]) + areai = getproperty(area_index, labels[i]) + areaip1 = getproperty(area_index, labels[ip1]) + + midpoint = hydraulics.h_stem / 2 + midpointip1 = hydraulics.h_stem + (hydraulics.h_leaf / 2) # Compute the flux*area between the current compartment `i` # and the compartment above. @. fa.:($$i) = PlantHydraulics.flux( - hydraulics.compartment_midpoints[i], - hydraulics.compartment_midpoints[ip1], + midpoint, + midpointip1, ψ.:($$i), ψ.:($$ip1), PlantHydraulics.hydraulic_conductivity( @@ -546,7 +548,7 @@ function ClimaLand.make_update_aux( T_canopy = canopy_temperature(canopy.energy, canopy, Y, p, t) # update moisture stress - i_end = n_stem + n_leaf + i_end = (hydraulics.h_stem > 0 ? 1 : 0) + 1 @. β = moisture_stress(ψ.:($$i_end) * ρ_l * grav, sc, pc) # Update Rd, An, Vcmax25 (if applicable to model) in place @@ -566,7 +568,7 @@ function ClimaLand.make_update_aux( @. GPP = compute_GPP(An, K, LAI, Ω) @. gs = medlyn_conductance(g0, Drel, medlyn_factor, An, c_co2_air) # update autotrophic respiration - h_canopy = hydraulics.compartment_surfaces[end] + h_canopy = hydraulics.h_stem + hydraulics.h_leaf @. Ra = compute_autrophic_respiration( canopy.autotrophic_respiration, Vcmax25, @@ -593,7 +595,7 @@ function make_compute_exp_tendency( <:Union{BeerLambertModel, TwoStreamModel}, <:Union{FarquharModel, OptimalityFarquharModel}, <:MedlynConductanceModel, - <:PlantHydraulicsModel, + <:BigLeafHydraulicsModel, <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel}, }, ) where {FT} diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl index 07de4f6909..c3b001935c 100644 --- a/src/standalone/Vegetation/PlantHydraulics.jl +++ b/src/standalone/Vegetation/PlantHydraulics.jl @@ -21,7 +21,7 @@ import ClimaLand: auxiliary_domain_names, prognostic_domain_names, name -export PlantHydraulicsModel, +export BigLeafHydraulicsModel, AbstractPlantHydraulicsModel, flux, effective_saturation, @@ -48,6 +48,9 @@ abstract type AbstractPlantHydraulicsModel{FT} <: AbstractCanopyComponent{FT} en ClimaLand.name(::AbstractPlantHydraulicsModel) = :hydraulics +# Hydraulics models should be broadcastable +Base.broadcastable(x::AbstractPlantHydraulicsModel) = tuple(x) + """ AbstractTranspiration{FT <: AbstractFloat} @@ -85,45 +88,6 @@ function PrescribedSiteAreaIndex{FT}( PrescribedSiteAreaIndex{FT, typeof(LAIfunction)}(LAIfunction, SAI, RAI) end -""" - lai_consistency_check( - n_stem::Int64, - n_leaf::Int64, - area_index::NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}}, - ) where {FT} - -Carries out consistency checks using the area indices supplied and the number of -stem elements being modeled. - -Note that it is possible to have a plant with no stem compartments -but with leaf compartments, and that a plant must have leaf compartments -(even if LAI = 0). - -Specifically, this checks that: -1. n_leaf > 0 -2. if LAI is nonzero or SAI is nonzero, RAI must be nonzero. -3. if SAI > 0, n_stem must be > 0 for consistency. If SAI == 0, n_stem must -be zero. -""" -function lai_consistency_check( - n_stem::Int64, - n_leaf::Int64, - area_index::NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}}, -) where {FT} - @assert n_leaf > 0 - if area_index.leaf > eps(FT) || area_index.stem > eps(FT) - @assert area_index.root > eps(FT) - end - # If there SAI is zero, there should be no stem compartment - if area_index.stem < eps(FT) - @assert n_stem == FT(0) - else - # if SAI is > 0, n_stem should be > 0 for consistency - @assert n_stem > 0 - end - -end - """ PlantHydraulicsParameters @@ -176,157 +140,6 @@ function PlantHydraulicsParameters(; ) end -""" - PlantHydraulicsModel{FT, PS, T, AA} <: AbstractPlantHydraulicsModel{FT} - -Defines, and constructs instances of, the PlantHydraulicsModel type, which is used -for simulation flux of water to/from soil, along roots of different depths, -along a stem, to a leaf, and ultimately being lost from the system by -transpiration. Note that the canopy height is specified as part of the -PlantHydraulicsModel, along with the area indices of the leaves, roots, and -stems. - -This model can also be combined with the soil model using ClimaLand, in which -case the prognostic soil water content is used to determine root extraction, and -the transpiration is also computed diagnostically. In global run with patches -of bare soil, you can "turn off" the canopy model (to get zero root extraction, zero absorption and -emission, zero transpiration and sensible heat flux from the canopy), by setting: -- n_leaf = 1 -- n_stem = 0 -- LAI = SAI = RAI = 0. - -A plant model can have leaves but no stem, but not vice versa. If n_stem = 0, SAI must be zero. - -Finally, the model can be used in Canopy standalone mode by prescribing -the soil matric potential at the root tips or flux in the roots. There is also the -option (intendend only for debugging) to use a prescribed transpiration rate. - -$(DocStringExtensions.FIELDS) -""" -struct PlantHydraulicsModel{FT, PS, T, AA <: AbstractArray{FT}} <: - AbstractPlantHydraulicsModel{FT} - "The number of stem compartments for the plant; can be zero" - n_stem::Int64 - "The number of leaf compartments for the plant; must be >=1" - n_leaf::Int64 - "The height of the center of each leaf compartment/stem compartment, in meters" - compartment_midpoints::AA - "The height of the compartments' top faces, in meters. The canopy height is the last element of the vector." - compartment_surfaces::AA - "The label (:stem or :leaf) of each compartment" - compartment_labels::Vector{Symbol} - "Parameters required by the Plant Hydraulics model" - parameters::PS - "The transpiration model, of type `AbstractTranspiration`" - transpiration::T -end - -function PlantHydraulicsModel{FT}(; - n_stem::Int64, - n_leaf::Int64, - compartment_midpoints::Vector{FT}, - compartment_surfaces::Vector{FT}, - parameters::PlantHydraulicsParameters{FT}, - transpiration::AbstractTranspiration{FT} = DiagnosticTranspiration{FT}(), -) where {FT} - args = (parameters, transpiration) - @assert (n_leaf + n_stem) == length(compartment_midpoints) - @assert (n_leaf + n_stem) + 1 == length(compartment_surfaces) - for i in 1:length(compartment_midpoints) - @assert compartment_midpoints[i] == - ((compartment_surfaces[i + 1] - compartment_surfaces[i]) / 2) + - compartment_surfaces[i] - end - compartment_labels = Vector{Symbol}(undef, n_stem + n_leaf) - for i in 1:(n_stem + n_leaf) - if i <= n_stem - compartment_labels[i] = :stem - else - compartment_labels[i] = :leaf - end - end - return PlantHydraulicsModel{ - FT, - typeof.(args)..., - typeof(compartment_midpoints), - }( - n_stem, - n_leaf, - compartment_midpoints, - compartment_surfaces, - compartment_labels, - args..., - ) -end - -""" - prognostic_vars(model::PlantHydraulicsModel) - -A function which returns the names of the prognostic -variables of the `PlantHydraulicsModel`. -""" -prognostic_vars(model::PlantHydraulicsModel) = (:ϑ_l,) - -""" - auxiliary_vars(model::PlantHydraulicsModel) - -A function which returns the names of the auxiliary -variables of the `PlantHydraulicsModel`, -the transpiration stress factor `β` (unitless), -the water potential `ψ` (m), the volume flux*cross section `fa` (1/s), -and the volume flux*root cross section in the roots `fa_roots` (1/s), -where the cross section can be represented by an area index. -""" -auxiliary_vars(model::PlantHydraulicsModel) = - (:β, :ψ, :fa, :fa_roots, :area_index) - -""" - ClimaLand.prognostic_types(model::PlantHydraulicsModel{FT}) where {FT} - -Defines the prognostic types for the PlantHydraulicsModel. -""" -ClimaLand.prognostic_types(model::PlantHydraulicsModel{FT}) where {FT} = - (NTuple{model.n_stem + model.n_leaf, FT},) -ClimaLand.prognostic_domain_names(::PlantHydraulicsModel) = (:surface,) - -""" - ClimaLand.auxiliary_types(model::PlantHydraulicsModel{FT}) where {FT} - -Defines the auxiliary types for the PlantHydraulicsModel. -""" -ClimaLand.auxiliary_types(model::PlantHydraulicsModel{FT}) where {FT} = ( - FT, - NTuple{model.n_stem + model.n_leaf, FT}, - NTuple{model.n_stem + model.n_leaf, FT}, - FT, - NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}}, -) -ClimaLand.auxiliary_domain_names(::PlantHydraulicsModel) = - (:surface, :surface, :surface, :surface, :surface) - - -""" - set_canopy_prescribed_field!(component::PlantHydraulics{FT}, - p, - t, - ) where {FT} - - -Sets the canopy prescribed fields pertaining to the PlantHydraulics -component (the area indices) with their values at time t. -""" -function ClimaLand.Canopy.set_canopy_prescribed_field!( - component::PlantHydraulicsModel{FT}, - p, - t, -) where {FT} - (; LAIfunction, SAI, RAI) = component.parameters.ai_parameterization - evaluate!(p.canopy.hydraulics.area_index.leaf, LAIfunction, t) - @. p.canopy.hydraulics.area_index.stem = SAI - @. p.canopy.hydraulics.area_index.root = RAI -end - - """ flux( z1, @@ -510,9 +323,167 @@ function effective_saturation(ν::FT, ϑ_l::FT) where {FT} end """ - make_compute_exp_tendency(model::PlantHydraulicsModel, _) + PrescribedTranspiration{FT, F <: Function} <: AbstractTranspiration{FT} + +A concrete type used for dispatch when computing the transpiration +from the leaves, in the case where transpiration is prescribed. +""" +struct PrescribedTranspiration{FT, F <: Function} <: AbstractTranspiration{FT} + T::F +end + +function PrescribedTranspiration{FT}(T::Function) where {FT <: AbstractFloat} + return PrescribedTranspiration{FT, typeof(T)}(T) +end + +""" + transpiration_per_ground_area( + transpiration::PrescribedTranspiration{FT}, + Y, + p, + t, + )::FT where {FT} + +A method which computes the transpiration in meters/sec between the leaf +and the atmosphere, +in the case of a standalone plant hydraulics model with prescribed +transpiration rate. + +Transpiration should be per unit ground area, not per leaf area. +""" +function transpiration_per_ground_area( + transpiration::PrescribedTranspiration{FT}, + _, + _, + t, +)::FT where {FT} + return FT(transpiration.T(t)) # (m/s) +end + +""" + DiagnosticTranspiration{FT} <: AbstractTranspiration{FT} + +A concrete type used for dispatch in the case where transpiration is computed +diagnostically, as a function of prognostic variables and parameters, +and stored in `p` during the `update_aux!` step. +""" +struct DiagnosticTranspiration{FT} <: AbstractTranspiration{FT} end -A function which creates the compute_exp_tendency! function for the PlantHydraulicsModel. +""" + transpiration_per_ground_area(transpiration::DiagnosticTranspiration, Y, p, t) + +Returns the transpiration computed diagnostically using local conditions. +In this case, it just returns the value which was computed and stored in +the `aux` state during the update_aux! step. + +Transpiration should be per unit ground area, not per leaf area. +""" +function transpiration_per_ground_area( + transpiration::DiagnosticTranspiration, + Y, + p, + t, +) + @inbounds return p.canopy.conductance.transpiration +end + +# BIGLEAF HYDRAULICS MODEL + +""" + BigLeafHydraulicsModel{FT, PS, T} <: AbstractPlantHydraulicsModel{FT} + +Defines, and constructs instances of, the BigLeafHydraulicsModel type which is +used for simulation flux of water to/from the soil for the case of a big leaf +model, meaning that the canopy is treated as a single stem and single leaf +compartment. Water flows along roots of different depths, along the stem, to +the leaf, and ultimately being lost from the system through transpiration. Note +that the canopy height is specified as part of the BigLeafHydraulicsModel, along +with the area indices of the leaves, roots, and stems. + +This model can also be combined with the soil model using ClimaLand, in which +case the prognostic soil water content is used to determine root extraction, and +the transpiration is also computed diagnostically. In a global run with patches +of bare soil, you can "turn off" the canopy model (to get zero root extraction, +zero absorption, and emission, zero transpiration, and sensible heat flux from +the canopy) by setting the canopy heights to zero. + +Finally, the model can be used in Canopy standalone mode by prescribing +the soil matric potential at the root tips or flux in the roots. There is also the +option (intendend only for debugging) to use a prescribed transpiration rate. + +$(DocStringExtensions.FIELDS) +""" +struct BigLeafHydraulicsModel{FT, PS, T} <: AbstractPlantHydraulicsModel{FT} + "The height of the stem compartment, in meters" + h_stem::FT + "The height of the leaf compartment, in meters" + h_leaf::FT + "Parameters required by the BigLeafHydraulics model" + parameters::PS + "The transpiration mode, of type `AbstractTranspiration`" + transpiration::T +end + +function BigLeafHydraulicsModel{FT}(; + h_stem::FT, + h_leaf::FT, + parameters::PlantHydraulicsParameters{FT}, + transpiration::AbstractTranspiration{FT} = DiagnosticTranspiration{FT}(), +) where {FT} + args = (parameters, transpiration) + @assert h_stem >= 0 + @assert h_leaf >= 0 + return BigLeafHydraulicsModel{FT, typeof.(args)...}(h_stem, h_leaf, args...) +end + +""" + prognostic_vars(model::BigLeafHydraulicsModel) + +A function which returns the names of the prognostic variables of the +`BigLeafHydraulicsModel`. +""" +prognostic_vars(model::BigLeafHydraulicsModel) = (:ϑ_l,) + +""" + auxiliary_vars(model::BigLeafHydraulicsModel) + +A function which returns the names of the auxiliary variables of the +`BigLeafHydraulicsModel`, the transpiration stress factor `β` (unitless), +the water potential `ψ` (m), the volume flux*cross section `fa` (1/s), +and the volume flux*root cross section in the roots `fa_roots` (1/s), +where the cross section can be represented by an area index. +""" +auxiliary_vars(model::BigLeafHydraulicsModel) = + (:β, :ψ, :fa, :fa_roots, :area_index) + +""" + ClimaLand.prognostic_types(model::BigLeafHydraulicsModel{FT}) where {FT} + +Defines the prognostic types for the BigLeafHydraulicsModel. +""" +ClimaLand.prognostic_types(model::BigLeafHydraulicsModel{FT}) where {FT} = + (NTuple{2, FT},) +ClimaLand.prognostic_domain_names(::BigLeafHydraulicsModel) = (:surface,) + +""" + ClimaLand.auxiliary_types(model::BigLeafHydraulicsModel{FT}) where {FT} + +Defines the auxiliary types for the BigLeafHydraulicsModel. +""" +ClimaLand.auxiliary_types(model::BigLeafHydraulicsModel{FT}) where {FT} = ( + FT, + NTuple{2, FT}, + NTuple{2, FT}, + FT, + NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}}, +) +ClimaLand.auxiliary_domain_names(::BigLeafHydraulicsModel) = + (:surface, :surface, :surface, :surface, :surface) + +""" + make_compute_exp_tendency(model::BigLeafHydraulicsModel, _) + +A function which creates the compute_exp_tendency! function for the BigLeafHydraulicsModel. The compute_exp_tendency! function must comply with a rhs function of SciMLBase.jl. Below, `fa` denotes a flux multiplied by the relevant cross section @@ -528,35 +499,24 @@ To prevent dividing by zero, we change AI/(AI x dz)" to "AI/max(AI x dz, eps(FT))" """ function make_compute_exp_tendency( - model::PlantHydraulicsModel{FT}, + model::BigLeafHydraulicsModel{FT}, canopy, ) where {FT} function compute_exp_tendency!(dY, Y, p, t) area_index = p.canopy.hydraulics.area_index - n_stem = model.n_stem - n_leaf = model.n_leaf fa = p.canopy.hydraulics.fa fa_roots = p.canopy.hydraulics.fa_roots + labels = model.h_stem > 0 ? [:stem, :leaf] : [:leaf] - # Inside of a loop, we need to use a single dollar sign - # for indexing into Fields of Tuples in non broadcasted - # expressions, and two dollar signs for - # for broadcasted expressions using the macro @. - # field.:($index) .= value # works - # @ field.:($$index) = value # works - @inbounds for i in 1:(n_stem + n_leaf) + @inbounds for i in eachindex(labels) im1 = i - 1 - ip1 = i + 1 + compartment_type = labels[i] + dz = + compartment_type == :stem ? model.h_stem : + model.h_leaf - model.h_stem # To prevent dividing by zero, change AI/(AI x dz)" to # "AI/max(AI x dz, eps(FT))" - AIdz = - max.( - getproperty(area_index, model.compartment_labels[i]) .* ( - model.compartment_surfaces[ip1] - - model.compartment_surfaces[i] - ), - eps(FT), - ) + AIdz = max.(getproperty(area_index, labels[i]) * (dz), eps(FT)) if i == 1 @inbounds @. dY.canopy.hydraulics.ϑ_l.:($$i) = 1 / AIdz * (fa_roots - fa.:($$i)) @@ -573,15 +533,15 @@ end root_water_flux_per_ground_area!( fa::ClimaCore.Fields.Field, s::PrescribedSoil, - model::PlantHydraulicsModel{FT}, + model::BigLeafHydraulicsModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, ) where {FT} -A method which computes the water flux between the soil and the stem, via the roots, -and multiplied by the RAI, in the case of a model running without an integrated -soil model. +A method which computes the water flux between the soil and the stem, via the +roots, and multiplied by the RAI, in the case of a big leaf model running +without an integrated soil model. The returned flux is per unit ground area. This assumes that the stem compartment is the first element of `Y.canopy.hydraulics.ϑ_l`. @@ -589,12 +549,11 @@ is the first element of `Y.canopy.hydraulics.ϑ_l`. function root_water_flux_per_ground_area!( fa::ClimaCore.Fields.Field, s::PrescribedSoil, - model::PlantHydraulicsModel{FT}, + model::BigLeafHydraulicsModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, ) where {FT} - (; conductivity_model, root_distribution) = model.parameters area_index = p.canopy.hydraulics.area_index # We can index into a field of Tuple{FT} to extract a field of FT @@ -604,15 +563,18 @@ function root_water_flux_per_ground_area!( n_root_layers = length(root_depths) ψ_soil::FT = s.ψ(t) fa .= FT(0.0) + labels = model.h_stem > 0 ? [:stem, :leaf] : [:leaf] @inbounds for i in 1:n_root_layers - above_ground_area_index = - getproperty(area_index, model.compartment_labels[1]) - + above_ground_area_index = getproperty(area_index, labels[1]) + compartment = labels[i] + midpoint = + compartment == :stem ? model.h_stem / 2 : + model.h_stem + (model.h_leaf / 2) if i != n_root_layers @. fa += flux( root_depths[i], - model.compartment_midpoints[1], + midpoint, ψ_soil, ψ_base, hydraulic_conductivity(conductivity_model, ψ_soil), @@ -625,82 +587,74 @@ function root_water_flux_per_ground_area!( @. fa += flux( root_depths[i], - model.compartment_midpoints[1], + midpoint, ψ_soil, ψ_base, hydraulic_conductivity(conductivity_model, ψ_soil), hydraulic_conductivity(conductivity_model, ψ_base), ) * root_distribution(root_depths[i]) * - (model.compartment_surfaces[1] - root_depths[n_root_layers]) * + (model.h_leaf - root_depths[n_root_layers]) * (area_index.root + above_ground_area_index) / 2 end end end -""" - PrescribedTranspiration{FT, F <: Function} <: AbstractTranspiration{FT} -A concrete type used for dispatch when computing the transpiration -from the leaves, in the case where transpiration is prescribed. """ -struct PrescribedTranspiration{FT, F <: Function} <: AbstractTranspiration{FT} - T::F -end - -function PrescribedTranspiration{FT}(T::Function) where {FT <: AbstractFloat} - return PrescribedTranspiration{FT, typeof(T)}(T) -end + lai_consistency_check( + model::BigLeafHydraulicsModel{FT}, + ) where {FT} -""" - transpiration_per_ground_area( - transpiration::PrescribedTranspiration{FT}, - Y, - p, - t, - )::FT where {FT} +Carries out consistency checks using the area indices supplied and the heights +of the stem and leaf compartments. -A method which computes the transpiration in meters/sec between the leaf -and the atmosphere, -in the case of a standalone plant hydraulics model with prescribed -transpiration rate. +Note that it is possible to have a plant with 0 height stem compartment +but with a leaf compartment, and that a plant must have leaf compartments +(even if LAI = 0). -Transpiration should be per unit ground area, not per leaf area. +Specifically, this checks that: +1. if LAI is nonzero or SAI is nonzero, RAI must be nonzero. +2. if SAI > 0, h_stem must be > 0 for consistency. If SAI == 0, h_stem must +be zero. """ -function transpiration_per_ground_area( - transpiration::PrescribedTranspiration{FT}, - _, - _, - t, -)::FT where {FT} - return FT(transpiration.T(t)) # (m/s) +function lai_consistency_check( + model::BigLeafHydraulicsModel{FT}, + area_index::NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}}, +) where {FT} + if area_index.leaf > eps(FT) || area_index.stem > eps(FT) + @assert area_index.root > eps(FT) + end + # If the SAI is 0, then the stem compartment should be height 0 + if area_index.stem < eps(FT) + @assert model.h_stem == FT(0) + else + # if SAI is > 0, h_stem should be > 0 for consistency + @assert model.h_stem > 0 + end end -""" - DiagnosticTranspiration{FT} <: AbstractTranspiration{FT} - -A concrete type used for dispatch in the case where transpiration is computed -diagnostically, as a function of prognostic variables and parameters, -and stored in `p` during the `update_aux!` step. -""" -struct DiagnosticTranspiration{FT} <: AbstractTranspiration{FT} end +# SHARED METHODS FOR BigLeafHydraulics AND PlantHydraulics """ - transpiration_per_ground_area(transpiration::DiagnosticTranspiration, Y, p, t) + set_canopy_prescribed_field!(component::PlantHydraulics{FT}, + p, + t, + ) where {FT} -Returns the transpiration computed diagnostically using local conditions. -In this case, it just returns the value which was computed and stored in -the `aux` state during the update_aux! step. -Transpiration should be per unit ground area, not per leaf area. +Sets the canopy prescribed fields pertaining to the PlantHydraulics +component (the area indices) with their values at time t. """ -function transpiration_per_ground_area( - transpiration::DiagnosticTranspiration, - Y, +function ClimaLand.Canopy.set_canopy_prescribed_field!( + component::BigLeafHydraulicsModel{FT}, p, t, -) - @inbounds return p.canopy.conductance.transpiration +) where {FT} + (; LAIfunction, SAI, RAI) = component.parameters.ai_parameterization + evaluate!(p.canopy.hydraulics.area_index.leaf, LAIfunction, t) + @. p.canopy.hydraulics.area_index.stem = SAI + @. p.canopy.hydraulics.area_index.root = RAI end end diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl index 5199989ddf..690e22aa9b 100644 --- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl +++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl @@ -16,7 +16,7 @@ model. See Cowan 1968; Brutsaert 1982, pp. 113–116; Campbell and Norman 1998, p. 71; Shuttleworth 2012, p. 343; Monteith and Unsworth 2013, p. 304. """ function ClimaLand.displacement_height(model::CanopyModel{FT}, Y, p) where {FT} - return FT(0.67) * model.hydraulics.compartment_surfaces[end] + return FT(0.67) * (model.hydraulics.h_stem + model.hydraulics.h_leaf) end """ @@ -68,8 +68,8 @@ end A helper function which returns the surface height for the canopy model, which is stored in the parameter struct. """ -function ClimaLand.surface_height(model::CanopyModel, _...) - return model.hydraulics.compartment_surfaces[1] +function ClimaLand.surface_height(model::CanopyModel{FT}, _...) where {FT} + return FT(0.0) end """ @@ -110,7 +110,7 @@ end <:Union{BeerLambertModel, TwoStreamModel}, <:Union{FarquharModel,OptimalityFarquharModel}, <:MedlynConductanceModel, - <:PlantHydraulicsModel, + <:BigLeafHydraulicsModel, <:Union{PrescribedCanopyTempModel,BigLeafEnergyModel} }, radiation::PrescribedRadiativeFluxes, @@ -141,7 +141,7 @@ function canopy_boundary_fluxes!( <:Union{BeerLambertModel, TwoStreamModel}, <:Union{FarquharModel, OptimalityFarquharModel}, <:MedlynConductanceModel, - <:PlantHydraulicsModel, + <:BigLeafHydraulicsModel, <:Union{PrescribedCanopyTempModel, BigLeafEnergyModel}, }, radiation::PrescribedRadiativeFluxes, @@ -157,7 +157,7 @@ function canopy_boundary_fluxes!( shf = p.canopy.energy.shf lhf = p.canopy.energy.lhf r_ae = p.canopy.energy.r_ae - i_end = canopy.hydraulics.n_stem + canopy.hydraulics.n_leaf + i_end = (canopy.hydraulics.h_stem > 0 ? 1 : 0) + 1 # Compute transpiration, SHF, LHF canopy_turbulent_fluxes = turbulent_fluxes(atmos, canopy, Y, p, t) diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl index bc2e8958d9..b0d9c2d065 100644 --- a/test/standalone/Vegetation/canopy_model.jl +++ b/test/standalone/Vegetation/canopy_model.jl @@ -143,38 +143,15 @@ for FT in (Float32, Float64) conductivity_model = conductivity_model, retention_model = retention_model, ) - Δz = FT(1.0) # height of compartments - n_stem = Int64(0) # number of stem elements - n_leaf = Int64(1) # number of leaf elements - compartment_centers = - FT.( - Vector( - range( - start = Δz / 2, - step = Δz, - stop = Δz * (n_stem + n_leaf) - (Δz / 2), - ), - ), - ) - compartment_faces = - FT.( - Vector( - range( - start = 0.0, - step = Δz, - stop = Δz * (n_stem + n_leaf), - ), - ) - ) + h_stem = FT(0.0) + h_leaf = FT(1.0) ψ_soil0 = FT(0.0) soil_driver = PrescribedSoil(FT) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + plant_hydraulics = PlantHydraulics.BigLeafHydraulicsModel{FT}(; parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, + h_stem = h_stem, + h_leaf = h_leaf, ) canopy = ClimaLand.Canopy.CanopyModel{FT}(; parameters = shared_params, @@ -332,7 +309,7 @@ for FT in (Float32, Float64) ) < 0.5 @test ClimaLand.surface_evaporative_scaling(canopy, Y, p) == FT(1.0) - @test ClimaLand.surface_height(canopy, Y, p) == compartment_faces[1] + @test ClimaLand.surface_height(canopy, Y, p) == h_stem T_sfc = FT.(T_atmos(t0)) @test Array(parent(ClimaLand.surface_temperature(canopy, Y, p, t0))) == [T_sfc] @@ -403,36 +380,13 @@ for FT in (Float32, Float64) conductivity_model = conductivity_model, retention_model = retention_model, ) - Δz = FT(1.0) # height of compartments - n_stem = Int64(0) # number of stem elements - n_leaf = Int64(1) # number of leaf elements - compartment_centers = - FT.( - Vector( - range( - start = Δz / 2, - step = Δz, - stop = Δz * (n_stem + n_leaf) - (Δz / 2), - ), - ), - ) - compartment_faces = - FT.( - Vector( - range( - start = 0.0, - step = Δz, - stop = Δz * (n_stem + n_leaf), - ), - ) - ) + h_stem = FT(0.0) + h_leaf = FT(1.0) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + plant_hydraulics = PlantHydraulics.BigLeafHydraulicsModel{FT}(; parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, + h_stem = h_stem, + h_leaf = h_leaf, ) t0 = FT(100) @@ -621,29 +575,8 @@ for FT in (Float32, Float64) conductivity_model = conductivity_model, retention_model = retention_model, ) - Δz = FT(1.0) # height of compartments - n_stem = Int64(0) # number of stem elements - n_leaf = Int64(1) # number of leaf elements - compartment_centers = - FT.( - Vector( - range( - start = Δz / 2, - step = Δz, - stop = Δz * (n_stem + n_leaf) - (Δz / 2), - ), - ), - ) - compartment_faces = - FT.( - Vector( - range( - start = 0.0, - step = Δz, - stop = Δz * (n_stem + n_leaf), - ), - ) - ) + h_stem = FT(0.0) + h_leaf = FT(1.0) ψ_soil0 = FT(0.0) T_soil0 = FT(290) @@ -656,12 +589,10 @@ for FT in (Float32, Float64) FT(0.98), ) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + plant_hydraulics = PlantHydraulics.BigLeafHydraulicsModel{FT}(; parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, + h_stem = h_stem, + h_leaf = h_leaf, ) autotrophic_parameters = AutotrophicRespirationParameters(FT) autotrophic_respiration_model = diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl index 341aff58d4..267669dce8 100644 --- a/test/standalone/Vegetation/plant_hydraulics_test.jl +++ b/test/standalone/Vegetation/plant_hydraulics_test.jl @@ -20,42 +20,42 @@ for FT in (Float32, Float64) LAI = FT.([0, 1]) SAI = FT.([0, 1]) RAI = FT.([0, 1]) - n_stem = [0, 1] - n_leaf = [1] + h_stem = FT.([0, 1]) + h_leaf = FT.([1]) for L in LAI for S in SAI for R in RAI - for n_s in n_stem - for n_l in n_leaf + for h_s in h_stem + for h_l in h_leaf area_index = (root = R, stem = S, leaf = L) if n_l == 0 @test_throws AssertionError ClimaLand.Canopy.PlantHydraulics.lai_consistency_check( - n_s, - n_l, + h_s, + h_l, area_index, ) end if (L > eps(FT) || S > eps(FT)) && R < eps(FT) @test_throws AssertionError ClimaLand.Canopy.PlantHydraulics.lai_consistency_check( - n_s, - n_l, + h_s, + h_l, area_index, ) end if S > FT(0) && n_s == 0 @test_throws AssertionError ClimaLand.Canopy.PlantHydraulics.lai_consistency_check( - n_s, - n_l, + h_s, + h_l, area_index, ) end if S < eps(FT) && n_s > 0 @test_throws AssertionError ClimaLand.Canopy.PlantHydraulics.lai_consistency_check( - n_s, - n_l, + h_s, + h_l, area_index, ) end @@ -101,309 +101,6 @@ for FT in (Float32, Float64) ) end - @testset "Plant hydraulics model integration tests, FT = $FT" begin - domains = [ - Point(; z_sfc = FT(0.0)), - Plane(; - xlim = FT.((0, 1)), - ylim = FT.((0, 1)), - nelements = (2, 2), - periodic = (true, true), - npolynomial = 1, - ), - ] - - AR_params = AutotrophicRespirationParameters(FT) - RTparams = BeerLambertParameters(FT) - photosynthesis_params = FarquharParameters(FT, C3()) - stomatal_g_params = MedlynConductanceParameters(FT) - - AR_model = AutotrophicRespirationModel{FT}(AR_params) - stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) - photosynthesis_model = FarquharModel{FT}(photosynthesis_params) - rt_model = BeerLambertModel{FT}(RTparams) - - earth_param_set = LP.LandParameters(FT) - LAI = (t) -> 1.0 # m2 [leaf] m-2 [ground] - z_0m = FT(2.0) # m, Roughness length for momentum - z_0b = FT(0.1) # m, Roughness length for scalars - h_canopy = FT(20.0) # m, canopy height - h_sfc = FT(20.0) # m, canopy height - h_int = FT(30.0) # m, "where measurements would be taken at a typical flux tower of a 20m canopy" - shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( - z_0m, - z_0b, - earth_param_set, - ) - lat = FT(0.0) # degree - long = FT(-180) # degree - - function zenith_angle( - t, - ref_time; - latitude = lat, - longitude = long, - insol_params = earth_param_set.insol_params, - ) - current_datetime = ref_time + Dates.Second(round(t)) - d, δ, η_UTC = - FT.( - Insolation.helper_instantaneous_zenith_angle( - current_datetime, - ref_time, - insol_params, - ) - ) - return Insolation.instantaneous_zenith_angle( - d, - δ, - η_UTC, - longitude, - latitude, - )[1] - end - - function shortwave_radiation( - t; - latitude = lat, - longitude = long, - insol_params = earth_param_set.insol_params, - ) - return 1000 # W/m^2 - end - - function longwave_radiation(t) - return 200 # W/m^2 - end - - u_atmos = t -> 10 #m.s-1 - - liquid_precip = (t) -> 0 # m - snow_precip = (t) -> 0 # m - T_atmos = t -> 290 # Kelvin - q_atmos = t -> 0.001 # kg/kg - P_atmos = t -> 1e5 # Pa - h_atmos = h_int # m - c_atmos = (t) -> 4.11e-4 # mol/mol - ref_time = DateTime(2005) - atmos = PrescribedAtmosphere( - TimeVaryingInput(liquid_precip), - TimeVaryingInput(snow_precip), - TimeVaryingInput(T_atmos), - TimeVaryingInput(u_atmos), - TimeVaryingInput(q_atmos), - TimeVaryingInput(P_atmos), - ref_time, - h_atmos, - earth_param_set; - c_co2 = TimeVaryingInput(c_atmos), - ) - radiation = PrescribedRadiativeFluxes( - FT, - TimeVaryingInput(shortwave_radiation), - TimeVaryingInput(longwave_radiation), - ref_time; - θs = zenith_angle, - ) - Δz = FT(1.0) # height of compartments - n_stem = Int64(5) # number of stem elements - n_leaf = Int64(6) # number of leaf elements - SAI = FT(1) # m2/m2 - RAI = FT(1) # m2/m2 - ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( - TimeVaryingInput(LAI), - SAI, - RAI, - ) - K_sat_plant = 1.8e-8 # m/s. - ψ63 = FT(-4 / 0.0098) # / MPa to m - Weibull_param = FT(4) # unitless - a = FT(0.05 * 0.0098) # 1/m - plant_ν = FT(0.7) # m3/m3 - plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m - conductivity_model = - PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) - retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) - root_depths = FT.(-Array{FT}(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth - function root_distribution(z::T) where {T} - return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m) - end - compartment_midpoints = Vector{FT}( - range( - start = Δz / 2, - step = Δz, - stop = Δz * (n_stem + n_leaf) - (Δz / 2), - ), - ) - - compartment_surfaces = Vector{FT}( - range(start = 0.0, step = Δz, stop = Δz * (n_stem + n_leaf)), - ) - - param_set = PlantHydraulics.PlantHydraulicsParameters(; - ai_parameterization = ai_parameterization, - ν = plant_ν, - S_s = plant_S_s, - root_distribution = root_distribution, - conductivity_model = conductivity_model, - retention_model = retention_model, - ) - - function leaf_transpiration(t) - T = FT(1e-8) # m/s - end - - ψ_soil0 = FT(0.0) - transpiration = PrescribedTranspiration{FT}(leaf_transpiration) - - soil_driver = PrescribedSoil(FT) - - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; - parameters = param_set, - transpiration = transpiration, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_surfaces, - compartment_midpoints = compartment_midpoints, - ) - autotrophic_parameters = AutotrophicRespirationParameters(FT) - autotrophic_respiration_model = - AutotrophicRespirationModel{FT}(autotrophic_parameters) - for domain in domains - model = ClimaLand.Canopy.CanopyModel{FT}(; - parameters = shared_params, - domain = domain, - autotrophic_respiration = autotrophic_respiration_model, - radiative_transfer = rt_model, - photosynthesis = photosynthesis_model, - conductance = stomatal_model, - hydraulics = plant_hydraulics, - soil_driver = soil_driver, - atmos = atmos, - radiation = radiation, - ) - # Set system to hydrostatic equilibrium state by setting fluxes to zero, and setting LHS of both ODEs to 0 - function initial_compute_exp_tendency!(F, Y) - AI = (; leaf = LAI(0.0), root = RAI, stem = SAI) - T0A = FT(1e-8) * AI[:leaf] - for i in 1:(n_leaf + n_stem) - if i == 1 - fa = - sum( - flux.( - root_depths, - plant_hydraulics.compartment_midpoints[i], - ψ_soil0, - Y[i], - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - ψ_soil0, - ), - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - Y[i], - ), - ) .* root_distribution.(root_depths) .* ( - vcat(root_depths, [0.0])[2:end] - - vcat(root_depths, [0.0])[1:(end - 1)] - ), - ) * AI[:root] - else - fa = - flux( - plant_hydraulics.compartment_midpoints[i - 1], - plant_hydraulics.compartment_midpoints[i], - Y[i - 1], - Y[i], - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - ψ_soil0, - ), - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - Y[i], - ), - ) * ( - AI[plant_hydraulics.compartment_labels[1]] + - AI[plant_hydraulics.compartment_labels[2]] - ) / 2 - end - F[i] = fa - T0A - end - end - - soln = nlsolve( - initial_compute_exp_tendency!, - Vector{FT}(-0.03:0.01:0.07); - ftol = 1e-11, - ) - - S_l = - inverse_water_retention_curve.( - retention_model, - soln.zero, - plant_ν, - plant_S_s, - ) - - ϑ_l_0 = augmented_liquid_fraction.(plant_ν, S_l) - - Y, p, coords = initialize(model) - if typeof(domain) <: ClimaLand.Domains.Point - @test propertynames(p) == (:canopy, :drivers) - elseif typeof(domain) <: ClimaLand.Domains.Plane - @test propertynames(p) == (:canopy, :dss_buffer_2d, :drivers) - @test typeof(p.dss_buffer_2d) == typeof( - ClimaCore.Spaces.create_dss_buffer( - ClimaCore.Fields.zeros(domain.space.surface), - ), - ) - end - - dY = similar(Y) - for i in 1:(n_stem + n_leaf) - Y.canopy.hydraulics.ϑ_l.:($i) .= ϑ_l_0[i] - p.canopy.hydraulics.ψ.:($i) .= NaN - p.canopy.hydraulics.fa.:($i) .= NaN - dY.canopy.hydraulics.ϑ_l.:($i) .= NaN - end - set_initial_cache! = make_set_initial_cache(model) - set_initial_cache!(p, Y, 0.0) - canopy_exp_tendency! = make_exp_tendency(model) - canopy_exp_tendency!(dY, Y, p, 0.0) - - m = similar(dY.canopy.hydraulics.ϑ_l.:1) - m .= FT(0) - for i in 1:(n_stem + n_leaf) - @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0) - end - @test maximum(parent(m)) < 200 * eps(FT) # starts in equilibrium - - - # repeat using the plant hydraulics model directly - # make sure it agrees with what we get when use the canopy model ODE - Y, p, coords = initialize(model) - standalone_dY = similar(Y) - for i in 1:(n_stem + n_leaf) - Y.canopy.hydraulics.ϑ_l.:($i) .= ϑ_l_0[i] - p.canopy.hydraulics.ψ.:($i) .= NaN - p.canopy.hydraulics.fa.:($i) .= NaN - standalone_dY.canopy.hydraulics.ϑ_l.:($i) .= NaN - end - set_initial_cache!(p, Y, 0.0) - standalone_exp_tendency! = - make_compute_exp_tendency(model.hydraulics, model) - standalone_exp_tendency!(standalone_dY, Y, p, 0.0) - - m = similar(dY.canopy.hydraulics.ϑ_l.:1) - m .= FT(0) - for i in 1:(n_stem + n_leaf) - @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0) - end - @test maximum(parent(m)) < 10^2.1 * eps(FT) # starts in equilibrium - end - end - @testset "No plant, FT = $FT" begin domain = Point(; z_sfc = FT(0.0)) @@ -500,8 +197,8 @@ for FT in (Float32, Float64) θs = zenith_angle, ) - n_stem = Int64(0) # number of stem elements - n_leaf = Int64(1) # number of leaf elements + h_stem = FT(0.0) # m + h_leaf = FT(0.0) # m SAI = FT(0) # m2/m2 RAI = FT(0) # m2/m2 lai_fun = t -> 0 @@ -523,8 +220,6 @@ for FT in (Float32, Float64) function root_distribution(z::T) where {T} return T(0) # (1/m) end - compartment_midpoints = [h_canopy] - compartment_surfaces = [FT(0.0), h_canopy] param_set = PlantHydraulics.PlantHydraulicsParameters(; ai_parameterization = ai_parameterization, @@ -537,13 +232,11 @@ for FT in (Float32, Float64) transpiration = DiagnosticTranspiration{FT}() soil_driver = PrescribedSoil(FT) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + plant_hydraulics = PlantHydraulics.BigLeafHydraulicsModel{FT}(; parameters = param_set, transpiration = transpiration, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_surfaces, - compartment_midpoints = compartment_midpoints, + h_stem = h_stem, + h_leaf = h_leaf, ) autotrophic_parameters = AutotrophicRespirationParameters(FT) @@ -565,7 +258,7 @@ for FT in (Float32, Float64) Y, p, coords = initialize(model) dY = similar(Y) - for i in 1:(n_stem + n_leaf) + for i in 1:(h_stem != 0 + 1) Y.canopy.hydraulics.ϑ_l.:($i) .= FT(0.1) p.canopy.hydraulics.ψ.:($i) .= NaN p.canopy.hydraulics.fa.:($i) .= NaN