From c76c154f3b2110139064aa10b432090d657988aa Mon Sep 17 00:00:00 2001 From: Alexis Renchon Date: Fri, 15 Dec 2023 08:35:02 -0800 Subject: [PATCH] first commit WIP WIP WIP buildkite test fix? buildkite error fix? mkdir figures WIP WIP fix attempt WIP WIP WIP fix rename fixes fixes added diurnals WIP small tweaks add timeseries swc, moisture stress, stomatal conductance WIP WIP tweaks format unitful SI or commonly used WIP WIP WIP WIP WIP WIP WIP WIP WIP WIP WIP WIP WIP up WIP WIP WIP precompile works rename WIP WIP WIP WIP WIP WIP WIP lib lib WIP WIP WIP WIP WIP WIP fluxnet_simulation works bug fix for Soil.k_dry Ozark only README update README tweaks path to configs as kwarg inputs and outputs df WIP WIP cleaning WIP update .buildkite/pipeline.yml climaformat remove ClimaLandDashboards from this PR compat .buildkite/pipeline.yml update export in files, not module update module format more compat entries dev in lib env instantiate buildkite fix attempt buildkite fix attempt 2 buildkite fix attempt 3 Buildkite: remove Manifest.toml buildkite fix attempt 4 add InverseFunc dep use Makie; loosen compats add MutableArith dep fix buildkite finger crossed pinpoint buildkite fix 1 pinpoint buildkite fix 2 pinpoint buildkite fix pinpoint buildkite fix 4 FLUXNET in fluxnet folder rename some files update data_tools.jl climaformat WIP WIP format --- .buildkite/pipeline.yml | 10 + lib/ClimaLandSimulations/LICENSE | 202 ++++++++ lib/ClimaLandSimulations/Project.toml | 66 +++ lib/ClimaLandSimulations/README.md | 16 + .../experiments/README.md | 38 ++ lib/ClimaLandSimulations/experiments/ozark.jl | 6 + .../src/ClimaLandSimulations.jl | 57 +++ .../src/fluxnet/US-Ha1/US-Ha1_parameters.jl | 125 +++++ .../src/fluxnet/US-Ha1/US-Ha1_simulation.jl | 26 ++ .../src/fluxnet/US-MOz/US-MOz_parameters.jl | 123 +++++ .../src/fluxnet/US-MOz/US-MOz_simulation.jl | 26 ++ .../src/fluxnet/US-NR1/US-NR1_parameters.jl | 128 +++++ .../src/fluxnet/US-NR1/US-NR1_simulation.jl | 26 ++ .../src/fluxnet/US-Var/US-Var_parameters.jl | 132 ++++++ .../src/fluxnet/US-Var/US-Var_simulation.jl | 26 ++ .../src/fluxnet_simulation.jl | 353 ++++++++++++++ .../utilities/climalsm_output_dataframe.jl | 65 +++ .../src/utilities/data_tools.jl | 174 +++++++ .../src/utilities/domain_setup.jl | 14 + .../src/utilities/inputs_dataframe.jl | 112 +++++ .../src/utilities/makie_plots.jl | 439 ++++++++++++++++++ .../src/utilities/met_drivers_FLUXNET.jl | 224 +++++++++ .../src/utilities/pull_MODIS.jl | 57 +++ .../src/utilities/timestepper_setup.jl | 7 + lib/ClimaLandSimulations/test/runtests.jl | 3 + lib/README.md | 4 + 26 files changed, 2459 insertions(+) create mode 100644 lib/ClimaLandSimulations/LICENSE create mode 100644 lib/ClimaLandSimulations/Project.toml create mode 100644 lib/ClimaLandSimulations/README.md create mode 100644 lib/ClimaLandSimulations/experiments/README.md create mode 100644 lib/ClimaLandSimulations/experiments/ozark.jl create mode 100644 lib/ClimaLandSimulations/src/ClimaLandSimulations.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_parameters.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_simulation.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_parameters.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_simulation.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_parameters.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_simulation.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_parameters.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_simulation.jl create mode 100644 lib/ClimaLandSimulations/src/fluxnet_simulation.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/climalsm_output_dataframe.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/data_tools.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/domain_setup.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/inputs_dataframe.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/makie_plots.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/met_drivers_FLUXNET.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/pull_MODIS.jl create mode 100644 lib/ClimaLandSimulations/src/utilities/timestepper_setup.jl create mode 100644 lib/ClimaLandSimulations/test/runtests.jl create mode 100644 lib/README.md diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 5307b76323..ec44e70b60 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -26,6 +26,10 @@ steps: - "julia --project=experiments -e 'using Pkg; Pkg.instantiate(;verbose=true)'" - "julia --project=experiments -e 'using Pkg; Pkg.status()'" + - echo "--- Instantiate lib/ClimaLandSimulations" + - "julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.develop(;path=\".\"); Pkg.instantiate(;verbose=true)'" + - "julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.status()'" + - echo "--- Instantiate test" - "julia --project=test -e 'using Pkg; Pkg.develop(;path=\".\"); Pkg.instantiate(;verbose=true)'" - "julia --project=test -e 'using Pkg; Pkg.status()'" @@ -74,6 +78,12 @@ steps: command: "julia --color=yes --project=experiments experiments/standalone/Soil/water_conservation.jl" artifact_paths: "experiments/standalone/Soil/water_conservation*png" + - group: "lib/ClimaLandSimulations" + steps: + - label: "Ozark figures Makie" + command: "julia --color=yes --project=lib/ClimaLandSimulations lib/ClimaLandSimulations/experiments/ozark.jl" + artifact_paths: "figures/*pdf" + - group: "GPU" steps: - label: "GPU runtests" diff --git a/lib/ClimaLandSimulations/LICENSE b/lib/ClimaLandSimulations/LICENSE new file mode 100644 index 0000000000..02e503f989 --- /dev/null +++ b/lib/ClimaLandSimulations/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright 2022 California Institute of Technology + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. \ No newline at end of file diff --git a/lib/ClimaLandSimulations/Project.toml b/lib/ClimaLandSimulations/Project.toml new file mode 100644 index 0000000000..ab3c228b9a --- /dev/null +++ b/lib/ClimaLandSimulations/Project.toml @@ -0,0 +1,66 @@ +name = "ClimaLandSimulations" +uuid = "348a0bd3-1299-4261-8002-d2cd97df6055" +authors = ["Clima Land Team"] +version = "0.1.0" + +[deps] +ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a" +CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" +ClimaLSM = "7884a58f-fab6-4fd0-82bb-ecfedb2d8430" +ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" +Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0" +HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" +Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" +RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f" +Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +UnitfulMoles = "999f2bd7-36bf-5ba7-9bc1-c9473aa75374" +cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" + +[compat] +ArtifactWrappers = "0.2" +CLIMAParameters = "0.7" +CairoMakie = "0.11" +ClimaCore = "0.11.7" +ClimaTimeSteppers = "0.7" +DataFrames = "1" +Dates = "1" +DelimitedFiles = "1" +Dierckx = "0.5" +Formatting = "0.4" +HTTP = "1" +Insolation = "0.8" +Interpolations = "0.15" +JSON = "0.21" +LaTeXStrings = "1" +NLsolve = "4" +PlotUtils = "1" +RootSolvers = "0.4" +SciMLBase = "2" +StaticArrays = "1" +Statistics = "1" +StatsBase = "0.34" +SurfaceFluxes = "0.8" +Thermodynamics = "0.11" +Unitful = "1" +UnitfulMoles = "0.1" +cuDNN = "1" +julia = "1.10" diff --git a/lib/ClimaLandSimulations/README.md b/lib/ClimaLandSimulations/README.md new file mode 100644 index 0000000000..d38b62edba --- /dev/null +++ b/lib/ClimaLandSimulations/README.md @@ -0,0 +1,16 @@ +# ClimaLandSimulations + +A library of methods to run ClimaLSM: +- at single sites, such as FLUXNET +- globally (WIP) + +## Installation + +```jl +julia> ] # to go into pkg mode +pkg> add ClimaLandSimulations +pkg> *backspace* # press backspace button to go back into julia mode +julia> using ClimaLandSimulations +``` + +For examples, see the [experiments](https://github.com/CliMA/ClimaLSM.jl/lib/ClimaLandSimulations/experiments) folder diff --git a/lib/ClimaLandSimulations/experiments/README.md b/lib/ClimaLandSimulations/experiments/README.md new file mode 100644 index 0000000000..4d294222cf --- /dev/null +++ b/lib/ClimaLandSimulations/experiments/README.md @@ -0,0 +1,38 @@ +# Running ClimaLand at FLUXNET sites + +To run a FLUXNET site, we just need the site ID, for example to run the Ozark site: + +```jl +julia> fluxnet_simulation("US-MOz") +``` + +This will run ClimaLand with default parameters and domain. Optionally, you can modify these: + +```jl +julia> fluxnet_simulation("US-MOz"; customparams = path_to_parameters, customdomain = path_to_domain) +``` + +where `customparams` and `customdomain` are modified from templates. Each site FLUXNET site has a default template. + +To get a list of ID of available FLUXNET sites, use: + +```jl +julia> list_fluxnet_ID() +``` + +To copy the parameters or domain template of a fluxnet site, use: + +```jl +julia> template_parameters(site_ID) +julia> template_domain(site_ID) +``` + +# Calibrating parameters + +The default parameters sites are already calibrated, but users can modify which parameter to calibrate, or +during which period, priors, etc. (see full documentation). + +```jl +julia> calibrate(site_ID, parameters_to_calibrate) +``` + diff --git a/lib/ClimaLandSimulations/experiments/ozark.jl b/lib/ClimaLandSimulations/experiments/ozark.jl new file mode 100644 index 0000000000..678930f67f --- /dev/null +++ b/lib/ClimaLandSimulations/experiments/ozark.jl @@ -0,0 +1,6 @@ +using ClimaLandSimulations + +sv, sol, Y, p, inputs, climalsm = fluxnet_simulation("US-MOz"); + +# Will save figures in current repo /figures +make_plots(inputs, climalsm) diff --git a/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl b/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl new file mode 100644 index 0000000000..1f1863d9f6 --- /dev/null +++ b/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl @@ -0,0 +1,57 @@ +module ClimaLandSimulations + +import SciMLBase +import ClimaTimeSteppers as CTS +using ClimaCore +import CLIMAParameters as CP +using Statistics +using Dates +using Insolation +using StatsBase +using Interpolations +using ClimaLSM +using ClimaLSM.Domains: Column +using ClimaLSM.Soil +using ClimaLSM.Soil.Biogeochemistry +using ClimaLSM.Canopy +using ClimaLSM.Canopy.PlantHydraulics +import ClimaLSM.Parameters as LSMP +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 HTTP +using JSON +using ArtifactWrappers +using DelimitedFiles +using Dierckx +using Thermodynamics +using Formatting +using CairoMakie +using LaTeXStrings +using PlotUtils: optimize_ticks +using DataFrames + +@compound CO₂ +const FT = Float64 + +climalsm_dir = pkgdir(ClimaLSM) +savedir = + joinpath(climalsm_dir, "experiments", "integrated", "fluxnet/figures/") +include(joinpath(climalsm_dir, "parameters", "create_parameters.jl")) +earth_param_set = create_lsm_parameters(FT) + +ClimaLandSimulations_dir = pkgdir(ClimaLandSimulations) + +include(joinpath("utilities", "makie_plots.jl")) +include(joinpath("utilities", "data_tools.jl")) +include(joinpath("utilities", "pull_MODIS.jl")) +include(joinpath("utilities", "met_drivers_FLUXNET.jl")) +include(joinpath("utilities", "climalsm_output_dataframe.jl")) +include(joinpath("utilities", "inputs_dataframe.jl")) +include("fluxnet_simulation.jl") + +function __init__() + Unitful.register(ClimaLandSimulations) +end + +end diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_parameters.jl b/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_parameters.jl new file mode 100644 index 0000000000..97bd225d59 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_parameters.jl @@ -0,0 +1,125 @@ +"""Site-specific model parameters for running CliMA LSM on the Harvard Forest +fluxtower site.""" + +# Data File download link +data_link = "https://caltech.box.com/shared/static/xixaod6511cutz51ag81k1mtvy05hbol.csv" + +# Timezone (offset from UTC in hrs) +time_offset = 5 + +# Site latitude and longitude +lat = FT(42.5378) # degree +long = FT(-72.1715) # degree + +# Height of sensor +atmos_h = FT(30) +# https://atmos.seas.harvard.edu/research-harvard_forest-instrumentation + +## LAI data from MODIS +# Heterotrophic respiration parameters +θ_a100 = FT(0.1816) +D_ref = FT(1.39e-5) +b = FT(4.547) +D_liq = FT(3.17) +# DAMM +α_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) + +# Autotrophic respiration parameters +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) + +# Soil parameters +soil_ν = FT(0.5) # m3/m3 +soil_K_sat = FT(4e-7) # m/s, matches Natan +soil_S_s = FT(1e-3) # 1/m, guess +soil_vg_n = FT(2.05) # unitless +soil_vg_α = FT(0.04) # inverse meters +θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + +# Soil heat transfer parameters; not needed for hydrology only test +ν_ss_quartz = FT(0.1) +ν_ss_om = FT(0.1) +ν_ss_gravel = FT(0.0); +κ_quartz = FT(7.7) # W/m/K +κ_minerals = FT(2.5) # W/m/K +κ_om = FT(0.25) # W/m/K +κ_liq = FT(0.57) # W/m/K +κ_ice = FT(2.29) # W/m/K +κ_air = FT(0.025); #W/m/K +ρp = FT(2700); # kg/m^3 +κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals) +κ_dry_val = Soil.κ_dry(ρp, soil_ν, κ_solid, κ_air) +κ_sat_frozen_val = Soil.κ_sat_frozen(κ_solid, soil_ν, κ_ice) +κ_sat_unfrozen_val = Soil.κ_sat_unfrozen(κ_solid, soil_ν, κ_liq); +ρc_ds = FT((1 - soil_ν) * 4e6); # J/m^3/K +z_0m_soil = FT(0.01) +z_0b_soil = FT(0.001) +soil_ϵ = FT(0.98) +soil_α_PAR = FT(0.2) +soil_α_NIR = FT(0.2) + +# TwoStreamModel parameters +Ω = FT(0.69) +ld = 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) + +# Energy Balance model +ac_canopy = FT(2.5e3) + +# Conductance Model +g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. +Drel = FT(1.6) +g0 = FT(1e-4) + +#Photosynthesis model +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) + +# Plant Hydraulics and general plant parameters +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 +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) +capacity = FT(10) # kg/m^2 +plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m +rooting_depth = FT(0.5) # from Natan +z0_m = FT(0.13) * h_canopy +z0_b = FT(0.1) * z0_m diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_simulation.jl b/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_simulation.jl new file mode 100644 index 0000000000..b0bff732b4 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-Ha1/US-Ha1_simulation.jl @@ -0,0 +1,26 @@ +"""This file contains simulation variables for running CliMA LSM on the US-Ha1 +fluxtower site. This includes both the domain variables and timestepping +variables for running the simulation.""" + +# DOMAIN SETUP: + +# Column dimensions - separation of layers at the top and bottom of the column: +dz_bottom = FT(1.5) +dz_top = FT(0.025) + +# Stem and leaf compartments and their heights: +n_stem = Int64(1) +n_leaf = Int64(1) +h_leaf = FT(12) # m +h_stem = FT(14) # m + +# TIME STEPPING: + +# Starting time: +t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow + +# Time step size: +dt = Float64(40) + +# Number of timesteps between saving output: +n = 45 diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_parameters.jl b/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_parameters.jl new file mode 100644 index 0000000000..e79edcfbab --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_parameters.jl @@ -0,0 +1,123 @@ +"""Site-specific model parameters for running CliMA LSM on the Ozark +fluxtower site.""" + +# Data File download link +data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv" + +# Timezone (offset from UTC in hrs) +time_offset = 7 + +# Height of sensor on flux tower +atmos_h = FT(32) + +# Site latitude and longitude +lat = FT(38.7441) # degree +long = FT(-92.2000) # degree + +# Heterotrophic respiration parameters +θ_a100 = FT(0.1816) +D_ref = FT(1.39e-5) +b = FT(4.547) +D_liq = FT(3.17) +# DAMM +α_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) + +# Autotrophic respiration parameters +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) + +# Soil parameters +soil_ν = FT(0.5) # m3/m3 +soil_K_sat = FT(4e-7) # m/s, matches Natan +soil_S_s = FT(1e-3) # 1/m, guess +soil_vg_n = FT(2.05) # unitless +soil_vg_α = FT(0.04) # inverse meters +θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + +# Soil heat transfer parameters; not needed for hydrology only test +ν_ss_quartz = FT(0.1) +ν_ss_om = FT(0.1) +ν_ss_gravel = FT(0.0); +κ_quartz = FT(7.7) # W/m/K +κ_minerals = FT(2.5) # W/m/K +κ_om = FT(0.25) # W/m/K +κ_liq = FT(0.57) # W/m/K +κ_ice = FT(2.29) # W/m/K +κ_air = FT(0.025); #W/m/K +ρp = FT(2700); # kg/m^3 +κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals) +κ_dry_val = Soil.κ_dry(ρp, soil_ν, κ_solid, κ_air) +κ_sat_frozen_val = Soil.κ_sat_frozen(κ_solid, soil_ν, κ_ice) +κ_sat_unfrozen_val = Soil.κ_sat_unfrozen(κ_solid, soil_ν, κ_liq); +ρc_ds = FT((1 - soil_ν) * 4e6); # J/m^3/K +z_0m_soil = FT(0.01) +z_0b_soil = FT(0.001) +soil_ϵ = FT(0.98) +soil_α_PAR = FT(0.2) +soil_α_NIR = FT(0.2) + +# TwoStreamModel parameters +Ω = FT(0.69) +ld = 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) + +# Energy Balance model +ac_canopy = FT(2.5e3) + +# Conductance Model +g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. +Drel = FT(1.6) +g0 = FT(1e-4) + +#Photosynthesis model +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) + +# Plant Hydraulics and general plant parameters +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 +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) +capacity = FT(10) # kg/m^2 +plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m +rooting_depth = FT(0.5) # from Natan +z0_m = FT(0.13) * h_canopy +z0_b = FT(0.1) * z0_m diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_simulation.jl b/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_simulation.jl new file mode 100644 index 0000000000..a17d64f1d1 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-MOz/US-MOz_simulation.jl @@ -0,0 +1,26 @@ +"""This file contains simulation variables for running CliMA LSM on the US-MOz +fluxtower site. This includes both the domain variables and timestepping +variables for running the simulation.""" + +# DOMAIN SETUP: + +# Column dimensions - separation of layers at the top and bottom of the column: +dz_bottom = FT(1.5) +dz_top = FT(0.025) + +# Stem and leaf compartments and their heights: +n_stem = Int64(1) +n_leaf = Int64(1) +h_stem = FT(9) # m +h_leaf = FT(9.5) # m + +# TIME STEPPING: + +# Starting time: +t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow + +# Time step size: +dt = Float64(120) + +# Number of timesteps between saving output: +n = 15 diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_parameters.jl b/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_parameters.jl new file mode 100644 index 0000000000..830355153f --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_parameters.jl @@ -0,0 +1,128 @@ +"""Site-specific model parameters for running CliMA LSM on the Niwot Ridge +fluxtower site.""" + +# Data download link +data_link = "https://caltech.box.com/shared/static/r6gvldgabk3mvtx53gevnlaq1ztsk41i.csv" + +# Timezone (offset from UTC in hrs) +time_offset = 7 + +# Site latitude and longitude +lat = FT(40.0329) # degree +long = FT(-105.5464) # degree + +# Height of sensor +atmos_h = FT(21.5) +# Metzger, Stefan & Burba, George & Burns, Sean & Blanken, Peter & Li, +# Jiahong & Luo, Hongyan & Zulueta, Rommel. (2016). Optimization of an enclosed +# gas analyzer sampling system for measuring eddy covariance fluxes of H2O and +# CO2. Atmospheric Measurement Techniques. 9. 1341-1359. 10.5194/amt-9-1341-2016. + +## LAI data from MODIS +# Heterotrophic respiration parameters +θ_a100 = FT(0.1816) +D_ref = FT(1.39e-5) +b = FT(4.547) +D_liq = FT(3.17) +# DAMM +α_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) + +# Autotrophic respiration parameters +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) + +# Soil parameters +soil_ν = FT(0.45) # m3/m3 +soil_K_sat = FT(4e-7) # m/s, matches Natan +soil_S_s = FT(1e-3) # 1/m, guess +soil_vg_n = FT(2.05) # unitless +soil_vg_α = FT(0.04) # inverse meters +θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021 + +# Soil heat transfer parameters; not needed for hydrology only test +ν_ss_quartz = FT(0.1) +ν_ss_om = FT(0.1) +ν_ss_gravel = FT(0.0); +κ_quartz = FT(7.7) # W/m/K +κ_minerals = FT(2.5) # W/m/K +κ_om = FT(0.25) # W/m/K +κ_liq = FT(0.57) # W/m/K +κ_ice = FT(2.29) # W/m/K +κ_air = FT(0.025); #W/m/K +ρp = FT(2700); # kg/m^3 +κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals) +κ_dry_val = Soil.κ_dry(ρp, soil_ν, κ_solid, κ_air) +κ_sat_frozen_val = Soil.κ_sat_frozen(κ_solid, soil_ν, κ_ice) +κ_sat_unfrozen_val = Soil.κ_sat_unfrozen(κ_solid, soil_ν, κ_liq); +ρc_ds = FT((1 - soil_ν) * 4e6); # J/m^3/K +z_0m_soil = FT(0.1) +z_0b_soil = FT(0.1) +soil_ϵ = FT(0.98) +soil_α_PAR = FT(0.2) +soil_α_NIR = FT(0.2) + +# TwoStreamModel parameters +Ω = FT(0.71) +ld = 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) + +# Energy Balance model +ac_canopy = FT(3e3) + +# Conductance Model +g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300. +Drel = FT(1.6) +g0 = FT(1e-4) + +#Photosynthesis model +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) + +# Plant Hydraulics and general plant parameters +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 +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) +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) +z0_m = FT(0.13) * h_canopy +z0_b = FT(0.1) * z0_m diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_simulation.jl b/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_simulation.jl new file mode 100644 index 0000000000..9fd32d130e --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-NR1/US-NR1_simulation.jl @@ -0,0 +1,26 @@ +"""This file contains simulation variables for running CliMA LSM on the US-NR1 +fluxtower site. This includes both the domain variables and timestepping +variables for running the simulation.""" + +# DOMAIN SETUP: + +# Column dimensions - separation of layers at the top and bottom of the column: +dz_bottom = FT(1.25) +dz_top = FT(0.05) + +# Stem and leaf compartments and their heights: +n_stem = Int64(1) +n_leaf = Int64(1) +h_leaf = FT(6.5) # m +h_stem = FT(7.5) # m + +# TIME STEPPING: + +# Starting time: +t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow + +# Time step size: +dt = Float64(40) + +# Number of timesteps between saving output: +n = 45 diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_parameters.jl b/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_parameters.jl new file mode 100644 index 0000000000..5679435b50 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_parameters.jl @@ -0,0 +1,132 @@ +"""Site-specific model parameters for running CliMA LSM on the Vaira Ranch +fluxtower site.""" + +## Some site parameters have been taken from +## Ma, S., Baldocchi, D. D., Xu, L., Hehn, T. (2007) +## Inter-Annual Variability In Carbon Dioxide Exchange Of An +## Oak/Grass Savanna And Open Grassland In California, Agricultural +## And Forest Meteorology, 147(3-4), 157-171. https://doi.org/10.1016/j.agrformet.2007.07.008 +## CLM 5.0 Tech Note: https://www2.cesm.ucar.edu/models/cesm2/land/CLM50_Tech_Note.pdf +# Bonan, G. Climate change and terrestrial ecosystem modeling. Cambridge University Press, 2019. + +# Data download link +data_link = "https://caltech.box.com/shared/static/dx0p5idbsbrdebsda10t9pfv2lbdaz95.csv" +LAI_link = "https://caltech.box.com/shared/static/y5vf8s9qkoogglc1bc2eyu1k95sbjsc3.csv" + +# Height of sensor +atmos_h = FT(2) + +# Timezone (offset from UTC in hrs) +time_offset = 8 + +# Site latitude and longitude +lat = FT(38.4133) # degree +long = FT(-120.9508) # degree + +# Heterotrophic respiration parameters +θ_a100 = FT(0.1816) +D_ref = FT(1.39e-5) +b = FT(4.547) +D_liq = FT(3.17) +# DAMM +α_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) + +# Autotrophic respiration parameters +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) + +# Soil parameters +soil_ν = FT(0.45) # m3/m3 +soil_K_sat = FT(0.45 / 3600 / 100) # m/s, +soil_S_s = FT(1e-3) # 1/m, guess +soil_vg_n = FT(2.0) # unitless +soil_vg_α = FT(2.0) # inverse meters +θ_r = FT(0.067) # m3/m3, + +# Soil heat transfer parameters; not needed for hydrology only test +ν_ss_quartz = FT(0.38) +ν_ss_om = FT(0.0) +ν_ss_gravel = FT(0.1); +κ_quartz = FT(7.7) # W/m/K +κ_minerals = FT(2.5) # W/m/K +κ_om = FT(0.25) # W/m/K +κ_liq = FT(0.57) # W/m/K +κ_ice = FT(2.29) # W/m/K +κ_air = FT(0.025); #W/m/K +ρp = FT(2700); # kg/m^3 +κ_solid = Soil.κ_solid(ν_ss_om, ν_ss_quartz, κ_om, κ_quartz, κ_minerals) +κ_dry_val = Soil.κ_dry(ρp, soil_ν, κ_solid, κ_air) +κ_sat_frozen_val = Soil.κ_sat_frozen(κ_solid, soil_ν, κ_ice) +κ_sat_unfrozen_val = Soil.κ_sat_unfrozen(κ_solid, soil_ν, κ_liq); +ρc_ds = FT((1 - soil_ν) * 4e6); # J/m^3/K +z_0m_soil = FT(0.01) +z_0b_soil = FT(0.001) +soil_ϵ = FT(0.98) +soil_α_PAR = FT(0.3) +soil_α_NIR = FT(0.4) + +# TwoStreamModel parameters +Ω = FT(1.0) +ld = FT(0.5) +α_PAR_leaf = FT(0.11) +λ_γ_PAR = FT(5e-7) +λ_γ_NIR = FT(1.65e-6) +τ_PAR_leaf = FT(0.05) +α_NIR_leaf = FT(0.35) +τ_NIR_leaf = FT(0.34) +n_layers = UInt64(20) +ϵ_canopy = FT(0.97) + +# Conductance Model +g1 = FT(166) # CLM C3 grass +Drel = FT(1.6) +g0 = FT(1e-4) + +#Photosynthesis model +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(4.225e-5) # CLM C3 grass, Slevin et al. 2015 +Γ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) + +# Energy Balance model +ac_canopy = FT(745) + +# Plant Hydraulics and general plant parameters +SAI = FT(0) +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 +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) +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 Bonan Table 2.3 +z0_m = FT(0.13) * h_canopy +z0_b = FT(0.1) * z0_m diff --git a/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_simulation.jl b/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_simulation.jl new file mode 100644 index 0000000000..186aac3eb8 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet/US-Var/US-Var_simulation.jl @@ -0,0 +1,26 @@ +"""This file contains simulation variables for running CliMA LSM on the US-Var +fluxtower site. This includes both the domain variables and timestepping +variables for running the simulation.""" + +# DOMAIN SETUP: + +# Column dimensions - separation of layers at the top and bottom of the column: +dz_bottom = FT(1.0) +dz_top = FT(0.05) + +# Stem and leaf compartments and their heights: +n_stem = Int64(0) +n_leaf = Int64(1) +h_leaf = FT(0.7) # m +h_stem = FT(0) # m + +# TIME STEPPING: + +# Starting time: +t0 = Float64(21 * 3600 * 24)# start day 21 of the year + +# Time step size: +dt = Float64(40) + +# Number of timesteps between saving output: +n = 45 diff --git a/lib/ClimaLandSimulations/src/fluxnet_simulation.jl b/lib/ClimaLandSimulations/src/fluxnet_simulation.jl new file mode 100644 index 0000000000..0b37dacd46 --- /dev/null +++ b/lib/ClimaLandSimulations/src/fluxnet_simulation.jl @@ -0,0 +1,353 @@ +export fluxnet_simulation + +""" + fluxnet_simulation(site_ID; kwargs) + +Run ClimaLSM at a fluxnet site. +""" +function fluxnet_simulation( + site_ID; + FT = Float64, + path_to_sim = joinpath( + ClimaLandSimulations_dir, + "src", + "fluxnet", + "$site_ID", + "$(site_ID)_simulation.jl", + ), + path_to_domain_setup = joinpath( + ClimaLandSimulations_dir, + "src", + "utilities", + "domain_setup.jl", + ), + path_to_params = joinpath( + ClimaLandSimulations_dir, + "src", + "fluxnet", + "$site_ID", + "$(site_ID)_parameters.jl", + ), + path_to_timestepper_setup = joinpath( + ClimaLandSimulations_dir, + "src", + "utilities", + "timestepper_setup.jl", + ), +) # why is there 2? + + include(path_to_sim) + include(path_to_domain_setup) + include(path_to_params) + LOCAL_DATETIME, + atmos_co2, + DATA_DT, + drivers, + atmos, + radiation, + LAIfunction, + maxLAI, + RAI, + plant_ν = setup_drivers(site_ID) + include(path_to_timestepper_setup) + + # Now we set up the model. For the soil model, we pick + # a model type and model args: + soil_domain = land_domain + soil_ps = Soil.EnergyHydrologyParameters{FT}(; + κ_dry = κ_dry_val, + κ_sat_frozen = κ_sat_frozen_val, + κ_sat_unfrozen = κ_sat_unfrozen_val, + ρc_ds = ρc_ds, + ν = soil_ν, + ν_ss_om = ν_ss_om, + ν_ss_quartz = ν_ss_quartz, + ν_ss_gravel = ν_ss_gravel, + hydrology_cm = vanGenuchten(; α = soil_vg_α, n = soil_vg_n), + K_sat = soil_K_sat, + S_s = soil_S_s, + θ_r = θ_r, + earth_param_set = earth_param_set, + z_0m = z_0m_soil, + z_0b = z_0b_soil, + emissivity = soil_ϵ, + PAR_albedo = soil_α_PAR, + NIR_albedo = soil_α_NIR, + ) + + 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}(; + ν = soil_ν, # same as soil + θ_a100 = θ_a100, + D_ref = D_ref, + b = b, + D_liq = D_liq, + # DAMM + α_sx = α_sx, + Ea_sx = Ea_sx, + kM_sx = kM_sx, + kM_o2 = kM_o2, + O2_a = O2_a, + D_oa = D_oa, + p_sx = p_sx, + earth_param_set = earth_param_set, + ) + + # soil microbes args + Csom = (z, t) -> eltype(z)(5.0) + + soilco2_top_bc = Soil.Biogeochemistry.SoilCO2StateBC((p, t) -> atmos_co2(t)) + soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux + soilco2_sources = (MicrobeProduction{FT}(),) + + soilco2_boundary_conditions = + (; top = (CO2 = soilco2_top_bc,), bottom = (CO2 = soilco2_bot_bc,)) + + soilco2_drivers = Soil.Biogeochemistry.SoilDrivers( + Soil.Biogeochemistry.PrognosticMet{FT}(), + Soil.Biogeochemistry.PrescribedSOC{FT}(Csom), + atmos, + ) + + soilco2_args = (; + boundary_conditions = soilco2_boundary_conditions, + sources = soilco2_sources, + domain = soil_domain, + parameters = soilco2_ps, + drivers = soilco2_drivers, + ) + + # Now we set up the canopy model, which we set up by component: + # Component Types + canopy_component_types = (; + autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT}, + radiative_transfer = Canopy.TwoStreamModel{FT}, + photosynthesis = Canopy.FarquharModel{FT}, + conductance = Canopy.MedlynConductanceModel{FT}, + hydraulics = Canopy.PlantHydraulicsModel{FT}, + energy = Canopy.BigLeafEnergyModel{FT}, + ) + # Individual Component arguments + # Set up autotrophic respiration + autotrophic_respiration_args = (; + parameters = AutotrophicRespirationParameters{FT}(; + ne = ne, + ηsl = ηsl, + σl = σl, + μr = μr, + μs = μs, + f1 = f1, + f2 = f2, + ) + ) + # Set up radiative transfer + radiative_transfer_args = (; + parameters = TwoStreamParameters{FT}(; + Ω = Ω, + ld = ld, + α_PAR_leaf = α_PAR_leaf, + λ_γ_PAR = λ_γ_PAR, + λ_γ_NIR = λ_γ_NIR, + τ_PAR_leaf = τ_PAR_leaf, + α_NIR_leaf = α_NIR_leaf, + τ_NIR_leaf = τ_NIR_leaf, + n_layers = n_layers, + ϵ_canopy = ϵ_canopy, + ) + ) + # Set up conductance + conductance_args = (; + parameters = MedlynConductanceParameters{FT}(; + g1 = g1, + Drel = Drel, + g0 = g0, + ) + ) + # Set up photosynthesis + photosynthesis_args = (; + parameters = FarquharParameters{FT}( + Canopy.C3(); + oi = oi, + ϕ = ϕ, + θj = θj, + f = f, + sc = sc, + pc = pc, + Vcmax25 = Vcmax25, + Γstar25 = Γstar25, + Kc25 = Kc25, + Ko25 = Ko25, + To = To, + ΔHkc = ΔHkc, + ΔHko = ΔHko, + ΔHVcmax = ΔHVcmax, + ΔHΓstar = ΔHΓstar, + ΔHJmax = ΔHJmax, + ΔHRd = ΔHRd, + ) + ) + # Set up plant hydraulics + ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) + + function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m + end + + plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, + ) + plant_hydraulics_args = ( + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_midpoints = compartment_midpoints, + compartment_surfaces = compartment_surfaces, + ) + + energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(ac_canopy),) + + # Canopy component args + canopy_component_args = (; + autotrophic_respiration = autotrophic_respiration_args, + radiative_transfer = radiative_transfer_args, + photosynthesis = photosynthesis_args, + conductance = conductance_args, + hydraulics = plant_hydraulics_args, + energy = energy_args, + ) + + # Other info needed + shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, + ) + + canopy_model_args = (; parameters = shared_params, domain = canopy_domain) + + # Integrated plant hydraulics and soil model + land_input = (atmos = atmos, radiation = radiation) + land = SoilCanopyModel{FT}(; + soilco2_type = soilco2_type, + soilco2_args = soilco2_args, + land_args = land_input, + soil_model_type = soil_model_type, + soil_args = soil_args, + canopy_component_types = canopy_component_types, + canopy_component_args = canopy_component_args, + canopy_model_args = canopy_model_args, + ) + Y, p, cds = initialize(land) + exp_tendency! = make_exp_tendency(land) + + #Initial conditions + Y.soil.ϑ_l = + drivers.SWC.status != absent ? + drivers.SWC.values[1 + Int(round(t0 / DATA_DT))] : 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.TS.status != absent ? + drivers.TS.values[1 + Int(round(t0 / DATA_DT))] : + drivers.TA.values[1 + Int(round(t0 / DATA_DT))] + 40# Get soil temperature at t0 + ρc_s = + volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + Ref(land.soil.parameters), + ) + Y.soil.ρe_int = + volumetric_internal_energy.( + Y.soil.θ_i, + ρc_s, + T_0, + Ref(land.soil.parameters), + ) + Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air + ψ_stem_0 = FT(-1e5 / 9800) # pressure in the leaf divided by rho_liquid*gravitational acceleration [m] + ψ_leaf_0 = FT(-2e5 / 9800) + ψ_comps = n_stem > 0 ? [ψ_stem_0, ψ_leaf_0] : ψ_leaf_0 + + S_l_ini = + inverse_water_retention_curve.( + retention_model, + ψ_comps, + plant_ν, + plant_S_s, + ) + + for i in 1:(n_stem + n_leaf) + Y.canopy.hydraulics.ϑ_l.:($i) .= + augmented_liquid_fraction.(plant_ν, S_l_ini[i]) + end + + Y.canopy.energy.T = drivers.TA.values[1 + Int(round(t0 / DATA_DT))] # Get atmos temperature at t0 + + set_initial_cache! = make_set_initial_cache(land) + set_initial_cache!(p, Y, t0) + + # Simulation + sv = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), + ) + saving_cb = ClimaLSM.NonInterpSavingCallback(sv, 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(t0:DATA_DT:tf) + updatefunc = ClimaLSM.make_update_drivers(atmos, radiation) + driver_cb = ClimaLSM.DriverUpdateCallback(updateat, updatefunc) + cb = SciMLBase.CallbackSet(driver_cb, saving_cb) + + prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction((T_exp!) = exp_tendency!), + Y, + (t0, tf), + p, + ) + sol = SciMLBase.solve( + prob, + ode_algo; + dt = dt, + callback = cb, + adaptive = false, + saveat = saveat, + ) + + inputs, inputs_SI, inputs_commonly_used = + make_inputs_df(LOCAL_DATETIME, drivers) # why does it needs args? + climalsm = make_output_df(sv, inputs) + + if isfile( + joinpath( + climalsm_dir, + "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", + ), + ) + rm( + joinpath( + climalsm_dir, + "lib/ClimaLandSimulations/src/utilities/$site_ID/Artifacts.toml", + ), + ) + else + nothing + end + + return sv, sol, Y, p, inputs, climalsm + +end + + + diff --git a/lib/ClimaLandSimulations/src/utilities/climalsm_output_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/climalsm_output_dataframe.jl new file mode 100644 index 0000000000..c8f406e5b0 --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/climalsm_output_dataframe.jl @@ -0,0 +1,65 @@ +export getoutput, make_output_df + +""" + getoutput(variable::Symbol, variables::Symbol...; result = sv.saveval, depth = 1) + +Return a vector of FT corresponding to the variable of interest at all times. +By default, get output from sv.saveval, but user can specify e.g., result = sol.u +By default, get surface value, but user can specify depth for e.g., soil temperature +""" +function getoutput( + sv, + variable::Symbol, + variables::Symbol...; + result = sv.saveval, + depth = 1, +) + for v in (variable, variables...) + result = getproperty.(result, v) + end + return [parent(r)[depth] for r in result] +end +# example: getoutput(:SW_out) +# example 2: getoutput(:canopy, :conductance, :gs) +# example 3: getoutput(:soil, :T; result = sol.u, depth = 2) + +""" + make_output_df(sv, inputs) + +Return a dataframe containing climalsm outputs +""" +function make_output_df(sv, inputs) + # List of output that we want + output_list = ( + (:SW_out,), + (:LW_out,), + (:canopy, :conductance, :gs), + (:canopy, :conductance, :transpiration), + (:canopy, :autotrophic_respiration, :Ra), + (:canopy, :photosynthesis, :GPP), + (:canopy, :hydraulics, :β), + (:canopy, :hydraulics, :area_index, :leaf), + (:canopy, :energy, :lhf), + (:soil, :turbulent_fluxes, :shf), + (:soil, :turbulent_fluxes, :lhf), + (:soil, :T), + (:soil, :θ_l), + (:soil, :turbulent_fluxes, :vapor_flux), + ) + + output_vectors = [getoutput(sv, args...) for args in output_list] + # Extract the last symbol from each tuple for column names + column_names = [names[end] for names in output_list] + # Create a dictionary with simplified column names and corresponding vectors + data_dict = Dict(zip(column_names, output_vectors)) + # Create a DataFrame from the dictionary + climalsm = DataFrame(data_dict) + # now if I want for example GPP, I can just do df.GPP + + index_t_start = 120 * 48 # we shouldn't hardcode that 120 in ozark_simulation.jl + index_t_end = 120 * 48 + (N_days - N_spinup_days) * 48 + model_dt = inputs.DateTime[index_t_start:index_t_end] + + insertcols!(climalsm, 1, :DateTime => model_dt) + return climalsm +end diff --git a/lib/ClimaLandSimulations/src/utilities/data_tools.jl b/lib/ClimaLandSimulations/src/utilities/data_tools.jl new file mode 100644 index 0000000000..5243362fbc --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/data_tools.jl @@ -0,0 +1,174 @@ +"""Data utilities for running fluxnet site experiments. Provides the DataColumn +representation of a data column read from the FLUXNET data which includes +a data status flag and units for the data. This enables the user to easily +convert units and check for missing data. Data is automatically filled when +read in based on the quality control flags provided in the FLUXNET data or +by simply filling with the mean value in the column.""" + +using Interpolations +using StatsBase + +# Define the valid data statuses +@enum DataStatus complete = 1 absent = 2 incomplete = 3 + +""" + replace_poor_quality_with_mean!(field, flag) + +Replace values indicated to be poor quality by a fluxnet +QC `flag` (Array) with the mean value +in the `field` (Array). + +This uses the Fluxnet convention of 0 = measured, 1 = good quality +gap fill, 2 = medium quality, and 3 = poor quality, and replaces +data with QC flag = 2,3. +""" +function replace_poor_quality_with_mean!(field, flag) + good_indices = (flag .== 0) .|| (flag .== 1) + fill_value = mean(field[good_indices]) + field[.~good_indices] .= fill_value + return field +end + +""" + replace_replace_missing_with_mean_by_value!(field) + +Replace values indicated to be missing + with the mean value in the `field` (Array). + +This uses the Fluxnet convention of a value of -9999 indicating +missing data. +""" +function replace_missing_with_mean_by_value!(field) + good_indices = .~(field .== -9999) + fill_value = mean(field[good_indices]) + field[.~good_indices] .= fill_value + return field +end + +""" + column_status(driver_data::Matrix, column_name::String) + +Checks that the `driver_data` matrix contains the column +specified by the string `column_name`; checks if that the column does not contain +missing values. Returns the corresponding DataStatus (absent if column name +is not present, incomplete if missing data is present, or complete if no missing +data is present). +""" +function column_status(driver_data::Matrix, column_name::String) + col_dat = [] + column_names = driver_data[1, :] + try + col_dat = driver_data[2:end, column_names .== column_name] + catch _ + return col_dat, absent + end + if isempty(col_dat) + return col_dat, absent + end + try + @assert(minimum(col_dat) > -9999) + catch _ + if maximum(col_dat) == -9999 + return col_dat, absent + end + return col_dat, incomplete + end + return col_dat, complete +end + +""" + DataColumn + +A struct for storing data columns along with thir units and +a status flag which indicates the data is complete (but not neccessarily QC +checked), absent entirely, or incomplete. +""" +struct DataColumn + "A nx1 matrix, where n is the number of data points, of data values" + values::Matrix + "The units of those values" + units::String + "The missing-data status of the column" + status::DataStatus +end + +""" + filter_column(driver_data::Matrix, column_name::String, units::String) + +Checks the `driver_data` for the present of `column_name`; returns a +DataColumn struct with cleaned data and a complete status if present, +returns an empty DataColumn with status absent if not. + +- Checks to see if a quality control flag +is available for the data column, and if so, replaces poor quality values with +the mean value per the QC flag. + - Checks for missing data and replaces with the mean value +- Returns an `absent` status if the column is not present. +""" +function filter_column(driver_data::Matrix, column_name::String, units::String) + # Check that the data exists and read it in if so + col_dat, status = column_status(driver_data, column_name) + # Set missing data threshold above which column is treated as absent + missing_threshold = 0.5 + # Set poor quality threshold above which the column undergoes no replacement using the quality flag + quality_threshold = 0.5 + # if it does not exist, exit + if status == absent + @info "Warning: Data for $column_name is absent" + return DataColumn(reshape([], 0, 0), "", status) + elseif status == complete # all data is present but some may be poor quality + # Check for a QC flag column and use it to clean data if possible + column_QC = column_name * "_QC" + QC_col, QCstatus = column_status(driver_data, column_QC) + if QCstatus != absent + # QC flag = [0,1] -> good data, [2,3] -> poorer data + num_poor = count(x -> x >= 2, QC_col) + percent_poor = 100.0 * num_poor / length(col_dat) + # If the percent_poor is below a threshold, fill any missing data with the mean value per the QC flag, + # otherwise return data column with a specific warning but use poorer quality data directly + if percent_poor < quality_threshold + replace_poor_quality_with_mean!(col_dat, QC_col) + @info "Warning: Data for $column_name $(round(percent_poor, sigdigits=3))% poor quality. Filled with mean value using QC flag" + return DataColumn(col_dat, units, status) + else + @info "Warning: Data for $column_name $(round(percent_poor, sigdigits=3))% poor quality. Returning with no replacement." + return DataColumn(col_dat, units, status) + end + else + @info "Information: Data for $column_name is complete and no QC flag present" + return DataColumn(col_dat, units, status) + end + elseif status == incomplete # data is missing + # Replace values of -9999 with mean + num_missing = count(x -> x == -9999, col_dat) + percent_missing = 100.0 * num_missing / length(col_dat) + if percent_missing > missing_threshold + @info "Warning: Data for $column_name $(round(percent_missing, + sigdigits=3))% has value of -9999. Treating as absent" + return DataColumn(col_dat, "", absent) + end + replace_missing_with_mean_by_value!(col_dat) + status = complete + @info "Warning: Data for $column_name $(round(percent_missing, + sigdigits=3))% has value of -9999. Filled with mean value" + return DataColumn(col_dat, units, status) + end +end + +""" + transform_column(column::DataColumn, + transform::Function, + units::String, + ) + +Returns a new DataColumn struct with the values of the input `column` +transformed via the function `transform`. Useful for unit conversions. Sets +the units of the column to the string units. +""" +function transform_column( + column::DataColumn, + transform::Function, + units::String, +) + return DataColumn(transform.(column.values), units, column.status) +end diff --git a/lib/ClimaLandSimulations/src/utilities/domain_setup.jl b/lib/ClimaLandSimulations/src/utilities/domain_setup.jl new file mode 100644 index 0000000000..cbda540545 --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/domain_setup.jl @@ -0,0 +1,14 @@ +nelements = 20 +zmin = FT(-10) +zmax = FT(0) +h_canopy = h_stem + h_leaf +compartment_midpoints = + n_stem > 0 ? [h_stem / 2, h_stem + h_leaf / 2] : [h_leaf / 2] +compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] + +land_domain = Column(; + zlim = (zmin, zmax), + nelements = nelements, + dz_tuple = (dz_bottom, dz_top), +) +canopy_domain = ClimaLSM.Domains.obtain_surface_domain(land_domain) diff --git a/lib/ClimaLandSimulations/src/utilities/inputs_dataframe.jl b/lib/ClimaLandSimulations/src/utilities/inputs_dataframe.jl new file mode 100644 index 0000000000..a37a41eb7b --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/inputs_dataframe.jl @@ -0,0 +1,112 @@ +export make_inputs_df + +""" + make_inputs_df(LOCAL_DATETIME, drivers) + +Returns DataFrames of inputs, one unitless, one in SI units, one in commonly used units. +""" +function make_inputs_df(LOCAL_DATETIME, drivers) + + variables_name = ( + "DateTime", + "RECO", + "TA", + "VPD", + "PA", + "P", + "WS", + "LW_IN", + "SW_IN", + "CO2", + "SWC", + "TS", + "GPP", + "LE", + "H", + "G", + "SW_OUT", + ) # I imagine this may not work for all fluxnet sites... + + variables = [ + vec(mat) for mat in ( + LOCAL_DATETIME, + drivers.RECO.values, + drivers.TA.values, + drivers.VPD.values, + drivers.PA.values, + drivers.P.values, + drivers.WS.values, + drivers.LW_IN.values, + drivers.SW_IN.values, + drivers.CO2.values, + drivers.SWC.values, + drivers.TS.values, + drivers.GPP.values, + drivers.LE.values, + drivers.H.values, + drivers.G.values, + drivers.SW_OUT.values, + ) + ] + + inputs = DataFrame([ + variables_name[i] => variables[i] for i in 1:length(variables_name) + ]) + + # make a Unitful dataframe with SI units and one with commonly used units + columns = variables_name[2:end] + units = [ + molCO₂ * m^-2 * s^-1, + K, + Pa, + Pa, + m * s^-1, + m * s^-1, + W * m^-2, + W * m^-2, + molCO₂, + m^3 * m^-3, + K, + molCO₂ * m^-2 * s^-1, + W * m^-2, + W * m^-2, + W * m^-2, + W * m^-2, + ] + inputs_SI = copy(inputs) + foreach( + (col, unit) -> inputs_SI[!, col] .= first([inputs_SI[!, col]]unit), + columns, + units, + ) + + units_to = [ + μmolCO₂ * m^-2 * s^-1, + °C, + kPa, + Pa, + mm * s^-1, + m * s^-1, + W * m^-2, + W * m^-2, + μmolCO₂, + m^3 * m^-3, + °C, + μmolCO₂ * m^-2 * s^-1, + W * m^-2, + W * m^-2, + W * m^-2, + W * m^-2, + ] + inputs_commonly_used = copy(inputs_SI) + foreach( + (col, unit_to) -> + inputs_commonly_used[!, col] = + uconvert.(unit_to, inputs_SI[!, col]), + columns, + units_to, + ) + + return (inputs, inputs_SI, inputs_commonly_used) + +end diff --git a/lib/ClimaLandSimulations/src/utilities/makie_plots.jl b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl new file mode 100644 index 0000000000..d285846352 --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/makie_plots.jl @@ -0,0 +1,439 @@ +export timeseries_fluxes_fig, + timeseries_H2O_fig, + fingerprint_fig, + diurnal, + diurnal_plot!, + diurnals_fig, + make_plots + +# Make all sort of plots with data and model output +# 1. Time series (e.g., C fluxes, h2o fluxes, energy fluxes, met drivers) +# 2. Seasonal pattern (i.e., monthly average and std) +# 3. Diurnal pattern (i.e., hourly average and std) +# 4. Response curves (e.g., NEE vs. PAR with VPD color and SWC brightness...) +# 5. Energy conservation (i.e., Rn - G vs. L + H) +# 6. Water budget (i.e., cumulative ET vs. P) +# 7. Data quality plots (e.g., NEE vs. u* by T and SWC bins) +# 8. Fingerprint plots (showing both seasonality and diurnal pattern) +# 9. Wavelet coherence +# to do in another script: animations + +# 1. Time series of GPP, ET and SW_OUT +function timeseries_fluxes_fig(inputs, climalsm, index_t_start, index_t_end) # will run for any inputs or climalsm output of FLUXNET sites + # create an empty figure + fig = Figure(size = (1000, 1000)) # note: do not load Plots.jl in this branch (it is loading in plot_utils) + fontsize_theme = Theme(fontsize = 20) + set_theme!(fontsize_theme) + + # create empty axis, with a specific layout + ax_C = Axis( + fig[1, 1], + ylabel = L"\text{GPP} \, (\mu\text{mol m}^{-2} \, \text{s}^{-1})", + ) # C fluxes + ax_W = Axis(fig[2, 1], ylabel = L"\text{ET (mm)}") # h2o fluxes + ax_SWOUT = Axis( + fig[3, 1], + ylabel = L"\text{SW OUT} \, (\text{W} \, \text{m}^{-2})", + ) # shortwave out + # ax_T = Axis(fig[4, 1]) # air, canopy, and soil temperature + + # for time series, CairoMakie should allow DateTime type soon (but not yet) + # so the 2 lines of code below are a trick to be able to use DateTime - will be removed later + dateticks = + optimize_ticks(climalsm.DateTime[1], climalsm.DateTime[end])[1][2:(end - 1)] # first and last are weirdly placed + + # add plots into axis ax_C + p_GPP_m = CairoMakie.scatter!( + ax_C, + datetime2unix.(climalsm.DateTime), + climalsm.GPP .* 1e6, + color = :blue, + ) + p_GPP_d = CairoMakie.scatter!( + ax_C, + datetime2unix.(inputs.DateTime[index_t_start:index_t_end]), + inputs.GPP[index_t_start:index_t_end] .* 1e6, + color = :black, + ) + + # ax_W + p_ET_m = CairoMakie.scatter!( + ax_W, + datetime2unix.(climalsm.DateTime), + (climalsm.vapor_flux .* 1e3 .* 24 .* 3600) .+ + (climalsm.transpiration .* 1e3 .* 24 .* 3600), + color = :blue, + ) # not sure about units + p_ET_d = CairoMakie.scatter!( + ax_W, + datetime2unix.(inputs.DateTime[index_t_start:index_t_end]), + inputs.LE[index_t_start:index_t_end] ./ + (LSMP.LH_v0(earth_param_set) * 1000) .* (1e3 * 24 * 3600), + color = :black, + ) # not sure units + + # ax_SW_OUT + p_SWOUT_m = CairoMakie.scatter!( + ax_SWOUT, + datetime2unix.(climalsm.DateTime), + climalsm.SW_out, + color = :blue, + ) + p_SWOUT_d = CairoMakie.scatter!( + ax_SWOUT, + datetime2unix.(inputs.DateTime[index_t_start:index_t_end]), + FT.(inputs.SW_OUT[index_t_start:index_t_end]), + color = :black, + ) + + # xticks + ax_C.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + ax_W.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + ax_SWOUT.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + + axislegend( + ax_C, + [p_GPP_d, p_GPP_m], + ["Observations", "ClimaLSM"], + "", + position = :rt, + orientation = :horizontal, + ) + + CairoMakie.ylims!(ax_C, (0, 40)) + CairoMakie.ylims!(ax_W, (0, 30)) + CairoMakie.ylims!(ax_SWOUT, (0, 200)) + + [ + CairoMakie.xlims!( + axes, + ( + datetime2unix(climalsm.DateTime[1]), + datetime2unix(climalsm.DateTime[end]), + ), + ) for axes in [ax_C, ax_W, ax_SWOUT] + ] + + fig + return fig +end + +# 2. Time series of SWC, Precip, moisture stress, stomatal conductance +function timeseries_H2O_fig(inputs, climalsm, index_t_start, index_t_end) # will run for any inputs or climalsm output of FLUXNET sites + # create an empty figure + fig = Figure(size = (1000, 1000)) # note: do not load Plots.jl in this branch (it is loading in plot_utils) + fontsize_theme = Theme(fontsize = 20) + set_theme!(fontsize_theme) + + # create empty axis, with a specific layout + ax_H2O = + Axis(fig[1, 1], ylabel = L"\theta \, (\text{m}^{3} \, \text{m}^{-3})") # soil moisture + ax_H2O_rain = Axis( + fig[1, 1], + ylabel = L"\text{Rainfall (mm)}", + yaxisposition = :right, + ygridvisible = false, + ylabelcolor = :blue, + yticklabelcolor = :blue, + ) + ax_MS = Axis(fig[2, 1], ylabel = L"\text{Moisture stress}") # moisture stress + ax_SC = Axis( + fig[3, 1], + ylabel = L"\text{Stomatal conductance} \, (\text{mol m}^{-2} \, \text{s}^{-1})", + ) # stomatal conductance + # ax_T = Axis(fig[4, 1]) # air, canopy, and soil temperature + + # for time series, CairoMakie should allow DateTime type soon (but not yet) + # so the 2 lines of code below are a trick to be able to use DateTime - will be removed later + dateticks = + optimize_ticks(climalsm.DateTime[1], climalsm.DateTime[end])[1][2:(end - 1)] # first and last are weirdly placed + + # add plots into axis ax_H2O + p_H2O_m = CairoMakie.scatter!( + ax_H2O, + datetime2unix.(climalsm.DateTime), + climalsm.θ_l, + color = :green, + ) + p_H2O_d = CairoMakie.scatter!( + ax_H2O, + datetime2unix.(inputs.DateTime[index_t_start:index_t_end]), + inputs.SWC[index_t_start:index_t_end], + color = :black, + ) + + # Rain on secondary axis + p_rain = barplot!( + ax_H2O_rain, + datetime2unix.(inputs.DateTime[index_t_start:index_t_end]), + inputs.P[index_t_start:index_t_end], + color = :blue, + ) + + # Moisture stress + p_MS = CairoMakie.scatter!( + ax_MS, + datetime2unix.(climalsm.DateTime), + climalsm.β, + color = :green, + ) # not sure about units + + # Stomatal conductance + p_SC = CairoMakie.scatter!( + ax_SC, + datetime2unix.(climalsm.DateTime), + climalsm.gs, + color = :green, + ) + + # xticks + ax_H2O.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + ax_H2O_rain.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + hidexdecorations!(ax_H2O_rain) + ax_MS.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + ax_SC.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + + axislegend( + ax_H2O, + [p_H2O_d, p_H2O_m], + ["Observations", "ClimaLSM"], + "", + position = :rt, + orientation = :horizontal, + ) + + #CairoMakie.ylims!(ax_C, (0, 40)) + #CairoMakie.ylims!(ax_W, (0, 30)) + #CairoMakie.ylims!(ax_SWOUT, (0, 200)) + + [ + CairoMakie.xlims!( + axes, + ( + datetime2unix(climalsm.DateTime[1]), + datetime2unix(climalsm.DateTime[end]), + ), + ) for axes in [ax_H2O, ax_H2O_rain, ax_MS, ax_SC] + ] + + fig + return fig +end + +# 3. Fingerprint plot +function fingerprint_fig(inputs, climalsm, index_t_start, index_t_end) + fig = Figure(size = (1000, 1000)) + fontsize_theme = Theme(fontsize = 20) + set_theme!(fontsize_theme) + + # create empty axis, with a specific layout + ax_C = Axis( + fig[1, 1], + ylabel = "Hour of the day", + xlabel = "Date", + title = L"\text{GPP} \, (\mu\text{mol m}^{-2} \, \text{s}^{-1})", + ) # C fluxes + ax_W = Axis(fig[2, 1]) # h2o fluxes + ax_P = Axis(fig[3, 1]) # Precip & soil moisture + ax_T = Axis(fig[4, 1]) # air, canopy, and soil temperature + + # for time series, CairoMakie should allow DateTime type soon (but not yet) + # so the 2 lines of code below are a trick to be able to use DateTime - will be removed later + dateticks = + optimize_ticks(inputs.DateTime[1], inputs.DateTime[end])[1][2:(end - 1)] # first and last are weirdly placed + + # Fingerprint plot + hm_GPP = CairoMakie.heatmap!( + ax_C, + datetime2unix.(DateTime.(Date.(inputs.DateTime))), + hour.(inputs.DateTime) .+ (minute.(inputs.DateTime) ./ 60), + inputs.GPP .* 1e6, + ) + Colorbar(fig[1, 2], hm_GPP) + + ax_C.xticks[] = + (datetime2unix.(dateticks), Dates.format.(dateticks, "mm/dd")) + + fig + return fig +end + +function diurnal(datetime, data) + hourlyquantile = [ + quantile(data[hour.(datetime) .== h], [0.25, 0.5, 0.75]) for h in 0:1:23 + ] + return hourlyquantile +end + +function diurnal_plot!( + fig, + ax, + datetime, + data, + color; + alpha = 0.3, + linestyle = :solid, +) + fig + hourlyquantile = diurnal(datetime, data) + diurnal_p = CairoMakie.lines!( + ax, + 0.5:1:23.5, + getindex.(hourlyquantile[1:24], 2), + color = color, + linestyle = linestyle, + ) + diurnal_q = CairoMakie.band!( + ax, + 0.5:1:23.5, + getindex.(hourlyquantile[1:24], 1), + getindex.(hourlyquantile[1:24], 3), + color = (color, alpha), + ) + return diurnal_p +end + +function diurnals_fig(inputs, climalsm, index_t_start, index_t_end) + fig = Figure(size = (1000, 1000)) + fontsize_theme = Theme(fontsize = 20) + set_theme!(fontsize_theme) + + ax_C = Axis( + fig[1, 1], + ylabel = L"\text{CO}_{2} \, (\mu\text{mol m}^{-2} \, \text{s}^{-1})", + ) # C fluxes + ax_W = Axis(fig[2, 1], ylabel = L"\text{H}_{2}\text{O} \, \text{(mm)}") # h2o fluxes + ax_E = Axis( + fig[3, 1], + ylabel = L"\text{Radiation} \, (\text{W} \, \text{m}^{-2})", + xlabel = L"\text{Hour of the day}", + xgridvisible = false, + ) # shortwave out + + # CO2 fluxes + # model + p_GPP_m = + diurnal_plot!(fig, ax_C, climalsm.DateTime, climalsm.GPP .* 1e6, :green) + diurnal_plot!(fig, ax_C, climalsm.DateTime, climalsm.Ra .* 1e6, :black) + # data + p_GPP_d = diurnal_plot!( + fig, + ax_C, + inputs.DateTime[index_t_start:index_t_end], + inputs.GPP[index_t_start:index_t_end] .* 1e6, + :green, + alpha = 0.1, + linestyle = :dot, + ) + + # H2O fluxes + # model + p_ET_m = diurnal_plot!( + fig, + ax_W, + climalsm.DateTime, + climalsm.transpiration .* 1e3 .* 24 .* 3600, + :blue, + ) + # data + p_ET_d = diurnal_plot!( + fig, + ax_W, + inputs.DateTime[index_t_start:index_t_end], + inputs.LE[index_t_start:index_t_end] ./ + (LSMP.LH_v0(earth_param_set) * 1000) .* (1e3 * 24 * 3600), + :blue, + alpha = 0.1, + linestyle = :dot, + ) + + # Energy fluxes + # model + # diurnal_plot!(fig, ax_E, climalsm.DateTime, climalsm.LW_out, :red) + p_SWout_m = + diurnal_plot!(fig, ax_E, climalsm.DateTime, climalsm.SW_out, :red) + # data + p_SWout_d = diurnal_plot!( + fig, + ax_E, + inputs.DateTime[index_t_start:index_t_end], + FT.(inputs.SW_OUT[index_t_start:index_t_end]), + :red, + alpha = 0.1, + linestyle = :dot, + ) + + [CairoMakie.xlims!(axes, (0, 24)) for axes in [ax_C, ax_W, ax_E]] + + axislegend( + ax_C, + [p_GPP_d, p_GPP_m], + ["Observations", "ClimaLSM"], + "", + position = :rt, + orientation = :horizontal, + ) + axislegend( + ax_W, + [p_ET_d, p_ET_m], + ["Observations", "ClimaLSM"], + "", + position = :rt, + orientation = :horizontal, + ) + axislegend( + ax_E, + [p_SWout_d, p_SWout_m], + ["Observations", "ClimaLSM"], + "", + position = :rt, + orientation = :horizontal, + ) + + hidexdecorations!(ax_C) + hidexdecorations!(ax_W) + + fig + return fig +end + +# 5. Cumulative P and ET + +# 6. Energy balance closure (L + H = Rn - G) + +function make_plots(inputs, climalsm) + # below shouldn't be hardcoded! + index_t_start = 120 * 48 # we shouldn't hardcode that 120 in ozark_simulation.jl + index_t_end = 120 * 48 + (60 - 30) * 48 + + if isdir("figures") + nothing + else + mkdir("figures") + end + + fig1 = timeseries_fluxes_fig(inputs, climalsm, index_t_start, index_t_end) + fig2 = timeseries_H2O_fig(inputs, climalsm, index_t_start, index_t_end) + fig3 = fingerprint_fig(inputs, climalsm, index_t_start, index_t_end) + fig4 = diurnals_fig(inputs, climalsm, index_t_start, index_t_end) + + names = [ + "timeseries_fluxes.pdf", + "timeseries_H2O.pdf", + "fingerprint.pdf", + "diurnals.pdf", + ] + + [ + save(joinpath("figures", name), fig) for + (name, fig) in zip(names, [fig1, fig2, fig3, fig4]) + ] + return nothing +end diff --git a/lib/ClimaLandSimulations/src/utilities/met_drivers_FLUXNET.jl b/lib/ClimaLandSimulations/src/utilities/met_drivers_FLUXNET.jl new file mode 100644 index 0000000000..769afa02f5 --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/met_drivers_FLUXNET.jl @@ -0,0 +1,224 @@ +export setup_drivers + +function setup_drivers(site_ID) + + af = ArtifactFile( + url = data_link, + filename = "AMF_$(site_ID)_FLUXNET_FULLSET.csv", + ) + dataset = ArtifactWrapper( + "$ClimaLandSimulations_dir/src/fluxnet/$site_ID", + "ameriflux_data", + ArtifactFile[af], + ) + dataset_path = get_data_folder(dataset) + data = joinpath(dataset_path, "AMF_$(site_ID)_FLUXNET_FULLSET.csv") + driver_data = readdlm(data, ',') + + LOCAL_DATETIME = DateTime.(format.(driver_data[2:end, 1]), "yyyymmddHHMM") + UTC_DATETIME = LOCAL_DATETIME .+ Dates.Hour(time_offset) + DATA_DT = Second(LOCAL_DATETIME[2] - LOCAL_DATETIME[1]).value # seconds + + # Label of every data column to be collected from the driver data file + labels = ( + :TA, + :VPD, + :PA, + :P, + :WS, + :LW_IN, + :SW_IN, + :G, + :GPP, + :LE, + :H, + :LW_OUT, + :SW_OUT, + :SWC, + :TS, + :CO2, + :RECO, + ) + + # For every data column to be collected, name of the column in the file, + # transformation function to desired unit, and string of final unit + collect_args = [ + ("TA_F", (x) -> x .+ 273.15, "K") + ("VPD_F", (x) -> x .* 100, "Pa") + ("PA_F", (x) -> x .* 1000, "Pa") + ("P_F", (x) -> x ./ (1000 * DATA_DT), "m/s") + ("WS_F", (x) -> x, "m/s") + ("LW_IN_F", (x) -> x, "W/m^2") + ("SW_IN_F", (x) -> x, "W/m^2") + ("G_F_MDS", (x) -> x, "W/m^2") + ("GPP_DT_VUT_REF", (x) -> x .* 1e-6, "mol/m^2/s") + ("LE_CORR", (x) -> x, "W/m^2") + ("H_CORR", (x) -> x, "W/m^2") + ("LW_OUT", (x) -> x, "W/m^2") + ("SW_OUT", (x) -> x, "W/m^2") + ("SWC_F_MDS_1", (x) -> x ./ 100, "m^3/m^3") + ("TS_F_MDS_1", (x) -> x .+ 273.15, "K") + ("CO2_F_MDS", (x) -> x .* 1e-6, "mol") + ("RECO_DT_VUT_REF", (x) -> x .* 1e-6, "mol/m^-2/s") + ] + + # Named tuple mapping every label to a DataColumn with the correct transformed + # unit and data status. Automatically fills gaps in the data + drivers = (; + zip( + labels, + [ + transform_column( + filter_column(driver_data, name, ""), + transform, + unit, + ) for (name, transform, unit) in collect_args + ], + )... + ) + # Check that all required driver data is present, otherwise throw error + required = [ + "TA" => drivers.TA, + "VPD" => drivers.VPD, + "PA" => drivers.PA, + "P" => drivers.P, + "WS" => drivers.WS, + "LW_IN" => drivers.LW_IN, + "SW_IN" => drivers.SW_IN, + ] + missing_drivers = [] + for (name, column) in required + if column.status == absent + push!(missing_drivers, name) + end + end + if length(missing_drivers) != 0 + error( + "Driver data missing for columns: $([missing_drivers[i] * " " for + i in 1:length(missing_drivers)]...)", + ) + end + + thermo_params = LSMP.thermodynamic_parameters(earth_param_set) + esat = + Thermodynamics.saturation_vapor_pressure.( + Ref(thermo_params), + drivers.TA.values, + Ref(Thermodynamics.Liquid()), + ) + e = @. esat - drivers.VPD.values + q = @. 0.622 * e ./ (drivers.PA.values - 0.378 * e) + + # Create splines for each atmospheric driver needed for PrescribedAtmosphere + # and for PrescribedRadiation + seconds = FT.(0:DATA_DT:((length(UTC_DATETIME) - 1) * DATA_DT)) + p_spline = Spline1D(seconds, -drivers.P.values[:]) # m/s + atmos_q = Spline1D(seconds, q[:]) + atmos_T = Spline1D(seconds, drivers.TA.values[:]) + atmos_p = Spline1D(seconds, drivers.PA.values[:]) + atmos_co2 = Spline1D(seconds, drivers.CO2.values[:]) + atmos_u = Spline1D(seconds, drivers.WS.values[:]) + LW_IN_spline = Spline1D(seconds, drivers.LW_IN.values[:]) + SW_IN_spline = Spline1D(seconds, drivers.SW_IN.values[:]) + precipitation_function(t::FT) where {FT} = + p_spline(t) < 0.0 ? p_spline(t) : 0.0 # m/s + snow_precip(t) = FT(0) # this is likely not correct + + # Construct the drivers. For the reference time we will use the UTC time at the + # start of the simulation + atmos = ClimaLSM.PrescribedAtmosphere( + precipitation_function, + snow_precip, + atmos_T, + atmos_u, + atmos_q, + atmos_p, + UTC_DATETIME[1], + atmos_h; + c_co2 = atmos_co2, + ) + + function zenith_angle( + t, + ref_time; + latitude = lat, + longitude = long, + insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params, + ) where {FT} + # This should be time in UTC + current_datetime = ref_time + Dates.Second(round(t)) + + # Orbital Data uses Float64, so we need to convert to our sim FT + d, δ, η_UTC = + FT.( + Insolation.helper_instantaneous_zenith_angle( + current_datetime, + ref_time, + insol_params, + ) + ) + + FT( + Insolation.instantaneous_zenith_angle( + d, + δ, + η_UTC, + longitude, + latitude, + )[1], + ) + end + + radiation = ClimaLSM.PrescribedRadiativeFluxes( + FT, + SW_IN_spline, + LW_IN_spline, + UTC_DATETIME[1]; + θs = zenith_angle, + ) + + # Start and end dates of data in MODIS format + modis_start_date = "A$(Dates.year(UTC_DATETIME[1]))$(lpad(Dates.dayofyear(UTC_DATETIME[1]), 3, "0"))" + modis_end_date = "A$(Dates.year(UTC_DATETIME[end]))$(lpad(Dates.dayofyear(UTC_DATETIME[end]), 3, "0"))" + + MODIS_LAI = single_col_data_matrix( + parse_response( + check_response( + send_get_subset( + "MCD15A2H", + modis_start_date, + modis_end_date, + site_ID, + band = "Lai_500m", + ), + ), + ), + ) + + LAI_dt = Second(MODIS_LAI[2, 1] - MODIS_LAI[1, 1]).value + LAI_seconds = FT.(0:LAI_dt:((length(MODIS_LAI[:, 1]) - 1) * LAI_dt)) + + # LAI spline for radiative transfer + LAIspline = Spline1D(LAI_seconds, MODIS_LAI[:, 2]) + LAIfunction = (t) -> eltype(t)(LAIspline(t)) + + # Necessary inputs from LAI Data + # Note that f_root_to_shoot, capacity, and h_leaf are site-specific parameters + # defined in the parameters file for each site + maxLAI = FT(maximum(MODIS_LAI[:, 2])) + RAI = maxLAI * f_root_to_shoot + plant_ν = capacity / (maxLAI * h_leaf) / FT(1000) + + return ( + LOCAL_DATETIME, + atmos_co2, + DATA_DT, + drivers, + atmos, + radiation, + LAIfunction, + maxLAI, + RAI, + plant_ν, + ) +end diff --git a/lib/ClimaLandSimulations/src/utilities/pull_MODIS.jl b/lib/ClimaLandSimulations/src/utilities/pull_MODIS.jl new file mode 100644 index 0000000000..111c311260 --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/pull_MODIS.jl @@ -0,0 +1,57 @@ +export send_get_subset, check_response, parse_response, single_col_data_matrix + +""" This file provides methods for interacting with the MODIS REST API at a site + in order to downloads specific subsets of MODIS products at a site level.""" + +function send_get_subset( + product::String, + start_date::String, + end_date::String, + site::String; + network::String = "AMERIFLUX", + band::String = "", +) + """Sends a GET request to the MODIS REST API for a specific product, time, + and site. Returns an HTTP.messages.response object.""" + url = "https://modis.ornl.gov/rst/api/v1/$product/$network/$site/subset?" + band = band == "" ? band : "band=$band" + start_date = "&startDate=$start_date" + end_date = "&endDate=$end_date" + url *= "$band$start_date$end_date" + return HTTP.get(url) +end + +function check_response(response::HTTP.Messages.Response) + """Checks the response from the MODIS REST API for errors. Returns the + response if no errors are found, otherwise throws an error.""" + if response.status != 200 + throw(ErrorException("Error in MODIS pull: $(response.status)")) + end + return response +end + +function parse_response(response::HTTP.Messages.Response) + """Parses the response from the MODIS REST API into a nested Julia + dictionary using the JSON module and returns the dictionary.""" + return JSON.parse(String(response.body))["subset"] +end + +function single_col_data_matrix(JSON_data::Vector) + """Takes in a vector of JSON data dicitonaries of a single column of chunked + MODIS data and assembles it into a data matrix with a time column and a + single data column. The data column gives the average of the grid cell + arround the site for each day, and the time column gives the Dates.DateTime + object for each day. In averaging, MODIS data marked as missing (>100) is + ignored.""" + cleansed_dat = [ + [ + chunk["data"][i] < 100 ? chunk["data"][i] : missing for + i in 1:length(chunk["data"]) + ] for chunk in JSON_data + ] + data_col = + [mean(filter(!ismissing, day_dat)) / 10 for day_dat in cleansed_dat] + time_col = + [DateTime(JSON_data[i]["calendar_date"]) for i in 1:length(JSON_data)] + return hcat(time_col, data_col) +end diff --git a/lib/ClimaLandSimulations/src/utilities/timestepper_setup.jl b/lib/ClimaLandSimulations/src/utilities/timestepper_setup.jl new file mode 100644 index 0000000000..398b82241d --- /dev/null +++ b/lib/ClimaLandSimulations/src/utilities/timestepper_setup.jl @@ -0,0 +1,7 @@ +N_spinup_days = 30 +N_days = N_spinup_days + 30 +tf = Float64(t0 + 3600 * 24 * N_days) +t_spinup = Float64(t0 + N_spinup_days * 3600 * 24) +saveat = Array(t_spinup:(n * dt):tf) +timestepper = CTS.RK4() +ode_algo = CTS.ExplicitAlgorithm(timestepper) diff --git a/lib/ClimaLandSimulations/test/runtests.jl b/lib/ClimaLandSimulations/test/runtests.jl new file mode 100644 index 0000000000..66ea8635e4 --- /dev/null +++ b/lib/ClimaLandSimulations/test/runtests.jl @@ -0,0 +1,3 @@ +using ClimaLSM # rename to ClimaLand WIP +using ClimaLandSite +using Test diff --git a/lib/README.md b/lib/README.md new file mode 100644 index 0000000000..d1a82e6d3a --- /dev/null +++ b/lib/README.md @@ -0,0 +1,4 @@ +# Packages that extend ClimaLand functionality + +- ClimaLandSimulations: ClimaLand single site or global simulations +- ClimaLandDashboards: ClimaLand web dashboards