From a4f4989580879f2cdf3b14f4ad401c324429415e Mon Sep 17 00:00:00 2001 From: Gabriele Bozzola Date: Mon, 19 Aug 2024 09:40:46 -0700 Subject: [PATCH 1/2] Switch to actual monthly means ClimaDiagnostics 0.2.4 introduces a new schedule to use actual months when computing reductions. This commit switches to using that when monthly reducitons are requested. --- Project.toml | 2 +- docs/src/diagnostics/users_diagnostics.md | 2 +- .../integrated/global/global_soil_canopy.jl | 3 +- experiments/long_runs/land.jl | 3 +- experiments/long_runs/soil.jl | 7 +- .../Bucket/global_bucket_function.jl | 7 +- src/diagnostics/Diagnostics.jl | 5 +- src/diagnostics/default_diagnostics.jl | 84 ++- .../standard_diagnostic_frequencies.jl | 534 +++++++++--------- 9 files changed, 354 insertions(+), 293 deletions(-) diff --git a/Project.toml b/Project.toml index c6abe0d563..4f6b76c4d8 100644 --- a/Project.toml +++ b/Project.toml @@ -42,7 +42,7 @@ CSV = "0.10" CUDA = "5.3" ClimaComms = "0.5.6, 0.6" ClimaCore = "0.13.2, 0.14" -ClimaDiagnostics = "0.2" +ClimaDiagnostics = "0.2.4" ClimaParams = "0.10.2" ClimaUtilities = "0.1.2" DataFrames = "1" diff --git a/docs/src/diagnostics/users_diagnostics.md b/docs/src/diagnostics/users_diagnostics.md index b474f4994a..3af6c1388f 100644 --- a/docs/src/diagnostics/users_diagnostics.md +++ b/docs/src/diagnostics/users_diagnostics.md @@ -46,7 +46,7 @@ providing the space and output_dir defined in steps 1. and 2. Now that you defined your model and your writter, you can create a callback function to be called when solving your model. For example: ``` -diags = ClimaLand.default_diagnostics(model, 1.0; output_writer = nc_writer) +diags = ClimaLand.default_diagnostics(model, 1.0, reference_date; output_writer = nc_writer) diagnostic_handler = ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl index 7d57c34470..d56e48d9d7 100644 --- a/experiments/integrated/global/global_soil_canopy.jl +++ b/experiments/integrated/global/global_soil_canopy.jl @@ -374,7 +374,8 @@ nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, output_dir) diags = ClimaLand.default_diagnostics( land, - t0; + t0, + ref_time; output_writer = nc_writer, average_period = :hourly, ) diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl index f9a79ce126..3fd0a5dabe 100644 --- a/experiments/long_runs/land.jl +++ b/experiments/long_runs/land.jl @@ -601,7 +601,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15)) diags = ClimaLand.default_diagnostics( land, - t0; + t0, + ref_time; output_writer = nc_writer, output_vars = :long, ) diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl index 36e162739d..e636e22d38 100644 --- a/experiments/long_runs/soil.jl +++ b/experiments/long_runs/soil.jl @@ -420,7 +420,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15)) nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, outdir) - diags = ClimaLand.default_diagnostics(soil, t0; output_writer = nc_writer) + diags = ClimaLand.default_diagnostics( + soil, + t0, + ref_time; + output_writer = nc_writer, + ) diagnostic_handler = ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) diff --git a/experiments/standalone/Bucket/global_bucket_function.jl b/experiments/standalone/Bucket/global_bucket_function.jl index 2572ba23c8..b13656107e 100644 --- a/experiments/standalone/Bucket/global_bucket_function.jl +++ b/experiments/standalone/Bucket/global_bucket_function.jl @@ -168,7 +168,12 @@ space = bucket_domain.space.subsurface nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir) -diags = ClimaLand.default_diagnostics(model, t0; output_writer = nc_writer) +diags = ClimaLand.default_diagnostics( + model, + t0, + ref_time; + output_writer = nc_writer, +) diagnostic_handler = ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index 81678eaee3..e41abb01a1 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -1,5 +1,7 @@ module Diagnostics +import Dates: Month, Period + import ClimaComms using ..Bucket: BucketModel @@ -13,7 +15,8 @@ import ..Domains: top_center_to_surface import ClimaDiagnostics: DiagnosticVariable, ScheduledDiagnostic, average_pre_output_hook! -import ClimaDiagnostics.Schedules: EveryStepSchedule, EveryDtSchedule +import ClimaDiagnostics.Schedules: + EveryStepSchedule, EveryDtSchedule, EveryCalendarDtSchedule import ClimaDiagnostics.Writers: HDF5Writer, NetCDFWriter diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl index f87471b284..e6e4065eb3 100644 --- a/src/diagnostics/default_diagnostics.jl +++ b/src/diagnostics/default_diagnostics.jl @@ -15,6 +15,7 @@ export default_diagnostics reduction, output_writer, t_start, + reference_date, short_names...; pre_output_hook! = nothing, ) @@ -26,25 +27,37 @@ function common_diagnostics( reduction, output_writer, t_start, + reference_date, short_names...; pre_output_hook! = nothing, ) - return [ - ScheduledDiagnostic( - variable = get_diagnostic_variable(short_name), - compute_schedule_func = EveryStepSchedule(), - output_schedule_func = EveryDtSchedule(period; t_start), - reduction_time_func = reduction, - output_writer = output_writer, - pre_output_hook! = pre_output_hook!, - ) for short_name in short_names - ] + return vcat( + map(short_names) do short_name + output_schedule_func = + period isa Period ? + EveryCalendarDtSchedule(period; t_start, reference_date) : + EveryDtSchedule(period; t_start) + return ScheduledDiagnostic( + variable = get_diagnostic_variable(short_name), + compute_schedule_func = EveryStepSchedule(), + output_schedule_func = output_schedule_func, + reduction_time_func = reduction, + output_writer = output_writer, + pre_output_hook! = pre_output_hook!, + ) + end..., + ) end include("standard_diagnostic_frequencies.jl") # Bucket -function default_diagnostics(land_model::BucketModel, t_start; output_writer) +function default_diagnostics( + land_model::BucketModel, + t_start, + reference_date; + output_writer, +) define_diagnostics!(land_model) @@ -64,15 +77,21 @@ function default_diagnostics(land_model::BucketModel, t_start; output_writer) "ssfc", ] - default_outputs = - hourly_averages(bucket_diagnostics...; output_writer, t_start) + default_outputs = hourly_averages( + bucket_diagnostics...; + output_writer, + t_start, + reference_date, + ) + return [default_outputs...] end # SoilCanopyModel function default_diagnostics( land_model::SoilCanopyModel, - t_start; + t_start, + reference_date; output_writer, output_vars = :long, average_period = :daily, @@ -143,14 +162,26 @@ function default_diagnostics( end if average_period == :hourly - default_outputs = - hourly_averages(soilcanopy_diagnostics...; output_writer, t_start) + default_outputs = hourly_averages( + soilcanopy_diagnostics...; + output_writer, + t_start, + reference_date, + ) elseif average_period == :daily - default_outputs = - daily_averages(soilcanopy_diagnostics...; output_writer, t_start) - elseif average_period == :monthly # !! this is currently 30 days, not exact month - default_outputs = - monthly_averages(soilcanopy_diagnostics...; output_writer, t_start) + default_outputs = daily_averages( + soilcanopy_diagnostics...; + output_writer, + t_start, + reference_date, + ) + elseif average_period == :monthly + default_outputs = monthly_averages( + soilcanopy_diagnostics...; + output_writer, + t_start, + reference_date, + ) end return [default_outputs...] @@ -160,7 +191,8 @@ end # SoilModel function default_diagnostics( land_model::EnergyHydrology, - t_start; + t_start, + reference_date; output_writer, ) @@ -168,7 +200,11 @@ function default_diagnostics( soil_diagnostics = ["swc", "si", "sie"] - default_outputs = - daily_averages(soil_diagnostics...; output_writer, t_start) + default_outputs = daily_averages( + soil_diagnostics...; + output_writer, + t_start, + reference_date, + ) return [default_outputs...] end diff --git a/src/diagnostics/standard_diagnostic_frequencies.jl b/src/diagnostics/standard_diagnostic_frequencies.jl index abf80f4d47..b68b1ca02a 100644 --- a/src/diagnostics/standard_diagnostic_frequencies.jl +++ b/src/diagnostics/standard_diagnostic_frequencies.jl @@ -1,262 +1,272 @@ -""" - monthly_maxs(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the monthly max for the given variables. - -A month is defined as 30 days. -""" -monthly_maxs(short_names...; output_writer, t_start) = common_diagnostics( - 30 * 24 * 60 * 60 * one(t_start), - max, - output_writer, - t_start, - short_names..., -) -""" - monthly_max(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the monthly max for the given variable. - -A month is defined as 30 days. -""" -monthly_max(short_names; output_writer, t_start) = - monthly_maxs(short_names; output_writer, t_start)[1] - -""" - monthly_mins(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the monthly min for the given variables. -""" -monthly_mins(short_names...; output_writer, t_start) = common_diagnostics( - 30 * 24 * 60 * 60 * one(t_start), - min, - output_writer, - t_start, - short_names..., -) -""" - monthly_min(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the monthly min for the given variable. - -A month is defined as 30 days. -""" -monthly_min(short_names; output_writer, t_start) = - monthly_mins(short_names; output_writer, t_start)[1] - -""" - monthly_averages(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the monthly average for the given variables. - -A month is defined as 30 days. -""" -# An average is just a sum with a normalization before output -monthly_averages(short_names...; output_writer, t_start) = common_diagnostics( - 30 * 24 * 60 * 60 * one(t_start), - (+), - output_writer, - t_start, - short_names...; - pre_output_hook! = average_pre_output_hook!, -) -""" - monthly_average(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that compute the monthly average for the given variable. - -A month is defined as 30 days. -""" -# An average is just a sum with a normalization before output -monthly_average(short_names; output_writer, t_start) = - monthly_averages(short_names; output_writer, t_start)[1] - -""" - tendaily_maxs(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the max over ten days for the given variables. -""" -tendaily_maxs(short_names...; output_writer, t_start) = common_diagnostics( - 10 * 24 * 60 * 60 * one(t_start), - max, - output_writer, - t_start, - short_names..., -) -""" - tendaily_max(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the max over ten days for the given variable. -""" -tendaily_max(short_names; output_writer, t_start) = - tendaily_maxs(short_names; output_writer, t_start)[1] - -""" - tendaily_mins(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the min over ten days for the given variables. -""" -tendaily_mins(short_names...; output_writer, t_start) = common_diagnostics( - 10 * 24 * 60 * 60 * one(t_start), - min, - output_writer, - t_start, - short_names..., -) -""" - tendaily_min(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the min over ten days for the given variable. -""" -tendaily_min(short_names; output_writer, t_start) = - tendaily_mins(short_names; output_writer, t_start)[1] - -""" - tendaily_averages(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the average over ten days for the given variables. -""" -# An average is just a sum with a normalization before output -tendaily_averages(short_names...; output_writer, t_start) = common_diagnostics( - 10 * 24 * 60 * 60 * one(t_start), - (+), - output_writer, - t_start, - short_names...; - pre_output_hook! = average_pre_output_hook!, -) -""" - tendaily_average(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that compute the average over ten days for the given variable. -""" -# An average is just a sum with a normalization before output -tendaily_average(short_names; output_writer, t_start) = - tendaily_averages(short_names; output_writer, t_start)[1] - -""" - daily_maxs(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the daily max for the given variables. -""" -daily_maxs(short_names...; output_writer, t_start) = common_diagnostics( - 24 * 60 * 60 * one(t_start), - max, - output_writer, - t_start, - short_names..., -) -""" - daily_max(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the daily max for the given variable. -""" -daily_max(short_names; output_writer, t_start) = - daily_maxs(short_names; output_writer, t_start)[1] - -""" - daily_mins(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the daily min for the given variables. -""" -daily_mins(short_names...; output_writer, t_start) = common_diagnostics( - 24 * 60 * 60 * one(t_start), - min, - output_writer, - t_start, - short_names..., -) -""" - daily_min(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the daily min for the given variable. -""" -daily_min(short_names; output_writer, t_start) = - daily_mins(short_names; output_writer, t_start)[1] - -""" - daily_averages(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the daily average for the given variables. -""" -# An average is just a sum with a normalization before output -daily_averages(short_names...; output_writer, t_start) = common_diagnostics( - 24 * 60 * 60 * one(t_start), - (+), - output_writer, - t_start, - short_names...; - pre_output_hook! = average_pre_output_hook!, -) -""" - daily_average(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that compute the daily average for the given variable. -""" -# An average is just a sum with a normalization before output -daily_average(short_names; output_writer, t_start) = - daily_averages(short_names; output_writer, t_start)[1] - -""" - hourly_maxs(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the hourly max for the given variables. -""" -hourly_maxs(short_names...; output_writer, t_start) = common_diagnostics( - 60 * 60 * one(t_start), - max, - output_writer, - t_start, - short_names..., -) - -""" - hourly_max(short_names; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the hourly max for the given variable. -""" -hourly_max(short_names; output_writer, t_start) = - hourly_maxs(short_names; output_writer, t_start)[1] - -""" - hourly_mins(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the hourly min for the given variables. -""" -hourly_mins(short_names...; output_writer, t_start) = common_diagnostics( - 60 * 60 * one(t_start), - min, - output_writer, - t_start, - short_names..., -) -""" - hourly_mins(short_names...; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the hourly min for the given variable. -""" -hourly_min(short_names; output_writer, t_start) = - hourly_mins(short_names; output_writer, t_start)[1] - -# An average is just a sum with a normalization before output -""" - hourly_averages(short_names...; output_writer, t_start) - -Return a list of `ScheduledDiagnostics` that compute the hourly average for the given variables. -""" -hourly_averages(short_names...; output_writer, t_start) = common_diagnostics( - 60 * 60 * one(t_start), - (+), - output_writer, - t_start, - short_names...; - pre_output_hook! = average_pre_output_hook!, -) - -""" - hourly_average(short_names...; output_writer, t_start) - -Return a `ScheduledDiagnostics` that computes the hourly average for the given variable. -""" -hourly_average(short_names; output_writer, t_start) = - hourly_averages(short_names; output_writer, t_start)[1] +""" + monthly_maxs(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the monthly max for the given variables. +""" +monthly_maxs(short_names...; output_writer, t_start, reference_date) = + common_diagnostics(Month(1), max, output_writer, t_start, short_names...) + +""" + monthly_max(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the monthly max for the given variable. +""" +monthly_max(short_names; output_writer, t_start, reference_date) = + monthly_maxs(short_names; output_writer, t_start, reference_date)[1] + +""" + monthly_mins(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the monthly min for the given variables. +""" +monthly_mins(short_names...; output_writer, t_start, reference_date) = + common_diagnostics(Month(1), min, output_writer, t_start, short_names...) + +""" + monthly_min(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the monthly min for the given variable. +""" +monthly_min(short_names; output_writer, t_start, reference_date) = + monthly_mins(short_names; output_writer, t_start, reference_date)[1] + +""" + monthly_averages(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the monthly average for the given variables. +""" +# An average is just a sum with a normalization before output +monthly_averages(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + Month(1), + (+), + output_writer, + t_start, + reference_date, + short_names...; + pre_output_hook! = average_pre_output_hook!, + ) + +""" + monthly_average(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that compute the monthly average for the given variable. +""" +# An average is just a sum with a normalization before output +monthly_average(short_names; output_writer, t_start, reference_date) = + monthly_averages(short_names; output_writer, t_start, reference_date)[1] + +""" + tendaily_maxs(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the max over ten days for the given variables. +""" +tendaily_maxs(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 10 * 24 * 60 * 60 * one(t_start), + max, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + tendaily_max(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the max over ten days for the given variable. +""" +tendaily_max(short_names; output_writer, t_start, reference_date) = + tendaily_maxs(short_names; output_writer, t_start, reference_date)[1] + +""" + tendaily_mins(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the min over ten days for the given variables. +""" +tendaily_mins(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 10 * 24 * 60 * 60 * one(t_start), + min, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + tendaily_min(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the min over ten days for the given variable. +""" +tendaily_min(short_names; output_writer, t_start, reference_date) = + tendaily_mins(short_names; output_writer, t_start, reference_date)[1] + +""" + tendaily_averages(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the average over ten days for the given variables. +""" +# An average is just a sum with a normalization before output +tendaily_averages(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 10 * 24 * 60 * 60 * one(t_start), + (+), + output_writer, + t_start, + reference_date, + short_names...; + pre_output_hook! = average_pre_output_hook!, + ) + +""" + tendaily_average(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that compute the average over ten days for the given variable. +""" +# An average is just a sum with a normalization before output +tendaily_average(short_names; output_writer, t_start, reference_date) = + tendaily_averages(short_names; output_writer, t_start, reference_date)[1] + +""" + daily_maxs(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the daily max for the given variables. +""" +daily_maxs(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 24 * 60 * 60 * one(t_start), + max, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + daily_max(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the daily max for the given variable. +""" +daily_max(short_names; output_writer, t_start, reference_date) = + daily_maxs(short_names; output_writer, t_start, reference_date)[1] + +""" + daily_mins(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the daily min for the given variables. +""" +daily_mins(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 24 * 60 * 60 * one(t_start), + min, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + daily_min(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the daily min for the given variable. +""" +daily_min(short_names; output_writer, t_start, reference_date) = + daily_mins(short_names; output_writer, t_start, reference_date)[1] + +""" + daily_averages(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the daily average for the given variables. +""" +# An average is just a sum with a normalization before output +daily_averages(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 24 * 60 * 60 * one(t_start), + (+), + output_writer, + t_start, + reference_date, + short_names...; + pre_output_hook! = average_pre_output_hook!, + ) + +""" + daily_average(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that compute the daily average for the given variable. +""" +# An average is just a sum with a normalization before output +daily_average(short_names; output_writer, t_start, reference_date) = + daily_averages(short_names; output_writer, t_start, reference_date)[1] + +""" + hourly_maxs(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the hourly max for the given variables. +""" +hourly_maxs(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 60 * 60 * one(t_start), + max, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + hourly_max(short_names; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the hourly max for the given variable. +""" +hourly_max(short_names; output_writer, t_start, reference_date) = + hourly_maxs(short_names; output_writer, t_start, reference_date)[1] + +""" + hourly_mins(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the hourly min for the given variables. +""" +hourly_mins(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 60 * 60 * one(t_start), + min, + output_writer, + t_start, + reference_date, + short_names..., + ) + +""" + hourly_mins(short_names...; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the hourly min for the given variable. +""" +hourly_min(short_names; output_writer, t_start, reference_date) = + hourly_mins(short_names; output_writer, t_start, reference_date)[1] + +# An average is just a sum with a normalization before output +""" + hourly_averages(short_names...; output_writer, t_start, reference_date) + +Return a list of `ScheduledDiagnostics` that compute the hourly average for the given variables. +""" +hourly_averages(short_names...; output_writer, t_start, reference_date) = + common_diagnostics( + 60 * 60 * one(t_start), + (+), + output_writer, + t_start, + reference_date, + short_names...; + pre_output_hook! = average_pre_output_hook!, + ) + +""" + hourly_average(short_names...; output_writer, t_start, reference_date) + +Return a `ScheduledDiagnostics` that computes the hourly average for the given variable. +""" +hourly_average(short_names; output_writer, t_start, reference_date) = + hourly_averages(short_names; output_writer, t_start, reference_date)[1] From 7ae418131a20e7b12d25bc6d3ee4272d5a7a803a Mon Sep 17 00:00:00 2001 From: Gabriele Bozzola Date: Mon, 19 Aug 2024 10:46:31 -0700 Subject: [PATCH 2/2] Remove stray ^M characters Some ^M characters infiltrated ClimaLand (as a result of editing files on Windows). This commit removes all of them. --- .github/workflows/ClimaLandSimulations.yml | 64 +- docs/list_diagnostics.jl | 10 +- docs/list_standalone_models.jl | 66 +- .../src/diagnostics/developers_diagnostics.md | 234 +++---- docs/src/diagnostics/users_diagnostics.md | 264 ++++--- docs/src/standalone/apps.jl | 94 +-- docs/src/standalone/apps/beer.jl | 70 +- docs/src/standalone/apps/hetero_resp.jl | 184 ++--- docs/src/standalone/apps/leaf_an.jl | 274 ++++---- docs/src/standalone/apps/leaf_an_ci.jl | 254 +++---- lib/ClimaLandSimulations/README.md | 2 +- .../src/Dashboards/Dashboards.jl | 24 +- .../src/Dashboards/fluxnet_dashboard.jl | 102 +-- .../src/Fluxnet/Fluxnet.jl | 118 ++-- .../fluxnet_sites/US-Ha1/US-Ha1_parameters.jl | 416 +++++------ .../fluxnet_sites/US-NR1/US-NR1_parameters.jl | 416 +++++------ .../fluxnet_sites/US-Var/US-Var_parameters.jl | 416 +++++------ .../src/Fluxnet/run_fluxnet.jl | 655 +++++++++--------- 18 files changed, 1821 insertions(+), 1842 deletions(-) diff --git a/.github/workflows/ClimaLandSimulations.yml b/.github/workflows/ClimaLandSimulations.yml index 5ae09e5b46..27d463b255 100644 --- a/.github/workflows/ClimaLandSimulations.yml +++ b/.github/workflows/ClimaLandSimulations.yml @@ -1,32 +1,32 @@ -name: ClimaLandSimulations CI -on: - pull_request: - push: - branches: - - main - -permissions: - actions: write - contents: read - -jobs: - lib-climalandsimulations: - runs-on: ubuntu-20.04 - timeout-minutes: 45 - steps: - - name: Checkout - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@latest - with: - version: '1.10' -# julia-actions/cache@v2 did not seem to invalidate cache correctly... -# - uses: julia-actions/cache@v2 - - name: Install Julia dependencies - run: > - julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.develop(path=".")' - - name: Run the tests - continue-on-error: true - env: - CI_OUTPUT_DIR: output - run: > - julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.test("ClimaLandSimulations")' +name: ClimaLandSimulations CI +on: + pull_request: + push: + branches: + - main + +permissions: + actions: write + contents: read + +jobs: + lib-climalandsimulations: + runs-on: ubuntu-20.04 + timeout-minutes: 45 + steps: + - name: Checkout + uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@latest + with: + version: '1.10' +# julia-actions/cache@v2 did not seem to invalidate cache correctly... +# - uses: julia-actions/cache@v2 + - name: Install Julia dependencies + run: > + julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.develop(path=".")' + - name: Run the tests + continue-on-error: true + env: + CI_OUTPUT_DIR: output + run: > + julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.test("ClimaLandSimulations")' diff --git a/docs/list_diagnostics.jl b/docs/list_diagnostics.jl index 4d76cb4f6f..b299d62a78 100644 --- a/docs/list_diagnostics.jl +++ b/docs/list_diagnostics.jl @@ -1,5 +1,5 @@ -diagnostics = [ - "Available diagnostics" => "diagnostics/available_diagnostics.md", - "For developers" => "diagnostics/developers_diagnostics.md", - "For users" => "diagnostics/users_diagnostics.md", -] +diagnostics = [ + "Available diagnostics" => "diagnostics/available_diagnostics.md", + "For developers" => "diagnostics/developers_diagnostics.md", + "For users" => "diagnostics/users_diagnostics.md", +] diff --git a/docs/list_standalone_models.jl b/docs/list_standalone_models.jl index 698a6b4fa9..d3f61407db 100644 --- a/docs/list_standalone_models.jl +++ b/docs/list_standalone_models.jl @@ -1,33 +1,33 @@ -standalone_models = [ - "Vegetation" => [ - "Radiative transfer" => [ - "Beer model" => "standalone/pages/vegetation/radiative_transfer/beer_model.md", - "Two-Stream model" => "standalone/pages/vegetation/radiative_transfer/twostream_model.md", - ] - "Photosynthesis" => [ - "Farquhar model" => "standalone/pages/vegetation/photosynthesis/farquhar_model.md", - "Optimality model" => "standalone/pages/vegetation/photosynthesis/optimality_model.md", - ] - "Stomatal conductance" => [ - "Medlyn model" => "standalone/pages/vegetation/stomatal_conductance/medlyn_model.md", - ] - "Plant hydraulics" => [ - "Van Genuchten model" => "standalone/pages/vegetation/plant_hydraulics/van_genuchten_model.md", - ] - ] - "Soil" => [ - "Biogeochemistry" => [ - "DAMM model" => "standalone/pages/soil/biogeochemistry/DAMM_model.md", - ] - "Hydrology" => [ - "Richards model" => "standalone/pages/soil/hydrology/richards_model.md", - ] - "Energy" => [ - "Energy model" => "standalone/pages/soil/energy/energy_model.md", - ] - ] - "Snow" => ["Snow model" => "standalone/pages/snow/snow_model.md"] - "SurfaceWater" => [ - "Surface water model" => "standalone/pages/surface_water/surface_water_model.md", - ] -] +standalone_models = [ + "Vegetation" => [ + "Radiative transfer" => [ + "Beer model" => "standalone/pages/vegetation/radiative_transfer/beer_model.md", + "Two-Stream model" => "standalone/pages/vegetation/radiative_transfer/twostream_model.md", + ] + "Photosynthesis" => [ + "Farquhar model" => "standalone/pages/vegetation/photosynthesis/farquhar_model.md", + "Optimality model" => "standalone/pages/vegetation/photosynthesis/optimality_model.md", + ] + "Stomatal conductance" => [ + "Medlyn model" => "standalone/pages/vegetation/stomatal_conductance/medlyn_model.md", + ] + "Plant hydraulics" => [ + "Van Genuchten model" => "standalone/pages/vegetation/plant_hydraulics/van_genuchten_model.md", + ] + ] + "Soil" => [ + "Biogeochemistry" => [ + "DAMM model" => "standalone/pages/soil/biogeochemistry/DAMM_model.md", + ] + "Hydrology" => [ + "Richards model" => "standalone/pages/soil/hydrology/richards_model.md", + ] + "Energy" => [ + "Energy model" => "standalone/pages/soil/energy/energy_model.md", + ] + ] + "Snow" => ["Snow model" => "standalone/pages/snow/snow_model.md"] + "SurfaceWater" => [ + "Surface water model" => "standalone/pages/surface_water/surface_water_model.md", + ] +] diff --git a/docs/src/diagnostics/developers_diagnostics.md b/docs/src/diagnostics/developers_diagnostics.md index 429973d632..2e1620c7d0 100644 --- a/docs/src/diagnostics/developers_diagnostics.md +++ b/docs/src/diagnostics/developers_diagnostics.md @@ -1,117 +1,117 @@ -# ClimaLand Diagnostics: why and how - -ClimaLand simulations generates variables in the integrator state and cache at each time step. -A user will need to use these variables in some form, i.e., access them from a file that contains variables at a given temporal and spatial resolution. -The user will also want to retrieve metadata about those variables, such as name and units. -This is where ClimaLand diagnostics comes in, it writes simulations variables (in a file, such as NetCDF or HDF5, or in Julia Dict), at a specified spatio-temporal reduction -(e.g., hourly averages, monthly max, instantaneous, integrated through soil depth...), along with metadata (e.g., soil temperature short name is t_soil, expressed in "K" units). -We want to provide users with default options, but also the possibility to define their own variables and reductions. - -Internally, this is done by using the [`ClimaDiagnostics.jl`](https://github.com/CliMA/ClimaDiagnostics.jl) package, that provides the functionality to produce a -[`ClimaLand.Diagnostics`](https://github.com/CliMA/ClimaLand.jl/tree/main/src/Diagnostics/Diagnostics.jl) module in the src/Diagnostics.jl folder. In this folder, - - `Diagnostics.jl` defines the module, - - `diagnostic.jl` defines `ALL_DIAGNOSTICS`, a Dict containing all diagnostics variables defined in `define_diagnostics.jl`, it also defines the function - `add_diagnostic_variable!` which defines a method to add diagnostic variables to ALL_DIAGNOSTICS, finally it contains a function `get_diagnostic_variable` which returns a - `DiagnosticVariable` from its `short_name`, if it exists. - - `define_diagnostics.jl`, mentioned above, contains a function `define_diagnostics!(land_model)` which contains all default diagnostic variables by calling. - `add_diagnostic_variable!`, and dispatch off the type of land_model to define how to compute a diagnostic (for example, surface temperature is computed in `p.bucket.T_sfc` in the bucket model). - - compute methods are defined in a separate file, for example, `bucket_compute_methods.jl`. - - `standard_diagnostic_frequencies.jl` defines standard functions to schedule diagnostics, for example, hourly average or monthly max, these functions are called on a list of diagnostic variables. As developers, we can add more standard functions that users may want to have access to easily in this file. - - `default_diagnostics.jl` defines default diagnostics functions to use on a model simulation. For example, `default_diagnostics(land_model::BucketModel, t_start; output_writer)`. - will return a `ScheduledDiagnostics` that computes hourly averages for all Bucket variables, along with their metadata, ready to be written on a NetCDF file when running a Bucket simulation. - -The following section give more details on these functions, along with examples. As developers, we want to extand these functionality as ClimaLand progresses. - -# Compute methods - -Each model defines all its compute methods in a file (bucket_compute_methods.jl for the bucket model, for example). -The structure of a diagnostic variable compute method is, for example: -``` -function compute_albedo!(out, Y, p, t, land_model::BucketModel) - if isnothing(out) - return copy(p.bucket.α_sfc) - else - out .= p.bucket.α_sfc - end -end -``` - -It defines how to access your diagnostic (here, p.bucket.α_sfc), in your model type (here, ::BucketModel). -Note that, as explained in the [ClimaDiagnostics.jl documentation](https://clima.github.io/ClimaDiagnostics.jl/dev/user_guide/), `out` will probably not be needed in the future. - -We also define helper functions returning error messages if a user tries to compute a diagnostic variable that doesn't exist in their model type. - -``` -error_diagnostic_variable(variable, land_model::T) where {T} = - error("Cannot compute $variable with model = $T") - -compute_albedo!(_, _, _, _, land_model) = - error_diagnostic_variable("albedo", land_model) -``` - -# Define diagnostics - -Once the compute functions have been defined, they are added to `define_diagnostics!(land_model)`, which adds diagnostics variables to ALL_DIAGNOSTICS dict, -defined in diagnostic.jl. In these functions, you also define a `short_name`, `long_name`, `standard_name`, `units` and `comment`. For example: - -``` -add_diagnostic_variable!( - short_name = "alpha", - long_name = "Albedo", - standard_name = "albedo", - units = "", - compute! = (out, Y, p, t) -> compute_albedo!(out, Y, p, t, land_model), - ) -``` - -# Default diagnostics - -For each model, we define a function `default_diagnostics` which will define what diagnostic variables to compute by default for a specific model, and -on what schedule (for example, hourly average). For example, -``` -function default_diagnostics(land_model::BucketModel, t_start; output_writer) - - define_diagnostics!(land_model) - - bucket_diagnostics = [ - "alpha", - "rn", - "tsfc", - "qsfc", - "lhf", - "rae", - "shf", - "vflux", - "rhosfc", - "t", - "w", - "ws", - "sigmas", - ] - - default_outputs = - hourly_averages(bucket_diagnostics...; output_writer, t_start) - return [default_outputs...] -end -``` - -is the default for the BucketModel, it will return hourly averages for the variables listed in `bucket_diagnostics` (which are all variables in the BucketModel). - -# Standard diagnostic frequencies - -We defined some functions of diagnostic schedule that may often be used in `standard_diagnostic_frequencies.jl`, for example - -``` -hourly_averages(short_names...; output_writer, t_start) = common_diagnostics( - 60 * 60 * one(t_start), - (+), - output_writer, - t_start, - short_names...; - pre_output_hook! = average_pre_output_hook!, -) -``` - -will return a list of `ScheduledDiagnostics` that compute the hourly average for the given variables listed in `short_names`. -We also, so far, provide functions for mins, maxs and averages aggregated monthly, over ten days, daily, and hourly. -As a developer, you may want to add more standard diagnostics here. +# ClimaLand Diagnostics: why and how + +ClimaLand simulations generates variables in the integrator state and cache at each time step. +A user will need to use these variables in some form, i.e., access them from a file that contains variables at a given temporal and spatial resolution. +The user will also want to retrieve metadata about those variables, such as name and units. +This is where ClimaLand diagnostics comes in, it writes simulations variables (in a file, such as NetCDF or HDF5, or in Julia Dict), at a specified spatio-temporal reduction +(e.g., hourly averages, monthly max, instantaneous, integrated through soil depth...), along with metadata (e.g., soil temperature short name is t_soil, expressed in "K" units). +We want to provide users with default options, but also the possibility to define their own variables and reductions. + +Internally, this is done by using the [`ClimaDiagnostics.jl`](https://github.com/CliMA/ClimaDiagnostics.jl) package, that provides the functionality to produce a +[`ClimaLand.Diagnostics`](https://github.com/CliMA/ClimaLand.jl/tree/main/src/Diagnostics/Diagnostics.jl) module in the src/Diagnostics.jl folder. In this folder, + - `Diagnostics.jl` defines the module, + - `diagnostic.jl` defines `ALL_DIAGNOSTICS`, a Dict containing all diagnostics variables defined in `define_diagnostics.jl`, it also defines the function + `add_diagnostic_variable!` which defines a method to add diagnostic variables to ALL_DIAGNOSTICS, finally it contains a function `get_diagnostic_variable` which returns a + `DiagnosticVariable` from its `short_name`, if it exists. + - `define_diagnostics.jl`, mentioned above, contains a function `define_diagnostics!(land_model)` which contains all default diagnostic variables by calling. + `add_diagnostic_variable!`, and dispatch off the type of land_model to define how to compute a diagnostic (for example, surface temperature is computed in `p.bucket.T_sfc` in the bucket model). + - compute methods are defined in a separate file, for example, `bucket_compute_methods.jl`. + - `standard_diagnostic_frequencies.jl` defines standard functions to schedule diagnostics, for example, hourly average or monthly max, these functions are called on a list of diagnostic variables. As developers, we can add more standard functions that users may want to have access to easily in this file. + - `default_diagnostics.jl` defines default diagnostics functions to use on a model simulation. For example, `default_diagnostics(land_model::BucketModel, t_start; output_writer)`. + will return a `ScheduledDiagnostics` that computes hourly averages for all Bucket variables, along with their metadata, ready to be written on a NetCDF file when running a Bucket simulation. + +The following section give more details on these functions, along with examples. As developers, we want to extand these functionality as ClimaLand progresses. + +# Compute methods + +Each model defines all its compute methods in a file (bucket_compute_methods.jl for the bucket model, for example). +The structure of a diagnostic variable compute method is, for example: +``` +function compute_albedo!(out, Y, p, t, land_model::BucketModel) + if isnothing(out) + return copy(p.bucket.α_sfc) + else + out .= p.bucket.α_sfc + end +end +``` + +It defines how to access your diagnostic (here, p.bucket.α_sfc), in your model type (here, ::BucketModel). +Note that, as explained in the [ClimaDiagnostics.jl documentation](https://clima.github.io/ClimaDiagnostics.jl/dev/user_guide/), `out` will probably not be needed in the future. + +We also define helper functions returning error messages if a user tries to compute a diagnostic variable that doesn't exist in their model type. + +``` +error_diagnostic_variable(variable, land_model::T) where {T} = + error("Cannot compute $variable with model = $T") + +compute_albedo!(_, _, _, _, land_model) = + error_diagnostic_variable("albedo", land_model) +``` + +# Define diagnostics + +Once the compute functions have been defined, they are added to `define_diagnostics!(land_model)`, which adds diagnostics variables to ALL_DIAGNOSTICS dict, +defined in diagnostic.jl. In these functions, you also define a `short_name`, `long_name`, `standard_name`, `units` and `comment`. For example: + +``` +add_diagnostic_variable!( + short_name = "alpha", + long_name = "Albedo", + standard_name = "albedo", + units = "", + compute! = (out, Y, p, t) -> compute_albedo!(out, Y, p, t, land_model), + ) +``` + +# Default diagnostics + +For each model, we define a function `default_diagnostics` which will define what diagnostic variables to compute by default for a specific model, and +on what schedule (for example, hourly average). For example, +``` +function default_diagnostics(land_model::BucketModel, t_start; output_writer) + + define_diagnostics!(land_model) + + bucket_diagnostics = [ + "alpha", + "rn", + "tsfc", + "qsfc", + "lhf", + "rae", + "shf", + "vflux", + "rhosfc", + "t", + "w", + "ws", + "sigmas", + ] + + default_outputs = + hourly_averages(bucket_diagnostics...; output_writer, t_start) + return [default_outputs...] +end +``` + +is the default for the BucketModel, it will return hourly averages for the variables listed in `bucket_diagnostics` (which are all variables in the BucketModel). + +# Standard diagnostic frequencies + +We defined some functions of diagnostic schedule that may often be used in `standard_diagnostic_frequencies.jl`, for example + +``` +hourly_averages(short_names...; output_writer, t_start) = common_diagnostics( + 60 * 60 * one(t_start), + (+), + output_writer, + t_start, + short_names...; + pre_output_hook! = average_pre_output_hook!, +) +``` + +will return a list of `ScheduledDiagnostics` that compute the hourly average for the given variables listed in `short_names`. +We also, so far, provide functions for mins, maxs and averages aggregated monthly, over ten days, daily, and hourly. +As a developer, you may want to add more standard diagnostics here. diff --git a/docs/src/diagnostics/users_diagnostics.md b/docs/src/diagnostics/users_diagnostics.md index 3af6c1388f..450b79c94e 100644 --- a/docs/src/diagnostics/users_diagnostics.md +++ b/docs/src/diagnostics/users_diagnostics.md @@ -1,141 +1,123 @@ -# Using ClimaLand Diagnostics when running a simulation - -When running a ClimaLand simulations, you have multiple options on how to write the outputs of that simulation. -You may want all variables, or just a selected few. -You may want instantaneous values, at the highest temporal and spatial resolution, or you may want to get averages at hourly or monthly time scale, and integrate in space -(for example soil moisture from 0 to 1 meter depth). -You may want to get more specific reductions, such as 10 days maximums, or compute a new variables that is a function of others. -You may want to get your outputs in memory in a Julia Dict, or write them in a NetCDF file. - -This is where ClimaLand Diagnostics comes in for users. - -In this documentation page, we first explain how to use default diagnostics and what are the defaults, and then explain how to define your own diagnostics for more advanced users. - -# Default Diagnostics - -Once you have defined your model and are ready to run a simulation, and after adding ClimaDiagnostics (using ClimaDiagnostics), - you can add default diagnostics to it by doing the following steps: - -1. define an output folder - -``` -output_dir = ClimaUtilities.OutputPathGenerator.generate_output_path("base_output_dir/") -``` - -2. define a space - -Your diagnostics will be written in time and space. These may be defined in your model, but usually land model space is a sphere with no vertical dimension. -You may have variables varying with soil depth, and so you will need: - -``` -space = bucket_domain.space.subsurface -``` - -3. define your writter - -Your diagnostics will be written in a Julia Dict or a netcdf file, for example. This is up to you. For a netcdf file, you define your writter like this: - -``` -nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir) -``` - -providing the space and output_dir defined in steps 1. and 2. - -4. make your diagnostics on your model, using your writter, and define a callback - -Now that you defined your model and your writter, you can create a callback function to be called when solving your model. For example: - -``` -diags = ClimaLand.default_diagnostics(model, 1.0, reference_date; output_writer = nc_writer) - -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) - -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) - -sol = SciMLBase.solve(prob, ode_algo; dt = Δt, callback = diag_cb) -``` - -Your diagnostics have now been written in netcdf files in your output folder. - -Note that by default, `default_diagnostics` assign two optional kwargs: `output_vars = :long` and `average_period` = :daily. -`output_vars = :long` will write all available diagnostics, whereas `output_vars = :short` will only write essentials diagnostics. -`average_period` defines the period over which diagnostics are averaged, it can be set to `:hourly`, `:daily` and `:monthly`. - -# Custom Diagnostics - -When defining a custom diagnostic, follow these steps: - - Define how to compute your diagnostic variable from your model state and cache. -For example, let's say you want the Bowen ratio (ratio between sensible heat and latent heat) in the Bucket model. -``` -function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel) - if isnothing(out) - return copy(p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf) - else - out .= p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf - end -end -``` -It is good practice to add error messages to inform other users that a diagnostic variable makes sense only with a -specific `land_model`. This can be accomplished by prepending the `@with_error` macro at the function declaration, -as in -``` -import ClimaLand.Diagnostics: @witherror - -@with_error function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel) - if isnothing(out) - return copy(p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf) - else - out .= p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf - end -end -``` -So, when someone tries outputting the Bowen ratio with a different model (e.g., `SnowModel`), `ClimaLand` will produce the following message: -``` -Cannot compute albedo with model = SnowModel -``` - - Add that diagnostic variable to your list of variables - ``` - add_diagnostic_variable!( - short_name = "bor", - long_name = "Bowen ratio", - standard_name = "bowen_ratio", - units = "", - compute! = (out, Y, p, t) -> compute_bowen_ratio!(out, Y, p, t, land_model), - ) - ``` - - Define how to schedule your variables. For example, you want the seasonal maximum of your variables, where season is defined as 90 days. - ``` - seasonal_maxs(short_names...; output_writer, t_start) = common_diagnostics( - 90 * 24 * 60 * 60 * one(t_start), - max, - output_writer, - t_start, - short_names..., -) - ``` - - Define a function to return your `ScheduledDiagnostics` - ``` - function default_diagnostics(land_model::BucketModel, t_start; output_writer) - - define_diagnostics!(land_model) - - add_diagnostic_variable!( - short_name = "bor", - long_name = "Bowen ratio", - standard_name = "bowen_ratio", - units = "", - compute! = (out, Y, p, t) -> compute_bowen_ratio!(out, Y, p, t, land_model), - ) - - my_custom_diagnostics = [ - "lhf", - "shf", - "bor", - ] - - my_custom_outputs = - seasonal_maxs(my_custom_diagnostics...; output_writer, t_start) - return [my_custom_outputs...] -end - ``` +# Using ClimaLand Diagnostics when running a simulation + +When running a ClimaLand simulations, you have multiple options on how to write the outputs of that simulation. +You may want all variables, or just a selected few. +You may want instantaneous values, at the highest temporal and spatial resolution, or you may want to get averages at hourly or monthly time scale, and integrate in space +(for example soil moisture from 0 to 1 meter depth). +You may want to get more specific reductions, such as 10 days maximums, or compute a new variables that is a function of others. +You may want to get your outputs in memory in a Julia Dict, or write them in a NetCDF file. + +This is where ClimaLand Diagnostics comes in for users. + +In this documentation page, we first explain how to use default diagnostics and what are the defaults, and then explain how to define your own diagnostics for more advanced users. + +# Default Diagnostics + +Once you have defined your model and are ready to run a simulation, and after adding ClimaDiagnostics (using ClimaDiagnostics), + you can add default diagnostics to it by doing the following steps: + +1. define an output folder + +``` +output_dir = ClimaUtilities.OutputPathGenerator.generate_output_path("base_output_dir/") +``` + +2. define a space + +Your diagnostics will be written in time and space. These may be defined in your model, but usually land model space is a sphere with no vertical dimension. +You may have variables varying with soil depth, and so you will need: + +``` +space = bucket_domain.space.subsurface +``` + +3. define your writter + +Your diagnostics will be written in a Julia Dict or a netcdf file, for example. This is up to you. For a netcdf file, you define your writter like this: + +``` +nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir) +``` + +providing the space and output_dir defined in steps 1. and 2. + +4. make your diagnostics on your model, using your writter, and define a callback + +Now that you defined your model and your writter, you can create a callback function to be called when solving your model. For example: + +``` +diags = ClimaLand.default_diagnostics(model, 1.0, reference_date; output_writer = nc_writer) + +diagnostic_handler = + ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) + +diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) + +sol = SciMLBase.solve(prob, ode_algo; dt = Δt, callback = diag_cb) +``` + +Your diagnostics have now been written in netcdf files in your output folder. + +Note that by default, `default_diagnostics` assign two optional kwargs: `output_vars = :long` and `average_period` = :daily. +`output_vars = :long` will write all available diagnostics, whereas `output_vars = :short` will only write essentials diagnostics. +`average_period` defines the period over which diagnostics are averaged, it can be set to `:hourly`, `:daily` and `:monthly`. + +# Custom Diagnostics + +When defining a custom diagnostic, follow these steps: + - Define how to compute your diagnostic variable from your model state and cache. + For example, let's say you want the bowen ratio (ratio between sensible heat and latent heat) in the Bucket model. + ``` + function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel) + if isnothing(out) + return copy(p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf) + else + out .= p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf + end +end + ``` + - Add that diagnostic variable to your list of variables + ``` + add_diagnostic_variable!( + short_name = "bor", + long_name = "Bowen ratio", + standard_name = "bowen_ratio", + units = "", + compute! = (out, Y, p, t) -> compute_bowen_ratio!(out, Y, p, t, land_model), + ) + ``` + - Define how to schedule your variables. For example, you want the seasonal maximum of your variables, where season is defined as 90 days. + ``` + seasonal_maxs(short_names...; output_writer, t_start) = common_diagnostics( + 90 * 24 * 60 * 60 * one(t_start), + max, + output_writer, + t_start, + short_names..., +) + ``` + - Define a function to return your `ScheduledDiagnostics` + ``` + function default_diagnostics(land_model::BucketModel, t_start; output_writer) + + define_diagnostics!(land_model) + + add_diagnostic_variable!( + short_name = "bor", + long_name = "Bowen ratio", + standard_name = "bowen_ratio", + units = "", + compute! = (out, Y, p, t) -> compute_bowen_ratio!(out, Y, p, t, land_model), + ) + + my_custom_diagnostics = [ + "lhf", + "shf", + "bor", + ] + + my_custom_outputs = + seasonal_maxs(my_custom_diagnostics...; output_writer, t_start) + return [my_custom_outputs...] +end + ``` diff --git a/docs/src/standalone/apps.jl b/docs/src/standalone/apps.jl index 0820734eaf..cdb1002497 100644 --- a/docs/src/standalone/apps.jl +++ b/docs/src/standalone/apps.jl @@ -1,47 +1,47 @@ -# Load packages -using ParamViz -using Bonito -using ClimaLand -using ClimaLand.Canopy -using ClimaLand.Soil.Biogeochemistry -import ClimaParams as CP -import ClimaLand.Parameters as LP -using ClimaLandSimulations -FT = Float64 -earth_param_set = LP.LandParameters(FT) -RTparams = BeerLambertParameters(FT) - -# Create server -IPa = "127.0.0.1" -port = 9385 -server = Server( - IPa, - port; - proxy_url = "https://clima.westus3.cloudapp.azure.com/jsserve/", -) - -# Load and create apps -include("apps/leaf_an.jl"); -An_app = An_app_f(); -An_app = An_app_f(); # need to run twice for unicode character... (bug) -include("apps/beer.jl"); -beer_app = Beer_app_f(); -beer_app = Beer_app_f(); -include("apps/hetero_resp.jl"); -Rh_app = Rh_app_f(); -Rh_app = Rh_app_f(); -include("apps/leaf_an_ci.jl"); -An_ci_app = An_ci_app_f(); -An_ci_app = An_ci_app_f(); -# Fluxnet DashBoard -fluxnet_app = ClimaLandSimulations.Dashboards.fluxnet_app(); -fluxnet_app = ClimaLandSimulations.Dashboards.fluxnet_app(); - -# Route apps -route!(server, "/leaf_An" => An_app) -route!(server, "/beer_APAR" => beer_app) -route!(server, "/leaf_An_ci" => An_ci_app) -route!(server, "/Rh" => Rh_app) -route!(server, "/fluxnet" => fluxnet_app) - -wait() +# Load packages +using ParamViz +using Bonito +using ClimaLand +using ClimaLand.Canopy +using ClimaLand.Soil.Biogeochemistry +import ClimaParams as CP +import ClimaLand.Parameters as LP +using ClimaLandSimulations +FT = Float64 +earth_param_set = LP.LandParameters(FT) +RTparams = BeerLambertParameters(FT) + +# Create server +IPa = "127.0.0.1" +port = 9385 +server = Server( + IPa, + port; + proxy_url = "https://clima.westus3.cloudapp.azure.com/jsserve/", +) + +# Load and create apps +include("apps/leaf_an.jl"); +An_app = An_app_f(); +An_app = An_app_f(); # need to run twice for unicode character... (bug) +include("apps/beer.jl"); +beer_app = Beer_app_f(); +beer_app = Beer_app_f(); +include("apps/hetero_resp.jl"); +Rh_app = Rh_app_f(); +Rh_app = Rh_app_f(); +include("apps/leaf_an_ci.jl"); +An_ci_app = An_ci_app_f(); +An_ci_app = An_ci_app_f(); +# Fluxnet DashBoard +fluxnet_app = ClimaLandSimulations.Dashboards.fluxnet_app(); +fluxnet_app = ClimaLandSimulations.Dashboards.fluxnet_app(); + +# Route apps +route!(server, "/leaf_An" => An_app) +route!(server, "/beer_APAR" => beer_app) +route!(server, "/leaf_An_ci" => An_ci_app) +route!(server, "/Rh" => Rh_app) +route!(server, "/fluxnet" => fluxnet_app) + +wait() diff --git a/docs/src/standalone/apps/beer.jl b/docs/src/standalone/apps/beer.jl index f16380b38e..266f3b0da7 100644 --- a/docs/src/standalone/apps/beer.jl +++ b/docs/src/standalone/apps/beer.jl @@ -1,35 +1,35 @@ -using Unitful: m, s, mol, μmol - -function ParamViz.parameterisation(PAR, LAI, ρ_leaf, K, Ω, a, b) - # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) - APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd - return APAR -end - -function Beer_app_f() - drivers = Drivers( - ("PAR (μmol m⁻² s⁻¹)", "LAI (m² m⁻²)"), - (FT.([0, 1500 * 1e-6]), FT.([0, 10])), - ((mol * m^-2 * s^-1, μmol * m^-2 * s^-1), (m^2 * m^-2, m^2 * m^-2)), - ) - - parameters = Parameters( - ( - "canopy reflectance, ρ_leaf", - "extinction coefficient, K", - "clumping index, Ω", - ), - (FT.([0, 1]), FT.([0, 1]), FT.([0, 1])), - ((m, m), (m, m), (m, m)), # dummy units, no conversion - ) - - constants = Constants(("a", "b"), (FT(1), FT(2))) # dummy constants - inputs = Inputs(drivers, parameters, constants) - output = Output( - "APAR (μmol m⁻² s⁻¹)", - [0, 1500 * 1e-6], - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), - ) - beer_app = webapp(ParamViz.parameterisation, inputs, output) - return beer_app -end +using Unitful: m, s, mol, μmol + +function ParamViz.parameterisation(PAR, LAI, ρ_leaf, K, Ω, a, b) + # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) + APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd + return APAR +end + +function Beer_app_f() + drivers = Drivers( + ("PAR (μmol m⁻² s⁻¹)", "LAI (m² m⁻²)"), + (FT.([0, 1500 * 1e-6]), FT.([0, 10])), + ((mol * m^-2 * s^-1, μmol * m^-2 * s^-1), (m^2 * m^-2, m^2 * m^-2)), + ) + + parameters = Parameters( + ( + "canopy reflectance, ρ_leaf", + "extinction coefficient, K", + "clumping index, Ω", + ), + (FT.([0, 1]), FT.([0, 1]), FT.([0, 1])), + ((m, m), (m, m), (m, m)), # dummy units, no conversion + ) + + constants = Constants(("a", "b"), (FT(1), FT(2))) # dummy constants + inputs = Inputs(drivers, parameters, constants) + output = Output( + "APAR (μmol m⁻² s⁻¹)", + [0, 1500 * 1e-6], + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), + ) + beer_app = webapp(ParamViz.parameterisation, inputs, output) + return beer_app +end diff --git a/docs/src/standalone/apps/hetero_resp.jl b/docs/src/standalone/apps/hetero_resp.jl index 134a5a7f00..0a6ef56fe6 100644 --- a/docs/src/standalone/apps/hetero_resp.jl +++ b/docs/src/standalone/apps/hetero_resp.jl @@ -1,92 +1,92 @@ -using Unitful: K, °C, mol, μmol, m, s, kg, J - -function ParamViz.parameterisation( - Tₛ, - θ, # drivers - ν, - α_sx, - Ea_sx, - kM_sx, - kM_O₂, - O₂_a, - p_sx, - Csom, # parameters - θ_a100, - D_liq, - D_oa, -) # constants - - params = SoilCO2ModelParameters(FT; ν) - # θ has to be lower than porosity - if θ > ν - θ = ν - end - - Rh = microbe_source(Tₛ, θ, Csom, params) - return Rh -end - -function Rh_app_f() - drivers = Drivers( - ("Tₛ (°C)", "θ (m³ m⁻³)"), # drivers.names - (FT.([273, 323]), FT.([0.0, 1.0])), # drivers.ranges - ((K, °C), (m^3 * m^-3, m^3 * m^-3)), - ) - - parameters = Parameters( - (# names - "Soil porosity, ν (m³ m⁻³)", - "Pre-exponential factor, α_sx (kg C m⁻³ s⁻¹)", - "Activation energy, Ea_sx (J mol⁻¹)", - "Michaelis constant for soil, kM_sx (kg C m⁻³)", - "Michaelis constant for O₂, kM_O₂ (m³ m⁻³)", - "Volumetric fraction of O₂ in the soil air, O₂_a (dimensionless)", - "Fraction of soil carbon that is considered soluble, p_sx (dimensionless)", - "Soil organic C, Csom (kg C m⁻³)", - ), - (# ranges - FT.([0.0, 1.0]), # porosity - FT.([100e3, 300e3]), # α_sx - FT.([50e3, 70e3]), # Ea_sx - FT.([1e-10, 0.1]), # kM_sx - FT.([1e-10, 0.1]), # kM_o2 - FT.([0.005, 0.5]), # O2_a - FT.([0.005, 0.5]), # p_sx - FT.([1.0, 10.0]), # Csom - ), - ( - (m^3 * m^-3, m^3 * m^-3), # porosity - (kg * m^-3 * s^-1, kg * m^-3 * s^-1), # α_sx - (J * mol^-1, J * mol^-1), # Ea_sx - (kg * m^-3, kg * m^-3), # kM_sx - (m^3 * m^-3, m^3 * m^-3), # kM_O2 - (m, m), # O2_a - (m, m), # p_sx - (kg * m^-3, kg * m^-3), # Csom - ), - ) - - constants = Constants( - (# names - "Air-filled porosity at soil water potential of -100 cm H₂O (~ 10 Pa)", - "Diffusivity of soil C substrate in liquid (unitless)", - "Diffusion coefficient of oxygen in air, dimensionless", - ), - (# values - FT(0.1816), # θ_a100 - FT(3.17), # D_liq - FT(1.67), # D_oa - ), - ) - - inputs = Inputs(drivers, parameters, constants) - - output = Output( - "Rh (mg C m⁻³ s⁻¹)", # name #NEED TO CONVERT THE UNIT - [0, 20 * 1e-6], # range - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), # unit from, unit to --> actually need to fix this, should be g to mg C m-3 s-1 - ) - - Rh_app = webapp(ParamViz.parameterisation, inputs, output) - return Rh_app -end +using Unitful: K, °C, mol, μmol, m, s, kg, J + +function ParamViz.parameterisation( + Tₛ, + θ, # drivers + ν, + α_sx, + Ea_sx, + kM_sx, + kM_O₂, + O₂_a, + p_sx, + Csom, # parameters + θ_a100, + D_liq, + D_oa, +) # constants + + params = SoilCO2ModelParameters(FT; ν) + # θ has to be lower than porosity + if θ > ν + θ = ν + end + + Rh = microbe_source(Tₛ, θ, Csom, params) + return Rh +end + +function Rh_app_f() + drivers = Drivers( + ("Tₛ (°C)", "θ (m³ m⁻³)"), # drivers.names + (FT.([273, 323]), FT.([0.0, 1.0])), # drivers.ranges + ((K, °C), (m^3 * m^-3, m^3 * m^-3)), + ) + + parameters = Parameters( + (# names + "Soil porosity, ν (m³ m⁻³)", + "Pre-exponential factor, α_sx (kg C m⁻³ s⁻¹)", + "Activation energy, Ea_sx (J mol⁻¹)", + "Michaelis constant for soil, kM_sx (kg C m⁻³)", + "Michaelis constant for O₂, kM_O₂ (m³ m⁻³)", + "Volumetric fraction of O₂ in the soil air, O₂_a (dimensionless)", + "Fraction of soil carbon that is considered soluble, p_sx (dimensionless)", + "Soil organic C, Csom (kg C m⁻³)", + ), + (# ranges + FT.([0.0, 1.0]), # porosity + FT.([100e3, 300e3]), # α_sx + FT.([50e3, 70e3]), # Ea_sx + FT.([1e-10, 0.1]), # kM_sx + FT.([1e-10, 0.1]), # kM_o2 + FT.([0.005, 0.5]), # O2_a + FT.([0.005, 0.5]), # p_sx + FT.([1.0, 10.0]), # Csom + ), + ( + (m^3 * m^-3, m^3 * m^-3), # porosity + (kg * m^-3 * s^-1, kg * m^-3 * s^-1), # α_sx + (J * mol^-1, J * mol^-1), # Ea_sx + (kg * m^-3, kg * m^-3), # kM_sx + (m^3 * m^-3, m^3 * m^-3), # kM_O2 + (m, m), # O2_a + (m, m), # p_sx + (kg * m^-3, kg * m^-3), # Csom + ), + ) + + constants = Constants( + (# names + "Air-filled porosity at soil water potential of -100 cm H₂O (~ 10 Pa)", + "Diffusivity of soil C substrate in liquid (unitless)", + "Diffusion coefficient of oxygen in air, dimensionless", + ), + (# values + FT(0.1816), # θ_a100 + FT(3.17), # D_liq + FT(1.67), # D_oa + ), + ) + + inputs = Inputs(drivers, parameters, constants) + + output = Output( + "Rh (mg C m⁻³ s⁻¹)", # name #NEED TO CONVERT THE UNIT + [0, 20 * 1e-6], # range + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), # unit from, unit to --> actually need to fix this, should be g to mg C m-3 s-1 + ) + + Rh_app = webapp(ParamViz.parameterisation, inputs, output) + return Rh_app +end diff --git a/docs/src/standalone/apps/leaf_an.jl b/docs/src/standalone/apps/leaf_an.jl index 43e27fd510..5960996a6c 100644 --- a/docs/src/standalone/apps/leaf_an.jl +++ b/docs/src/standalone/apps/leaf_an.jl @@ -1,137 +1,137 @@ -using Unitful: K, °C, mol, μmol, m, s, Pa, kPa - -function ParamViz.parameterisation( - PAR, - T, - β, - LAI, - ca, - VPD, - θs, - ld, - ρ_leaf, - Ω, - Γstar25, - ΔHΓstar, - To, - R, - Vcmax25, - ΔHJmax, - θj, - ϕ, - ΔHVcmax, - g1, - Kc25, - ΔHkc, - Ko25, - ΔHko, - oi, - f, - ΔHRd, -) - K = extinction_coeff(RTparams.G_Function, θs) - # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) - APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd - Jmax = max_electron_transport(Vcmax25, ΔHJmax, T, To, R) - J = electron_transport(APAR, Jmax, θj, ϕ) - Vcmax = compute_Vcmax(Vcmax25, T, To, R, ΔHVcmax) - Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R) - medlynt = 1 + g1 / sqrt(VPD) - ci = intercellular_co2(ca, Γstar, medlynt) - Aj = light_assimilation(Canopy.C3(), J, ci, Γstar) - Kc = MM_Kc(Kc25, ΔHkc, T, To, R) - Ko = MM_Ko(Ko25, ΔHko, T, To, R) - Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi) - Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R) - An = net_photosynthesis(Ac, Aj, Rd, β) - return An -end - -function An_app_f() - drivers = Drivers( - ("PAR (μmol m⁻² s⁻¹)", "T (°C)"), - (FT.([0, 1500 * 1e-6]), FT.([273, 323])), - ((mol * m^-2 * s^-1, μmol * m^-2 * s^-1), (K, °C)), - ) - - parameters = Parameters( - ( - "Moisture stress, β", - "Leaf area index, LAI (m² m⁻²)", - "CO2 concentration, ca (ppm)", - "Vapor pressure deficit, VPD (Pa)", - ), - ( - FT.([0, 1]), # β - FT.([1, 10]), # LAI, m2 m-2 - FT.([300 * 1e-6, 500 * 1e-6]), # ca - FT.([500, 10000]), # VPD, Pa - ), - ( - (m, m), # dummy units, no conversion - (m^2 * m^-2, m^2 * m^-2), - (mol / mol, μmol / mol), - (Pa, kPa), - ), - ) - - constants = Constants( - ( - "θs", # Sun zenith angle - "ld", # Leaf angle distribution - "ρ_leaf", # PAR canopy reflectance - "Ω", # Clumping index - "Γstar25", # co2 compensation at 25c - "ΔHΓstar", # a constant energy of activation - "To", # a standard temperature - "R", # the universal gas constant - "Vcmax25", # the maximum rate of carboxylation of Rubisco - "ΔHJmax", # a constant - "θj", # an empirical "curvature parameter" - "ϕ", # the quantum yield of photosystem II - "ΔHVcmax", # a constant - "g1", # a constant - "Kc25", # a constant - "ΔHkc", # a constant - "Ko25", # a constant - "ΔHko", # a constant - "oi", # an empirical parameter - "f", # an empirical factor - "ΔHRd", - ), - ( - FT(0.6), # θs - is that a good value? in radian, right? -------- - FT(0.5), # ld - ozark val - FT(0.1), # ρ_leaf - ozark val - FT(0.69), # Ω - ozark val - FT(4.275e-5), # Γstar25 - ozark val - FT(37830), # ΔHΓstar - ozark val - FT(298.15), # To - ozark vak - FT(8.314), # R, J/mol - FT(LSMP.gas_constant(earth_param_set)) - FT(5e-5), # Vcmax25 - ozark model - FT(43540), # ΔHJmax - ozark model - FT(0.9), # θj - ozark model - FT(0.6), # ϕ - ozark model - FT(58520), # ΔHVcmax - ozark model - FT(141), # g1 - ozark model - FT(4.049e-4), # Kc25 - ozark model - FT(79430), # ΔHkc - ozark model - FT(0.2874), # Ko25 - ozark model - FT(36380), # ΔHko - ozark model - FT(0.209), # oi - ozark model - FT(0.015), # f - ozark model - FT(43390), # ΔHRd - ozark model - ), - ) - - inputs = Inputs(drivers, parameters, constants) - - output = Output( - "An (μmol m⁻² s⁻¹)", - [0, 20 * 1e-6], - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), - ) - - An_app = webapp(ParamViz.parameterisation, inputs, output) - return An_app -end +using Unitful: K, °C, mol, μmol, m, s, Pa, kPa + +function ParamViz.parameterisation( + PAR, + T, + β, + LAI, + ca, + VPD, + θs, + ld, + ρ_leaf, + Ω, + Γstar25, + ΔHΓstar, + To, + R, + Vcmax25, + ΔHJmax, + θj, + ϕ, + ΔHVcmax, + g1, + Kc25, + ΔHkc, + Ko25, + ΔHko, + oi, + f, + ΔHRd, +) + K = extinction_coeff(RTparams.G_Function, θs) + # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) + APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd + Jmax = max_electron_transport(Vcmax25, ΔHJmax, T, To, R) + J = electron_transport(APAR, Jmax, θj, ϕ) + Vcmax = compute_Vcmax(Vcmax25, T, To, R, ΔHVcmax) + Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R) + medlynt = 1 + g1 / sqrt(VPD) + ci = intercellular_co2(ca, Γstar, medlynt) + Aj = light_assimilation(Canopy.C3(), J, ci, Γstar) + Kc = MM_Kc(Kc25, ΔHkc, T, To, R) + Ko = MM_Ko(Ko25, ΔHko, T, To, R) + Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi) + Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R) + An = net_photosynthesis(Ac, Aj, Rd, β) + return An +end + +function An_app_f() + drivers = Drivers( + ("PAR (μmol m⁻² s⁻¹)", "T (°C)"), + (FT.([0, 1500 * 1e-6]), FT.([273, 323])), + ((mol * m^-2 * s^-1, μmol * m^-2 * s^-1), (K, °C)), + ) + + parameters = Parameters( + ( + "Moisture stress, β", + "Leaf area index, LAI (m² m⁻²)", + "CO2 concentration, ca (ppm)", + "Vapor pressure deficit, VPD (Pa)", + ), + ( + FT.([0, 1]), # β + FT.([1, 10]), # LAI, m2 m-2 + FT.([300 * 1e-6, 500 * 1e-6]), # ca + FT.([500, 10000]), # VPD, Pa + ), + ( + (m, m), # dummy units, no conversion + (m^2 * m^-2, m^2 * m^-2), + (mol / mol, μmol / mol), + (Pa, kPa), + ), + ) + + constants = Constants( + ( + "θs", # Sun zenith angle + "ld", # Leaf angle distribution + "ρ_leaf", # PAR canopy reflectance + "Ω", # Clumping index + "Γstar25", # co2 compensation at 25c + "ΔHΓstar", # a constant energy of activation + "To", # a standard temperature + "R", # the universal gas constant + "Vcmax25", # the maximum rate of carboxylation of Rubisco + "ΔHJmax", # a constant + "θj", # an empirical "curvature parameter" + "ϕ", # the quantum yield of photosystem II + "ΔHVcmax", # a constant + "g1", # a constant + "Kc25", # a constant + "ΔHkc", # a constant + "Ko25", # a constant + "ΔHko", # a constant + "oi", # an empirical parameter + "f", # an empirical factor + "ΔHRd", + ), + ( + FT(0.6), # θs - is that a good value? in radian, right? -------- + FT(0.5), # ld - ozark val + FT(0.1), # ρ_leaf - ozark val + FT(0.69), # Ω - ozark val + FT(4.275e-5), # Γstar25 - ozark val + FT(37830), # ΔHΓstar - ozark val + FT(298.15), # To - ozark vak + FT(8.314), # R, J/mol - FT(LSMP.gas_constant(earth_param_set)) + FT(5e-5), # Vcmax25 - ozark model + FT(43540), # ΔHJmax - ozark model + FT(0.9), # θj - ozark model + FT(0.6), # ϕ - ozark model + FT(58520), # ΔHVcmax - ozark model + FT(141), # g1 - ozark model + FT(4.049e-4), # Kc25 - ozark model + FT(79430), # ΔHkc - ozark model + FT(0.2874), # Ko25 - ozark model + FT(36380), # ΔHko - ozark model + FT(0.209), # oi - ozark model + FT(0.015), # f - ozark model + FT(43390), # ΔHRd - ozark model + ), + ) + + inputs = Inputs(drivers, parameters, constants) + + output = Output( + "An (μmol m⁻² s⁻¹)", + [0, 20 * 1e-6], + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), + ) + + An_app = webapp(ParamViz.parameterisation, inputs, output) + return An_app +end diff --git a/docs/src/standalone/apps/leaf_an_ci.jl b/docs/src/standalone/apps/leaf_an_ci.jl index 00acc1f35a..351338a131 100644 --- a/docs/src/standalone/apps/leaf_an_ci.jl +++ b/docs/src/standalone/apps/leaf_an_ci.jl @@ -1,127 +1,127 @@ -using Unitful: K, °C, mol, μmol, m, s - -function ParamViz.parameterisation( - ci, - T, - β, - LAI, - PAR, - Vcmax25, - θs, - ld, - ρ_leaf, - Ω, - Γstar25, - ΔHΓstar, - To, - R, - ΔHJmax, - θj, - ϕ, - ΔHVcmax, - Kc25, - ΔHkc, - Ko25, - ΔHko, - oi, - f, - ΔHRd, -) - K = extinction_coeff(RTparams.G_Function, θs) - # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) - APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd - Jmax = max_electron_transport(Vcmax25, ΔHJmax, T, To, R) - J = electron_transport(APAR, Jmax, θj, ϕ) - Vcmax = compute_Vcmax(Vcmax25, T, To, R, ΔHVcmax) - Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R) - Kc = MM_Kc(Kc25, ΔHkc, T, To, R) - Ko = MM_Ko(Ko25, ΔHko, T, To, R) - Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi) - Aj = light_assimilation(Canopy.C3(), J, ci, Γstar) - Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R) - An = net_photosynthesis(Ac, Aj, Rd, β) - return An -end - -function An_ci_app_f() - drivers = Drivers( - ("ci (ppm)", "T (°C)"), # names - (FT.([0.0001, 0.001]), FT.([273, 323])), # ranges - ((mol / mol, μmol / mol), (K, °C)), # units from, unit to - ) - - parameters = Parameters( - ( - "Moisture stress, β", - "Leaf area index, LAI (m² m⁻²)", - "Photosynthetic radiation, PAR (μmol m⁻² s⁻¹)", - "Vcmax25, (μmol m⁻² s⁻¹)", - ), - ( - FT.([0, 1]), # β - FT.([1, 10]), # LAI, m2 m-2 - FT.([0, 1500 * 1e-6]), # PAR - FT.([1e-5, 1e-4]), # Vcmax25 - ), - ( - (m, m), # dummy unit, no conversion - (m^2 * m^-2, m^2 * m^-2), - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), - ), - ) - - constants = Constants( - ( - "θs", # Sun zenith angle - "ld", # Leaf angle distribution - "ρ_leaf", # PAR canopy reflectance - "Ω", # Clumping index - "Γstar25", # co2 compensation at 25c - "ΔHΓstar", # a constant energy of activation - "To", # a standard temperature - "R", # the universal gas constant - "ΔHJmax", # a constant - "θj", # an empirical "curvature parameter" - "ϕ", # the quantum yield of photosystem II - "ΔHVcmax", # a constant - "Kc25", # a constant - "ΔHkc", # a constant - "Ko25", # a constant - "ΔHko", # a constant - "oi", # an empirical parameter - "f", # an empirical factor - "ΔHRd",# a constant - ), - ( - FT(0.6), # θs - is that a good value? in radian, right? -------- - FT(0.5), # ld - ozark val - FT(0.1), # ρ_leaf - ozark val - FT(0.69), # Ω - ozark val - FT(4.275e-5), # Γstar25 - ozark val - FT(37830), # ΔHΓstar - ozark val - FT(298.15), # To - ozark vak - FT(8.314), # R, J/mol - FT(LSMP.gas_constant(earth_param_set)) - FT(43540), # ΔHJmax - ozark model - FT(0.9), # θj - ozark model - FT(0.6), # ϕ - ozark model - FT(58520), # ΔHVcmax - ozark model - FT(4.049e-4), # Kc25 - ozark model - FT(79430), # ΔHkc - ozark model - FT(0.2874), # Ko25 - ozark model - FT(36380), # ΔHko - ozark model - FT(0.209), # oi - ozark model - FT(0.015), # f - ozark model - FT(43390), # ΔHRd - ozark model - ), - ) - - inputs = Inputs(drivers, parameters, constants) - output = Output( - "An (μmol m⁻² s⁻¹)", - [0, 20 * 1e-6], - (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), - ) - An_ci_app = webapp(ParamViz.parameterisation, inputs, output) - return An_ci_app -end +using Unitful: K, °C, mol, μmol, m, s + +function ParamViz.parameterisation( + ci, + T, + β, + LAI, + PAR, + Vcmax25, + θs, + ld, + ρ_leaf, + Ω, + Γstar25, + ΔHΓstar, + To, + R, + ΔHJmax, + θj, + ϕ, + ΔHVcmax, + Kc25, + ΔHkc, + Ko25, + ΔHko, + oi, + f, + ΔHRd, +) + K = extinction_coeff(RTparams.G_Function, θs) + # APAR = plant_absorbed_ppfd(PAR, ρ_leaf, K, LAI, Ω) + APAR = PAR * (1 - ρ_leaf) * (1 - exp(-K * LAI * Ω)) # not sure how to call new plant_absorbed_pfd + Jmax = max_electron_transport(Vcmax25, ΔHJmax, T, To, R) + J = electron_transport(APAR, Jmax, θj, ϕ) + Vcmax = compute_Vcmax(Vcmax25, T, To, R, ΔHVcmax) + Γstar = co2_compensation(Γstar25, ΔHΓstar, T, To, R) + Kc = MM_Kc(Kc25, ΔHkc, T, To, R) + Ko = MM_Ko(Ko25, ΔHko, T, To, R) + Ac = rubisco_assimilation(Canopy.C3(), Vcmax, ci, Γstar, Kc, Ko, oi) + Aj = light_assimilation(Canopy.C3(), J, ci, Γstar) + Rd = dark_respiration(Vcmax25, β, f, ΔHRd, T, To, R) + An = net_photosynthesis(Ac, Aj, Rd, β) + return An +end + +function An_ci_app_f() + drivers = Drivers( + ("ci (ppm)", "T (°C)"), # names + (FT.([0.0001, 0.001]), FT.([273, 323])), # ranges + ((mol / mol, μmol / mol), (K, °C)), # units from, unit to + ) + + parameters = Parameters( + ( + "Moisture stress, β", + "Leaf area index, LAI (m² m⁻²)", + "Photosynthetic radiation, PAR (μmol m⁻² s⁻¹)", + "Vcmax25, (μmol m⁻² s⁻¹)", + ), + ( + FT.([0, 1]), # β + FT.([1, 10]), # LAI, m2 m-2 + FT.([0, 1500 * 1e-6]), # PAR + FT.([1e-5, 1e-4]), # Vcmax25 + ), + ( + (m, m), # dummy unit, no conversion + (m^2 * m^-2, m^2 * m^-2), + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), + ), + ) + + constants = Constants( + ( + "θs", # Sun zenith angle + "ld", # Leaf angle distribution + "ρ_leaf", # PAR canopy reflectance + "Ω", # Clumping index + "Γstar25", # co2 compensation at 25c + "ΔHΓstar", # a constant energy of activation + "To", # a standard temperature + "R", # the universal gas constant + "ΔHJmax", # a constant + "θj", # an empirical "curvature parameter" + "ϕ", # the quantum yield of photosystem II + "ΔHVcmax", # a constant + "Kc25", # a constant + "ΔHkc", # a constant + "Ko25", # a constant + "ΔHko", # a constant + "oi", # an empirical parameter + "f", # an empirical factor + "ΔHRd",# a constant + ), + ( + FT(0.6), # θs - is that a good value? in radian, right? -------- + FT(0.5), # ld - ozark val + FT(0.1), # ρ_leaf - ozark val + FT(0.69), # Ω - ozark val + FT(4.275e-5), # Γstar25 - ozark val + FT(37830), # ΔHΓstar - ozark val + FT(298.15), # To - ozark vak + FT(8.314), # R, J/mol - FT(LSMP.gas_constant(earth_param_set)) + FT(43540), # ΔHJmax - ozark model + FT(0.9), # θj - ozark model + FT(0.6), # ϕ - ozark model + FT(58520), # ΔHVcmax - ozark model + FT(4.049e-4), # Kc25 - ozark model + FT(79430), # ΔHkc - ozark model + FT(0.2874), # Ko25 - ozark model + FT(36380), # ΔHko - ozark model + FT(0.209), # oi - ozark model + FT(0.015), # f - ozark model + FT(43390), # ΔHRd - ozark model + ), + ) + + inputs = Inputs(drivers, parameters, constants) + output = Output( + "An (μmol m⁻² s⁻¹)", + [0, 20 * 1e-6], + (mol * m^-2 * s^-1, μmol * m^-2 * s^-1), + ) + An_ci_app = webapp(ParamViz.parameterisation, inputs, output) + return An_ci_app +end diff --git a/lib/ClimaLandSimulations/README.md b/lib/ClimaLandSimulations/README.md index 9bc8eaf0fa..448368d0dc 100644 --- a/lib/ClimaLandSimulations/README.md +++ b/lib/ClimaLandSimulations/README.md @@ -15,7 +15,7 @@ julia> using ClimaLandSimulations For examples, see the [experiments](https://github.com/CliMA/ClimaLand.jl/lib/ClimaLandSimulations/experiments) folder -### Development of the `ClimaLandSimulations` subpackage +### Development of the `ClimaLandSimulations` subpackage cd ClimaLand/lib/ClimaLandSimulations diff --git a/lib/ClimaLandSimulations/src/Dashboards/Dashboards.jl b/lib/ClimaLandSimulations/src/Dashboards/Dashboards.jl index 1356a25da7..57d8521e6e 100644 --- a/lib/ClimaLandSimulations/src/Dashboards/Dashboards.jl +++ b/lib/ClimaLandSimulations/src/Dashboards/Dashboards.jl @@ -1,12 +1,12 @@ -module Dashboards - -using ClimaLandSimulations -using ClimaLandSimulations.Fluxnet -using Bonito -using WGLMakie -import ClimaLand.Parameters as LP -FT = Float64 - -include("fluxnet_dashboard.jl") - -end +module Dashboards + +using ClimaLandSimulations +using ClimaLandSimulations.Fluxnet +using Bonito +using WGLMakie +import ClimaLand.Parameters as LP +FT = Float64 + +include("fluxnet_dashboard.jl") + +end diff --git a/lib/ClimaLandSimulations/src/Dashboards/fluxnet_dashboard.jl b/lib/ClimaLandSimulations/src/Dashboards/fluxnet_dashboard.jl index d33d55d1e0..1cbb3f8d7f 100644 --- a/lib/ClimaLandSimulations/src/Dashboards/fluxnet_dashboard.jl +++ b/lib/ClimaLandSimulations/src/Dashboards/fluxnet_dashboard.jl @@ -1,51 +1,51 @@ -export fluxnet_dashboard, fluxnet_app - -""" - fluxnet_dashboard(menu, menu2) - -Return selected figure from selected site. -""" -function fluxnet_dashboard(menu, menu2) - display_fig = map(menu.value) do value - @info value - sv = run_fluxnet(value)[1] - simulation_inputs = make_inputs_df(value)[1] - simulation_outputs = make_output_df(value, sv, simulation_inputs) - figs = make_plots( - simulation_inputs, - simulation_outputs; - save_fig = false, - dashboard = true, - ) - return Dict( - "timeseries" => figs.timeseries, - "water" => figs.water, - "fingerprint" => figs.fingerprint, - "diurnals" => figs.diurnals, - ) - end - - fig = map(menu2.value) do value - @info value - return get(display_fig[], value, nothing) - end - return fig -end - -""" - fluxnet_app() - -Make a web dashboard to interact with fluxnet_dashboard(menu, menu2) -""" -function fluxnet_app() - ClimaLand_dashboard = App() do - sites_ID = ["US-MOz", "US-Ha1", "US-NR1", "US-Var"] - menu = Dropdown(sites_ID) - displayed_figure = - ["timeseries", "water", "fingerprint", "diurnals"] - menu2 = Dropdown(displayed_figure) - fig = fluxnet_dashboard(menu, menu2) - return DOM.div(menu, menu2, fig) - end - return ClimaLand_dashboard -end +export fluxnet_dashboard, fluxnet_app + +""" + fluxnet_dashboard(menu, menu2) + +Return selected figure from selected site. +""" +function fluxnet_dashboard(menu, menu2) + display_fig = map(menu.value) do value + @info value + sv = run_fluxnet(value)[1] + simulation_inputs = make_inputs_df(value)[1] + simulation_outputs = make_output_df(value, sv, simulation_inputs) + figs = make_plots( + simulation_inputs, + simulation_outputs; + save_fig = false, + dashboard = true, + ) + return Dict( + "timeseries" => figs.timeseries, + "water" => figs.water, + "fingerprint" => figs.fingerprint, + "diurnals" => figs.diurnals, + ) + end + + fig = map(menu2.value) do value + @info value + return get(display_fig[], value, nothing) + end + return fig +end + +""" + fluxnet_app() + +Make a web dashboard to interact with fluxnet_dashboard(menu, menu2) +""" +function fluxnet_app() + ClimaLand_dashboard = App() do + sites_ID = ["US-MOz", "US-Ha1", "US-NR1", "US-Var"] + menu = Dropdown(sites_ID) + displayed_figure = + ["timeseries", "water", "fingerprint", "diurnals"] + menu2 = Dropdown(displayed_figure) + fig = fluxnet_dashboard(menu, menu2) + return DOM.div(menu, menu2, fig) + end + return ClimaLand_dashboard +end diff --git a/lib/ClimaLandSimulations/src/Fluxnet/Fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/Fluxnet.jl index c1b7f74b0c..117e7c32d5 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/Fluxnet.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/Fluxnet.jl @@ -1,59 +1,59 @@ -module Fluxnet - -using ClimaLandSimulations -import SciMLBase -import ClimaTimeSteppers as CTS -using ClimaCore -import ClimaParams as CP -using Statistics -using Dates -using Insolation -using StatsBase -using Interpolations -using ClimaLand -using ClimaLand.Domains: Column -using ClimaLand.Soil -using ClimaLand.Soil.Biogeochemistry -using ClimaLand.Canopy -using ClimaLand.Canopy.PlantHydraulics -import ClimaLand.Parameters as LP -import ClimaComms -using Unitful: R, L, mol, K, kJ, °C, m, g, cm, hr, mg, s, μmol, Pa, W, mm, kPa -using UnitfulMoles: molC -using Unitful, UnitfulMoles -using ArtifactWrappers -using DelimitedFiles -using Thermodynamics -using Formatting -using DataFrames - -# Directory paths -climalandsimulations_dir = pkgdir(ClimaLandSimulations) -climaland_dir = pkgdir(ClimaLand) - -# ClimaParams -const FT = Float64 -earth_param_set = LP.LandParameters(FT) - -include("run_fluxnet.jl") - -sites = ["US-MOz", "US-Ha1", "US-NR1", "US-Var"] -for site in sites - include(joinpath("fluxnet_sites", site, "$(site)_simulation.jl")) - include(joinpath("fluxnet_sites", site, "$(site)_config.jl")) - include(joinpath("fluxnet_sites", site, "$(site)_parameters.jl")) -end - -include(joinpath("fluxnet_utilities", "make_config.jl")) -include(joinpath("fluxnet_utilities", "make_setup.jl")) -include(joinpath("fluxnet_utilities", "make_parameters.jl")) -include(joinpath("fluxnet_utilities", "data_tools.jl")) -include(joinpath("fluxnet_utilities", "make_drivers.jl")) -include(joinpath("fluxnet_utilities", "inputs_dataframe.jl")) - -# Register Fluxnet with Unitful.jl -function __init__() - Unitful.register(Fluxnet) -end - -end +module Fluxnet + +using ClimaLandSimulations +import SciMLBase +import ClimaTimeSteppers as CTS +using ClimaCore +import ClimaParams as CP +using Statistics +using Dates +using Insolation +using StatsBase +using Interpolations +using ClimaLand +using ClimaLand.Domains: Column +using ClimaLand.Soil +using ClimaLand.Soil.Biogeochemistry +using ClimaLand.Canopy +using ClimaLand.Canopy.PlantHydraulics +import ClimaLand.Parameters as LP +import ClimaComms +using Unitful: R, L, mol, K, kJ, °C, m, g, cm, hr, mg, s, μmol, Pa, W, mm, kPa +using UnitfulMoles: molC +using Unitful, UnitfulMoles +using ArtifactWrappers +using DelimitedFiles +using Thermodynamics +using Formatting +using DataFrames + +# Directory paths +climalandsimulations_dir = pkgdir(ClimaLandSimulations) +climaland_dir = pkgdir(ClimaLand) + +# ClimaParams +const FT = Float64 +earth_param_set = LP.LandParameters(FT) + +include("run_fluxnet.jl") + +sites = ["US-MOz", "US-Ha1", "US-NR1", "US-Var"] +for site in sites + include(joinpath("fluxnet_sites", site, "$(site)_simulation.jl")) + include(joinpath("fluxnet_sites", site, "$(site)_config.jl")) + include(joinpath("fluxnet_sites", site, "$(site)_parameters.jl")) +end + +include(joinpath("fluxnet_utilities", "make_config.jl")) +include(joinpath("fluxnet_utilities", "make_setup.jl")) +include(joinpath("fluxnet_utilities", "make_parameters.jl")) +include(joinpath("fluxnet_utilities", "data_tools.jl")) +include(joinpath("fluxnet_utilities", "make_drivers.jl")) +include(joinpath("fluxnet_utilities", "inputs_dataframe.jl")) + +# Register Fluxnet with Unitful.jl +function __init__() + Unitful.register(Fluxnet) +end + +end diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl index a7fc8a915e..10c61b4d0b 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Ha1/US-Ha1_parameters.jl @@ -1,208 +1,208 @@ -export harvard_default_params, - hetero_resp_harvard, - auto_resp_harvard, - soil_harvard, - radiative_transfer_harvard, - canopy_energy_balance_harvard, - conductance_harvard, - photosynthesis_harvard, - plant_hydraulics_harvard - -function harvard_default_params(; - hetero_resp = hetero_resp_harvard(), - auto_resp = auto_resp_harvard(), - soil = soil_harvard(), - radiative_transfer = radiative_transfer_harvard(), - canopy_energy_balance = canopy_energy_balance_harvard(), - conductance = conductance_harvard(), - photosynthesis = photosynthesis_harvard(), - plant_hydraulics = plant_hydraulics_harvard(), -) - return SiteParams( - hetero_resp, - auto_resp, - soil, - radiative_transfer, - canopy_energy_balance, - conductance, - photosynthesis, - plant_hydraulics, - ) -end - -function hetero_resp_harvard(; - D_ref = FT(1.39e-5), - D_liq = FT(3.17), - α_sx = FT(194e3), - Ea_sx = FT(61e3), - kM_sx = FT(5e-3), - kM_o2 = FT(0.004), - O2_a = FT(0.209), - D_oa = FT(1.67), - p_sx = FT(0.024), -) - return HeteroRespP( - D_ref, - D_liq, - α_sx, - Ea_sx, - kM_sx, - kM_o2, - O2_a, - D_oa, - p_sx, - ) -end - -function auto_resp_harvard(; - ne = FT(8 * 1e-4), - ηsl = FT(0.01), - σl = FT(0.05), - μr = FT(1.0), - μs = FT(0.1), - f1 = FT(0.012), - f2 = FT(0.25), -) - return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) -end - -function soil_harvard(; - vg_α = FT(0.04), - vg_n = FT(2.05), - ν = FT(0.5), # m3/m3 - ν_ss_quartz = FT(0.1), - ν_ss_om = FT(0.1), - K_sat = FT(4e-7), # m/s, matches Natan - S_s = FT(1e-3), # 1/m, guess - θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 - ν_ss_gravel = FT(0.0), - z_0m = FT(0.01), - z_0b = FT(0.001), - emissivity = FT(0.98), - PAR_albedo = FT(0.2), - NIR_albedo = FT(0.2), -) - hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) - return EnergyHydrologyParameters( - FT; - ν, - ν_ss_om, - ν_ss_quartz, - ν_ss_gravel, - hydrology_cm, - K_sat, - S_s, - θ_r, - PAR_albedo, - NIR_albedo, - emissivity, - z_0m, - z_0b, - ) -end - -function radiative_transfer_harvard(; - Ω = FT(0.69), - ld = ConstantGFunction(FT(0.5)), - α_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), - n_layers = UInt64(20), - ϵ_canopy = FT(0.97), -) - return TwoStreamParameters( - G_Function = ld, - α_PAR_leaf = α_PAR_leaf, - τ_PAR_leaf = τ_PAR_leaf, - α_NIR_leaf = α_NIR_leaf, - τ_NIR_leaf = τ_NIR_leaf, - ϵ_canopy = ϵ_canopy, - Ω = Ω, - λ_γ_PAR = λ_γ_PAR, - λ_γ_NIR = λ_γ_NIR, - n_layers = n_layers, - ) -end - -function canopy_energy_balance_harvard(; ac_canopy = FT(2.5e3)) - return BigLeafEnergyParameters(ac_canopy) -end - -function conductance_harvard(; - g1 = FT(141), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. - Drel = FT(1.6), - g0 = FT(1e-4), -) - return MedlynConductanceParameters(Drel, g0, g1) -end - -function photosynthesis_harvard(; - mechanism = Canopy.C3(), - oi = FT(0.209), - ϕ = FT(0.6), - θj = FT(0.9), - f = FT(0.015), - sc = FT(2e-6), # Bonan's book: range of 2-5e-6 - pc = FT(-2e6), # Bonan's book: -2e6 - Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 - Γstar25 = FT(4.275e-5), - Kc25 = FT(4.049e-4), - Ko25 = FT(0.2874), - To = FT(298.15), - ΔHkc = FT(79430), - ΔHko = FT(36380), - ΔHVcmax = FT(58520), - ΔHΓstar = FT(37830), - ΔHJmax = FT(43540), - ΔHRd = FT(46390), -) - return FarquharParameters( - Vcmax25, - Γstar25, - Kc25, - Ko25, - ΔHkc, - ΔHko, - ΔHVcmax, - ΔHΓstar, - ΔHJmax, - ΔHRd, - To, - oi, - ϕ, - θj, - f, - sc, - pc, - mechanism, - ) -end - -function plant_hydraulics_harvard(; - SAI = FT(1.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? - f_root_to_shoot = FT(3.5), - 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 - capacity = FT(10), # kg/m^2 - S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m - rooting_depth = FT(0.5), # from Natan -) - conductivity_model = - PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) - retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) - - return PlantHydraulicsP( - SAI, - f_root_to_shoot, - capacity, - S_s, - rooting_depth, - conductivity_model, - retention_model, - ) -end +export harvard_default_params, + hetero_resp_harvard, + auto_resp_harvard, + soil_harvard, + radiative_transfer_harvard, + canopy_energy_balance_harvard, + conductance_harvard, + photosynthesis_harvard, + plant_hydraulics_harvard + +function harvard_default_params(; + hetero_resp = hetero_resp_harvard(), + auto_resp = auto_resp_harvard(), + soil = soil_harvard(), + radiative_transfer = radiative_transfer_harvard(), + canopy_energy_balance = canopy_energy_balance_harvard(), + conductance = conductance_harvard(), + photosynthesis = photosynthesis_harvard(), + plant_hydraulics = plant_hydraulics_harvard(), +) + return SiteParams( + hetero_resp, + auto_resp, + soil, + radiative_transfer, + canopy_energy_balance, + conductance, + photosynthesis, + plant_hydraulics, + ) +end + +function hetero_resp_harvard(; + D_ref = FT(1.39e-5), + D_liq = FT(3.17), + α_sx = FT(194e3), + Ea_sx = FT(61e3), + kM_sx = FT(5e-3), + kM_o2 = FT(0.004), + O2_a = FT(0.209), + D_oa = FT(1.67), + p_sx = FT(0.024), +) + return HeteroRespP( + D_ref, + D_liq, + α_sx, + Ea_sx, + kM_sx, + kM_o2, + O2_a, + D_oa, + p_sx, + ) +end + +function auto_resp_harvard(; + ne = FT(8 * 1e-4), + ηsl = FT(0.01), + σl = FT(0.05), + μr = FT(1.0), + μs = FT(0.1), + f1 = FT(0.012), + f2 = FT(0.25), +) + return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) +end + +function soil_harvard(; + vg_α = FT(0.04), + vg_n = FT(2.05), + ν = FT(0.5), # m3/m3 + ν_ss_quartz = FT(0.1), + ν_ss_om = FT(0.1), + K_sat = FT(4e-7), # m/s, matches Natan + S_s = FT(1e-3), # 1/m, guess + θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + ν_ss_gravel = FT(0.0), + z_0m = FT(0.01), + z_0b = FT(0.001), + emissivity = FT(0.98), + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), +) + hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) + return EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo, + NIR_albedo, + emissivity, + z_0m, + z_0b, + ) +end + +function radiative_transfer_harvard(; + Ω = FT(0.69), + ld = ConstantGFunction(FT(0.5)), + α_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), + n_layers = UInt64(20), + ϵ_canopy = FT(0.97), +) + return TwoStreamParameters( + G_Function = ld, + α_PAR_leaf = α_PAR_leaf, + τ_PAR_leaf = τ_PAR_leaf, + α_NIR_leaf = α_NIR_leaf, + τ_NIR_leaf = τ_NIR_leaf, + ϵ_canopy = ϵ_canopy, + Ω = Ω, + λ_γ_PAR = λ_γ_PAR, + λ_γ_NIR = λ_γ_NIR, + n_layers = n_layers, + ) +end + +function canopy_energy_balance_harvard(; ac_canopy = FT(2.5e3)) + return BigLeafEnergyParameters(ac_canopy) +end + +function conductance_harvard(; + g1 = FT(141), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + Drel = FT(1.6), + g0 = FT(1e-4), +) + return MedlynConductanceParameters(Drel, g0, g1) +end + +function photosynthesis_harvard(; + mechanism = Canopy.C3(), + oi = FT(0.209), + ϕ = FT(0.6), + θj = FT(0.9), + f = FT(0.015), + sc = FT(2e-6), # Bonan's book: range of 2-5e-6 + pc = FT(-2e6), # Bonan's book: -2e6 + Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 + Γstar25 = FT(4.275e-5), + Kc25 = FT(4.049e-4), + Ko25 = FT(0.2874), + To = FT(298.15), + ΔHkc = FT(79430), + ΔHko = FT(36380), + ΔHVcmax = FT(58520), + ΔHΓstar = FT(37830), + ΔHJmax = FT(43540), + ΔHRd = FT(46390), +) + return FarquharParameters( + Vcmax25, + Γstar25, + Kc25, + Ko25, + ΔHkc, + ΔHko, + ΔHVcmax, + ΔHΓstar, + ΔHJmax, + ΔHRd, + To, + oi, + ϕ, + θj, + f, + sc, + pc, + mechanism, + ) +end + +function plant_hydraulics_harvard(; + SAI = FT(1.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? + f_root_to_shoot = FT(3.5), + 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 + capacity = FT(10), # kg/m^2 + S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m + rooting_depth = FT(0.5), # from Natan +) + conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) + + return PlantHydraulicsP( + SAI, + f_root_to_shoot, + capacity, + S_s, + rooting_depth, + conductivity_model, + retention_model, + ) +end diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl index 311552a439..b798ffd179 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-NR1/US-NR1_parameters.jl @@ -1,208 +1,208 @@ -export niwotridge_default_params, - hetero_resp_niwotridge, - auto_resp_niwotridge, - soil_niwotridge, - radiative_transfer_niwotridge, - canopy_energy_balance_niwotridge, - conductance_niwotridge, - photosynthesis_niwotridge, - plant_hydraulics_niwotridge - -function niwotridge_default_params(; - hetero_resp = hetero_resp_niwotridge(), - auto_resp = auto_resp_niwotridge(), - soil = soil_niwotridge(), - radiative_transfer = radiative_transfer_niwotridge(), - canopy_energy_balance = canopy_energy_balance_niwotridge(), - conductance = conductance_niwotridge(), - photosynthesis = photosynthesis_niwotridge(), - plant_hydraulics = plant_hydraulics_niwotridge(), -) - return SiteParams( - hetero_resp, - auto_resp, - soil, - radiative_transfer, - canopy_energy_balance, - conductance, - photosynthesis, - plant_hydraulics, - ) -end - -function hetero_resp_niwotridge(; - D_ref = FT(1.39e-5), - D_liq = FT(3.17), - α_sx = FT(194e3), - Ea_sx = FT(61e3), - kM_sx = FT(5e-3), - kM_o2 = FT(0.004), - O2_a = FT(0.209), - D_oa = FT(1.67), - p_sx = FT(0.024), -) - return HeteroRespP( - D_ref, - D_liq, - α_sx, - Ea_sx, - kM_sx, - kM_o2, - O2_a, - D_oa, - p_sx, - ) -end - -function auto_resp_niwotridge(; - ne = FT(8 * 1e-4), - ηsl = FT(0.01), - σl = FT(0.05), - μr = FT(1.0), - μs = FT(0.1), - f1 = FT(0.012), - f2 = FT(0.25), -) - return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) -end - -function soil_niwotridge(; - vg_α = FT(0.04), - vg_n = FT(2.05), - ν = FT(0.45), # m3/m3 - ν_ss_quartz = FT(0.1), - ν_ss_om = FT(0.1), - K_sat = FT(4e-7), # m/s, matches Natan - S_s = FT(1e-3), # 1/m, guess - θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 - ν_ss_gravel = FT(0.0), - z_0m = FT(0.1), - z_0b = FT(0.1), - emissivity = FT(0.98), - PAR_albedo = FT(0.2), - NIR_albedo = FT(0.2), -) - hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) - return EnergyHydrologyParameters( - FT; - ν, - ν_ss_om, - ν_ss_quartz, - ν_ss_gravel, - hydrology_cm, - K_sat, - S_s, - θ_r, - PAR_albedo, - NIR_albedo, - emissivity, - z_0m, - z_0b, - ) -end - -function radiative_transfer_niwotridge(; - Ω = FT(0.69), - ld = ConstantGFunction(FT(0.5)), - α_PAR_leaf = FT(0.1), - λ_γ_PAR = FT(5e-7), - λ_γ_NIR = FT(1.65e-6), - τ_PAR_leaf = FT(0.05), - α_NIR_leaf = FT(0.35), - τ_NIR_leaf = FT(0.25), - n_layers = UInt64(20), - ϵ_canopy = FT(0.97), -) - return TwoStreamParameters( - G_Function = ld, - α_PAR_leaf = α_PAR_leaf, - τ_PAR_leaf = τ_PAR_leaf, - α_NIR_leaf = α_NIR_leaf, - τ_NIR_leaf = τ_NIR_leaf, - ϵ_canopy = ϵ_canopy, - Ω = Ω, - λ_γ_PAR = λ_γ_PAR, - λ_γ_NIR = λ_γ_NIR, - n_layers = n_layers, - ) -end - -function canopy_energy_balance_niwotridge(; ac_canopy = FT(3e3)) - return BigLeafEnergyParameters(ac_canopy) -end - -function conductance_niwotridge(; - g1 = FT(141), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. - Drel = FT(1.6), - g0 = FT(1e-4), -) - return MedlynConductanceParameters(Drel, g0, g1) -end - -function photosynthesis_niwotridge(; - mechanism = Canopy.C3(), - oi = FT(0.209), - ϕ = FT(0.6), - θj = FT(0.9), - f = FT(0.015), - sc = FT(2e-6), # Bonan's book: range of 2-5e-6 - pc = FT(-2e6), # Bonan's book: -2e6 - Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 - Γstar25 = FT(4.275e-5), - Kc25 = FT(4.049e-4), - Ko25 = FT(0.2874), - To = FT(298.15), - ΔHkc = FT(79430), - ΔHko = FT(36380), - ΔHVcmax = FT(58520), - ΔHΓstar = FT(37830), - ΔHJmax = FT(43540), - ΔHRd = FT(46390), -) - return FarquharParameters( - Vcmax25, - Γstar25, - Kc25, - Ko25, - ΔHkc, - ΔHko, - ΔHVcmax, - ΔHΓstar, - ΔHJmax, - ΔHRd, - To, - oi, - ϕ, - θj, - f, - sc, - pc, - mechanism, - ) -end - -function plant_hydraulics_niwotridge(; - SAI = FT(1.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? - f_root_to_shoot = FT(3.5), - 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 - capacity = FT(10), # kg/m^2 - plant_S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m - rooting_depth = FT(1.0), # from Natan -) - conductivity_model = - PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) - retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) - - return PlantHydraulicsP( - SAI, - f_root_to_shoot, - capacity, - plant_S_s, - rooting_depth, - conductivity_model, - retention_model, - ) -end +export niwotridge_default_params, + hetero_resp_niwotridge, + auto_resp_niwotridge, + soil_niwotridge, + radiative_transfer_niwotridge, + canopy_energy_balance_niwotridge, + conductance_niwotridge, + photosynthesis_niwotridge, + plant_hydraulics_niwotridge + +function niwotridge_default_params(; + hetero_resp = hetero_resp_niwotridge(), + auto_resp = auto_resp_niwotridge(), + soil = soil_niwotridge(), + radiative_transfer = radiative_transfer_niwotridge(), + canopy_energy_balance = canopy_energy_balance_niwotridge(), + conductance = conductance_niwotridge(), + photosynthesis = photosynthesis_niwotridge(), + plant_hydraulics = plant_hydraulics_niwotridge(), +) + return SiteParams( + hetero_resp, + auto_resp, + soil, + radiative_transfer, + canopy_energy_balance, + conductance, + photosynthesis, + plant_hydraulics, + ) +end + +function hetero_resp_niwotridge(; + D_ref = FT(1.39e-5), + D_liq = FT(3.17), + α_sx = FT(194e3), + Ea_sx = FT(61e3), + kM_sx = FT(5e-3), + kM_o2 = FT(0.004), + O2_a = FT(0.209), + D_oa = FT(1.67), + p_sx = FT(0.024), +) + return HeteroRespP( + D_ref, + D_liq, + α_sx, + Ea_sx, + kM_sx, + kM_o2, + O2_a, + D_oa, + p_sx, + ) +end + +function auto_resp_niwotridge(; + ne = FT(8 * 1e-4), + ηsl = FT(0.01), + σl = FT(0.05), + μr = FT(1.0), + μs = FT(0.1), + f1 = FT(0.012), + f2 = FT(0.25), +) + return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) +end + +function soil_niwotridge(; + vg_α = FT(0.04), + vg_n = FT(2.05), + ν = FT(0.45), # m3/m3 + ν_ss_quartz = FT(0.1), + ν_ss_om = FT(0.1), + K_sat = FT(4e-7), # m/s, matches Natan + S_s = FT(1e-3), # 1/m, guess + θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + ν_ss_gravel = FT(0.0), + z_0m = FT(0.1), + z_0b = FT(0.1), + emissivity = FT(0.98), + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), +) + hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) + return EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo, + NIR_albedo, + emissivity, + z_0m, + z_0b, + ) +end + +function radiative_transfer_niwotridge(; + Ω = FT(0.69), + ld = ConstantGFunction(FT(0.5)), + α_PAR_leaf = FT(0.1), + λ_γ_PAR = FT(5e-7), + λ_γ_NIR = FT(1.65e-6), + τ_PAR_leaf = FT(0.05), + α_NIR_leaf = FT(0.35), + τ_NIR_leaf = FT(0.25), + n_layers = UInt64(20), + ϵ_canopy = FT(0.97), +) + return TwoStreamParameters( + G_Function = ld, + α_PAR_leaf = α_PAR_leaf, + τ_PAR_leaf = τ_PAR_leaf, + α_NIR_leaf = α_NIR_leaf, + τ_NIR_leaf = τ_NIR_leaf, + ϵ_canopy = ϵ_canopy, + Ω = Ω, + λ_γ_PAR = λ_γ_PAR, + λ_γ_NIR = λ_γ_NIR, + n_layers = n_layers, + ) +end + +function canopy_energy_balance_niwotridge(; ac_canopy = FT(3e3)) + return BigLeafEnergyParameters(ac_canopy) +end + +function conductance_niwotridge(; + g1 = FT(141), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + Drel = FT(1.6), + g0 = FT(1e-4), +) + return MedlynConductanceParameters(Drel, g0, g1) +end + +function photosynthesis_niwotridge(; + mechanism = Canopy.C3(), + oi = FT(0.209), + ϕ = FT(0.6), + θj = FT(0.9), + f = FT(0.015), + sc = FT(2e-6), # Bonan's book: range of 2-5e-6 + pc = FT(-2e6), # Bonan's book: -2e6 + Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 + Γstar25 = FT(4.275e-5), + Kc25 = FT(4.049e-4), + Ko25 = FT(0.2874), + To = FT(298.15), + ΔHkc = FT(79430), + ΔHko = FT(36380), + ΔHVcmax = FT(58520), + ΔHΓstar = FT(37830), + ΔHJmax = FT(43540), + ΔHRd = FT(46390), +) + return FarquharParameters( + Vcmax25, + Γstar25, + Kc25, + Ko25, + ΔHkc, + ΔHko, + ΔHVcmax, + ΔHΓstar, + ΔHJmax, + ΔHRd, + To, + oi, + ϕ, + θj, + f, + sc, + pc, + mechanism, + ) +end + +function plant_hydraulics_niwotridge(; + SAI = FT(1.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? + f_root_to_shoot = FT(3.5), + 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 + capacity = FT(10), # kg/m^2 + plant_S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m + rooting_depth = FT(1.0), # from Natan +) + conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) + + return PlantHydraulicsP( + SAI, + f_root_to_shoot, + capacity, + plant_S_s, + rooting_depth, + conductivity_model, + retention_model, + ) +end diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl index 7228ccab2a..ac423ba8cf 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_sites/US-Var/US-Var_parameters.jl @@ -1,208 +1,208 @@ -export vairaranch_default_params, - hetero_resp_vairaranch, - auto_resp_vairaranch, - soil_vairaranch, - radiative_transfer_vairaranch, - canopy_energy_balance_vairaranch, - conductance_vairaranch, - photosynthesis_vairaranch, - plant_hydraulics_vairaranch - -function vairaranch_default_params(; - hetero_resp = hetero_resp_vairaranch(), - auto_resp = auto_resp_vairaranch(), - soil = soil_vairaranch(), - radiative_transfer = radiative_transfer_vairaranch(), - canopy_energy_balance = canopy_energy_balance_vairaranch(), - conductance = conductance_vairaranch(), - photosynthesis = photosynthesis_vairaranch(), - plant_hydraulics = plant_hydraulics_vairaranch(), -) - return SiteParams( - hetero_resp, - auto_resp, - soil, - radiative_transfer, - canopy_energy_balance, - conductance, - photosynthesis, - plant_hydraulics, - ) -end - -function hetero_resp_vairaranch(; - D_ref = FT(1.39e-5), - D_liq = FT(3.17), - α_sx = FT(194e3), - Ea_sx = FT(61e3), - kM_sx = FT(5e-3), - kM_o2 = FT(0.004), - O2_a = FT(0.209), - D_oa = FT(1.67), - p_sx = FT(0.024), -) - return HeteroRespP( - D_ref, - D_liq, - α_sx, - Ea_sx, - kM_sx, - kM_o2, - O2_a, - D_oa, - p_sx, - ) -end - -function auto_resp_vairaranch(; - ne = FT(8 * 1e-4), - ηsl = FT(0.01), - σl = FT(0.05), - μr = FT(1.0), - μs = FT(0.1), - f1 = FT(0.012), - f2 = FT(0.25), -) - return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) -end - -function soil_vairaranch(; - vg_n = FT(2.00), # unitless - vg_α = FT(2.00), # inverse meters - ν = FT(0.45), # m3/m3 - ν_ss_quartz = FT(0.1), - ν_ss_om = FT(0.1), - K_sat = FT(4e-7), # m/s, matches Natan - S_s = FT(1e-3), # 1/m, guess - θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 - ν_ss_gravel = FT(0.0), - z_0m = FT(0.01), - z_0b = FT(0.001), - emissivity = FT(0.98), - PAR_albedo = FT(0.2), - NIR_albedo = FT(0.2), -) - hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) - return EnergyHydrologyParameters( - FT; - ν, - ν_ss_om, - ν_ss_quartz, - ν_ss_gravel, - hydrology_cm, - K_sat, - S_s, - θ_r, - PAR_albedo, - NIR_albedo, - emissivity, - z_0m, - z_0b, - ) -end - -function radiative_transfer_vairaranch(; - Ω = FT(0.69), - ld = ConstantGFunction(FT(0.5)), - α_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), - n_layers = UInt64(20), - ϵ_canopy = FT(0.97), -) - return TwoStreamParameters( - G_Function = ld, - α_PAR_leaf = α_PAR_leaf, - τ_PAR_leaf = τ_PAR_leaf, - α_NIR_leaf = α_NIR_leaf, - τ_NIR_leaf = τ_NIR_leaf, - ϵ_canopy = ϵ_canopy, - Ω = Ω, - λ_γ_PAR = λ_γ_PAR, - λ_γ_NIR = λ_γ_NIR, - n_layers = n_layers, - ) -end - -function canopy_energy_balance_vairaranch(; ac_canopy = FT(745)) - return BigLeafEnergyParameters(ac_canopy) -end - -function conductance_vairaranch(; - g1 = FT(166), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. - Drel = FT(1.6), - g0 = FT(1e-4), -) - return MedlynConductanceParameters(Drel, g0, g1) -end - -function photosynthesis_vairaranch(; - mechanism = Canopy.C3(), - oi = FT(0.209), - ϕ = FT(0.6), - θj = FT(0.9), - f = FT(0.015), - sc = FT(2e-6), # Bonan's book: range of 2-5e-6 - pc = FT(-2e6), # Bonan's book: -2e6 - Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 - Γstar25 = FT(4.275e-5), - Kc25 = FT(4.049e-4), - Ko25 = FT(0.2874), - To = FT(298.15), - ΔHkc = FT(79430), - ΔHko = FT(36380), - ΔHVcmax = FT(58520), - ΔHΓstar = FT(37830), - ΔHJmax = FT(43540), - ΔHRd = FT(46390), -) - return FarquharParameters( - Vcmax25, - Γstar25, - Kc25, - Ko25, - ΔHkc, - ΔHko, - ΔHVcmax, - ΔHΓstar, - ΔHJmax, - ΔHRd, - To, - oi, - ϕ, - θj, - f, - sc, - pc, - mechanism, - ) -end - -function plant_hydraulics_vairaranch(; - SAI = FT(0.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? - f_root_to_shoot = FT(1.0), - K_sat_plant = 5e-9, # m/s # seems much too small? - ψ63 = FT(-2.7 / 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 - capacity = FT(2.0), # kg/m^2 - plant_S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m - rooting_depth = FT(2.6), # from Natan -) - conductivity_model = - PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) - retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) - - return PlantHydraulicsP( - SAI, - f_root_to_shoot, - capacity, - plant_S_s, - rooting_depth, - conductivity_model, - retention_model, - ) -end +export vairaranch_default_params, + hetero_resp_vairaranch, + auto_resp_vairaranch, + soil_vairaranch, + radiative_transfer_vairaranch, + canopy_energy_balance_vairaranch, + conductance_vairaranch, + photosynthesis_vairaranch, + plant_hydraulics_vairaranch + +function vairaranch_default_params(; + hetero_resp = hetero_resp_vairaranch(), + auto_resp = auto_resp_vairaranch(), + soil = soil_vairaranch(), + radiative_transfer = radiative_transfer_vairaranch(), + canopy_energy_balance = canopy_energy_balance_vairaranch(), + conductance = conductance_vairaranch(), + photosynthesis = photosynthesis_vairaranch(), + plant_hydraulics = plant_hydraulics_vairaranch(), +) + return SiteParams( + hetero_resp, + auto_resp, + soil, + radiative_transfer, + canopy_energy_balance, + conductance, + photosynthesis, + plant_hydraulics, + ) +end + +function hetero_resp_vairaranch(; + D_ref = FT(1.39e-5), + D_liq = FT(3.17), + α_sx = FT(194e3), + Ea_sx = FT(61e3), + kM_sx = FT(5e-3), + kM_o2 = FT(0.004), + O2_a = FT(0.209), + D_oa = FT(1.67), + p_sx = FT(0.024), +) + return HeteroRespP( + D_ref, + D_liq, + α_sx, + Ea_sx, + kM_sx, + kM_o2, + O2_a, + D_oa, + p_sx, + ) +end + +function auto_resp_vairaranch(; + ne = FT(8 * 1e-4), + ηsl = FT(0.01), + σl = FT(0.05), + μr = FT(1.0), + μs = FT(0.1), + f1 = FT(0.012), + f2 = FT(0.25), +) + return AutotrophicRespirationParameters(ne, ηsl, σl, μr, μs, f1, f2) +end + +function soil_vairaranch(; + vg_n = FT(2.00), # unitless + vg_α = FT(2.00), # inverse meters + ν = FT(0.45), # m3/m3 + ν_ss_quartz = FT(0.1), + ν_ss_om = FT(0.1), + K_sat = FT(4e-7), # m/s, matches Natan + S_s = FT(1e-3), # 1/m, guess + θ_r = FT(0.067), # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + ν_ss_gravel = FT(0.0), + z_0m = FT(0.01), + z_0b = FT(0.001), + emissivity = FT(0.98), + PAR_albedo = FT(0.2), + NIR_albedo = FT(0.2), +) + hydrology_cm = vanGenuchten{FT}(; α = vg_α, n = vg_n) + return EnergyHydrologyParameters( + FT; + ν, + ν_ss_om, + ν_ss_quartz, + ν_ss_gravel, + hydrology_cm, + K_sat, + S_s, + θ_r, + PAR_albedo, + NIR_albedo, + emissivity, + z_0m, + z_0b, + ) +end + +function radiative_transfer_vairaranch(; + Ω = FT(0.69), + ld = ConstantGFunction(FT(0.5)), + α_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), + n_layers = UInt64(20), + ϵ_canopy = FT(0.97), +) + return TwoStreamParameters( + G_Function = ld, + α_PAR_leaf = α_PAR_leaf, + τ_PAR_leaf = τ_PAR_leaf, + α_NIR_leaf = α_NIR_leaf, + τ_NIR_leaf = τ_NIR_leaf, + ϵ_canopy = ϵ_canopy, + Ω = Ω, + λ_γ_PAR = λ_γ_PAR, + λ_γ_NIR = λ_γ_NIR, + n_layers = n_layers, + ) +end + +function canopy_energy_balance_vairaranch(; ac_canopy = FT(745)) + return BigLeafEnergyParameters(ac_canopy) +end + +function conductance_vairaranch(; + g1 = FT(166), # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. + Drel = FT(1.6), + g0 = FT(1e-4), +) + return MedlynConductanceParameters(Drel, g0, g1) +end + +function photosynthesis_vairaranch(; + mechanism = Canopy.C3(), + oi = FT(0.209), + ϕ = FT(0.6), + θj = FT(0.9), + f = FT(0.015), + sc = FT(2e-6), # Bonan's book: range of 2-5e-6 + pc = FT(-2e6), # Bonan's book: -2e6 + Vcmax25 = FT(9e-5), # from Yujie's paper 4.5e-5 , Natan used 9e-5 + Γstar25 = FT(4.275e-5), + Kc25 = FT(4.049e-4), + Ko25 = FT(0.2874), + To = FT(298.15), + ΔHkc = FT(79430), + ΔHko = FT(36380), + ΔHVcmax = FT(58520), + ΔHΓstar = FT(37830), + ΔHJmax = FT(43540), + ΔHRd = FT(46390), +) + return FarquharParameters( + Vcmax25, + Γstar25, + Kc25, + Ko25, + ΔHkc, + ΔHko, + ΔHVcmax, + ΔHΓstar, + ΔHJmax, + ΔHRd, + To, + oi, + ϕ, + θj, + f, + sc, + pc, + mechanism, + ) +end + +function plant_hydraulics_vairaranch(; + SAI = FT(0.0), # m2/m2 or: estimated from Wang et al, FT(0.00242) ? + f_root_to_shoot = FT(1.0), + K_sat_plant = 5e-9, # m/s # seems much too small? + ψ63 = FT(-2.7 / 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 + capacity = FT(2.0), # kg/m^2 + plant_S_s = FT(1e-2 * 0.0098), # m3/m3/MPa to m3/m3/m + rooting_depth = FT(2.6), # from Natan +) + conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) + + return PlantHydraulicsP( + SAI, + f_root_to_shoot, + capacity, + plant_S_s, + rooting_depth, + conductivity_model, + retention_model, + ) +end diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl index 97907d5fa3..4cafbd28b2 100644 --- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl +++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl @@ -1,329 +1,326 @@ -export run_fluxnet - -""" - fluxnet_simulation(site_ID; kwargs) - -Run ClimaLand at a fluxnet site. -""" -function run_fluxnet( - site_ID; - FT = Float64, - context = ClimaComms.SingletonCommsContext(), - setup = make_setup(site_ID), - domain = make_domain(setup, FT), - config = make_config(site_ID), - params = make_parameters(site_ID), - drivers = make_drivers(site_ID, setup, config, params, context), - timestepper = make_timestepper(setup), -) - #earth_param_set = create_lsm_parameters(FT) - - # Now we set up the model. For the soil model, we pick - # a model type and model args: - soil_domain = domain.land_domain - soil_ps = Soil.EnergyHydrologyParameters( - FT; - ν = params.soil.ν, - ν_ss_om = params.soil.ν_ss_om, - ν_ss_quartz = params.soil.ν_ss_quartz, - ν_ss_gravel = params.soil.ν_ss_gravel, - hydrology_cm = params.soil.hydrology_cm, - K_sat = params.soil.K_sat, - S_s = params.soil.S_s, - θ_r = params.soil.θ_r, - z_0m = params.soil.z_0m, - z_0b = params.soil.z_0b, - emissivity = params.soil.emissivity, - PAR_albedo = params.soil.PAR_albedo, - NIR_albedo = params.soil.NIR_albedo, - ) - - soil_args = (domain = soil_domain, parameters = soil_ps) - soil_model_type = Soil.EnergyHydrology{FT} - - # Soil microbes model - soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT} - - soilco2_ps = SoilCO2ModelParameters( - FT; - D_ref = params.hetero_resp.D_ref, - D_liq = params.hetero_resp.D_liq, - # DAMM - α_sx = params.hetero_resp.α_sx, - Ea_sx = params.hetero_resp.Ea_sx, - kM_sx = params.hetero_resp.kM_sx, - kM_o2 = params.hetero_resp.kM_o2, - O2_a = params.hetero_resp.O2_a, - D_oa = params.hetero_resp.D_oa, - p_sx = params.hetero_resp.p_sx, - earth_param_set = earth_param_set, - ) - - # soil microbes args - Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5)) - - soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC() - soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux - soilco2_sources = (MicrobeProduction{FT}(),) - - soilco2_boundary_conditions = - (; top = soilco2_top_bc, bottom = CO2 = soilco2_bot_bc) - - soilco2_args = (; - boundary_conditions = soilco2_boundary_conditions, - sources = soilco2_sources, - domain = soil_domain, - parameters = soilco2_ps, - ) - - # Now we set up the canopy model, which we set up by component: - # Component Types - canopy_component_types = (; - autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, - radiative_transfer = Canopy.TwoStreamModel{FT}, - photosynthesis = Canopy.FarquharModel{FT}, - conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, - energy = Canopy.BigLeafEnergyModel{FT}, - ) - # Individual Component arguments - # Set up autotrophic respiration - autotrophic_respiration_args = (; - parameters = AutotrophicRespirationParameters{FT}(; - ne = params.auto_resp.ne, - ηsl = params.auto_resp.ηsl, - σl = params.auto_resp.σl, - μr = params.auto_resp.μr, - μs = params.auto_resp.μs, - f1 = params.auto_resp.f1, - f2 = params.auto_resp.f2, - ) - ) - # Set up radiative transfer - radiative_transfer_args = (; - parameters = TwoStreamParameters( - FT; - Ω = params.radiative_transfer.Ω, - G_Function = params.radiative_transfer.G_Function, - α_PAR_leaf = params.radiative_transfer.α_PAR_leaf, - λ_γ_PAR = params.radiative_transfer.λ_γ_PAR, - λ_γ_NIR = params.radiative_transfer.λ_γ_NIR, - τ_PAR_leaf = params.radiative_transfer.τ_PAR_leaf, - α_NIR_leaf = params.radiative_transfer.α_NIR_leaf, - τ_NIR_leaf = params.radiative_transfer.τ_NIR_leaf, - n_layers = params.radiative_transfer.n_layers, - ϵ_canopy = params.radiative_transfer.ϵ_canopy, - ) - ) - # Set up conductance - conductance_args = (; - parameters = MedlynConductanceParameters( - FT; - g1 = params.conductance.g1, - Drel = params.conductance.Drel, - g0 = params.conductance.g0, - ) - ) - # Set up photosynthesis - photosynthesis_args = (; - parameters = FarquharParameters( - FT, - Canopy.C3(); - oi = params.photosynthesis.oi, - ϕ = params.photosynthesis.ϕ, - θj = params.photosynthesis.θj, - f = params.photosynthesis.f, - sc = params.photosynthesis.sc, - pc = params.photosynthesis.pc, - Vcmax25 = params.photosynthesis.Vcmax25, - Γstar25 = params.photosynthesis.Γstar25, - Kc25 = params.photosynthesis.Kc25, - Ko25 = params.photosynthesis.Ko25, - To = params.photosynthesis.To, - ΔHkc = params.photosynthesis.ΔHkc, - ΔHko = params.photosynthesis.ΔHko, - ΔHVcmax = params.photosynthesis.ΔHVcmax, - ΔHΓstar = params.photosynthesis.ΔHΓstar, - ΔHJmax = params.photosynthesis.ΔHJmax, - ΔHRd = params.photosynthesis.ΔHRd, - ) - ) - # Set up plant hydraulics - ai_parameterization = PrescribedSiteAreaIndex{FT}( - drivers.LAIfunction, - params.plant_hydraulics.SAI, - drivers.RAI, - ) - - function root_distribution( - z::T; - rooting_depth = params.plant_hydraulics.rooting_depth, - ) where {T} - return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m - end - - plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; - ai_parameterization = ai_parameterization, - ν = drivers.plant_ν, - S_s = params.plant_hydraulics.S_s, - root_distribution = root_distribution, - conductivity_model = params.plant_hydraulics.conductivity_model, - retention_model = params.plant_hydraulics.retention_model, - ) - plant_hydraulics_args = ( - parameters = plant_hydraulics_ps, - n_stem = setup.n_stem, - n_leaf = setup.n_leaf, - compartment_midpoints = domain.compartment_midpoints, - compartment_surfaces = domain.compartment_surfaces, - ) - - energy_args = ( - parameters = Canopy.BigLeafEnergyParameters{FT}( - params.canopy_energy_balance.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, - ) - - z0_m = FT(0.13) * domain.h_canopy - z0_b = FT(0.1) * z0_m - - # Other info needed - shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( - z0_m, - z0_b, - earth_param_set, - ) - - canopy_model_args = - (; parameters = shared_params, domain = domain.canopy_domain) - - # Integrated plant hydraulics and soil model - land_input = ( - atmos = drivers.atmos, - radiation = drivers.radiation, - soil_organic_carbon = Csom, - ) - 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) - exp_tendency! = make_exp_tendency(land) - - #Initial conditions - Y.soil.ϑ_l = - drivers.drivers.SWC.status != absent ? - drivers.drivers.SWC.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] : - params.soil.soil_ν / 2 # Get soil water content at t0 - # Both data and simulation are reference to 2005-01-01-00 (LOCAL) - # or 2005-01-01-06 (UTC) - Y.soil.θ_i = FT(0.0) - T_0 = - drivers.drivers.TS.status != absent ? - drivers.drivers.TS.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] : - drivers.drivers.TA.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] + - 40# Get soil temperature at t0 - ρc_s = - volumetric_heat_capacity.( - Y.soil.ϑ_l, - Y.soil.θ_i, - 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, - 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] - ψ_leaf_0 = FT(-2e5 / 9800) - ψ_comps = setup.n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 - - S_l_ini = - inverse_water_retention_curve.( - params.plant_hydraulics.retention_model, - ψ_comps, - drivers.plant_ν, - params.plant_hydraulics.S_s, - ) - - for i in 1:(setup.n_stem + setup.n_leaf) - Y.canopy.hydraulics.ϑ_l.:($i) .= - augmented_liquid_fraction.(drivers.plant_ν, S_l_ini[i]) - end - - Y.canopy.energy.T = - drivers.drivers.TA.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] # Get atmos temperature at t0 - - set_initial_cache! = make_set_initial_cache(land) - set_initial_cache!(p, Y, setup.t0) - - # Simulation - sv = (; - t = Array{Float64}(undef, length(timestepper.saveat)), - saveval = Array{NamedTuple}(undef, length(timestepper.saveat)), - ) - saving_cb = ClimaLand.NonInterpSavingCallback(sv, timestepper.saveat) - ## How often we want to update the drivers. Note that this uses the defined `t0` and `tf` - ## defined in the simulatons file - updateat = Array((setup.t0):(drivers.DATA_DT):(timestepper.tf)) - land_drivers = ClimaLand.get_drivers(land) - updatefunc = ClimaLand.make_update_drivers(land_drivers) - driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - cb = SciMLBase.CallbackSet(driver_cb, saving_cb) - - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction((T_exp!) = exp_tendency!), - Y, - (setup.t0, timestepper.tf), - p, - ) - sol = SciMLBase.solve( - prob, - timestepper.ode_algo; - dt = setup.dt, - callback = cb, - adaptive = false, - saveat = timestepper.saveat, - ) - - if isfile( - joinpath( - climaland_dir, - "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", - ), - ) - rm( - joinpath( - climaland_dir, - "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", - ), - ) - else - nothing - end - - return sv, sol, Y, p -end - - - +export run_fluxnet + +""" + fluxnet_simulation(site_ID; kwargs) + +Run ClimaLand at a fluxnet site. +""" +function run_fluxnet( + site_ID; + FT = Float64, + context = ClimaComms.SingletonCommsContext(), + setup = make_setup(site_ID), + domain = make_domain(setup, FT), + config = make_config(site_ID), + params = make_parameters(site_ID), + drivers = make_drivers(site_ID, setup, config, params, context), + timestepper = make_timestepper(setup), +) + #earth_param_set = create_lsm_parameters(FT) + + # Now we set up the model. For the soil model, we pick + # a model type and model args: + soil_domain = domain.land_domain + soil_ps = Soil.EnergyHydrologyParameters( + FT; + ν = params.soil.ν, + ν_ss_om = params.soil.ν_ss_om, + ν_ss_quartz = params.soil.ν_ss_quartz, + ν_ss_gravel = params.soil.ν_ss_gravel, + hydrology_cm = params.soil.hydrology_cm, + K_sat = params.soil.K_sat, + S_s = params.soil.S_s, + θ_r = params.soil.θ_r, + z_0m = params.soil.z_0m, + z_0b = params.soil.z_0b, + emissivity = params.soil.emissivity, + PAR_albedo = params.soil.PAR_albedo, + NIR_albedo = params.soil.NIR_albedo, + ) + + soil_args = (domain = soil_domain, parameters = soil_ps) + soil_model_type = Soil.EnergyHydrology{FT} + + # Soil microbes model + soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT} + + soilco2_ps = SoilCO2ModelParameters( + FT; + D_ref = params.hetero_resp.D_ref, + D_liq = params.hetero_resp.D_liq, + # DAMM + α_sx = params.hetero_resp.α_sx, + Ea_sx = params.hetero_resp.Ea_sx, + kM_sx = params.hetero_resp.kM_sx, + kM_o2 = params.hetero_resp.kM_o2, + O2_a = params.hetero_resp.O2_a, + D_oa = params.hetero_resp.D_oa, + p_sx = params.hetero_resp.p_sx, + earth_param_set = earth_param_set, + ) + + # soil microbes args + Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(TimeVaryingInput((t) -> 5)) + + soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC() + soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux + soilco2_sources = (MicrobeProduction{FT}(),) + + soilco2_boundary_conditions = + (; top = soilco2_top_bc, bottom = CO2 = soilco2_bot_bc) + + soilco2_args = (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = soil_domain, + parameters = soilco2_ps, + ) + + # Now we set up the canopy model, which we set up by component: + # Component Types + canopy_component_types = (; + autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = Canopy.TwoStreamModel{FT}, + photosynthesis = Canopy.FarquharModel{FT}, + conductance = Canopy.MedlynConductanceModel{FT}, + hydraulics = Canopy.PlantHydraulicsModel{FT}, + energy = Canopy.BigLeafEnergyModel{FT}, + ) + # Individual Component arguments + # Set up autotrophic respiration + autotrophic_respiration_args = (; + parameters = AutotrophicRespirationParameters{FT}(; + ne = params.auto_resp.ne, + ηsl = params.auto_resp.ηsl, + σl = params.auto_resp.σl, + μr = params.auto_resp.μr, + μs = params.auto_resp.μs, + f1 = params.auto_resp.f1, + f2 = params.auto_resp.f2, + ) + ) + # Set up radiative transfer + radiative_transfer_args = (; + parameters = TwoStreamParameters( + FT; + Ω = params.radiative_transfer.Ω, + G_Function = params.radiative_transfer.G_Function, + α_PAR_leaf = params.radiative_transfer.α_PAR_leaf, + λ_γ_PAR = params.radiative_transfer.λ_γ_PAR, + λ_γ_NIR = params.radiative_transfer.λ_γ_NIR, + τ_PAR_leaf = params.radiative_transfer.τ_PAR_leaf, + α_NIR_leaf = params.radiative_transfer.α_NIR_leaf, + τ_NIR_leaf = params.radiative_transfer.τ_NIR_leaf, + n_layers = params.radiative_transfer.n_layers, + ϵ_canopy = params.radiative_transfer.ϵ_canopy, + ) + ) + # Set up conductance + conductance_args = (; + parameters = MedlynConductanceParameters( + FT; + g1 = params.conductance.g1, + Drel = params.conductance.Drel, + g0 = params.conductance.g0, + ) + ) + # Set up photosynthesis + photosynthesis_args = (; + parameters = FarquharParameters( + FT, + Canopy.C3(); + oi = params.photosynthesis.oi, + ϕ = params.photosynthesis.ϕ, + θj = params.photosynthesis.θj, + f = params.photosynthesis.f, + sc = params.photosynthesis.sc, + pc = params.photosynthesis.pc, + Vcmax25 = params.photosynthesis.Vcmax25, + Γstar25 = params.photosynthesis.Γstar25, + Kc25 = params.photosynthesis.Kc25, + Ko25 = params.photosynthesis.Ko25, + To = params.photosynthesis.To, + ΔHkc = params.photosynthesis.ΔHkc, + ΔHko = params.photosynthesis.ΔHko, + ΔHVcmax = params.photosynthesis.ΔHVcmax, + ΔHΓstar = params.photosynthesis.ΔHΓstar, + ΔHJmax = params.photosynthesis.ΔHJmax, + ΔHRd = params.photosynthesis.ΔHRd, + ) + ) + # Set up plant hydraulics + ai_parameterization = PrescribedSiteAreaIndex{FT}( + drivers.LAIfunction, + params.plant_hydraulics.SAI, + drivers.RAI, + ) + + function root_distribution( + z::T; + rooting_depth = params.plant_hydraulics.rooting_depth, + ) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m + end + + plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = drivers.plant_ν, + S_s = params.plant_hydraulics.S_s, + root_distribution = root_distribution, + conductivity_model = params.plant_hydraulics.conductivity_model, + retention_model = params.plant_hydraulics.retention_model, + ) + plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = setup.n_stem, + n_leaf = setup.n_leaf, + compartment_midpoints = domain.compartment_midpoints, + compartment_surfaces = domain.compartment_surfaces, + ) + + energy_args = ( + parameters = Canopy.BigLeafEnergyParameters{FT}( + params.canopy_energy_balance.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, + ) + + z0_m = FT(0.13) * domain.h_canopy + z0_b = FT(0.1) * z0_m + + # Other info needed + shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, + ) + + canopy_model_args = + (; parameters = shared_params, domain = domain.canopy_domain) + + # Integrated plant hydraulics and soil model + land_input = ( + atmos = drivers.atmos, + radiation = drivers.radiation, + soil_organic_carbon = Csom, + ) + 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) + exp_tendency! = make_exp_tendency(land) + + #Initial conditions + Y.soil.ϑ_l = + drivers.drivers.SWC.status != absent ? + drivers.drivers.SWC.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] : + params.soil.soil_ν / 2 # Get soil water content at t0 + # Both data and simulation are reference to 2005-01-01-00 (LOCAL) + # or 2005-01-01-06 (UTC) + Y.soil.θ_i = FT(0.0) + T_0 = + drivers.drivers.TS.status != absent ? + drivers.drivers.TS.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] : + drivers.drivers.TA.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] + + 40# Get soil temperature at t0 + ρc_s = + volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + 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, + 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] + ψ_leaf_0 = FT(-2e5 / 9800) + ψ_comps = setup.n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 + + S_l_ini = + inverse_water_retention_curve.( + params.plant_hydraulics.retention_model, + ψ_comps, + drivers.plant_ν, + params.plant_hydraulics.S_s, + ) + + for i in 1:(setup.n_stem + setup.n_leaf) + Y.canopy.hydraulics.ϑ_l.:($i) .= + augmented_liquid_fraction.(drivers.plant_ν, S_l_ini[i]) + end + + Y.canopy.energy.T = + drivers.drivers.TA.values[1 + Int(round(setup.t0 / drivers.DATA_DT))] # Get atmos temperature at t0 + + set_initial_cache! = make_set_initial_cache(land) + set_initial_cache!(p, Y, setup.t0) + + # Simulation + sv = (; + t = Array{Float64}(undef, length(timestepper.saveat)), + saveval = Array{NamedTuple}(undef, length(timestepper.saveat)), + ) + saving_cb = ClimaLand.NonInterpSavingCallback(sv, timestepper.saveat) + ## How often we want to update the drivers. Note that this uses the defined `t0` and `tf` + ## defined in the simulatons file + updateat = Array((setup.t0):(drivers.DATA_DT):(timestepper.tf)) + land_drivers = ClimaLand.get_drivers(land) + updatefunc = ClimaLand.make_update_drivers(land_drivers) + driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) + cb = SciMLBase.CallbackSet(driver_cb, saving_cb) + + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction((T_exp!) = exp_tendency!), + Y, + (setup.t0, timestepper.tf), + p, + ) + sol = SciMLBase.solve( + prob, + timestepper.ode_algo; + dt = setup.dt, + callback = cb, + adaptive = false, + saveat = timestepper.saveat, + ) + + if isfile( + joinpath( + climaland_dir, + "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", + ), + ) + rm( + joinpath( + climaland_dir, + "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", + ), + ) + else + nothing + end + + return sv, sol, Y, p +end