diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl index 5afcd666af..4a27d246dd 100644 --- a/experiments/benchmarks/land.jl +++ b/experiments/benchmarks/land.jl @@ -365,11 +365,21 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15)) # Energy Balance model ac_canopy = FT(2.5e3) + #clm_data is used for g1 and vcmax maps + clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context) + # Conductance Model - g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa) + g1 = SpaceVaryingInput( + joinpath(clm_artifact_path, "vegetation_properties_map.nc"), + "medlynslope", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),), + ) #Photosynthesis model - clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context) # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s Vcmax25 = SpaceVaryingInput( joinpath(clm_artifact_path, "vegetation_properties_map.nc"), diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl index 760c1dadcd..46163bf2ee 100644 --- a/experiments/long_runs/land.jl +++ b/experiments/long_runs/land.jl @@ -365,11 +365,21 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15)) # Energy Balance model ac_canopy = FT(2.5e4) # this will likely be 10x smaller! + #clm_data is used for g1 and vcmax maps + clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context) + # Conductance Model - g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + # g1 is read in units of sqrt(kPa) and then converted to sqrt(Pa) + g1 = SpaceVaryingInput( + joinpath(clm_artifact_path, "vegetation_properties_map.nc"), + "medlynslope", + surface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data * 10^(3 / 2),), + ) #Photosynthesis model - clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context) # vcmax is read in units of umol CO2/m^2/s and then converted to mol CO2/m^2/s Vcmax25 = SpaceVaryingInput( joinpath(clm_artifact_path, "vegetation_properties_map.nc"), diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl index 18520b6901..f9abc5359b 100644 --- a/ext/CreateParametersExt.jl +++ b/ext/CreateParametersExt.jl @@ -535,7 +535,9 @@ function MedlynConductanceParameters( parameters = CP.get_parameter_values(toml_dict, name_map, "Land") FT = CP.float_type(toml_dict) - return MedlynConductanceParameters{FT}(; g1, parameters..., kwargs...) + g1 = FT.(g1) + G1 = typeof(g1) + return MedlynConductanceParameters{FT, G1}(; g1, parameters..., kwargs...) end """ diff --git a/src/standalone/Vegetation/stomatalconductance.jl b/src/standalone/Vegetation/stomatalconductance.jl index 7e3ba1eb27..d177f34bc5 100644 --- a/src/standalone/Vegetation/stomatalconductance.jl +++ b/src/standalone/Vegetation/stomatalconductance.jl @@ -9,13 +9,16 @@ abstract type AbstractStomatalConductanceModel{FT} <: The required parameters for the Medlyn stomatal conductance model. $(DocStringExtensions.FIELDS) """ -Base.@kwdef struct MedlynConductanceParameters{FT <: AbstractFloat} +Base.@kwdef struct MedlynConductanceParameters{ + FT <: AbstractFloat, + G1 <: Union{FT, ClimaCore.Fields.Field}, +} "Relative diffusivity of water vapor (unitless)" Drel::FT "Minimum stomatal conductance mol/m^2/s" g0::FT "Slope parameter, inversely proportional to the square root of marginal water use efficiency (Pa^{1/2})" - g1::FT + g1::G1 end Base.eltype(::MedlynConductanceParameters{FT}) where {FT} = FT diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl index daf4034d80..3a79dcf5d1 100644 --- a/test/standalone/Vegetation/canopy_model.jl +++ b/test/standalone/Vegetation/canopy_model.jl @@ -22,355 +22,379 @@ import ClimaParams @testset "Canopy software pipes" begin for FT in (Float32, Float64) domain = Point(; z_sfc = FT(0.0)) - # create new field with constant value everywhere - Vcmax25 = fill(FT(9e-5), domain.space.surface) - AR_params = AutotrophicRespirationParameters(FT) - RTparams = BeerLambertParameters(FT) - photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25 = Vcmax25) - 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 = FT(8.0) # m2 [leaf] m-2 [ground] - z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT - z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT - 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, + for g1 in (FT(790), fill(FT(790), domain.space.surface)) + # create new field with constant value everywhere + Vcmax25 = fill(FT(9e-5), domain.space.surface) + AR_params = AutotrophicRespirationParameters(FT) + RTparams = BeerLambertParameters(FT) + photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25) + stomatal_g_params = MedlynConductanceParameters(FT; g1) + + 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 = FT(8.0) # m2 [leaf] m-2 [ground] + z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT + z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT + 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 + 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 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 + 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, - ) + 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, + ) - # Plant Hydraulics - RAI = FT(1) - SAI = FT(0) - lai_fun = t -> LAI - ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( - TimeVaryingInput(lai_fun), - SAI, - RAI, - ) - K_sat_plant = FT(1.8e-8) # m/s - ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value - 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 - 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 = SVector{10, 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 - 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, - ) - Δ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), - ), - ), + # Plant Hydraulics + RAI = FT(1) + SAI = FT(0) + lai_fun = t -> LAI + ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( + TimeVaryingInput(lai_fun), + SAI, + RAI, ) - compartment_faces = - FT.( - Vector( - range( - start = 0.0, - step = Δz, - stop = Δz * (n_stem + n_leaf), + K_sat_plant = FT(1.8e-8) # m/s + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value + 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 + 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 = + SVector{10, 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 + 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, + ) + Δ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), + ), + ) + ) - ψ_soil0 = FT(0.0) - soil_driver = PrescribedSoil(FT) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; - parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, - ) - canopy = ClimaLand.Canopy.CanopyModel{FT}(; - parameters = shared_params, - domain = domain, - autotrophic_respiration = AR_model, - radiative_transfer = rt_model, - photosynthesis = photosynthesis_model, - conductance = stomatal_model, - hydraulics = plant_hydraulics, - soil_driver = soil_driver, - atmos = atmos, - radiation = radiation, - ) - drivers = ClimaLand.get_drivers(canopy) - @test drivers == (atmos, radiation) - Y, p, coords = ClimaLand.initialize(canopy) - @test propertynames(p.drivers) == ( - :P_liq, - :P_snow, - :T, - :P, - :u, - :q, - :c_co2, - :thermal_state, - :SW_d, - :LW_d, - :θs, - ) - # Check that structure of Y is value (will error if not) - @test !isnothing(zero(Y)) - @test typeof(canopy.energy) == PrescribedCanopyTempModel{FT} - @test propertynames(p) == (:canopy, :drivers) - for component in ClimaLand.Canopy.canopy_components(canopy) - # Only hydraulics has a prognostic variable - if component == :hydraulics - @test propertynames(getproperty(Y.canopy, component)) == - ClimaLand.prognostic_vars(getproperty(canopy, component)) + ψ_soil0 = FT(0.0) + soil_driver = PrescribedSoil(FT) + plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = param_set, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_faces, + compartment_midpoints = compartment_centers, + ) + canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = domain, + autotrophic_respiration = AR_model, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, + ) + drivers = ClimaLand.get_drivers(canopy) + @test drivers == (atmos, radiation) + Y, p, coords = ClimaLand.initialize(canopy) + @test propertynames(p.drivers) == ( + :P_liq, + :P_snow, + :T, + :P, + :u, + :q, + :c_co2, + :thermal_state, + :SW_d, + :LW_d, + :θs, + ) + # Check that structure of Y is value (will error if not) + @test !isnothing(zero(Y)) + @test typeof(canopy.energy) == PrescribedCanopyTempModel{FT} + @test propertynames(p) == (:canopy, :drivers) + for component in ClimaLand.Canopy.canopy_components(canopy) + # Only hydraulics has a prognostic variable + if component == :hydraulics + @test propertynames(getproperty(Y.canopy, component)) == + ClimaLand.prognostic_vars( + getproperty(canopy, component), + ) + end + @test propertynames(getproperty(p.canopy, component)) == + ClimaLand.auxiliary_vars(getproperty(canopy, component)) + @test getproperty(auxiliary_types(canopy), component) == + auxiliary_types(getproperty(canopy, component)) + @test getproperty(auxiliary_vars(canopy), component) == + auxiliary_vars(getproperty(canopy, component)) + @test getproperty(prognostic_types(canopy), component) == + prognostic_types(getproperty(canopy, component)) + @test getproperty(prognostic_types(canopy), component) == + prognostic_types(getproperty(canopy, component)) end - @test propertynames(getproperty(p.canopy, component)) == - ClimaLand.auxiliary_vars(getproperty(canopy, component)) - @test getproperty(auxiliary_types(canopy), component) == - auxiliary_types(getproperty(canopy, component)) - @test getproperty(auxiliary_vars(canopy), component) == - auxiliary_vars(getproperty(canopy, component)) - @test getproperty(prognostic_types(canopy), component) == - prognostic_types(getproperty(canopy, component)) - @test getproperty(prognostic_types(canopy), component) == - prognostic_types(getproperty(canopy, component)) - end - parent(Y.canopy.hydraulics.ϑ_l) .= plant_ν - set_initial_cache! = make_set_initial_cache(canopy) - exp_tendency! = make_exp_tendency(canopy) - t0 = FT(0.0) - dY = similar(Y) - set_initial_cache!(p, Y, t0) - # check that this is updated correctly: - # @test p.canopy.autotrophic_respiration.Ra == - exp_tendency!(dY, Y, p, t0) - turb_fluxes = ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) - - @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux - @test p.canopy.energy.shf == turb_fluxes.shf - @test p.canopy.energy.lhf == turb_fluxes.lhf - @test p.canopy.conductance.transpiration == turb_fluxes.vapor_flux - c = FT(LP.light_speed(earth_param_set)) - h = FT(LP.planck_constant(earth_param_set)) - N_a = FT(LP.avogadro_constant(earth_param_set)) - _σ = FT(LP.Stefan(earth_param_set)) - (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR) = canopy.radiative_transfer.parameters - APAR = p.canopy.radiative_transfer.par.abs - ANIR = p.canopy.radiative_transfer.nir.abs - energy_per_photon_PAR = h * c / λ_γ_PAR - energy_per_photon_NIR = h * c / λ_γ_NIR - - @test p.canopy.radiative_transfer.SW_n == - @. (energy_per_photon_PAR * N_a * APAR) + - (energy_per_photon_NIR * N_a * ANIR) - ϵ_canopy = p.canopy.radiative_transfer.ϵ - T_canopy = FT.(T_atmos(t0)) - T_soil = FT.(soil_driver.T(t0)) - ϵ_soil = FT.(soil_driver.ϵ) - LW_d = FT.(longwave_radiation(t0)) - LW_d_canopy = @. (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 - LW_u_soil = @. ϵ_soil * _σ * T_soil^4 + (1 - ϵ_soil) * LW_d_canopy - @test p.canopy.radiative_transfer.LW_n == @. ( - ϵ_canopy * LW_d - 2 * ϵ_canopy * _σ * T_canopy^4 + - ϵ_canopy * LW_u_soil - ) + parent(Y.canopy.hydraulics.ϑ_l) .= plant_ν + set_initial_cache! = make_set_initial_cache(canopy) + exp_tendency! = make_exp_tendency(canopy) + t0 = FT(0.0) + dY = similar(Y) + set_initial_cache!(p, Y, t0) + # check that this is updated correctly: + # @test p.canopy.autotrophic_respiration.Ra == + exp_tendency!(dY, Y, p, t0) + turb_fluxes = + ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) + + @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux + @test p.canopy.energy.shf == turb_fluxes.shf + @test p.canopy.energy.lhf == turb_fluxes.lhf + @test p.canopy.conductance.transpiration == turb_fluxes.vapor_flux + c = FT(LP.light_speed(earth_param_set)) + h = FT(LP.planck_constant(earth_param_set)) + N_a = FT(LP.avogadro_constant(earth_param_set)) + _σ = FT(LP.Stefan(earth_param_set)) + (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR) = + canopy.radiative_transfer.parameters + APAR = p.canopy.radiative_transfer.par.abs + ANIR = p.canopy.radiative_transfer.nir.abs + energy_per_photon_PAR = h * c / λ_γ_PAR + energy_per_photon_NIR = h * c / λ_γ_NIR + + @test p.canopy.radiative_transfer.SW_n == + @. (energy_per_photon_PAR * N_a * APAR) + + (energy_per_photon_NIR * N_a * ANIR) + ϵ_canopy = p.canopy.radiative_transfer.ϵ + T_canopy = FT.(T_atmos(t0)) + T_soil = FT.(soil_driver.T(t0)) + ϵ_soil = FT.(soil_driver.ϵ) + LW_d = FT.(longwave_radiation(t0)) + LW_d_canopy = @. (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 + LW_u_soil = @. ϵ_soil * _σ * T_soil^4 + (1 - ϵ_soil) * LW_d_canopy + @test p.canopy.radiative_transfer.LW_n == @. ( + ϵ_canopy * LW_d - 2 * ϵ_canopy * _σ * T_canopy^4 + + ϵ_canopy * LW_u_soil + ) - # Penman-monteith - Δ = FT(100 * (0.444017302 + (290 - 273.15) * 0.0286064092)) - Rn = FT(shortwave_radiation(t0)) - G = FT(0.0) - thermo_params = canopy.parameters.earth_param_set.thermo_params - ts_in = - Thermodynamics.PhaseEquil_pTq.( - thermo_params, - p.drivers.P, - p.drivers.T, - p.drivers.q, - ) - ρa = Thermodynamics.air_density.(thermo_params, ts_in) - cp = - FT(cp_m(thermo_params, Thermodynamics.PhasePartition.(q_atmos(t0)))) - - es = - Thermodynamics.saturation_vapor_pressure.( - Ref(thermo_params), - FT.(T_atmos(t0)), - Ref(Thermodynamics.Liquid()), + # Penman-monteith + Δ = FT(100 * (0.444017302 + (290 - 273.15) * 0.0286064092)) + Rn = FT(shortwave_radiation(t0)) + G = FT(0.0) + thermo_params = canopy.parameters.earth_param_set.thermo_params + ts_in = + Thermodynamics.PhaseEquil_pTq.( + thermo_params, + p.drivers.P, + p.drivers.T, + p.drivers.q, + ) + ρa = Thermodynamics.air_density.(thermo_params, ts_in) + cp = FT( + cp_m( + thermo_params, + Thermodynamics.PhasePartition.(q_atmos(t0)), + ), ) - ea = - Thermodynamics.partial_pressure_vapor.( - thermo_params, - FT(P_atmos(t0)), - Thermodynamics.PhasePartition.(FT.(q_atmos(t0))), - ) - - VPD = es .- ea - - conditions = turbulent_fluxes(atmos, canopy, Y, p, t0) #Per unit m^2 of leaf - r_ae = Array(parent(conditions.r_ae))[1] # s/m - ga = 1 / r_ae - γ = FT(66) - R = FT(LP.gas_constant(earth_param_set)) - gs = Array( - parent( - ClimaLand.Canopy.upscale_leaf_conductance.( - p.canopy.conductance.gs, - LAI, + + es = + Thermodynamics.saturation_vapor_pressure.( + Ref(thermo_params), FT.(T_atmos(t0)), - R, - FT.(P_atmos(t0)), + Ref(Thermodynamics.Liquid()), + ) + ea = + Thermodynamics.partial_pressure_vapor.( + thermo_params, + FT(P_atmos(t0)), + Thermodynamics.PhasePartition.(FT.(q_atmos(t0))), + ) + + VPD = es .- ea + + conditions = turbulent_fluxes(atmos, canopy, Y, p, t0) #Per unit m^2 of leaf + r_ae = Array(parent(conditions.r_ae))[1] # s/m + ga = 1 / r_ae + γ = FT(66) + R = FT(LP.gas_constant(earth_param_set)) + gs = Array( + parent( + ClimaLand.Canopy.upscale_leaf_conductance.( + p.canopy.conductance.gs, + LAI, + FT.(T_atmos(t0)), + R, + FT.(P_atmos(t0)), + ), ), - ), - )[1] - Lv = FT(2453e6) #J/m^3 - - ET = penman_monteith.( - Δ, # Rate of change of saturation specific humidity with air temperature. (Pa K−1) - Rn, # Net irradiance (W m−2), the external source of energy flux - G, # Ground heat flux (W m−2) - ρa, # Dry air density (kg m−3) - cp, # Specific heat capacity of air (J kg−1 K−1) - VPD, # vapor pressure deficit (Pa) - ga, # Conductivity of air, atmospheric conductance (m s−1) - γ, # Psychrometric constant (γ ≈ 66 Pa K−1) - gs, # Conductivity of stoma, surface or stomatal conductance (m s−1) - Lv, # Volumetric latent heat of vaporization. Energy required per water volume vaporized. (Lv = 2453 MJ m−3) - ) + )[1] + Lv = FT(2453e6) #J/m^3 + + ET = penman_monteith.( + Δ, # Rate of change of saturation specific humidity with air temperature. (Pa K−1) + Rn, # Net irradiance (W m−2), the external source of energy flux + G, # Ground heat flux (W m−2) + ρa, # Dry air density (kg m−3) + cp, # Specific heat capacity of air (J kg−1 K−1) + VPD, # vapor pressure deficit (Pa) + ga, # Conductivity of air, atmospheric conductance (m s−1) + γ, # Psychrometric constant (γ ≈ 66 Pa K−1) + gs, # Conductivity of stoma, surface or stomatal conductance (m s−1) + Lv, # Volumetric latent heat of vaporization. Energy required per water volume vaporized. (Lv = 2453 MJ m−3) + ) - @test abs( - (Array(parent(turb_fluxes.vapor_flux .- ET))[1]) / - Array(parent(turb_fluxes.vapor_flux))[1], - ) < 0.5 - - @test ClimaLand.surface_evaporative_scaling(canopy, Y, p) == FT(1.0) - @test ClimaLand.surface_height(canopy, Y, p) == compartment_faces[1] - T_sfc = FT.(T_atmos(t0)) - @test Array(parent(ClimaLand.surface_temperature(canopy, Y, p, t0))) == - [T_sfc] - @test ClimaLand.surface_temperature(canopy, Y, p, t0) isa - ClimaCore.Fields.Field - @test Array( - parent( - ClimaLand.Canopy.canopy_temperature( - canopy.energy, - canopy, - Y, - p, - t0, + @test abs( + (Array(parent(turb_fluxes.vapor_flux .- ET))[1]) / + Array(parent(turb_fluxes.vapor_flux))[1], + ) < 0.5 + + @test ClimaLand.surface_evaporative_scaling(canopy, Y, p) == FT(1.0) + @test ClimaLand.surface_height(canopy, Y, p) == compartment_faces[1] + T_sfc = FT.(T_atmos(t0)) + @test Array( + parent(ClimaLand.surface_temperature(canopy, Y, p, t0)), + ) == [T_sfc] + @test ClimaLand.surface_temperature(canopy, Y, p, t0) isa + ClimaCore.Fields.Field + @test Array( + parent( + ClimaLand.Canopy.canopy_temperature( + canopy.energy, + canopy, + Y, + p, + t0, + ), ), - ), - ) == [T_sfc] - @test ClimaLand.Canopy.canopy_temperature( - canopy.energy, - canopy, - Y, - p, - t0, - ) isa ClimaCore.Fields.Field - - ρ_sfc = - ClimaLand.surface_air_density(canopy.atmos, canopy, Y, p, t0, T_sfc) - @test ClimaLand.surface_specific_humidity(canopy, Y, p, T_sfc, ρ_sfc) == - Thermodynamics.q_vap_saturation_generic.( - Ref(thermo_params), - T_sfc, - ρ_sfc, - Ref(Thermodynamics.Liquid()), - ) - @test ρ_sfc == compute_ρ_sfc.(thermo_params, ts_in, T_sfc) + ) == [T_sfc] + @test ClimaLand.Canopy.canopy_temperature( + canopy.energy, + canopy, + Y, + p, + t0, + ) isa ClimaCore.Fields.Field + + ρ_sfc = ClimaLand.surface_air_density( + canopy.atmos, + canopy, + Y, + p, + t0, + T_sfc, + ) + @test ClimaLand.surface_specific_humidity( + canopy, + Y, + p, + T_sfc, + ρ_sfc, + ) == + Thermodynamics.q_vap_saturation_generic.( + Ref(thermo_params), + T_sfc, + ρ_sfc, + Ref(Thermodynamics.Liquid()), + ) + @test ρ_sfc == compute_ρ_sfc.(thermo_params, ts_in, T_sfc) + end end end @@ -510,441 +534,175 @@ end @testset "Canopy software pipes with energy model" begin for FT in (Float32, Float64) domain = Point(; z_sfc = FT(0.0)) - # create new field with constant value everywhere - Vcmax25 = fill(FT(9e-5), domain.space.surface) - RTparams = BeerLambertParameters(FT) - photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25 = Vcmax25) - stomatal_g_params = MedlynConductanceParameters(FT) - - stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) - photosynthesis_model = FarquharModel{FT}(photosynthesis_params) - rt_model = BeerLambertModel{FT}(RTparams) - energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}()) - earth_param_set = LP.LandParameters(FT) - LAI = FT(8.0) # m2 [leaf] m-2 [ground] - z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT - z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT - 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, + for g1 in (FT(790), fill(FT(790), domain.space.surface)) + # create new field with constant value everywhere + Vcmax25 = fill(FT(9e-5), domain.space.surface) + RTparams = BeerLambertParameters(FT) + photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25) + stomatal_g_params = MedlynConductanceParameters(FT; g1) + + stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) + photosynthesis_model = FarquharModel{FT}(photosynthesis_params) + rt_model = BeerLambertModel{FT}(RTparams) + energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}()) + earth_param_set = LP.LandParameters(FT) + LAI = FT(8.0) # m2 [leaf] m-2 [ground] + z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT + z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT + 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 + return Insolation.instantaneous_zenith_angle( + d, + δ, + η_UTC, + longitude, + latitude, + )[1] + end - function longwave_radiation(t) - return 200 # W/m^2 - end + function shortwave_radiation( + t; + latitude = lat, + longitude = long, + insol_params = earth_param_set.insol_params, + ) + return 1000 # 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, - ) + function longwave_radiation(t) + return 200 # W/m^2 + end - # Plant Hydraulics - RAI = FT(1) - SAI = FT(0) - lai_fun = t -> LAI - ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( - TimeVaryingInput(lai_fun), - SAI, - RAI, - ) - K_sat_plant = FT(1.8e-8) # m/s - ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value - 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 - 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 = SVector{10, 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 - 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, - ) - Δ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), - ), - ), + 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), ) - compartment_faces = - FT.( - Vector( - range( - start = 0.0, - step = Δz, - stop = Δz * (n_stem + n_leaf), - ), - ) + radiation = PrescribedRadiativeFluxes( + FT, + TimeVaryingInput(shortwave_radiation), + TimeVaryingInput(longwave_radiation), + ref_time; + θs = zenith_angle, ) - ψ_soil0 = FT(0.0) - T_soil0 = FT(290) - soil_driver = PrescribedSoil( - root_depths, - (t) -> ψ_soil0, - (t) -> T_soil0, - FT(0.2), - FT(0.4), - FT(0.98), - ) - - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; - parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, - ) - autotrophic_parameters = AutotrophicRespirationParameters(FT) - autotrophic_respiration_model = - AutotrophicRespirationModel{FT}(autotrophic_parameters) - - canopy = ClimaLand.Canopy.CanopyModel{FT}(; - parameters = shared_params, - domain = domain, - radiative_transfer = rt_model, - photosynthesis = photosynthesis_model, - conductance = stomatal_model, - autotrophic_respiration = autotrophic_respiration_model, - energy = energy_model, - hydraulics = plant_hydraulics, - soil_driver = soil_driver, - atmos = atmos, - radiation = radiation, - ) - # This test needs to match ClimaParams `canopy_emissivity` - @test canopy.radiative_transfer.parameters.ϵ_canopy == FT(0.97) - @test canopy.energy.parameters.ac_canopy == FT(2.0e3) - Y, p, coords = ClimaLand.initialize(canopy) - - # Check that structure of Y is value (will error if not) - @test !isnothing(zero(Y)) - @test propertynames(p) == (:canopy, :drivers) - for component in ClimaLand.Canopy.canopy_components(canopy) - # Only hydraulics has a prognostic variable - if component == :hydraulics - @test propertynames(getproperty(Y.canopy, component)) == - ClimaLand.prognostic_vars(getproperty(canopy, component)) + # Plant Hydraulics + RAI = FT(1) + SAI = FT(0) + lai_fun = t -> LAI + ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( + TimeVaryingInput(lai_fun), + SAI, + RAI, + ) + K_sat_plant = FT(1.8e-8) # m/s + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value + 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 + 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 = + SVector{10, 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 - @test propertynames(getproperty(p.canopy, component)) == - ClimaLand.auxiliary_vars(getproperty(canopy, component)) - - @test getproperty(auxiliary_types(canopy), component) == - auxiliary_types(getproperty(canopy, component)) - @test getproperty(auxiliary_vars(canopy), component) == - auxiliary_vars(getproperty(canopy, component)) - @test getproperty(prognostic_types(canopy), component) == - prognostic_types(getproperty(canopy, component)) - @test getproperty(prognostic_types(canopy), component) == - prognostic_types(getproperty(canopy, component)) - end - Y.canopy.hydraulics .= plant_ν - Y.canopy.energy.T = FT(289) - - set_initial_cache! = make_set_initial_cache(canopy) - t0 = FT(0.0) - set_initial_cache!(p, Y, t0) - exp_tendency! = make_exp_tendency(canopy) - dY = similar(Y) - exp_tendency!(dY, Y, p, t0) - turb_fluxes = ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) - @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux - @test p.canopy.energy.lhf == turb_fluxes.lhf - @test p.canopy.energy.shf == turb_fluxes.shf - @test all(Array(parent(p.canopy.energy.fa_energy_roots)) .== FT(0)) - - @test all( - Array(parent(ClimaLand.surface_temperature(canopy, Y, p, t0))) .== - FT(289), - ) - @test all( - Array( - parent( - ClimaLand.Canopy.canopy_temperature( - canopy.energy, - canopy, - Y, - p, - t0, + 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, + ) + Δ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), + ), ), - ), - ) .== FT(289), - ) - end -end - - -@testset "Zero LAI;" begin - for FT in (Float32, Float64) - domain = Point(; z_sfc = FT(0.0)) - # create new field with constant value everywhere - Vcmax25 = fill(FT(9e-5), domain.space.surface) - BeerLambertparams = BeerLambertParameters(FT) - # TwoStreamModel parameters - Ω = FT(0.69) - χl = FT(0.1) - G_Function = CLMGFunction(χl) - α_PAR_leaf = FT(0.1) - λ_γ_PAR = FT(5e-7) - λ_γ_NIR = FT(1.65e-6) - τ_PAR_leaf = FT(0.05) - α_NIR_leaf = FT(0.45) - τ_NIR_leaf = FT(0.25) - ϵ_canopy = FT(0.97) - TwoStreamparams = TwoStreamParameters( - FT; - Ω, - α_PAR_leaf, - τ_PAR_leaf, - α_NIR_leaf, - τ_NIR_leaf, - G_Function, - ) - photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25 = Vcmax25) - stomatal_g_params = MedlynConductanceParameters(FT) - - stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) - photosynthesis_model = FarquharModel{FT}(photosynthesis_params) - rt_models = ( - BeerLambertModel{FT}(BeerLambertparams), - TwoStreamModel{FT}(TwoStreamparams), - ) - energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}()) - earth_param_set = LP.LandParameters(FT) - LAI = FT(0.0) # m2 [leaf] m-2 [ground] - z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT - z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT - 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 = + ) + compartment_faces = FT.( - Insolation.helper_instantaneous_zenith_angle( - current_datetime, - ref_time, - insol_params, + Vector( + range( + start = 0.0, + step = Δz, + stop = Δz * (n_stem + n_leaf), + ), ) ) - 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, - ) - # Plant Hydraulics - RAI = FT(0) - SAI = FT(0) - lai_fun = t -> LAI - ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( - TimeVaryingInput(lai_fun), - SAI, - RAI, - ) - K_sat_plant = FT(1.8e-8) # m/s - ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value - 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 - 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 = SVector{10, 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 - 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, - ) - Δ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), - ), - ) + ψ_soil0 = FT(0.0) + T_soil0 = FT(290) + soil_driver = PrescribedSoil( + root_depths, + (t) -> ψ_soil0, + (t) -> T_soil0, + FT(0.2), + FT(0.4), + FT(0.98), ) - ψ_soil0 = FT(0.0) - T_soil0 = FT(290) - soil_driver = PrescribedSoil( - root_depths, - (t) -> ψ_soil0, - (t) -> T_soil0, - FT(0.2), - FT(0.4), - FT(0.98), - ) + plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = param_set, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_faces, + compartment_midpoints = compartment_centers, + ) + autotrophic_parameters = AutotrophicRespirationParameters(FT) + autotrophic_respiration_model = + AutotrophicRespirationModel{FT}(autotrophic_parameters) - plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; - parameters = param_set, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_surfaces = compartment_faces, - compartment_midpoints = compartment_centers, - ) - autotrophic_parameters = AutotrophicRespirationParameters(FT) - autotrophic_respiration_model = - AutotrophicRespirationModel{FT}(autotrophic_parameters) - for rt_model in rt_models canopy = ClimaLand.Canopy.CanopyModel{FT}(; parameters = shared_params, domain = domain, @@ -958,25 +716,302 @@ end atmos = atmos, radiation = radiation, ) - + # This test needs to match ClimaParams `canopy_emissivity` + @test canopy.radiative_transfer.parameters.ϵ_canopy == FT(0.97) + @test canopy.energy.parameters.ac_canopy == FT(2.0e3) Y, p, coords = ClimaLand.initialize(canopy) + + # Check that structure of Y is value (will error if not) + @test !isnothing(zero(Y)) + @test propertynames(p) == (:canopy, :drivers) + for component in ClimaLand.Canopy.canopy_components(canopy) + # Only hydraulics has a prognostic variable + if component == :hydraulics + @test propertynames(getproperty(Y.canopy, component)) == + ClimaLand.prognostic_vars( + getproperty(canopy, component), + ) + end + @test propertynames(getproperty(p.canopy, component)) == + ClimaLand.auxiliary_vars(getproperty(canopy, component)) + + @test getproperty(auxiliary_types(canopy), component) == + auxiliary_types(getproperty(canopy, component)) + @test getproperty(auxiliary_vars(canopy), component) == + auxiliary_vars(getproperty(canopy, component)) + @test getproperty(prognostic_types(canopy), component) == + prognostic_types(getproperty(canopy, component)) + @test getproperty(prognostic_types(canopy), component) == + prognostic_types(getproperty(canopy, component)) + end Y.canopy.hydraulics .= plant_ν Y.canopy.energy.T = FT(289) + set_initial_cache! = make_set_initial_cache(canopy) t0 = FT(0.0) set_initial_cache!(p, Y, t0) + exp_tendency! = make_exp_tendency(canopy) + dY = similar(Y) + exp_tendency!(dY, Y, p, t0) + turb_fluxes = + ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) + + @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux + @test p.canopy.energy.lhf == turb_fluxes.lhf + @test p.canopy.energy.shf == turb_fluxes.shf + @test all(Array(parent(p.canopy.energy.fa_energy_roots)) .== FT(0)) + + @test all( + Array( + parent(ClimaLand.surface_temperature(canopy, Y, p, t0)), + ) .== FT(289), + ) + @test all( + Array( + parent( + ClimaLand.Canopy.canopy_temperature( + canopy.energy, + canopy, + Y, + p, + t0, + ), + ), + ) .== FT(289), + ) + end + end +end + + +@testset "Zero LAI;" begin + for FT in (Float32, Float64) + domain = Point(; z_sfc = FT(0.0)) + for g1 in (FT(790), fill(FT(790), domain.space.surface)) + # create new field with constant value everywhere + Vcmax25 = fill(FT(9e-5), domain.space.surface) + BeerLambertparams = BeerLambertParameters(FT) + # TwoStreamModel parameters + Ω = FT(0.69) + χl = FT(0.1) + G_Function = CLMGFunction(χl) + α_PAR_leaf = FT(0.1) + λ_γ_PAR = FT(5e-7) + λ_γ_NIR = FT(1.65e-6) + τ_PAR_leaf = FT(0.05) + α_NIR_leaf = FT(0.45) + τ_NIR_leaf = FT(0.25) + ϵ_canopy = FT(0.97) + TwoStreamparams = TwoStreamParameters( + FT; + Ω, + α_PAR_leaf, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + G_Function, + ) + photosynthesis_params = FarquharParameters(FT, C3(); Vcmax25) + stomatal_g_params = MedlynConductanceParameters(FT; g1) + + stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) + photosynthesis_model = FarquharModel{FT}(photosynthesis_params) + rt_models = ( + BeerLambertModel{FT}(BeerLambertparams), + TwoStreamModel{FT}(TwoStreamparams), + ) + energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}()) + earth_param_set = LP.LandParameters(FT) + LAI = FT(0.0) # m2 [leaf] m-2 [ground] + z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT + z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT + 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, + ) + + # Plant Hydraulics + RAI = FT(0) + SAI = FT(0) + lai_fun = t -> LAI + ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( + TimeVaryingInput(lai_fun), + SAI, + RAI, + ) + K_sat_plant = FT(1.8e-8) # m/s + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value + 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 + 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 = + SVector{10, 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 + 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, + ) + Δ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), + ), + ) + ) - @test all(parent(p.canopy.hydraulics.fa.:1) .== FT(0)) - @test all(parent(p.canopy.energy.lhf) .== FT(0)) - @test all(parent(p.canopy.energy.shf) .== FT(0)) - @test all(parent(p.canopy.energy.fa_energy_roots) .== FT(0)) - @test all(parent(p.canopy.hydraulics.fa_roots) .== FT(0)) - @test all(parent(p.canopy.conductance.transpiration) .== FT(0)) - @test all(parent(p.canopy.radiative_transfer.LW_n) .== FT(0)) - @test all(parent(p.canopy.radiative_transfer.SW_n) .== FT(0)) - @test all(parent(p.canopy.radiative_transfer.par.abs) .== FT(0)) - @test all(parent(p.canopy.radiative_transfer.nir.abs) .== FT(0)) - @test all(parent(p.canopy.autotrophic_respiration.Ra) .== FT(0)) + ψ_soil0 = FT(0.0) + T_soil0 = FT(290) + soil_driver = PrescribedSoil( + root_depths, + (t) -> ψ_soil0, + (t) -> T_soil0, + FT(0.2), + FT(0.4), + FT(0.98), + ) + + plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = param_set, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_faces, + compartment_midpoints = compartment_centers, + ) + autotrophic_parameters = AutotrophicRespirationParameters(FT) + autotrophic_respiration_model = + AutotrophicRespirationModel{FT}(autotrophic_parameters) + for rt_model in rt_models + canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = domain, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + autotrophic_respiration = autotrophic_respiration_model, + energy = energy_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, + ) + + Y, p, coords = ClimaLand.initialize(canopy) + Y.canopy.hydraulics .= plant_ν + Y.canopy.energy.T = FT(289) + set_initial_cache! = make_set_initial_cache(canopy) + t0 = FT(0.0) + set_initial_cache!(p, Y, t0) + + @test all(parent(p.canopy.hydraulics.fa.:1) .== FT(0)) + @test all(parent(p.canopy.energy.lhf) .== FT(0)) + @test all(parent(p.canopy.energy.shf) .== FT(0)) + @test all(parent(p.canopy.energy.fa_energy_roots) .== FT(0)) + @test all(parent(p.canopy.hydraulics.fa_roots) .== FT(0)) + @test all(parent(p.canopy.conductance.transpiration) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.LW_n) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.SW_n) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.par.abs) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.nir.abs) .== FT(0)) + @test all(parent(p.canopy.autotrophic_respiration.Ra) .== FT(0)) + end end end end