From c96affcab4ab0989de59092f35cc43d1187f5c38 Mon Sep 17 00:00:00 2001 From: "Katherine M. Deck" Date: Wed, 3 Apr 2024 11:33:42 -0700 Subject: [PATCH] global run of soil/canopy model --- .buildkite/pipeline.yml | 6 + .buildkite/target/pipeline.yml | 10 + .../integrated/soil_canopy_tutorial.jl | 12 +- docs/tutorials/standalone/Soil/evaporation.jl | 9 +- .../Soil/evaporation_gilat_loess.jl | 14 +- .../standalone/Soil/freezing_front.jl | 14 +- .../standalone/Soil/soil_energy_hydrology.jl | 15 +- docs/tutorials/standalone/Soil/sublimation.jl | 9 +- experiments/benchmarks/land.jl | 676 ++++++++++++++++++ experiments/integrated/fluxnet/ozark_pft.jl | 13 +- experiments/integrated/fluxnet/run_fluxnet.jl | 12 +- .../integrated/global/global_parameters.jl | 211 ++++++ .../integrated/global/global_soil_canopy.jl | 444 ++++++++++++ .../conservation/ozark_conservation_setup.jl | 12 +- .../performance/profile_allocations.jl | 5 +- .../standalone/Biogeochemistry/experiment.jl | 14 +- ext/CreateParametersExt.jl | 22 +- .../src/Fluxnet/run_fluxnet.jl | 5 +- src/integrated/soil_canopy_model.jl | 28 +- src/standalone/Soil/energy_hydrology.jl | 112 ++- .../Soil/soil_heat_parameterizations.jl | 92 +-- .../Soil/soil_hydrology_parameterizations.jl | 21 +- src/standalone/Vegetation/PlantHydraulics.jl | 2 +- src/standalone/Vegetation/radiation.jl | 2 - .../soil_energy_hydrology_biogeochemistry.jl | 26 +- test/standalone/Soil/climate_drivers.jl | 9 +- test/standalone/Soil/soil_bc.jl | 6 +- .../standalone/Soil/soil_parameterizations.jl | 79 +- test/standalone/Soil/soil_test_3d.jl | 9 +- test/standalone/Soil/soiltest.jl | 40 +- 30 files changed, 1752 insertions(+), 177 deletions(-) create mode 100644 experiments/benchmarks/land.jl create mode 100644 experiments/integrated/global/global_parameters.jl create mode 100644 experiments/integrated/global/global_soil_canopy.jl diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 05b5b4495a..efee397207 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -109,6 +109,12 @@ steps: agents: slurm_mem: 16G + - label: "Global Run CPU" + command: "julia --color=yes --project=.buildkite experiments/integrated/global/global_soil_canopy.jl" + artifact_paths: "experiments/integrated/global/plots/*png" + agents: + slurm_mem: 16G + - group: "Experiments on GPU" steps: - label: "Richards Runoff GPU" diff --git a/.buildkite/target/pipeline.yml b/.buildkite/target/pipeline.yml index 7c2cf49512..28adb2ae93 100644 --- a/.buildkite/target/pipeline.yml +++ b/.buildkite/target/pipeline.yml @@ -45,3 +45,13 @@ steps: agents: slurm_mem: 8GB slurm_gpus: 1 + + - label: "Soil/Canopy" + command: "julia --color=yes --project=.buildkite experiments/benchmarks/land.jl" + artifact_paths: + - "land_benchmark_gpu/*html" + env: + CLIMACOMMS_DEVICE: CUDA + agents: + slurm_mem: 8GB + slurm_gpus: 1 diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl index 339ebb253e..1b0cb32710 100644 --- a/docs/tutorials/integrated/soil_canopy_tutorial.jl +++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl @@ -339,14 +339,14 @@ Y.soil.ϑ_l = FT(0.4) Y.soil.θ_i = FT(0.0) T_0 = FT(288.7) ρc_s = - volumetric_heat_capacity.(Y.soil.ϑ_l, Y.soil.θ_i, Ref(land.soil.parameters)) -Y.soil.ρe_int = - volumetric_internal_energy.( + volumetric_heat_capacity.( + Y.soil.ϑ_l, Y.soil.θ_i, - ρc_s, - T_0, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + earth_param_set, ) +Y.soil.ρe_int = + volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air diff --git a/docs/tutorials/standalone/Soil/evaporation.jl b/docs/tutorials/standalone/Soil/evaporation.jl index dd2053cf0b..4156c5d062 100644 --- a/docs/tutorials/standalone/Soil/evaporation.jl +++ b/docs/tutorials/standalone/Soil/evaporation.jl @@ -158,9 +158,14 @@ function init_soil!(Y, z, params) Y.soil.ϑ_l .= hydrostatic_equilibrium.(z, FT(-0.001), params) Y.soil.θ_i .= 0 T = FT(296.15) - ρc_s = @. Soil.volumetric_heat_capacity(Y.soil.ϑ_l, FT(0), params) + ρc_s = @. Soil.volumetric_heat_capacity( + Y.soil.ϑ_l, + FT(0), + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int = - Soil.volumetric_internal_energy.(FT(0), ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.(FT(0), ρc_s, T, params.earth_param_set) end init_soil!(Y, z, soil.parameters); diff --git a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl index 1605b02635..da8a0512b0 100644 --- a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl +++ b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl @@ -146,9 +146,19 @@ function init_soil!(Y, z, params) Y.soil.ϑ_l .= estimated_ic.(z) Y.soil.θ_i .= 0 T = FT(294.15) - ρc_s = @. Soil.volumetric_heat_capacity(Y.soil.ϑ_l, Y.soil.θ_i, params) + ρc_s = @. Soil.volumetric_heat_capacity( + Y.soil.ϑ_l, + Y.soil.θ_i, + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int = - Soil.volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T, + params.earth_param_set, + ) end init_soil!(Y, z, soil.parameters) diff --git a/docs/tutorials/standalone/Soil/freezing_front.jl b/docs/tutorials/standalone/Soil/freezing_front.jl index 4c27cb0bbd..987d60805d 100644 --- a/docs/tutorials/standalone/Soil/freezing_front.jl +++ b/docs/tutorials/standalone/Soil/freezing_front.jl @@ -183,9 +183,19 @@ function init_soil!(Ysoil, z, params) Ysoil.soil.ϑ_l .= FT(0.33) Ysoil.soil.θ_i .= FT(0.0) T = FT(279.85) - ρc_s = Soil.volumetric_heat_capacity(FT(0.33), FT(0.0), params) + ρc_s = Soil.volumetric_heat_capacity( + FT(0.33), + FT(0.0), + params.ρc_ds, + params.earth_param_set, + ) Ysoil.soil.ρe_int .= - Soil.volumetric_internal_energy.(FT(0.0), ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.( + FT(0.0), + ρc_s, + T, + params.earth_param_set, + ) end init_soil!(Y, coords.subsurface.z, soil.parameters); diff --git a/docs/tutorials/standalone/Soil/soil_energy_hydrology.jl b/docs/tutorials/standalone/Soil/soil_energy_hydrology.jl index ea322381c4..3a15500f9e 100644 --- a/docs/tutorials/standalone/Soil/soil_energy_hydrology.jl +++ b/docs/tutorials/standalone/Soil/soil_energy_hydrology.jl @@ -221,9 +221,20 @@ function init_soil!(Y, z, params) T = @.(T_min + (T_max - T_min) * exp(-(z - zmax) / (zmin - zmax) * c)) θ_l = Soil.volumetric_liquid_fraction.(Y.soil.ϑ_l, ν, θ_r) - ρc_s = Soil.volumetric_heat_capacity.(θ_l, Y.soil.θ_i, Ref(params)) + ρc_s = + Soil.volumetric_heat_capacity.( + θ_l, + Y.soil.θ_i, + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int .= - Soil.volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T, + params.earth_param_set, + ) end init_soil!(Y, coords.subsurface.z, soil.parameters); diff --git a/docs/tutorials/standalone/Soil/sublimation.jl b/docs/tutorials/standalone/Soil/sublimation.jl index 414f027a00..ad77993bbb 100644 --- a/docs/tutorials/standalone/Soil/sublimation.jl +++ b/docs/tutorials/standalone/Soil/sublimation.jl @@ -144,9 +144,14 @@ function init_soil!(Y, z, params) Y.soil.ϑ_l .= hydrostatic_equilibrium.(z, FT(-0.1), params) Y.soil.θ_i .= 0 T = FT(275.0) - ρc_s = @. Soil.volumetric_heat_capacity(Y.soil.ϑ_l, FT(0), params) + ρc_s = @. Soil.volumetric_heat_capacity( + Y.soil.ϑ_l, + FT(0), + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int = - Soil.volumetric_internal_energy.(FT(0), ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.(FT(0), ρc_s, T, params.earth_param_set) end init_soil!(Y, z, soil.parameters); diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl new file mode 100644 index 0000000000..5514a6a40f --- /dev/null +++ b/experiments/benchmarks/land.jl @@ -0,0 +1,676 @@ +# # Global run of land model + +# The code sets up and runs the soil/canopy model for 6 hours on a spherical domain, +# using ERA5 data. In this simulation, we have +# turned lateral flow off because horizontal boundary conditions and the +# land/sea mask are not yet supported by ClimaCore. + +# This code runs the model multiple times and collects statistics for execution time and +# allocations +# +# When run with buildkite on clima, this code also compares with the previous best time +# saved at the bottom of this file + +# Simulation Setup +# Number of spatial elements: 101 in horizontal, 15 in vertical +# Soil depth: 50 m +# Simulation duration: 6 hours +# Timestep: 180 s +# Timestepper: ARS343 +# Maximum iterations: 1 +# Convergence criterion: 1e-8 +# Jacobian update: Every Newton iteration +# Precipitation data update: every timestep +import SciMLBase +import ClimaComms +@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends +import ClimaTimeSteppers as CTS +using ClimaCore +using ClimaUtilities.ClimaArtifacts +import Interpolations +using Insolation + +import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput +import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput +import ClimaUtilities.Regridders: InterpolationsRegridder +import ClimaUtilities.ClimaArtifacts: @clima_artifact +import ClimaParams as CP + +using ClimaLand +using ClimaLand.Soil +using ClimaLand.Canopy +import ClimaLand +import ClimaLand.Parameters as LP + +using Statistics +using Dates +import NCDatasets +import Profile, ProfileCanvas + +const FT = Float64; + +regridder_type = :InterpolationsRegridder +context = ClimaComms.context() +device = ClimaComms.device() +device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu" + +outdir = "land_benchmark_$(device_suffix)" +!ispath(outdir) && mkpath(outdir) + +function setup_prob(t0, tf, Δt; nelements = (101, 15)) + + earth_param_set = LP.LandParameters(FT) + radius = FT(6378.1e3) + depth = FT(50) + domain = ClimaLand.Domains.SphericalShell(; + radius = radius, + depth = depth, + nelements = nelements, + npolynomial = 1, + dz_tuple = FT.((10.0, 0.1)), + ) + surface_space = domain.space.surface + subsurface_space = domain.space.subsurface + + ref_time = DateTime(2021) + t_start = 0.0 + + # Forcing data + era5_artifact_path = + ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context) # Precipitation: + precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "rf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), + ) + + snow_precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), + ) + + u_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "ws", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + q_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "q", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + P_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sp", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + + T_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "t2m", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + ) + h_atmos = FT(10) + + atmos = PrescribedAtmosphere( + precip, + snow_precip, + T_atmos, + u_atmos, + q_atmos, + P_atmos, + ref_time, + h_atmos, + earth_param_set, + ) + + # Prescribed radiation + SW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "ssrd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), + ) + LW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "strd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), + ) + + function zenith_angle( + t, + ref_time; + latitude = ClimaCore.Fields.coordinate_field(surface_space).lat, + longitude = ClimaCore.Fields.coordinate_field(surface_space).long, + insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params, + ) where {FT} + # This should be time in UTC + current_datetime = ref_time + Dates.Second(round(t)) + + # Orbital Data uses Float64, so we need to convert to our sim FT + d, δ, η_UTC = + FT.( + Insolation.helper_instantaneous_zenith_angle( + current_datetime, + ref_time, + insol_params, + ) + ) + + Insolation.instantaneous_zenith_angle.( + d, + δ, + η_UTC, + longitude, + latitude, + ).:1 + end + radiation = + PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle) + + soil_params_artifact_path = + ClimaLand.Artifacts.soil_params_artifact_folder_path(; context) + extrapolation_bc = ( + Interpolations.Periodic(), + Interpolations.Flat(), + Interpolations.Flat(), + ) + soil_params_mask = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data > 0,), + ) + oceans_to_value(field, mask, value) = + mask == 1.0 ? field : eltype(field)(value) + + vg_α = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + vg_n = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGn_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "n", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + x = parent(vg_α) + μ = mean(log10.(x[x .> 0])) + vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ) + + x = parent(vg_n) + μ = mean(x[x .> 0]) + vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ) + + vg_fields_to_hcm_field(α::FT, n::FT) where {FT} = + ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...) + hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n) + + θ_r = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "residual_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "θ_r", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + + ν = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "porosity_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "ν", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + K_sat = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "ksat_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "Ksat", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + ) + + x = parent(K_sat) + μ = mean(log10.(x[x .> 0])) + K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ) + + ν .= oceans_to_value.(ν, soil_params_mask, 1) + + θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0) + + + S_s = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3), + soil_params_mask, + 1, + ) + ν_ss_om = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + ν_ss_quartz = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + ν_ss_gravel = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) + PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) + NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) + soil_params = Soil.EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo = PAR_albedo, + NIR_albedo = NIR_albedo, + ) + + + # TwoStreamModel parameters + Ω = FT(0.69) + ld = FT(0.5) + α_PAR_leaf = FT(0.1) + τ_PAR_leaf = FT(0.05) + α_NIR_leaf = FT(0.45) + τ_NIR_leaf = FT(0.25) + + # Energy Balance model + ac_canopy = FT(2.5e5) + + # Conductance Model + g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + + #Photosynthesis model + Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5 + + # Plant Hydraulics and general plant parameters + SAI = FT(0.0) # m2/m2 + f_root_to_shoot = FT(3.5) + RAI = FT(1.0) + K_sat_plant = 5e-9 # m/s # seems much too small? + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa + Weibull_param = FT(4) # unitless, Holtzman's original c param value + a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity + conductivity_model = + Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a) + plant_ν = FT(1.44e-4) + plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m + rooting_depth = FT(0.5) # from Natan + n_stem = 0 + n_leaf = 1 + h_stem = 0.0 + h_leaf = 1.0 + zmax = 0.0 + h_canopy = h_stem + h_leaf + compartment_midpoints = + n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2] + compartment_surfaces = + n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] + + z0_m = FT(0.13) * h_canopy + z0_b = FT(0.1) * z0_m + + + soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters( + FT; + ν = 1.0,# INCORRECT! + ) + + soil_args = (domain = domain, parameters = soil_params) + soil_model_type = Soil.EnergyHydrology{FT} + + # Soil microbes model + soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT} + + # soil microbes args + Csom = (z, t) -> eltype(z)(5.0) + + # Set the soil CO2 BC to being atmospheric CO2 + soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC() + soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux + soilco2_sources = (Soil.Biogeochemistry.MicrobeProduction{FT}(),) + + soilco2_boundary_conditions = + (; top = soilco2_top_bc, bottom = soilco2_bot_bc) + + soilco2_drivers = Soil.Biogeochemistry.SoilDrivers( + Soil.Biogeochemistry.PrognosticMet{FT}(), + Soil.Biogeochemistry.PrescribedSOC{FT}(Csom), + atmos, + ) + + soilco2_args = (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = domain, + parameters = soilco2_ps, + drivers = soilco2_drivers, + ) + + # Now we set up the canopy model, which we set up by component: + # Component Types + canopy_component_types = (; + autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = Canopy.TwoStreamModel{FT}, + photosynthesis = Canopy.FarquharModel{FT}, + conductance = Canopy.MedlynConductanceModel{FT}, + hydraulics = Canopy.PlantHydraulicsModel{FT}, + energy = Canopy.BigLeafEnergyModel{FT}, + ) + # Individual Component arguments + # Set up autotrophic respiration + autotrophic_respiration_args = + (; parameters = Canopy.AutotrophicRespirationParameters(FT)) + # Set up radiative transfer + radiative_transfer_args = (; + parameters = Canopy.TwoStreamParameters( + FT; + Ω, + α_PAR_leaf, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ) + ) + # Set up conductance + conductance_args = + (; parameters = Canopy.MedlynConductanceParameters(FT; g1)) + # Set up photosynthesis + photosynthesis_args = (; + parameters = Canopy.FarquharParameters( + FT, + Canopy.C3(); + Vcmax25 = Vcmax25, + ) + ) + # Set up plant hydraulics + # Note that we clip all values of LAI below 0.05 to zero. + # This is because we currently run into issues when LAI is + # of order eps(FT) in the SW radiation code. + # Please see Issue #644 + # or PR #645 for details. + # For now, this clipping is similar to what CLM does. + LAIfunction = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "lai", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; + preprocess_func = (data) -> data > 0.05 ? data : 0.0, + ), + ) + ai_parameterization = + Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) + + function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m + end + + plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, + ) + plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_midpoints = compartment_midpoints, + compartment_surfaces = compartment_surfaces, + ) + + energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) + + # Canopy component args + canopy_component_args = (; + autotrophic_respiration = autotrophic_respiration_args, + radiative_transfer = radiative_transfer_args, + photosynthesis = photosynthesis_args, + conductance = conductance_args, + hydraulics = plant_hydraulics_args, + energy = energy_args, + ) + + # Other info needed + shared_params = Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, + ) + + canopy_model_args = (; + parameters = shared_params, + domain = ClimaLand.obtain_surface_domain(domain), + ) + + # Integrated plant hydraulics and soil model + land_input = (atmos = atmos, radiation = radiation) + land = SoilCanopyModel{FT}(; + soilco2_type = soilco2_type, + soilco2_args = soilco2_args, + land_args = land_input, + soil_model_type = soil_model_type, + soil_args = soil_args, + canopy_component_types = canopy_component_types, + canopy_component_args = canopy_component_args, + canopy_model_args = canopy_model_args, + ) + + Y, p, cds = initialize(land) + + init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2 + Y.soil.ϑ_l .= init_soil.(ν, θ_r) + Y.soil.θ_i .= FT(0.0) + T = FT(276.85) + ρc_s = + Soil.volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + soil_params.ρc_ds, + soil_params.earth_param_set, + ) + Y.soil.ρe_int .= + Soil.volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T, + soil_params.earth_param_set, + ) + Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air + Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν + evaluate!(Y.canopy.energy.T, atmos.T, t0) + + set_initial_cache! = make_set_initial_cache(land) + exp_tendency! = make_exp_tendency(land) + imp_tendency! = ClimaLand.make_imp_tendency(land) + tendency_jacobian! = ClimaLand.make_tendency_jacobian(land) + set_initial_cache!(p, Y, t0) + + # set up jacobian info + jac_kwargs = (; + jac_prototype = ImplicitEquationJacobian(Y), + Wfact = tendency_jacobian!, + ) + + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction( + T_exp! = exp_tendency!, + T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), + dss! = ClimaLand.dss!, + ), + Y, + (t0, tf), + p, + ) + + updateat = Array(t0:(3600 * 3):tf) + updatefunc = ClimaLand.make_update_drivers(atmos, radiation) + cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) + return prob, cb +end + +function setup_and_solve_problem(; greet = false) + # We profile the setup phase as well here. This is not intended, but it is the easiest + # to set up for both CPU/GPU at the same time + t0 = 0.0 + tf = 60 * 60.0 * 6 + Δt = 180.0 + nelements = (101, 15) + if greet + @info "Run: Global RichardsModel" + @info "Resolution: $nelements" + @info "Timestep: $Δt s" + @info "Duration: $(tf - t0) s" + end + + prob, cb = setup_prob(t0, tf, Δt; nelements) + + # Define timestepper and ODE algorithm + stepper = CTS.ARS343() + norm_condition = CTS.MaximumError(FT(1e-8)) + conv_checker = CTS.ConvergenceChecker(; norm_condition = norm_condition) + ode_algo = CTS.IMEXAlgorithm( + stepper, + CTS.NewtonsMethod( + max_iters = 1, + update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), + convergence_checker = conv_checker, + ), + ) + + SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb) + return nothing +end + +# Warm up and greet +setup_and_solve_problem(; greet = true) + +@info "Starting profiling" +# Stop when we profile for MAX_PROFILING_TIME_SECONDS or MAX_PROFILING_SAMPLES +MAX_PROFILING_TIME_SECONDS = 1000 +MAX_PROFILING_SAMPLES = 100 +time_now = time() +timings_s = Float64[] +while (time() - time_now) < MAX_PROFILING_TIME_SECONDS && + length(timings_s) < MAX_PROFILING_SAMPLES + push!(timings_s, ClimaComms.@elapsed device setup_and_solve_problem()) +end +num_samples = length(timings_s) +average_timing_s = round(sum(timings_s) / num_samples, sigdigits = 3) +max_timing_s = round(maximum(timings_s), sigdigits = 3) +min_timing_s = round(minimum(timings_s), sigdigits = 3) +std_timing_s = round( + sum(((timings_s .- average_timing_s) .^ 2) / num_samples), + sigdigits = 3, +) +@info "Num samples: $num_samples" +@info "Average time: $(average_timing_s) s" +@info "Max time: $(max_timing_s) s" +@info "Min time: $(min_timing_s) s" +@info "Standard deviation time: $(std_timing_s) s" +@info "Done profiling" + +Profile.@profile setup_and_solve_problem() +results = Profile.fetch() +flame_file = joinpath(outdir, "flame_$device_suffix.html") +ProfileCanvas.html_file(flame_file, results) +@info "Saved compute flame to $flame_file" + +Profile.Allocs.@profile sample_rate = 0.01 setup_and_solve_problem() +results = Profile.Allocs.fetch() +profile = ProfileCanvas.view_allocs(results) +alloc_flame_file = joinpath(outdir, "alloc_flame_$device_suffix.html") +ProfileCanvas.html_file(alloc_flame_file, profile) +@info "Saved allocation flame to $alloc_flame_file" + +#= +if ClimaComms.device() isa ClimaComms.CUDADevice + import CUDA + CUDA.@profile setup_and_solve_problem() +end + +if get(ENV, "BUILDKITE_PIPELINE_SLUG", nothing) == "climaland-benchmark" + PREVIOUS_BEST_TIME = 9.3 + if average_timing_s > PREVIOUS_BEST_TIME + std_timing_s + @info "Possible performance regression, previous average time was $(PREVIOUS_BEST_TIME)" + elseif average_timing_s < PREVIOUS_BEST_TIME - std_timing_s + @info "Possible significant performance improvement, please update PREVIOUS_BEST_TIME in $(@__DIR__)" + end + @testset "Performance" begin + @test PREVIOUS_BEST_TIME - std_timing_s <= + average_timing_s ≤ + PREVIOUS_BEST_TIME + std_timing_s + end +end +=# diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl index aca2e043ea..1d707910cc 100644 --- a/experiments/integrated/fluxnet/ozark_pft.jl +++ b/experiments/integrated/fluxnet/ozark_pft.jl @@ -271,14 +271,15 @@ T_0 = drivers.TS.values[1 + Int(round(t0 / DATA_DT))] : drivers.TA.values[1 + Int(round(t0 / DATA_DT))] + 40# Get soil temperature at t0 ρc_s = - volumetric_heat_capacity.(Y.soil.ϑ_l, Y.soil.θ_i, Ref(land.soil.parameters)) -Y.soil.ρe_int = - volumetric_internal_energy.( + volumetric_heat_capacity.( + Y.soil.ϑ_l, Y.soil.θ_i, - ρc_s, - T_0, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + earth_param_set, ) +Y.soil.ρe_int = + volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) + 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) diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index ffa74b5059..7f4b586869 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -231,14 +231,14 @@ T_0 = drivers.TS.values[1 + Int(round(t0 / DATA_DT))] : drivers.TA.values[1 + Int(round(t0 / DATA_DT))] + 40# Get soil temperature at t0 ρc_s = - volumetric_heat_capacity.(Y.soil.ϑ_l, Y.soil.θ_i, Ref(land.soil.parameters)) -Y.soil.ρe_int = - volumetric_internal_energy.( + volumetric_heat_capacity.( + Y.soil.ϑ_l, Y.soil.θ_i, - ρc_s, - T_0, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + earth_param_set, ) +Y.soil.ρe_int = + volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) 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) diff --git a/experiments/integrated/global/global_parameters.jl b/experiments/integrated/global/global_parameters.jl new file mode 100644 index 0000000000..2d02fa982b --- /dev/null +++ b/experiments/integrated/global/global_parameters.jl @@ -0,0 +1,211 @@ +using ArtifactWrappers +#currently not being used +#= +# Read in f_max data and land sea mask +topmodel_dataset = ArtifactWrapper( + @__DIR__, + "processed_topographic_index 2.5 deg", + ArtifactFile[ArtifactFile( + url = "https://caltech.box.com/shared/static/dwa7g0uzhxd50a2z3thbx3a12n0r887s.nc", + filename = "means_2.5_new.nc", + ),], +) +infile_path = joinpath(get_data_folder(topmodel_dataset), "means_2.5_new.nc") + +f_max = SpaceVaryingInput( + infile_path, + "fmax", + surface_space; + regridder_type +) +mask = SpaceVaryingInput( + infile_path, + "landsea_mask", + surface_space; + regridder_type +) + +oceans_to_zero(field, mask) = mask > 0.5 ? field : eltype(field)(0) +f_over = FT(3.28) # 1/m +R_sb = FT(1.484e-4 / 1000) # m/s +runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(; + f_over = f_over, + f_max = f_max, + R_sb = R_sb, +) +=# +soil_params_artifact_path = + ClimaLand.Artifacts.soil_params_artifact_folder_path(; context) +extrapolation_bc = + (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat()) +soil_params_mask = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), + file_reader_kwargs = (; preprocess_func = (data) -> data > 0,), +) +oceans_to_value(field, mask, value) = mask == 1.0 ? field : eltype(field)(value) + +vg_α = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +vg_n = SpaceVaryingInput( + joinpath(soil_params_artifact_path, "vGn_map_gupta_etal2020_2.5x2.5x4.nc"), + "n", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +x = parent(vg_α) +μ = mean(log10.(x[x .> 0])) +vg_α .= oceans_to_value.(vg_α, soil_params_mask, 10.0^μ) + +x = parent(vg_n) +μ = mean(x[x .> 0]) +vg_n .= oceans_to_value.(vg_n, soil_params_mask, μ) + +vg_fields_to_hcm_field(α::FT, n::FT) where {FT} = + ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...) +hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n) + +θ_r = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "residual_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "θ_r", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) + +ν = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "porosity_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "ν", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +K_sat = SpaceVaryingInput( + joinpath(soil_params_artifact_path, "ksat_map_gupta_etal2020_2.5x2.5x4.nc"), + "Ksat", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) + +x = parent(K_sat) +μ = mean(log10.(x[x .> 0])) +K_sat .= oceans_to_value.(K_sat, soil_params_mask, 10.0^μ) + +ν .= oceans_to_value.(ν, soil_params_mask, 1) + +θ_r .= oceans_to_value.(θ_r, soil_params_mask, 0) + + +S_s = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3), + soil_params_mask, + 1, + ) +ν_ss_om = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) +ν_ss_quartz = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) +ν_ss_gravel = + oceans_to_value.( + ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.1), + soil_params_mask, + 0, + ) +PAR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) +NIR_albedo = ClimaCore.Fields.zeros(surface_space) .+ FT(0.2) +soil_params = Soil.EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo = PAR_albedo, + NIR_albedo = NIR_albedo, +); + + +# TwoStreamModel parameters +Ω = FT(0.69) +ld = FT(0.5) +α_PAR_leaf = FT(0.1) +τ_PAR_leaf = FT(0.05) +α_NIR_leaf = FT(0.45) +τ_NIR_leaf = FT(0.25) + +# Energy Balance model +ac_canopy = FT(2.5e4) + +# Conductance Model +g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + +#Photosynthesis model +Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5 + +# Plant Hydraulics and general plant parameters +SAI = FT(0.0) # m2/m2 +f_root_to_shoot = FT(3.5) +RAI = FT(1.0) +K_sat_plant = 5e-9 # m/s # seems much too small? +ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa +Weibull_param = FT(4) # unitless, Holtzman's original c param value +a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity +conductivity_model = + Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) +retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a) +plant_ν = FT(1.44e-4) +plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m +rooting_depth = FT(0.5) # from Natan +n_stem = 0 +n_leaf = 1 +h_stem = 0.0 +h_leaf = 1.0 +zmax = 0.0 +h_canopy = h_stem + h_leaf +compartment_midpoints = + n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2] +compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] + +z0_m = FT(0.13) * h_canopy +z0_b = FT(0.1) * z0_m + + +soilco2_ps = Soil.Biogeochemistry.SoilCO2ModelParameters( + FT; + ν = 1.0,# INCORRECT! This should be the same as the soil porosity, but + # currently, SoilCO2 does not support spatially varying parameters +) diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl new file mode 100644 index 0000000000..ff5b65182c --- /dev/null +++ b/experiments/integrated/global/global_soil_canopy.jl @@ -0,0 +1,444 @@ +import SciMLBase +import ClimaComms +@static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends +import ClimaTimeSteppers as CTS +using ClimaCore +using ClimaUtilities.ClimaArtifacts +import Interpolations +using Insolation + +import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput +import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput +import ClimaUtilities.Regridders: InterpolationsRegridder +import ClimaUtilities.ClimaArtifacts: @clima_artifact +import ClimaParams as CP + +using ClimaLand +using ClimaLand.Soil +using ClimaLand.Canopy +import ClimaLand +import ClimaLand.Parameters as LP + +using CairoMakie +using Statistics +using Dates +import NCDatasets + +regridder_type = :InterpolationsRegridder +extrapolation_bc = + (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat()) +context = ClimaComms.context() +outdir = joinpath(pkgdir(ClimaLand), "experiments/integrated/global/plots") +!ispath(outdir) && mkpath(outdir) + +device_suffix = + typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? + "cpu" : "gpu" + +FT = Float64 +earth_param_set = LP.LandParameters(FT) +radius = FT(6378.1e3); +depth = FT(50) +nelements = (50, 10) +domain = ClimaLand.Domains.SphericalShell(; + radius = radius, + depth = depth, + nelements = nelements, + npolynomial = 1, + dz_tuple = FT.((10.0, 0.1)), +); +surface_space = domain.space.surface +subsurface_space = domain.space.subsurface + +ref_time = DateTime(2021); +t_start = 0.0 + +# Forcing data +era5_artifact_path = + ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context) +precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "rf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), +) + +snow_precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), +) + +u_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "ws", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +q_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "q", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +P_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sp", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) + +T_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "t2m", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +h_atmos = FT(10); + +atmos = PrescribedAtmosphere( + precip, + snow_precip, + T_atmos, + u_atmos, + q_atmos, + P_atmos, + ref_time, + h_atmos, + earth_param_set, +); + +# Prescribed radiation +SW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "ssrd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), +) +LW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "strd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), +) + +function zenith_angle( + t, + ref_time; + latitude = ClimaCore.Fields.coordinate_field(surface_space).lat, + longitude = ClimaCore.Fields.coordinate_field(surface_space).long, + insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params, +) where {FT} + # This should be time in UTC + current_datetime = ref_time + Dates.Second(round(t)) + + # Orbital Data uses Float64, so we need to convert to our sim FT + d, δ, η_UTC = + FT.( + Insolation.helper_instantaneous_zenith_angle( + current_datetime, + ref_time, + insol_params, + ) + ) + + Insolation.instantaneous_zenith_angle.(d, δ, η_UTC, longitude, latitude).:1 +end +radiation = + PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle); + +include( + joinpath( + pkgdir(ClimaLand), + "experiments/integrated/global/global_parameters.jl", + ), +) +soil_args = (domain = domain, parameters = soil_params) +soil_model_type = Soil.EnergyHydrology{FT} + +# Soil microbes model +soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT} + +# soil microbes args +Csom = (z, t) -> eltype(z)(5.0) + +# Set the soil CO2 BC to being atmospheric CO2 +soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC() +soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux +soilco2_sources = (Soil.Biogeochemistry.MicrobeProduction{FT}(),) + +soilco2_boundary_conditions = (; top = soilco2_top_bc, bottom = soilco2_bot_bc) + +soilco2_drivers = Soil.Biogeochemistry.SoilDrivers( + Soil.Biogeochemistry.PrognosticMet{FT}(), + Soil.Biogeochemistry.PrescribedSOC{FT}(Csom), + atmos, +) + +soilco2_args = (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = domain, + parameters = soilco2_ps, + drivers = soilco2_drivers, +) + +# Now we set up the canopy model, which we set up by component: +# Component Types +canopy_component_types = (; + autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = Canopy.TwoStreamModel{FT}, + photosynthesis = Canopy.FarquharModel{FT}, + conductance = Canopy.MedlynConductanceModel{FT}, + hydraulics = Canopy.PlantHydraulicsModel{FT}, + energy = Canopy.BigLeafEnergyModel{FT}, +) +# Individual Component arguments +# Set up autotrophic respiration +autotrophic_respiration_args = + (; parameters = Canopy.AutotrophicRespirationParameters(FT)) +# Set up radiative transfer +radiative_transfer_args = (; + parameters = Canopy.TwoStreamParameters( + FT; + Ω, + α_PAR_leaf, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + ) +) +# Set up conductance +conductance_args = (; parameters = Canopy.MedlynConductanceParameters(FT; g1)) +# Set up photosynthesis +photosynthesis_args = (; + parameters = Canopy.FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25) +) +# Set up plant hydraulics +# Not ideal +LAIfunction = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "lai", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; + preprocess_func = (data) -> data > 0.05 ? data : 0.0, + ), +) +ai_parameterization = Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) + +function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m +end + +plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, +) +plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_midpoints = compartment_midpoints, + compartment_surfaces = compartment_surfaces, +) + +energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) + +# Canopy component args +canopy_component_args = (; + autotrophic_respiration = autotrophic_respiration_args, + radiative_transfer = radiative_transfer_args, + photosynthesis = photosynthesis_args, + conductance = conductance_args, + hydraulics = plant_hydraulics_args, + energy = energy_args, +) + +# Other info needed +shared_params = Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, +) + +canopy_model_args = (; + parameters = shared_params, + domain = ClimaLand.obtain_surface_domain(domain), +) + +# Integrated plant hydraulics and soil model +land_input = (atmos = atmos, radiation = radiation) +land = SoilCanopyModel{FT}(; + soilco2_type = soilco2_type, + soilco2_args = soilco2_args, + land_args = land_input, + soil_model_type = soil_model_type, + soil_args = soil_args, + canopy_component_types = canopy_component_types, + canopy_component_args = canopy_component_args, + canopy_model_args = canopy_model_args, +) + +Y, p, cds = initialize(land) + +t0 = 0.0 +dt = 180.0 +tf = 3600 + +init_soil(ν, θ_r) = θ_r + (ν - θ_r) / 2 +Y.soil.ϑ_l .= init_soil.(ν, θ_r) +Y.soil.θ_i .= FT(0.0) +T = FT(276.85) +ρc_s = + Soil.volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + soil_params.ρc_ds, + soil_params.earth_param_set, + ) +Y.soil.ρe_int .= + Soil.volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T, + soil_params.earth_param_set, + ) +Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air +# On the Caltech cluster, this is not permitted, so we allowscalar. +# This operation is fine on clima cluster +Y.canopy.hydraulics.ϑ_l.:1 .= plant_ν +evaluate!(Y.canopy.energy.T, atmos.T, t0) + +set_initial_cache! = make_set_initial_cache(land) +exp_tendency! = make_exp_tendency(land); +imp_tendency! = ClimaLand.make_imp_tendency(land); +tendency_jacobian! = ClimaLand.make_tendency_jacobian(land); +set_initial_cache!(p, Y, t0) +stepper = CTS.ARS343() +norm_condition = CTS.MaximumError(FT(1e-8)) +conv_checker = CTS.ConvergenceChecker(; norm_condition = norm_condition) +ode_algo = CTS.IMEXAlgorithm( + stepper, + CTS.NewtonsMethod( + max_iters = 1, + update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), + convergence_checker = conv_checker, + ), +) + +# set up jacobian info +jac_kwargs = + (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = tendency_jacobian!) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction( + T_exp! = exp_tendency!, + T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), + dss! = ClimaLand.dss!, + ), + Y, + (t0, tf), + p, +) + +saveat = [t0, tf] +sv = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), +) +saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) +updateat = Array(t0:(3600 * 3):tf) +updatefunc = ClimaLand.make_update_drivers(atmos, radiation) +driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) +cb = SciMLBase.CallbackSet(driver_cb, saving_cb) +@time sol = SciMLBase.solve( + prob, + ode_algo; + dt = dt, + callback = cb, + adaptive = false, + saveat = saveat, +) + +if device_suffix == "cpu" + longpts = range(-180.0, 180.0, 101) + latpts = range(-90.0, 90.0, 101) + hcoords = [ + ClimaCore.Geometry.LatLongPoint(lat, long) for long in longpts, + lat in reverse(latpts) + ] + remapper = ClimaCore.Remapping.Remapper(surface_space, hcoords) + S = ClimaLand.Domains.top_center_to_surface( + (sol.u[end].soil.ϑ_l .- θ_r) ./ (ν .- θ_r), + ) + S_ice = ClimaLand.Domains.top_center_to_surface(sol.u[end].soil.θ_i ./ ν) + T_soil = ClimaLand.Domains.top_center_to_surface(sv.saveval[end].soil.T) + SW = sv.saveval[end].drivers.SW_d + + GPP = sv.saveval[end].canopy.photosynthesis.GPP .* 1e6 + T_canopy = sol.u[end].canopy.energy.T + fields = [S, S_ice, T_soil, GPP, T_canopy, SW] + titles = [ + "Effective saturation", + "Effective ice saturation", + "Temperature (K) - Soil", + "GPP", + "Temperature (K) - Canopy", + "Incident SW", + ] + plotnames = ["S", "Sice", "temp", "gpp", "temp_canopy", "sw"] + mask_top = ClimaLand.Domains.top_center_to_surface(soil_params_mask) + mask_remap = ClimaCore.Remapping.interpolate(remapper, mask_top) + for (id, x) in enumerate(fields) + title = titles[id] + plotname = plotnames[id] + x_remap = ClimaCore.Remapping.interpolate(remapper, x) + + fig = Figure(size = (600, 400)) + ax = Axis( + fig[1, 1], + xlabel = "Longitude", + ylabel = "Latitude", + title = title, + ) + clims = extrema(x_remap) + CairoMakie.heatmap!( + ax, + longpts, + latpts, + oceans_to_value.(x_remap, mask_remap, 0.0), + colorrange = clims, + ) + Colorbar(fig[:, end + 1], colorrange = clims) + outfile = joinpath(outdir, "$plotname.png") + CairoMakie.save(outfile, fig) + end +end diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl index 5a91676078..6604d3f7e8 100644 --- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl +++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl @@ -209,14 +209,14 @@ Y.soil.ϑ_l = drivers.SWC.values[1 + Int(round(t0 / 1800))] # Get soil water con Y.soil.θ_i = FT(0.0) T_0 = FT(drivers.TS.values[1 + Int(round(t0 / 1800))]) # Get soil temperature at t0 ρc_s = - volumetric_heat_capacity.(Y.soil.ϑ_l, Y.soil.θ_i, Ref(land.soil.parameters)) -Y.soil.ρe_int = - volumetric_internal_energy.( + volumetric_heat_capacity.( + Y.soil.ϑ_l, Y.soil.θ_i, - ρc_s, - T_0, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + earth_param_set, ) +Y.soil.ρe_int = + volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl index 638f5b17ea..388a789274 100644 --- a/experiments/integrated/performance/profile_allocations.jl +++ b/experiments/integrated/performance/profile_allocations.jl @@ -37,14 +37,15 @@ function set_initial_conditions(land, t0) volumetric_heat_capacity.( Y.soil.ϑ_l, Y.soil.θ_i, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + land.soil.parameters.earth_param_set, ) Y.soil.ρe_int = volumetric_internal_energy.( Y.soil.θ_i, ρc_s, T_0, - Ref(land.soil.parameters), + land.soil.parameters.earth_param_set, ) Y.soilco2.C = FT(0.000412) # set to atmospheric co2, mol co2 per mol air diff --git a/experiments/standalone/Biogeochemistry/experiment.jl b/experiments/standalone/Biogeochemistry/experiment.jl index 28374bb089..44903485a7 100644 --- a/experiments/standalone/Biogeochemistry/experiment.jl +++ b/experiments/standalone/Biogeochemistry/experiment.jl @@ -132,9 +132,19 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf)) Y.soil.ϑ_l .= FT(0.33) Y.soil.θ_i .= FT(0.1) T = FT(279.85) - ρc_s = Soil.volumetric_heat_capacity(FT(0.33), FT(0.1), params) + ρc_s = Soil.volumetric_heat_capacity( + FT(0.33), + FT(0.1), + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int .= - Soil.volumetric_internal_energy.(FT(0.0), ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.( + FT(0.0), + ρc_s, + T, + params.earth_param_set, + ) end function init_co2!(Y, z) diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl index 1e7a4935e3..4a2770ae6b 100644 --- a/ext/CreateParametersExt.jl +++ b/ext/CreateParametersExt.jl @@ -66,6 +66,7 @@ function LP.LandParameters(toml_dict::CP.AbstractTOMLDict) insol_params, ) end +Base.broadcastable(ps::LP.LandParameters) = tuple(ps) """ AutotrophicRespirationParameters(FT; kwargs...) @@ -259,7 +260,14 @@ end EnergyHydrologyParameters has two constructors: float-type and toml dict based. Additional parameters must be added manually: ν, ν_ss_om, ν_ss_quartz, ν_ss_gravel, hydrology_cm, K_sat, S_s, and θ_r -All parameters can be manually overriden via keyword arguments. +All parameters can be manually overriden via keyword arguments. Note, however, +that certain parameters must have the same type (e.g, if a field is +supplied for porosity, it must be supplied for all other parameters +defined in the interior of the domain). Some parameters are defined only +on the surface of the domain (e.g albedo), while other are defined everywhere +(e.g. porosity). These are indicated with types `F` and `SF`. + +Please see the EnergyHydrologyParameters documentation for a complete list. """ function EnergyHydrologyParameters( ::Type{FT}; @@ -297,10 +305,14 @@ function EnergyHydrologyParameters( K_sat::F, S_s::F, θ_r::F, - PAR_albedo = 0.2, - NIR_albedo = 0.4, + PAR_albedo::SF = 0.2, + NIR_albedo::SF = 0.4, kwargs..., -) where {F <: Union{<:AbstractFloat, ClimaCore.Fields.FieldVector}, C} +) where { + F <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, + SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, + C, +} earth_param_set = LP.LandParameters(toml_dict) # Obtain parameters needed to calculate the derived parameters @@ -360,7 +372,7 @@ function EnergyHydrologyParameters( parameters = CP.get_parameter_values(toml_dict, name_map, "Land") PSE = typeof(earth_param_set) FT = CP.float_type(toml_dict) - EnergyHydrologyParameters{FT, F, C, PSE}(; + EnergyHydrologyParameters{FT, F, SF, C, PSE}(; PAR_albedo, NIR_albedo, ν, diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl index 0c4a36fa11..72ad5333ce 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl @@ -249,14 +249,15 @@ function run_fluxnet( volumetric_heat_capacity.( Y.soil.ϑ_l, Y.soil.θ_i, - Ref(land.soil.parameters), + land.soil.parameters.ρc_ds, + land.soil.parameters.earth_param_set, ) Y.soil.ρe_int = volumetric_internal_energy.( Y.soil.θ_i, ρc_s, T_0, - Ref(land.soil.parameters), + land.soil.parameters.earth_param_set, ) 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] diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl index 67cd03d063..e66a9c91aa 100644 --- a/src/integrated/soil_canopy_model.jl +++ b/src/integrated/soil_canopy_model.jl @@ -92,7 +92,7 @@ function SoilCanopyModel{FT}(; (; atmos, radiation) = land_args # These should always be set by the constructor. sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}()) - # add heat BC + # How to add TOPMODEL?? top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(atmos, CanopyRadiativeFluxes{FT}()) zero_flux = Soil.HeatFluxBC((p, t) -> 0.0) @@ -110,7 +110,7 @@ function SoilCanopyModel{FT}(; ) transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}() - soil_driver = PrognosticSoil{FT}( + soil_driver = PrognosticSoil{typeof(soil.parameters.PAR_albedo)}( soil.parameters.PAR_albedo, soil.parameters.NIR_albedo, ) @@ -276,7 +276,7 @@ function make_update_boundary_fluxes( @. p.root_energy_extraction = p.root_extraction * ClimaLand.Soil.volumetric_internal_energy_liq( p.soil.T, - land.soil.parameters, + land.soil.parameters.earth_param_set, ) # Radiation @@ -352,7 +352,6 @@ function lsm_radiant_energy_fluxes!( R_net_soil = p.soil.R_n LW_out = p.LW_out SW_out = p.SW_out - # in total: INC - OUT = CANOPY_ABS + (1-α_soil)*CANOPY_TRANS # SW out = reflected par + reflected nir @. SW_out = @@ -431,11 +430,12 @@ Concrete type of AbstractSoilDriver used for dispatch in cases where both a canopy model and soil model are run. $(DocStringExtensions.FIELDS) """ -struct PrognosticSoil{FT} <: AbstractSoilDriver +struct PrognosticSoil{F <: Union{AbstractFloat, ClimaCore.Fields.Field}} <: + AbstractSoilDriver "Soil albedo for PAR" - α_PAR::FT + α_PAR::F "Soil albedo for NIR" - α_NIR::FT + α_NIR::F end function Canopy.ground_albedo_PAR(soil_driver::PrognosticSoil, Y, p, t) @@ -482,12 +482,12 @@ end """ root_energy_flux_per_ground_area!( fa_energy::ClimaCore.Fields.Field, - s::PrognosticSoil{FT}, + s::PrognosticSoil{F}, model::Canopy.AbstractCanopyEnergyModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, - ) where {FT} + ) where {FT, F} A method computing the energy flux associated with the root-soil @@ -503,12 +503,12 @@ must account for it as well. """ function Canopy.root_energy_flux_per_ground_area!( fa_energy::ClimaCore.Fields.Field, - s::PrognosticSoil{FT}, + s::PrognosticSoil{F}, model::Canopy.AbstractCanopyEnergyModel{FT}, Y::ClimaCore.Fields.FieldVector, p::NamedTuple, t, -) where {FT} +) where {FT, F} ClimaCore.Operators.column_integral_definite!( fa_energy, p.root_energy_extraction, @@ -555,7 +555,7 @@ end """ Canopy.canopy_radiant_energy_fluxes!(p::NamedTuple, - s::PrognosticSoil{FT}, + s::PrognosticSoil{F}, canopy, radiation::PrescribedRadiativeFluxes, earth_param_set::PSE, @@ -575,13 +575,13 @@ and `p.canopy.radiative_transfer.SW_n`. """ function Canopy.canopy_radiant_energy_fluxes!( p::NamedTuple, - s::PrognosticSoil{FT}, + s::PrognosticSoil{F}, canopy, radiation::PrescribedRadiativeFluxes, earth_param_set::PSE, Y::ClimaCore.Fields.FieldVector, t, -) where {FT, PSE} +) where {F, PSE} nothing end diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl index 6f148f995f..f5cf74d331 100644 --- a/src/standalone/Soil/energy_hydrology.jl +++ b/src/standalone/Soil/energy_hydrology.jl @@ -1,17 +1,30 @@ """ - EnergyHydrologyParameters{FT <: AbstractFloat, - F <: Union{<:AbstractFloat, - ClimaCore.Fields.Field}, - C, - PSE,} + EnergyHydrologyParameters{ + FT <: AbstractFloat, + F <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, + SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, + C, + PSE, + } A parameter structure for the integrated soil water and energy equation system. + +Note that we require two different parameter types F and SF; these +are for parameters that are defined on the surface only +and those defined in the interior of the soil domain: + +- Surface parameters: albedo in each wavelength band (SF) +- Scalar parameters: emissivity, α, β, γ, γT_ref, Ω, + roughness lengths z_0, d_ds ) (FT) +- Parameters defined in the interior: all else (F) + $(DocStringExtensions.FIELDS) """ Base.@kwdef struct EnergyHydrologyParameters{ FT <: AbstractFloat, F <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, + SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field}, C, PSE, } @@ -50,9 +63,9 @@ Base.@kwdef struct EnergyHydrologyParameters{ "Reference temperature for the viscosity factor" γT_ref::FT "Soil PAR Albedo" - PAR_albedo::F + PAR_albedo::SF "Soil NIR Albedo" - NIR_albedo::F + NIR_albedo::SF "Soil Emissivity" emissivity::FT "Roughness length for momentum" @@ -196,13 +209,17 @@ function ClimaLand.make_compute_exp_tendency( top = Operators.SetValue(Geometry.WVector.(heat_top_flux_bc)), bottom = Operators.SetValue(Geometry.WVector.(heat_bottom_flux_bc)), ) - ρe_int_l = volumetric_internal_energy_liq.(p.soil.T, model.parameters) # GradC2F returns a Covariant3Vector, so no need to convert. @. dY.soil.ρe_int = -divf2c_heat( -interpc2f(p.soil.κ) * gradc2f(p.soil.T) - - interpc2f(ρe_int_l * p.soil.K) * gradc2f(p.soil.ψ + z), + interpc2f( + volumetric_internal_energy_liq( + p.soil.T, + model.parameters.earth_param_set, + ) * p.soil.K, + ) * gradc2f(p.soil.ψ + z), ) # Horizontal contributions horizontal_components!( @@ -383,7 +400,10 @@ function horizontal_components!( -hdiv( -p.soil.κ * hgrad(p.soil.T) - p.soil.K * - volumetric_internal_energy_liq(p.soil.T, model.parameters) * + volumetric_internal_energy_liq( + p.soil.T, + model.parameters.earth_param_set, + ) * hgrad(p.soil.ψ + z), ) end @@ -485,9 +505,16 @@ function ClimaLand.make_update_aux(model::EnergyHydrology) θ_r, Ω, γ, + α, + β, + ν_ss_om, + ν_ss_gravel, + ν_ss_quartz, γT_ref, κ_sat_frozen, κ_sat_unfrozen, + ρc_ds, + earth_param_set, ) = model.parameters @. p.soil.θ_l = @@ -498,7 +525,11 @@ function ClimaLand.make_update_aux(model::EnergyHydrology) kersten_number( Y.soil.θ_i, relative_saturation(p.soil.θ_l, Y.soil.θ_i, ν), - model.parameters, + α, + β, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, ), κ_sat(p.soil.θ_l, Y.soil.θ_i, κ_sat_unfrozen, κ_sat_frozen), ) @@ -506,8 +537,13 @@ function ClimaLand.make_update_aux(model::EnergyHydrology) @. p.soil.T = temperature_from_ρe_int( Y.soil.ρe_int, Y.soil.θ_i, - volumetric_heat_capacity(p.soil.θ_l, Y.soil.θ_i, model.parameters), - model.parameters, + volumetric_heat_capacity( + p.soil.θ_l, + Y.soil.θ_i, + ρc_ds, + earth_param_set, + ), + earth_param_set, ) @. p.soil.K = @@ -551,7 +587,7 @@ function ClimaLand.source!( model, ) where {FT} params = model.parameters - (; ν, earth_param_set) = params + (; ν, ρc_ds, θ_r, hydrology_cm, earth_param_set) = params _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set)) _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set)) Δz_top = model.domain.fields.Δz_top # center face distance @@ -561,11 +597,19 @@ function ClimaLand.source!( Y.soil.θ_i, p.soil.T, thermal_time( - volumetric_heat_capacity(p.soil.θ_l, Y.soil.θ_i, params), + volumetric_heat_capacity( + p.soil.θ_l, + Y.soil.θ_i, + ρc_ds, + earth_param_set, + ), 2 * Δz_top, # the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance p.soil.κ, ), - params, + ν, + θ_r, + hydrology_cm, + earth_param_set, ) @. dY.soil.θ_i += (_ρ_l / _ρ_i) * phase_change_source( @@ -573,11 +617,19 @@ function ClimaLand.source!( Y.soil.θ_i, p.soil.T, thermal_time( - volumetric_heat_capacity(p.soil.θ_l, Y.soil.θ_i, params), + volumetric_heat_capacity( + p.soil.θ_l, + Y.soil.θ_i, + ρc_ds, + earth_param_set, + ), 2 * Δz_top, #the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance p.soil.κ, ), - params, + ν, + θ_r, + hydrology_cm, + earth_param_set, ) end @@ -666,7 +718,7 @@ function ClimaLand.surface_resistance( p, t, ) where {FT} - (; ν, θ_r, hydrology_cm) = model.parameters + (; ν, θ_r, d_ds, earth_param_set, hydrology_cm) = model.parameters θ_l_sfc = p.soil.sfc_scratch ClimaLand.Domains.linear_interpolation_to_surface!( θ_l_sfc, @@ -674,8 +726,20 @@ function ClimaLand.surface_resistance( model.domain.fields.z, model.domain.fields.Δz_top, ) + # These are non-allocating θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) - return ClimaLand.Soil.soil_resistance.(θ_l_sfc, θ_i_sfc, model.parameters) + hydrology_cm_sfc = ClimaLand.Domains.top_center_to_surface(hydrology_cm) + ν_sfc = ClimaLand.Domains.top_center_to_surface(ν) + θ_r_sfc = ClimaLand.Domains.top_center_to_surface(θ_r) + return ClimaLand.Soil.soil_resistance.( + θ_l_sfc, + θ_i_sfc, + hydrology_cm_sfc, + ν_sfc, + θ_r_sfc, + d_ds, + earth_param_set, + ) end """ @@ -773,9 +837,11 @@ function ClimaLand.surface_specific_humidity( Thermodynamics.Liquid(), ) .* (@. exp(g * ψ_sfc * M_w / (R * T_sfc))) ice_frac = p.soil.ice_frac - ν_sfc = model.parameters.ν - θ_r_sfc = model.parameters.θ_r - S_c_sfc = model.parameters.hydrology_cm.S_c + surface_space = model.domain.space.surface + ν_sfc = get_top_surface_field(model.parameters.ν, surface_space) + θ_r_sfc = get_top_surface_field(model.parameters.θ_r, surface_space) + S_c_sfc = + get_top_surface_field(model.parameters.hydrology_cm.S_c, surface_space) θ_l_sfc = ClimaLand.Domains.top_center_to_surface(p.soil.θ_l) θ_i_sfc = ClimaLand.Domains.top_center_to_surface(Y.soil.θ_i) @. ice_frac = ice_fraction(θ_l_sfc, θ_i_sfc, ν_sfc, θ_r_sfc) diff --git a/src/standalone/Soil/soil_heat_parameterizations.jl b/src/standalone/Soil/soil_heat_parameterizations.jl index 152289c7d9..dc6fb35427 100644 --- a/src/standalone/Soil/soil_heat_parameterizations.jl +++ b/src/standalone/Soil/soil_heat_parameterizations.jl @@ -29,9 +29,11 @@ end θ_i::FT, T::FT, τ::FT, - params::EnergyHydrologyParameters{FT}, - ) where {FT} - + ν::FT, + θ_r::FT, + hydrology_cm::C, + earth_param_set::EP, + ) where {FT, EP, C} Returns the source term (1/s) used for converting liquid water and ice into each other during phase changes. Note that there are unitless prefactors multiplying this term in the @@ -46,9 +48,11 @@ function phase_change_source( θ_i::FT, T::FT, τ::FT, - params::EnergyHydrologyParameters{FT}, -) where {FT} - (; ν, θ_r, hydrology_cm, earth_param_set) = params + ν::FT, + θ_r::FT, + hydrology_cm::C, + earth_param_set::EP, +) where {FT, EP, C} _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set)) _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set)) _LH_f0 = FT(LP.LH_f0(earth_param_set)) @@ -72,28 +76,30 @@ end volumetric_heat_capacity( θ_l::FT, θ_i::FT, - parameters::EnergyHydrologyParameters{FT}, - ) where {FT} + ρc_ds::FT, + earth_param_set::EP, + ) where {FT,EP} Compute the expression for volumetric heat capacity. """ function volumetric_heat_capacity( θ_l::FT, θ_i::FT, - parameters::EnergyHydrologyParameters{FT}, -) where {FT} - _ρ_i = FT(LP.ρ_cloud_ice(parameters.earth_param_set)) - ρcp_i = FT(LP.cp_i(parameters.earth_param_set) * _ρ_i) + ρc_ds::FT, + earth_param_set::EP, +) where {FT, EP} + _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set)) + ρcp_i = FT(LP.cp_i(earth_param_set) * _ρ_i) - _ρ_l = FT(LP.ρ_cloud_liq(parameters.earth_param_set)) - ρcp_l = FT(LP.cp_l(parameters.earth_param_set) * _ρ_l) + _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set)) + ρcp_l = FT(LP.cp_l(earth_param_set) * _ρ_l) - ρc_s = parameters.ρc_ds + θ_l * ρcp_l + θ_i * ρcp_i + ρc_s = ρc_ds + θ_l * ρcp_l + θ_i * ρcp_i return ρc_s end """ temperature_from_ρe_int(ρe_int::FT, θ_i::FT, ρc_s::FT - parameters::EnergyHydrologyParameters{FT}) where {FT} + earth_param_set::EP) where {FT, EP} A pointwise function for computing the temperature from the volumetric internal energy, volumetric ice content, and volumetric heat capacity of @@ -103,19 +109,19 @@ function temperature_from_ρe_int( ρe_int::FT, θ_i::FT, ρc_s::FT, - parameters::EnergyHydrologyParameters{FT}, -) where {FT} + earth_param_set::EP, +) where {FT, EP} - _ρ_i = FT(LP.ρ_cloud_ice(parameters.earth_param_set)) - _T_ref = FT(LP.T_0(parameters.earth_param_set)) - _LH_f0 = FT(LP.LH_f0(parameters.earth_param_set)) + _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set)) + _T_ref = FT(LP.T_0(earth_param_set)) + _LH_f0 = FT(LP.LH_f0(earth_param_set)) T = _T_ref + (ρe_int + θ_i * _ρ_i * _LH_f0) / ρc_s return T end """ volumetric_internal_energy(θ_i::FT, ρc_s::FT, T::FT, - parameters::EnergyHydrologyParameters{FT}) where {FT} + earth_param_set::EP) where {FT, EP} A pointwise function for computing the volumetric internal energy of the soil, given the volumetric ice content, volumetric heat capacity, and temperature. @@ -124,17 +130,17 @@ function volumetric_internal_energy( θ_i::FT, ρc_s::FT, T::FT, - parameters::EnergyHydrologyParameters{FT}, -) where {FT} - _ρ_i = FT(LP.ρ_cloud_ice(parameters.earth_param_set)) - _LH_f0 = FT(LP.LH_f0(parameters.earth_param_set)) - _T_ref = FT(LP.T_0(parameters.earth_param_set)) + earth_param_set::EP, +) where {FT, EP} + _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set)) + _LH_f0 = FT(LP.LH_f0(earth_param_set)) + _T_ref = FT(LP.T_0(earth_param_set)) ρe_int = ρc_s * (T - _T_ref) - θ_i * _ρ_i * _LH_f0 return ρe_int end """ - volumetric_internal_energy_liq(T::FT, parameters::EnergyHydrologyParameters{FT}) where {FT} + volumetric_internal_energy_liq(T::FT, earth_param_set::EP) where {FT, EP} A pointwise function for computing the volumetric internal energy @@ -142,12 +148,12 @@ of the liquid water in the soil, given the temperature T. """ function volumetric_internal_energy_liq( T::FT, - parameters::EnergyHydrologyParameters{FT}, -) where {FT} + earth_param_set::EP, +) where {FT, EP} - _T_ref = FT(LP.T_0(parameters.earth_param_set)) - _ρ_l = FT(LP.ρ_cloud_liq(parameters.earth_param_set)) - ρcp_l = FT(LP.cp_l(parameters.earth_param_set) * _ρ_l) + _T_ref = FT(LP.T_0(earth_param_set)) + _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set)) + ρcp_l = FT(LP.cp_l(earth_param_set) * _ρ_l) ρe_int_l = ρcp_l * (T - _T_ref) return ρe_int_l end @@ -210,8 +216,12 @@ end kersten_number( θ_i::FT, S_r::FT, - parameters::EnergyHydrologyParameters{FT}, - ) where {FT} + α::FT, + β::FT, + ν_ss_om::FT, + ν_ss_quartz::FT, + ν_ss_gravel::FT, + ) where {FT} Compute the expression for the Kersten number, using the Balland and Arp model. @@ -219,14 +229,12 @@ and Arp model. function kersten_number( θ_i::FT, S_r::FT, - parameters::EnergyHydrologyParameters{FT}, + α::FT, + β::FT, + ν_ss_om::FT, + ν_ss_quartz::FT, + ν_ss_gravel::FT, ) where {FT} - α = parameters.α - β = parameters.β - ν_ss_om = parameters.ν_ss_om - ν_ss_quartz = parameters.ν_ss_quartz - ν_ss_gravel = parameters.ν_ss_gravel - if θ_i < eps(FT) K_e = S_r^((FT(1) + ν_ss_om - α * ν_ss_quartz - ν_ss_gravel) / FT(2)) * diff --git a/src/standalone/Soil/soil_hydrology_parameterizations.jl b/src/standalone/Soil/soil_hydrology_parameterizations.jl index ef113f7890..43ccc68590 100644 --- a/src/standalone/Soil/soil_hydrology_parameterizations.jl +++ b/src/standalone/Soil/soil_hydrology_parameterizations.jl @@ -89,7 +89,8 @@ function pressure_head( ν_eff::FT, S_s::FT, ) where {FT} - S_l_eff = effective_saturation(ν_eff, ϑ_l, θ_r) + ϑ_l_safe = max(ϑ_l, θ_r + eps(FT)) + S_l_eff = effective_saturation(ν_eff, ϑ_l_safe, θ_r) if S_l_eff <= FT(1.0) ψ = matric_potential(cm, S_l_eff) else @@ -308,7 +309,14 @@ function soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} end """ - soil_resistance(θ_l, θ_i, parameters::EnergyHydrologyParameters::FT) + soil_resistance(θ_l::FT, + θ_i::FT, + hydrology_cm::C, + ν::FT, + θ_r::FT, + d_ds::FT, + earth_param_set::EP, + ) where {FT, EP, C} Computes the resistance of the top of the soil column to water vapor diffusion, as a function of the surface @@ -319,9 +327,12 @@ fraction `θ_i`, and other soil parameters. function soil_resistance( θ_l::FT, θ_i::FT, - parameters::EnergyHydrologyParameters{FT}, -) where {FT} - (; ν, hydrology_cm, θ_r, d_ds, earth_param_set) = parameters + hydrology_cm::C, + ν::FT, + θ_r::FT, + d_ds::FT, + earth_param_set::EP, +) where {FT, EP, C} (; S_c) = hydrology_cm _D_vapor = FT(LP.D_vapor(earth_param_set)) S_w = effective_saturation(ν, θ_l + θ_i, θ_r) diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl index e55e04b0d7..26727dc5d8 100644 --- a/src/standalone/Vegetation/PlantHydraulics.jl +++ b/src/standalone/Vegetation/PlantHydraulics.jl @@ -321,7 +321,7 @@ function ClimaLand.Canopy.set_canopy_prescribed_field!( t, ) where {FT} (; LAIfunction, SAI, RAI) = component.parameters.ai_parameterization - evaluate!(p.canopy.hydraulics.area_index.leaf, LAIfunction, t) + evaluate!(p.canopy.hydraulics.area_index.leaf, LAIfunction, floor(t)) @. p.canopy.hydraulics.area_index.stem = SAI @. p.canopy.hydraulics.area_index.root = RAI end diff --git a/src/standalone/Vegetation/radiation.jl b/src/standalone/Vegetation/radiation.jl index 7b0b7b46ef..d39f7da9e1 100644 --- a/src/standalone/Vegetation/radiation.jl +++ b/src/standalone/Vegetation/radiation.jl @@ -236,8 +236,6 @@ function canopy_radiant_energy_fluxes!( @. p.canopy.radiative_transfer.SW_n = (energy_per_photon_PAR * N_a * APAR) + (energy_per_photon_NIR * N_a * ANIR) - LAI = p.canopy.hydraulics.area_index.leaf - SAI = p.canopy.hydraulics.area_index.stem ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI # Long wave: use soil conditions from the PrescribedSoil driver T_soil::FT = s.T(t) diff --git a/test/integrated/soil_energy_hydrology_biogeochemistry.jl b/test/integrated/soil_energy_hydrology_biogeochemistry.jl index 5f0b95342d..9025ff6240 100644 --- a/test/integrated/soil_energy_hydrology_biogeochemistry.jl +++ b/test/integrated/soil_energy_hydrology_biogeochemistry.jl @@ -124,9 +124,22 @@ for FT in (Float32, Float64) Y.soil.ϑ_l .= FT(0.33) Y.soil.θ_i .= FT(0.0) T = FT(279.85) - ρc_s = FT.(Soil.volumetric_heat_capacity(FT(0.33), FT(0.0), params)) + ρc_s = + FT.( + Soil.volumetric_heat_capacity( + FT(0.33), + FT(0.0), + params.ρc_ds, + params.earth_param_set, + ) + ) Y.soil.ρe_int .= - Soil.volumetric_internal_energy.(FT(0.0), ρc_s, T, Ref(params)) + Soil.volumetric_internal_energy.( + FT(0.0), + ρc_s, + T, + params.earth_param_set, + ) end function init_co2!(Y, C_0) @@ -195,13 +208,18 @@ for FT in (Float32, Float64) Y.soil.ϑ_l .= FT(0.33) Y.soil.θ_i .= FT(0.0) T = FT(279.85) - ρc_s = Soil.volumetric_heat_capacity(FT(0.33), FT(0.0), params) + ρc_s = Soil.volumetric_heat_capacity( + FT(0.33), + FT(0.0), + params.ρc_ds, + params.earth_params_set, + ) Y.soil.ρe_int .= Soil.volumetric_internal_energy.( FT(0.0), ρc_s, T, - Ref(params), + params.earth_param_set, ) end diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl index 271828c6c4..124669791f 100644 --- a/test/standalone/Soil/climate_drivers.jl +++ b/test/standalone/Soil/climate_drivers.jl @@ -151,13 +151,18 @@ for FT in (Float32, Float64) Y.soil.ϑ_l .= ν / 2 Y.soil.θ_i .= 0 T = FT(280) - ρc_s = Soil.volumetric_heat_capacity(ν / 2, FT(0), params) + ρc_s = Soil.volumetric_heat_capacity( + ν / 2, + FT(0), + params.ρc_ds, + params.earth_param_set, + ) Y.soil.ρe_int = Soil.volumetric_internal_energy.( FT(0), ρc_s, T, - Ref(params), + params.earth_param_set, ) end diff --git a/test/standalone/Soil/soil_bc.jl b/test/standalone/Soil/soil_bc.jl index f83411309f..549408169a 100644 --- a/test/standalone/Soil/soil_bc.jl +++ b/test/standalone/Soil/soil_bc.jl @@ -174,7 +174,11 @@ for FT in (Float32, Float64) kersten_number.( θ_i, relative_saturation.(ϑ_c, θ_i, ν), - Ref(parameters), + parameters.α, + parameters.β, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, ), κ_sat.( ϑ_c, diff --git a/test/standalone/Soil/soil_parameterizations.jl b/test/standalone/Soil/soil_parameterizations.jl index 713a6b08c1..b9d150c44e 100644 --- a/test/standalone/Soil/soil_parameterizations.jl +++ b/test/standalone/Soil/soil_parameterizations.jl @@ -104,17 +104,21 @@ for FT in (Float32, Float64) FT(5.4e7), FT(0.05), FT(2.1415e6), - parameters, + param_set, ) == FT(_T_ref + (5.4e7 + 0.05 * _ρ_i * _LH_f0) / 2.1415e6) - @test volumetric_heat_capacity(FT(0.25), FT(0.05), parameters) == - FT(ρc_ds + 0.25 * _ρcp_l + 0.05 * _ρcp_i) + @test volumetric_heat_capacity( + FT(0.25), + FT(0.05), + FT(ρc_ds), + param_set, + ) == FT(ρc_ds + 0.25 * _ρcp_l + 0.05 * _ρcp_i) @test volumetric_internal_energy( FT(0.05), FT(2.1415e6), FT(300), - parameters, + param_set, ) == FT(2.1415e6 * (300.0 - _T_ref) - 0.05 * _ρ_i * _LH_f0) @test κ_sat(FT(0.25), FT(0.05), FT(0.57), FT(2.29)) ≈ @@ -126,7 +130,15 @@ for FT in (Float32, Float64) FT((0.25 + 0.05) / 0.4) # ice fraction = 0 - @test kersten_number(FT(0), FT(0.75), parameters) ≈ FT( + @test kersten_number( + FT(0), + FT(0.75), + parameters.α, + parameters.β, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + ) ≈ FT( 0.75^((FT(1) + 0.1 - 0.24 * 0.1 - 0.1) / FT(2)) * ( (FT(1) + exp(-18.3 * 0.75))^(-FT(3)) - @@ -135,13 +147,20 @@ for FT in (Float32, Float64) ) # ice fraction ~= 0 - @test kersten_number(FT(0.05), FT(0.75), parameters) == - FT(0.75^(FT(1) + 0.1)) + @test kersten_number( + FT(0.05), + FT(0.75), + parameters.α, + parameters.β, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + ) == FT(0.75^(FT(1) + 0.1)) @test thermal_conductivity(FT(1.5), FT(0.7287), FT(0.7187)) == FT(0.7287 * 0.7187 + (FT(1) - 0.7287) * 1.5) - @test volumetric_internal_energy_liq(FT(300), parameters) == + @test volumetric_internal_energy_liq(FT(300), param_set) == FT(_ρcp_l * (300.0 - _T_ref)) @test Soil.κ_solid(FT(0.5), FT(0.25), FT(2.0), FT(3.0), FT(2.0)) ≈ @@ -171,12 +190,20 @@ for FT in (Float32, Float64) x = [0.35, 0.25, 0.0, 0.1] y = [0.0, 0.0, 0.1, 0.1 * 0.15 / 0.25] @test dry_soil_layer_thickness.(x, 0.25, 0.1) ≈ y - @test soil_resistance(θ_r, θ_r + eps(FT), parameters) ≈ + @test soil_resistance( + θ_r + eps(FT), + θ_r, + parameters.hydrology_cm, + parameters.ν, + parameters.θ_r, + parameters.d_ds, + param_set, + ) ≈ dry_soil_layer_thickness( Soil.effective_saturation(ν, 2 * θ_r + eps(FT), θ_r), hcm.S_c, parameters.d_ds, - ) / FT(LP.D_vapor(param_set)) / soil_tortuosity(θ_r, θ_r + eps(FT), ν) + ) / FT(LP.D_vapor(param_set)) / soil_tortuosity(θ_r + eps(FT), θ_r, ν) end @testset "Brooks and Corey closure, FT = $FT" begin @@ -317,8 +344,8 @@ for FT in (Float32, Float64) ) Δz = FT(1.0) - @test thermal_time(parameters.ρc_ds, Δz, parameters.κ_dry) == - parameters.ρc_ds * Δz^2 / parameters.κ_dry + τ = thermal_time(parameters.ρc_ds, Δz, parameters.κ_dry) + @test τ == parameters.ρc_ds * Δz^2 / parameters.κ_dry θ_l = FT.([0.11, 0.15, ν]) θ_i = FT(0.0) T = FT(273) @@ -326,12 +353,14 @@ for FT in (Float32, Float64) ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r)) ψT = @.(_LH_f0 / _T_freeze / _grav * (T - _T_freeze)) θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r + ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set) + τ = thermal_time.(ρc_s, Δz, parameters.κ_dry) @test ( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) ≈ - θ_l .- θ_star + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) ≈ + (θ_l .- θ_star) ./ τ ) @test sum( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) .> 0.0, + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) .> 0.0, ) == 3 # try θ_l = 0.1 @@ -342,8 +371,10 @@ for FT in (Float32, Float64) ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r)) ψT = FT(0.0) θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r + ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set) + τ = thermal_time.(ρc_s, Δz, parameters.κ_dry) @test ( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) ≈ + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) ≈ zeros(FT, 3) ) @test (θ_star ≈ θ_l) @@ -356,12 +387,14 @@ for FT in (Float32, Float64) ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r)) ψT = FT(0.0) θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r + ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set) + τ = thermal_time.(ρc_s, Δz, parameters.κ_dry) @test ( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) ≈ - θ_l .- θ_star + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) ≈ + (θ_l .- θ_star) ./ τ ) @test sum( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) .< 0.0, + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) .< 0.0, ) == 2 @@ -372,12 +405,14 @@ for FT in (Float32, Float64) ψ0 = @. matric_potential(hcm, Soil.effective_saturation(ν, θtot, θ_r)) ψT = @.(_LH_f0 / _T_freeze / _grav * (T - _T_freeze)) θ_star = @. inverse_matric_potential(hcm, ψ0 + ψT) * (ν - θ_r) + θ_r + ρc_s = volumetric_heat_capacity.(θ_l, θ_i, parameters.ρc_ds, param_set) + τ = thermal_time.(ρc_s, Δz, parameters.κ_dry) @test ( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) ≈ - θ_l .- θ_star + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) ≈ + (θ_l .- θ_star) ./ τ ) @test sum( - phase_change_source.(θ_l, θ_i, T, FT(1.0), Ref(parameters)) .> 0.0, + phase_change_source.(θ_l, θ_i, T, τ, ν, θ_r, hcm, param_set) .> 0.0, ) == 2 end end diff --git a/test/standalone/Soil/soil_test_3d.jl b/test/standalone/Soil/soil_test_3d.jl index 2ce6af730e..a58a0924a9 100644 --- a/test/standalone/Soil/soil_test_3d.jl +++ b/test/standalone/Soil/soil_test_3d.jl @@ -316,14 +316,19 @@ for FT in (Float32, Float64) θ_l = Soil.volumetric_liquid_fraction.(Y.soil.ϑ_l, ν, θ_r) volumetric_heat_capacity = - Soil.volumetric_heat_capacity.(θ_l, Y.soil.θ_i, Ref(parameters)) + Soil.volumetric_heat_capacity.( + θ_l, + Y.soil.θ_i, + parameters.ρc_ds, + parameters.earth_param_set, + ) Y.soil.ρe_int .= ClimaCore.Fields.zeros(FT, axes(Y.soil.ρe_int)) .+ volumetric_internal_energy.( FT(0), volumetric_heat_capacity, FT(280), - Ref(parameters), + parameters.earth_param_set, ) diff --git a/test/standalone/Soil/soiltest.jl b/test/standalone/Soil/soiltest.jl index 32e146e5d3..300db2c9ac 100644 --- a/test/standalone/Soil/soiltest.jl +++ b/test/standalone/Soil/soiltest.jl @@ -122,6 +122,7 @@ for FT in (Float32, Float64) ### + K_sat = FT(0) hyd_off_en_on = Soil.EnergyHydrologyParameters( FT; ν, @@ -129,7 +130,7 @@ for FT in (Float32, Float64) ν_ss_quartz, ν_ss_gravel, hydrology_cm, - K_sat = FT(0), + K_sat, S_s, θ_r, ) @@ -152,13 +153,14 @@ for FT in (Float32, Float64) ρc_s = @. Soil.volumetric_heat_capacity( Ysoil.soil.ϑ_l, Ysoil.soil.θ_i, - params, + params.ρc_ds, + params.earth_param_set, ) Ysoil.soil.ρe_int = @. Soil.volumetric_internal_energy.( Ysoil.soil.θ_i, ρc_s, T, - params, + params.earth_param_set, ) end @@ -198,6 +200,7 @@ for FT in (Float32, Float64) ### the thermal conductivities to zero, but at the expense of being unreadable. ### Because this is only for a test, we tolerate it :) hyd_on_en_off = Soil.EnergyHydrologyParameters{ + FT, FT, FT, typeof(hydrology_cm), @@ -246,13 +249,14 @@ for FT in (Float32, Float64) ρc_s = @. Soil.volumetric_heat_capacity.( Y.soil.ϑ_l, Ysoil.soil.θ_i, - params, + params.ρc_ds, + params.earth_param_set, ) Ysoil.soil.ρe_int = @. volumetric_internal_energy.( Ysoil.soil.θ_i, ρc_s, FT(288), - params, + params.earth_param_set, ) end @@ -351,7 +355,7 @@ for FT in (Float32, Float64) parent( Soil.volumetric_internal_energy_liq.( p.soil.T, - Ref(soil_water_on.parameters), + Ref(soil_water_on.parameters.earth_param_set), ), ), ) @@ -393,6 +397,7 @@ for FT in (Float32, Float64) ### the thermal conductivities to zero, but at the expense of being unreadable. ### Because this is only for a test, we tolerate it :) hyd_off_en_off = Soil.EnergyHydrologyParameters{ + FT, FT, FT, typeof(hydrology_cm), @@ -443,13 +448,14 @@ for FT in (Float32, Float64) ρc_s = @. Soil.volumetric_heat_capacity( Ysoil.soil.ϑ_l, Ysoil.soil.θ_i, - params, + params.ρc_ds, + params.earth_param_set, ) Ysoil.soil.ρe_int = @. Soil.volumetric_internal_energy.( Ysoil.soil.θ_i, ρc_s, T, - params, + params.earth_param_set, ) end @@ -508,10 +514,15 @@ for FT in (Float32, Float64) ρc_s = @. Soil.volumetric_heat_capacity.( Y.soil.ϑ_l, Ysoil.soil.θ_i, - params, + params.ρc_ds, + params.earth_param_set, + ) + Ysoil.soil.ρe_int = @. volumetric_internal_energy.( + Ysoil.soil.θ_i, + ρc_s, + T, + params.earth_param_set, ) - Ysoil.soil.ρe_int = - @. volumetric_internal_energy.(Ysoil.soil.θ_i, ρc_s, T, params) end t0 = FT(0) @@ -533,7 +544,7 @@ for FT in (Float32, Float64) parent( Soil.volumetric_internal_energy_liq.( p.soil.T, - Ref(soil_both_on.parameters), + Ref(soil_both_on.parameters.earth_param_set), ), ), ) @@ -681,13 +692,14 @@ for FT in (Float32, Float64) ρc_s = @. Soil.volumetric_heat_capacity( Ysoil.soil.ϑ_l, Ysoil.soil.θ_i, - params, + params.ρc_ds, + params.earth_param_set, ) Ysoil.soil.ρe_int = @. Soil.volumetric_internal_energy.( Ysoil.soil.θ_i, ρc_s, T, - params, + params.earth_param_set, ) end