diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml
index da552fa2f8..ff971690de 100644
--- a/.buildkite/Manifest.toml
+++ b/.buildkite/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "dae5b8a0ef6cad3a789e38991e9bafc7e71f7ab5"
+project_hash = "8f674e311ce137cbeb26c8e00409aca502afec54"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "fa0822e5baee6e23081c2685ae27265dabee23d8"
@@ -231,12 +231,6 @@ git-tree-sha1 = "9a9610fbe5779636f75229e423e367124034af41"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 version = "0.16.43"
 
-[[deps.Blosc_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "19b98ee7e3db3b4eff74c5c9c72bf32144e24f10"
-uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
-version = "1.21.5+0"
-
 [[deps.Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd"
@@ -357,6 +351,20 @@ weakdeps = ["SparseArrays"]
     [deps.ChainRulesCore.extensions]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
+[[deps.ClimaAnalysis]]
+deps = ["NCDatasets", "OrderedCollections", "Statistics"]
+git-tree-sha1 = "af5e012e13bc9609f712dea0434cb4ce2b3e709e"
+uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
+version = "0.5.2"
+
+    [deps.ClimaAnalysis.extensions]
+    CairoMakieExt = "CairoMakie"
+    GeoMakieExt = "GeoMakie"
+
+    [deps.ClimaAnalysis.weakdeps]
+    CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+    GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
+
 [[deps.ClimaComms]]
 git-tree-sha1 = "d76590e99fa942e07f1e992a4d4a5e25121a26d6"
 uuid = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
@@ -378,8 +386,14 @@ weakdeps = ["CUDA", "Krylov"]
     ClimaCoreCUDAExt = "CUDA"
     KrylovExt = "Krylov"
 
+[[deps.ClimaDiagnostics]]
+deps = ["Accessors", "ClimaComms", "ClimaCore", "Dates", "NCDatasets", "SciMLBase"]
+git-tree-sha1 = "aff194804df0fcfcf69a80c58978b84777272619"
+uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
+version = "0.2.2"
+
 [[deps.ClimaLand]]
-deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
+deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.12.4"
@@ -830,10 +844,10 @@ uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
 version = "0.4.1"
 
 [[deps.FFMPEG_jll]]
-deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e"
+deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Pkg", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
+git-tree-sha1 = "74faea50c1d007c85837327f6775bea60b5492dd"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "4.4.4+1"
+version = "4.4.2+2"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
@@ -1033,11 +1047,6 @@ git-tree-sha1 = "ff38ba61beff76b8f4acad8ab0c97ef73bb670cb"
 uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
 version = "3.3.9+0"
 
-[[deps.GMP_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
-version = "6.2.1+6"
-
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]
 git-tree-sha1 = "c154546e322a9c73364e8a60430b0f79b812d320"
@@ -1103,12 +1112,6 @@ git-tree-sha1 = "7c82e6a6cd34e9d935e9aa4051b66c6ff3af59ba"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
 version = "2.80.2+0"
 
-[[deps.GnuTLS_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Nettle_jll", "P11Kit_jll", "Zlib_jll"]
-git-tree-sha1 = "383db7d3f900f4c1f47a8a04115b053c095e48d3"
-uuid = "0951126a-58fd-58f1-b5b3-b08c7c4a876d"
-version = "3.8.4+0"
-
 [[deps.Graphics]]
 deps = ["Colors", "LinearAlgebra", "NaNMath"]
 git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd"
@@ -1129,9 +1132,9 @@ version = "1.11.1"
 
 [[deps.GridLayoutBase]]
 deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
-git-tree-sha1 = "fc713f007cff99ff9e50accba6373624ddd33588"
+git-tree-sha1 = "6f93a83ca11346771a93bbde2bdad2f65b61498f"
 uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
-version = "0.11.0"
+version = "0.10.2"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -1149,10 +1152,10 @@ weakdeps = ["MPI"]
     MPIExt = "MPI"
 
 [[deps.HDF5_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
-git-tree-sha1 = "82a471768b513dc39e471540fdadc84ff80ff997"
+deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "4cc2bb72df6ff40b055295fdef6d92955f9dede8"
 uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
-version = "1.14.3+3"
+version = "1.12.2+2"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
@@ -1778,9 +1781,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[deps.MathTeXEngine]]
 deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"]
-git-tree-sha1 = "1865d0b8a2d91477c8b16b49152a32764c7b1f5f"
+git-tree-sha1 = "96ca8a313eb6437db5ffe946c457a401bbb8ce1d"
 uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
-version = "0.6.0"
+version = "0.5.7"
 
 [[deps.MatrixFactorizations]]
 deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"]
@@ -1916,10 +1919,10 @@ uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391"
 version = "0.1.5"
 
 [[deps.NetCDF_jll]]
-deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML", "XML2_jll", "Zlib_jll", "Zstd_jll", "libzip_jll"]
-git-tree-sha1 = "4686378c4ae1d1948cfbe46c002a11a4265dcb07"
+deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg", "XML2_jll", "Zlib_jll"]
+git-tree-sha1 = "072f8371f74c3b9e1b26679de7fbf059d45ea221"
 uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3"
-version = "400.902.211+1"
+version = "400.902.5+1"
 
 [[deps.Netpbm]]
 deps = ["FileIO", "ImageCore", "ImageMetadata"]
@@ -1927,12 +1930,6 @@ git-tree-sha1 = "d92b107dbb887293622df7697a2223f9f8176fcd"
 uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
 version = "1.1.1"
 
-[[deps.Nettle_jll]]
-deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
-uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
-version = "3.7.2+0"
-
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
@@ -2058,12 +2055,6 @@ git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.6.3"
 
-[[deps.P11Kit_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "2cd396108e178f3ae8dedbd8e938a18726ab2fbf"
-uuid = "c2071276-7c44-58a7-b746-946036e04d0a"
-version = "0.24.1+0"
-
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
@@ -3215,12 +3206,6 @@ git-tree-sha1 = "51b5eeb3f98367157a7a12a1fb0aa5328946c03c"
 uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4"
 version = "0.2.3+0"
 
-[[deps.libaec_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "46bf7be2917b59b761247be3f317ddf75e50e997"
-uuid = "477f73a3-ac25-53e9-8cc3-50b2fa2566f0"
-version = "1.1.2+0"
-
 [[deps.libaom_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
 git-tree-sha1 = "1827acba325fdcdf1d2647fc8d5301dd9ba43a9d"
@@ -3274,12 +3259,6 @@ git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
 uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
 version = "1.3.7+1"
 
-[[deps.libzip_jll]]
-deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "3282b7d16ae7ac3e57ec2f3fa8fafb564d8f9f7f"
-uuid = "337d8026-41b4-5cde-a456-74a10e5b31d1"
-version = "1.10.1+0"
-
 [[deps.mtdev_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "814e154bdb7be91d78b6802843f76b6ece642f11"
diff --git a/.buildkite/Project.toml b/.buildkite/Project.toml
index 04bc1a43be..3cee10ee64 100644
--- a/.buildkite/Project.toml
+++ b/.buildkite/Project.toml
@@ -5,8 +5,10 @@ BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
+ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 ClimaLand = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index efee397207..0e0497ffce 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -137,7 +137,9 @@ steps:
       - label: "Global Bucket on CPU (functional albedo)"
         key: "global_bucket_function_cpu"
         command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/global_bucket_function.jl"
-        artifact_paths: "experiments/standalone/Bucket/artifacts/*cpu*"
+        artifact_paths:
+          - "experiments/standalone/Bucket/artifacts/*cpu*"
+          - "experiments/global_bucket_function/output_active/*.png"
 
       - label: "Global Bucket on CPU (static map albedo)"
         key: "global_bucket_staticmap_cpu"
@@ -196,20 +198,12 @@ steps:
 
   - group: "CPU/GPU comparisons"
     steps:
-      - label: "Compare GPU bucket with CPU bucket (functional albedo)"
+      - label: "Compare GPU bucket with CPU bucket"
         command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/compare_gpu_cpu_output.jl"
         depends_on:
           - "global_bucket_function_cpu"
           - "global_bucket_function_gpu"
-
-      - label: "Compare GPU bucket with CPU bucket (static map albedo)"
-        command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/compare_gpu_cpu_output.jl"
-        depends_on:
           - "global_bucket_staticmap_cpu"
           - "global_bucket_staticmap_gpu"
-
-      - label: "Compare GPU bucket with CPU bucket (temporal map albedo)"
-        command: "julia --color=yes --project=.buildkite experiments/standalone/Bucket/compare_gpu_cpu_output.jl"
-        depends_on:
           - "global_bucket_temporalmap_cpu"
           - "global_bucket_temporalmap_gpu"
diff --git a/.github/workflows/ClimaLandSimulations.yml b/.github/workflows/ClimaLandSimulations.yml
index 76cdb84487..5ae09e5b46 100644
--- a/.github/workflows/ClimaLandSimulations.yml
+++ b/.github/workflows/ClimaLandSimulations.yml
@@ -19,13 +19,14 @@ jobs:
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1.10'
-      - uses: julia-actions/cache@v2
+#  julia-actions/cache@v2 did not seem to invalidate cache correctly...
+#    - uses: julia-actions/cache@v2
       - name: Install Julia dependencies
         run: >
-          julia --project= -e 'using Pkg; Pkg.develop(path="$(pwd())"); Pkg.develop(path="$(pwd())/lib/ClimaLandSimulations")'
+          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= -e 'using Pkg; Pkg.test("ClimaLandSimulations")'
+          julia --project=lib/ClimaLandSimulations -e 'using Pkg; Pkg.test("ClimaLandSimulations")'
diff --git a/Artifacts.toml b/Artifacts.toml
index fdabc938f9..0e9a3ff72d 100644
--- a/Artifacts.toml
+++ b/Artifacts.toml
@@ -1,9 +1,9 @@
+[era5_land_forcing_data2021]
+git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a"
+
 [soil_params_Gupta2020_2022]
 git-tree-sha1 = "8e28b4274b10020b6cdd54b8e7585221379d9d33"
 
     [[soil_params_Gupta2020_2022.download]]
     sha256 = "97dcf1158cba03b1fd397262bdfaf85a523f57038c337bcce163e32664d3616b"
     url = "https://caltech.box.com/shared/static/f2y23qx0lggjskftzgh7ht7fsbh36gmm.gz"
-
-[era5_land_forcing_data2021]
-git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a"
diff --git a/Project.toml b/Project.toml
index e990493153..eba8a74e3f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
+ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -44,6 +45,7 @@ CSV = "0.10"
 CUDA = "5.3"
 ClimaComms = "0.5.6, 0.6"
 ClimaCore = "0.13.2, 0.14"
+ClimaDiagnostics = "0.2"
 ClimaParams = "0.10.2"
 ClimaUtilities = "0.1.2"
 DataFrames = "1"
diff --git a/README.md b/README.md
index b4325e0912..3125c38bf3 100644
--- a/README.md
+++ b/README.md
@@ -8,6 +8,7 @@
 component) or standalone (single component) modes.  </strong>
 </p>
 
+[![Downloads](https://img.shields.io/badge/dynamic/json?url=http%3A%2F%2Fjuliapkgstats.com%2Fapi%2Fv1%2Ftotal_downloads%2FClimaLand&query=total_requests&suffix=%2Ftotal&label=Downloads)](http://juliapkgstats.com/pkg/ClimaLand)
 
 Recommended Julia Version: Stable release v1.10.0. CI no longer tests earlier
 versions of Julia.
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 738768a9bf..3c9ec87cc5 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "bf43a7452b446a9ba967dab8eeb0bc21ba10184f"
+project_hash = "6c3f79ccae513840cc13b9ac947575c2305b507f"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "fa0822e5baee6e23081c2685ae27265dabee23d8"
@@ -375,8 +375,14 @@ weakdeps = ["CUDA", "Krylov"]
     ClimaCoreCUDAExt = "CUDA"
     KrylovExt = "Krylov"
 
+[[deps.ClimaDiagnostics]]
+deps = ["Accessors", "ClimaComms", "ClimaCore", "Dates", "NCDatasets", "SciMLBase"]
+git-tree-sha1 = "aff194804df0fcfcf69a80c58978b84777272619"
+uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
+version = "0.2.2"
+
 [[deps.ClimaLand]]
-deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
+deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 version = "0.12.4"
diff --git a/docs/Project.toml b/docs/Project.toml
index fcacd6e6a2..e7f6f974ed 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -20,6 +20,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
diff --git a/docs/list_diagnostics.jl b/docs/list_diagnostics.jl
new file mode 100644
index 0000000000..4d76cb4f6f
--- /dev/null
+++ b/docs/list_diagnostics.jl
@@ -0,0 +1,5 @@
+diagnostics = [
+    "Available diagnostics" => "diagnostics/available_diagnostics.md",
+    "For developers" => "diagnostics/developers_diagnostics.md",
+    "For users" => "diagnostics/users_diagnostics.md",
+]
diff --git a/docs/make.jl b/docs/make.jl
index f1355b2e21..985750b47b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -84,6 +84,7 @@ ext_jl2md(x) = joinpath(basename(GENERATED_DIR), replace(x, ".jl" => ".md"))
 tutorials = transform_second(x -> ext_jl2md(x), tutorials)
 include("list_of_apis.jl")
 include("list_standalone_models.jl")
+include("list_diagnostics.jl")
 pages = Any[
     "Home" => "index.md",
     "APIs" => apis,
@@ -91,6 +92,7 @@ pages = Any[
     "Tutorials" => tutorials,
     "Repository structure" => "folderstructure.md",
     "Standalone models" => standalone_models,
+    "Diagnostics" => diagnostics,
 ]
 
 
diff --git a/docs/src/diagnostics/available_diagnostics.md b/docs/src/diagnostics/available_diagnostics.md
new file mode 100644
index 0000000000..3642907644
--- /dev/null
+++ b/docs/src/diagnostics/available_diagnostics.md
@@ -0,0 +1,6 @@
+# Available diagnostic variables
+
+Autogenerate table of available diagnostics:
+```@example
+include("make_diagnostic_table.jl")
+```
diff --git a/docs/src/diagnostics/developers_diagnostics.md b/docs/src/diagnostics/developers_diagnostics.md
new file mode 100644
index 0000000000..429973d632
--- /dev/null
+++ b/docs/src/diagnostics/developers_diagnostics.md
@@ -0,0 +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.
diff --git a/docs/src/diagnostics/make_diagnostic_table.jl b/docs/src/diagnostics/make_diagnostic_table.jl
new file mode 100644
index 0000000000..e7b7cedab6
--- /dev/null
+++ b/docs/src/diagnostics/make_diagnostic_table.jl
@@ -0,0 +1,28 @@
+import ClimaLand as CL
+using PrettyTables
+
+# Print all available diagnostics to an ASCII table
+
+CL.Diagnostics.define_diagnostics!(nothing)
+short_names = []
+long_names = []
+units = []
+comments = []
+standard_names = []
+for d in values(CL.Diagnostics.ALL_DIAGNOSTICS)
+    push!(short_names, d.short_name)
+    push!(long_names, d.long_name)
+    push!(units, d.units)
+    push!(comments, d.comments)
+    push!(standard_names, d.standard_name)
+end
+data = hcat(short_names, long_names, units, comments, standard_names)
+pretty_table(
+    data;
+    autowrap = true,
+    linebreaks = true,
+    columns_width = [10, 15, 8, 32, 15],  # Width = 80
+    body_hlines = collect(1:size(data)[1]),
+    header = ["Short name", "Long name", "Units", "Comments", "Standard name"],
+    alignment = :l,
+)
diff --git a/docs/src/diagnostics/users_diagnostics.md b/docs/src/diagnostics/users_diagnostics.md
new file mode 100644
index 0000000000..0c9728b579
--- /dev/null
+++ b/docs/src/diagnostics/users_diagnostics.md
@@ -0,0 +1,119 @@
+# 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.CLD.default_diagnostics(model, 1.0; 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.
+
+# 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/experiments/Manifest.toml b/experiments/Manifest.toml
index c977b0bf30..7524dfe66a 100644
--- a/experiments/Manifest.toml
+++ b/experiments/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.4"
 manifest_format = "2.0"
-project_hash = "973443a0a609649a9289253b6f5bc3de181d9006"
+project_hash = "306f35a5dfb80400ec575833e6ee1eaa425bdc2f"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "daf26bbdec60d9ca1c0003b70f389d821ddb4224"
@@ -256,6 +256,20 @@ weakdeps = ["SparseArrays"]
     [deps.ChainRulesCore.extensions]
     ChainRulesCoreSparseArraysExt = "SparseArrays"
 
+[[deps.ClimaAnalysis]]
+deps = ["NCDatasets", "OrderedCollections", "Statistics"]
+git-tree-sha1 = "c2e1c0d5c30a2519a4282988037b255dbc9aee00"
+uuid = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
+version = "0.5.3"
+
+    [deps.ClimaAnalysis.extensions]
+    CairoMakieExt = "CairoMakie"
+    GeoMakieExt = "GeoMakie"
+
+    [deps.ClimaAnalysis.weakdeps]
+    CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+    GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
+
 [[deps.ClimaComms]]
 git-tree-sha1 = "f19a9d42ef27affa8d57a5702dcb8a318077aa86"
 uuid = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
@@ -283,11 +297,17 @@ version = "0.14.5"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 
+[[deps.ClimaDiagnostics]]
+deps = ["Accessors", "ClimaComms", "ClimaCore", "Dates", "NCDatasets", "SciMLBase"]
+git-tree-sha1 = "aff194804df0fcfcf69a80c58978b84777272619"
+uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
+version = "0.2.2"
+
 [[deps.ClimaLand]]
-deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
+deps = ["Adapt", "ArtifactWrappers", "ClimaComms", "ClimaCore", "ClimaDiagnostics", "ClimaUtilities", "DataFrames", "Dates", "DocStringExtensions", "Insolation", "Interpolations", "IntervalSets", "LazyArtifacts", "LinearAlgebra", "NCDatasets", "SciMLBase", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnrolledUtilities"]
 path = ".."
 uuid = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
-version = "0.12.2"
+version = "0.12.3"
 
     [deps.ClimaLand.extensions]
     CreateParametersExt = "ClimaParams"
diff --git a/experiments/Project.toml b/experiments/Project.toml
index cbf600f9d1..6ccb4f3696 100644
--- a/experiments/Project.toml
+++ b/experiments/Project.toml
@@ -1,7 +1,9 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
+ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 ClimaLand = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
diff --git a/experiments/standalone/Bucket/compare_gpu_cpu_output.jl b/experiments/standalone/Bucket/compare_gpu_cpu_output.jl
index f51b89a679..c49162124d 100644
--- a/experiments/standalone/Bucket/compare_gpu_cpu_output.jl
+++ b/experiments/standalone/Bucket/compare_gpu_cpu_output.jl
@@ -1,10 +1,19 @@
 using DelimitedFiles
 using Statistics
 import ClimaLand
-outdir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Bucket/artifacts")
-cpu_state = readdlm(joinpath(outdir, "tf_state_cpu.txt"), ',')
-gpu_state = readdlm(joinpath(outdir, "tf_state_gpu.txt"), ',')
-@show abs(maximum(cpu_state .- gpu_state))
-@show abs(median(cpu_state .- gpu_state))
-@show abs(mean(cpu_state .- gpu_state))
-@assert isapprox(cpu_state, gpu_state)
+function check(job)
+    outdir = joinpath(
+        pkgdir(ClimaLand),
+        "experiments/standalone/Bucket/artifacts_$job",
+    )
+    cpu_state = readdlm(joinpath(outdir, "tf_state_cpu_$job.txt"), ',')
+    gpu_state = readdlm(joinpath(outdir, "tf_state_gpu_$job.txt"), ',')
+    @show abs(maximum(cpu_state .- gpu_state))
+    @show abs(median(cpu_state .- gpu_state))
+    @show abs(mean(cpu_state .- gpu_state))
+    @assert isapprox(cpu_state, gpu_state)
+end
+
+check("function")
+check("staticmap")
+check("temporalmap")
diff --git a/experiments/standalone/Bucket/global_bucket_function.jl b/experiments/standalone/Bucket/global_bucket_function.jl
index fae2f9111b..45e03aaa1c 100644
--- a/experiments/standalone/Bucket/global_bucket_function.jl
+++ b/experiments/standalone/Bucket/global_bucket_function.jl
@@ -41,6 +41,11 @@ using ClimaLand:
     PrescribedAtmosphere,
     PrescribedRadiativeFluxes
 
+using ClimaDiagnostics
+using ClimaAnalysis
+import ClimaAnalysis.Visualize as viz
+using ClimaUtilities
+
 """
    compute_extrema(v)
 
@@ -58,7 +63,10 @@ anim_plots = false
 FT = Float64;
 context = ClimaComms.context()
 earth_param_set = LP.LandParameters(FT);
-outdir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Bucket/artifacts")
+outdir = joinpath(
+    pkgdir(ClimaLand),
+    "experiments/standalone/Bucket/artifacts_function",
+)
 !ispath(outdir) && mkpath(outdir)
 
 # Construct simulation domain
@@ -153,24 +161,55 @@ prob = SciMLBase.ODEProblem(
     p,
 );
 
-saveat = collect(t0:(Δt * 3):tf);
-saved_values = (;
-    t = Array{Float64}(undef, length(saveat)),
-    saveval = Array{NamedTuple}(undef, length(saveat)),
-);
-saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat);
-updateat = copy(saveat)
-updatefunc = ClimaLand.make_update_drivers(bucket_atmos, bucket_rad)
-driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
-cb = SciMLBase.CallbackSet(driver_cb, saving_cb)
+# ClimaDiagnostics
+base_output_dir = "global_bucket_function/"
+output_dir =
+    ClimaUtilities.OutputPathGenerator.generate_output_path(base_output_dir)
+
+space = bucket_domain.space.subsurface
+
+nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir)
+
+diags = ClimaLand.CLD.default_diagnostics(model, t0; output_writer = nc_writer)
+
+diagnostic_handler =
+    ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt)
+
+diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler)
 
 sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve(
     prob,
     ode_algo;
     dt = Δt,
-    saveat = saveat,
-    callback = cb,
-);
+    callback = diag_cb,
+)
+
+#### ClimaAnalysis ####
+
+# all
+simdir = ClimaAnalysis.SimDir(output_dir)
+short_names_2D = [
+    "alpha",
+    "rn",
+    "tsfc",
+    "qsfc",
+    "lhf",
+    "rae",
+    "shf",
+    "vflux",
+    "rhosfc",
+    "wsoil",
+    "wsfc",
+    "ssfc",
+]
+short_names_3D = ["tsoil"]
+for short_name in vcat(short_names_2D..., short_names_3D...)
+    var = get(simdir; short_name)
+    fig = CairoMakie.Figure(size = (800, 600))
+    kwargs = short_name in short_names_2D ? Dict() : Dict(:z => 1)
+    viz.plot!(fig, var, lat = 0; kwargs...)
+    CairoMakie.save(joinpath(output_dir, "$short_name.png"), fig)
+end
 
 # Interpolate to grid
 space = axes(coords.surface)
@@ -179,92 +218,17 @@ latpts = range(-90.0, 90.0, 21)
 hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
 remapper = Remapping.Remapper(space, hcoords)
 
-W = [
-    Array(Remapping.interpolate(remapper, sol.u[k].bucket.W)) for
-    k in 1:length(sol.t)
-];
-Ws = [
-    Array(Remapping.interpolate(remapper, sol.u[k].bucket.Ws)) for
-    k in 1:length(sol.t)
-];
-σS = [
-    Array(Remapping.interpolate(remapper, sol.u[k].bucket.σS)) for
-    k in 1:length(sol.t)
-];
-T_sfc = [
-    Array(
-        Remapping.interpolate(remapper, saved_values.saveval[k].bucket.T_sfc),
-    ) for k in 1:length(sol.t)
-];
-evaporation = [
-    Array(
-        Remapping.interpolate(
-            remapper,
-            saved_values.saveval[k].bucket.turbulent_fluxes.vapor_flux,
-        ),
-    ) for k in 1:length(sol.t)
-];
-F_sfc = [
-    Array(
-        Remapping.interpolate(
-            remapper,
-            saved_values.saveval[k].bucket.R_n .+
-            saved_values.saveval[k].bucket.turbulent_fluxes.lhf .+
-            saved_values.saveval[k].bucket.turbulent_fluxes.shf,
-        ),
-    ) for k in 1:length(sol.t)
-];
+W = Array(Remapping.interpolate(remapper, sol.u[end].bucket.W))
+Ws = Array(Remapping.interpolate(remapper, sol.u[end].bucket.Ws))
+σS = Array(Remapping.interpolate(remapper, sol.u[end].bucket.σS))
+T_sfc = Array(Remapping.interpolate(remapper, prob.p.bucket.T_sfc))
+
 
 # save prognostic state to CSV - for comparison between
 # GPU and CPU output
 device_suffix =
     typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ?
     "cpu" : "gpu"
-open(joinpath(outdir, "tf_state_$device_suffix.txt"), "w") do io
-    writedlm(io, hcat(T_sfc[end][:], W[end][:], Ws[end][:], σS[end][:]), ',')
+open(joinpath(outdir, "tf_state_$(device_suffix)_function.txt"), "w") do io
+    writedlm(io, hcat(T_sfc[:], W[:], Ws[:], σS[:]), ',')
 end;
-# animation settings
-nframes = length(W)
-framerate = 2
-fig_ts = Figure(size = (1000, 1000))
-for (i, (field_ts, field_name)) in enumerate(
-    zip(
-        [W, σS, T_sfc, evaporation, F_sfc],
-        ["W", "σS", "T_sfc", "evaporation", "F_sfc"],
-    ),
-)
-    if anim_plots
-        fig = Figure(size = (1000, 1000))
-        ax = Axis(
-            fig[1, 1],
-            xlabel = "Longitude",
-            ylabel = "Latitude",
-            title = field_name,
-        )
-        clims = compute_extrema(field_ts)
-        CairoMakie.Colorbar(fig[:, end + 1], colorrange = clims)
-        outfile = joinpath(
-            outdir,
-            string("anim_$(device_suffix)_", field_name, ".mp4"),
-        )
-        record(fig, outfile, 1:nframes; framerate = framerate) do frame
-            CairoMakie.heatmap!(
-                longpts,
-                latpts,
-                field_ts[frame],
-                colorrange = clims,
-            )
-        end
-    end
-    # Plot the timeseries of the mean value as well.
-    xlabel = i == 5 ? "Time (days)" : ""
-    ax2 = Axis(
-        fig_ts[i, 1],
-        xlabel = xlabel,
-        ylabel = field_name,
-        title = "Global bucket with analytic albedo",
-    )
-    CairoMakie.lines!(ax2, sol.t ./ 3600 ./ 24, [mean(x) for x in field_ts])
-end
-outfile = joinpath(outdir, string("ts_$device_suffix.png"))
-CairoMakie.save(outfile, fig_ts)
diff --git a/experiments/standalone/Bucket/global_bucket_staticmap.jl b/experiments/standalone/Bucket/global_bucket_staticmap.jl
index cbf07ace4c..71fed1d816 100644
--- a/experiments/standalone/Bucket/global_bucket_staticmap.jl
+++ b/experiments/standalone/Bucket/global_bucket_staticmap.jl
@@ -305,7 +305,7 @@ sw_forcing = [
 ];
 
 # save prognostic state to CSV - for comparison between GPU and CPU output
-open(joinpath(outdir, "tf_state_$device_suffix.txt"), "w") do io
+open(joinpath(outdir, "tf_state_$(device_suffix)_staticmap.txt"), "w") do io
     writedlm(io, hcat(T_sfc[end][:], W[end][:], Ws[end][:], σS[end][:]), ',')
 end;
 # animation settings
diff --git a/experiments/standalone/Bucket/global_bucket_temporalmap.jl b/experiments/standalone/Bucket/global_bucket_temporalmap.jl
index 05cee96a01..b2cab6d207 100644
--- a/experiments/standalone/Bucket/global_bucket_temporalmap.jl
+++ b/experiments/standalone/Bucket/global_bucket_temporalmap.jl
@@ -277,7 +277,7 @@ F_sfc = [
 ];
 
 # save prognostic state to CSV - for comparison between GPU and CPU output
-open(joinpath(outdir, "tf_state_$device_suffix.txt"), "w") do io
+open(joinpath(outdir, "tf_state_$(device_suffix)_temporalmap.txt"), "w") do io
     writedlm(io, hcat(T_sfc[end][:], W[end][:], Ws[end][:], σS[end][:]), ',')
 end;
 # animation settings
diff --git a/lib/ClimaLandSimulations/Project.toml b/lib/ClimaLandSimulations/Project.toml
index 812b3fe89b..225e6abd77 100644
--- a/lib/ClimaLandSimulations/Project.toml
+++ b/lib/ClimaLandSimulations/Project.toml
@@ -9,6 +9,7 @@ Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
+ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 ClimaLand = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
@@ -42,9 +43,12 @@ cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd"
 ArtifactWrappers = "0.2"
 Bonito = "3"
 CairoMakie = "0.11"
-ClimaComms = "0.5.6"
-ClimaParams = "0.10"
+ClimaComms = "0.5.6, 0.6"
+ClimaCore = "0.13.2, 0.14"
+ClimaDiagnostics = "0.2"
+ClimaParams = "0.10.2"
 ClimaTimeSteppers = "0.7"
+ClimaUtilities = "0.1.2"
 DataFrames = "1"
 Dates = "1"
 DelimitedFiles = "1"
diff --git a/lib/ClimaLandSimulations/README.md b/lib/ClimaLandSimulations/README.md
index d24b377315..9bc8eaf0fa 100644
--- a/lib/ClimaLandSimulations/README.md
+++ b/lib/ClimaLandSimulations/README.md
@@ -17,11 +17,11 @@ For examples, see the [experiments](https://github.com/CliMA/ClimaLand.jl/lib/Cl
 
 ### Development of the `ClimaLandSimulations` subpackage
 
-    cd ClimaCore/lib/ClimaLandSimulations
+    cd ClimaLand/lib/ClimaLandSimulations
 
-    # Add ClimaCore to subpackage environment
+    # Add ClimaLand to subpackage environment
     julia --project -e 'using Pkg; Pkg.develop(path="../../")'
 
-    # Instantiate ClimaCoreMakie project environment
+    # Instantiate ClimaLandSimulations project environment
     julia --project -e 'using Pkg; Pkg.instantiate()'
     julia --project -e 'using Pkg; Pkg.test()'
diff --git a/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl b/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl
index 0bea798b49..bdc60e6e9d 100644
--- a/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl
+++ b/lib/ClimaLandSimulations/src/ClimaLandSimulations.jl
@@ -1,6 +1,7 @@
 module ClimaLandSimulations
 
 using CairoMakie
+using ClimaDiagnostics
 using WGLMakie
 using DataFrames
 using LaTeXStrings
diff --git a/src/ClimaLand.jl b/src/ClimaLand.jl
index 1c68029a1d..c93930ffbf 100644
--- a/src/ClimaLand.jl
+++ b/src/ClimaLand.jl
@@ -318,4 +318,8 @@ include("integrated/soil_energy_hydrology_biogeochemistry.jl")
 include("integrated/pond_soil_model.jl")
 include("integrated/soil_canopy_model.jl")
 
+# Diagnostics
+include(joinpath("diagnostics", "Diagnostics.jl"))
+import .Diagnostics as CLD # ClimaLand Diagnostics
+
 end
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
new file mode 100644
index 0000000000..ba163b8d44
--- /dev/null
+++ b/src/diagnostics/Diagnostics.jl
@@ -0,0 +1,18 @@
+module Diagnostics
+
+import ClimaComms
+
+using ..Bucket: BucketModel
+
+import ..SoilCanopyModel
+
+import ClimaDiagnostics:
+    DiagnosticVariable, ScheduledDiagnostic, average_pre_output_hook!
+
+import ClimaDiagnostics.Schedules: EveryStepSchedule, EveryDtSchedule
+
+import ClimaDiagnostics.Writers: HDF5Writer, NetCDFWriter
+
+include("diagnostic.jl")
+
+end
diff --git a/src/diagnostics/bucket_compute_methods.jl b/src/diagnostics/bucket_compute_methods.jl
new file mode 100644
index 0000000000..9ab65958f6
--- /dev/null
+++ b/src/diagnostics/bucket_compute_methods.jl
@@ -0,0 +1,119 @@
+# stored in p
+
+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
+
+function compute_net_radiation!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.R_n)
+    else
+        out .= p.bucket.R_n
+    end
+end
+
+function compute_surface_temperature!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.T_sfc)
+    else
+        out .= p.bucket.T_sfc
+    end
+end
+
+function compute_surface_specific_humidity!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
+    if isnothing(out)
+        return copy(p.bucket.q_sfc)
+    else
+        out .= p.bucket.q_sfc
+    end
+end
+
+function compute_latent_heat_flux!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.turbulent_fluxes.lhf)
+    else
+        out .= p.bucket.turbulent_fluxes.lhf
+    end
+end
+
+function compute_aerodynamic_resistance!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.turbulent_fluxes.r_ae)
+    else
+        out .= p.bucket.turbulent_fluxes.r_ae
+    end
+end
+
+function compute_sensible_heat_flux!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.turbulent_fluxes.shf)
+    else
+        out .= p.bucket.turbulent_fluxes.shf
+    end
+end
+
+function compute_vapor_flux!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.turbulent_fluxes.vapor_flux)
+    else
+        out .= p.bucket.turbulent_fluxes.vapor_flux
+    end
+end
+
+function compute_surface_air_density!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.ρ_sfc)
+    else
+        out .= p.bucket.ρ_sfc
+    end
+end
+
+# stored in Y
+
+function compute_soil_temperature!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(Y.bucket.T)
+    else
+        out .= Y.bucket.T
+    end
+end
+
+function compute_subsurface_water_storage!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
+    if isnothing(out)
+        return copy(Y.bucket.W)
+    else
+        out .= Y.bucket.W
+    end
+end
+
+function compute_surface_water_content!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(Y.bucket.Ws)
+    else
+        out .= Y.bucket.Ws
+    end
+end
+
+function compute_snow_water_equivalent!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(Y.bucket.σS)
+    else
+        out .= Y.bucket.σS
+    end
+end
diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
new file mode 100644
index 0000000000..e255d0f800
--- /dev/null
+++ b/src/diagnostics/default_diagnostics.jl
@@ -0,0 +1,117 @@
+export default_diagnostics
+
+# This file is included by Diagnostics.jl and defines all the defaults for
+# various models (e.g., Bucket, SoilCanopyModel). A model here is either a
+# standalone (e.g., Bucket) or integrated (e.g., SoilCanopy) model.
+#
+# If you are developing new models, add your defaults here. If you want to add
+# more high level interfaces, add them here. Feel free to include extra files.
+
+# Bucket model
+
+"""
+    function common_diagnostics(
+                                period,
+                                reduction,
+                                output_writer,
+                                t_start,
+                                short_names...;
+                                pre_output_hook! = nothing,
+                               )
+
+Helper function to define functions like `daily_max`.
+"""
+function common_diagnostics(
+    period,
+    reduction,
+    output_writer,
+    t_start,
+    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
+    ]
+end
+
+include("standard_diagnostic_frequencies.jl")
+
+# Bucket
+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",
+        "tsoil",
+        "wsoil",
+        "wsfc",
+        "ssfc",
+    ] # TO DO: would it be helpful to return this list?
+
+    default_outputs =
+        hourly_averages(bucket_diagnostics...; output_writer, t_start)
+    return [default_outputs...]
+end
+
+# SoilCanopyModel
+function default_diagnostics(
+    land_model::SoilCanopyModel,
+    t_start;
+    output_writer,
+)
+
+    define_diagnostics!(land_model)
+
+    soilcanopy_diagnostics = [
+        "rn",
+        "lhf",
+        "rae",
+        "shf",
+        "vflux",
+        "tsoil",
+        "slw",
+        "infil",
+        "scd",
+        "scms",
+        "gs",
+        "mt",
+        "trans",
+        "rain", # do we want?
+        "an",
+        "gpp",
+        "rd",
+        "vcmax25",
+        "par",
+        "apar",
+        "rpar",
+        "tpar",
+        "nir",
+        "anir",
+        "rnir",
+        "tnir",
+        "swn",
+        "lwn",
+        "ra",
+        "soilco2",
+    ]
+
+    default_outputs =
+        hourly_averages(soilcanopy_diagnostics...; output_writer, t_start)
+    return [default_outputs...]
+end
diff --git a/src/diagnostics/define_diagnostics.jl b/src/diagnostics/define_diagnostics.jl
new file mode 100644
index 0000000000..cb14a796a1
--- /dev/null
+++ b/src/diagnostics/define_diagnostics.jl
@@ -0,0 +1,645 @@
+# General helper functions for undefined diagnostics for a particular model
+error_diagnostic_variable(variable, land_model::T) where {T} =
+    error("Cannot compute $variable with model = $T")
+
+# generate_error_functions is helper macro that generates the error message
+# when the user tries calling something that is incompatible with the model
+macro generate_error_functions(variable_names...)
+    functions = Expr[]
+    for variable in variable_names
+        function_name_sym = Symbol("compute_", variable, "!")
+        body = esc(quote
+            function $function_name_sym(_, _, _, _, land_model)
+                error_diagnostic_variable($variable, land_model)
+            end
+        end)
+        push!(functions, body)
+    end
+    return quote
+        $(functions...)
+    end
+end
+
+# TODO: Automatically generate this list from the names of the diagnostics
+@generate_error_functions "soil_net_radiation" "soil_latent_heat_flux" "soil_aerodynamic_resistance" "soil_sensible_heat_flux" "vapor_flux" "soil_temperature" "soil_water_liquid" "infiltration" "soilco2_diffusivity" "soilco2_source_microbe" "stomatal_conductance" "medlyn_term" "canopy_transpiration" "rainfall" "snowfall" "pressure" "wind_speed" "specific_humidity" "air_co2" "radiation_shortwave_down" "radiation_longwave_down" "photosynthesis_net_leaf" "photosynthesis_net_canopy" "respiration_leaf" "vcmax25" "photosynthetically_active_radiation" "photosynthetically_active_radiation_absorbed" "photosynthetically_active_radiation_reflected" "photosynthetically_active_radiation_transmitted" "near_infrared_radiation" "near_infrared_radiation_absorbed" "near_infrared_radiation_reflected" "near_infrared_radiation_transmitted" "radiation_shortwave_net" "radiation_longwave_net" "autotrophic_respiration" "soilco2" "heterotrophic_respiration" "soil_hydraulic_conductivity" "soil_water_potential" "soil_thermal_conductivity" "solar_zenith_angle" "moisture_stress_factor" "canopy_water_potential" "cross_section" "cross_section_roots" "area_index" "canopy_latent_heat_flux" "canopy_sensible_heat_flux" "canopy_aerodynamic_resistance" "canopy_temperature" "soil_ice"
+
+"""
+    define_diagnostics!(land_model)
+
+Calls `add_diagnostic_variable!` for all available variables specializing the
+compute function for `land_model`.
+"""
+function define_diagnostics!(land_model)
+
+    # Stored in p
+
+    # Albedo
+    add_diagnostic_variable!(
+        short_name = "alpha",
+        long_name = "Albedo",
+        standard_name = "albedo",
+        units = "",
+        comments = "The fraction of incoming radiation reflected by the land surface.",
+        compute! = (out, Y, p, t) -> compute_albedo!(out, Y, p, t, land_model),
+    )
+
+    # Net radiation
+    add_diagnostic_variable!(
+        short_name = "rn",
+        long_name = "Net Radiation",
+        standard_name = "net_radiation",
+        units = "W m^-2",
+        comments = "Difference between incoming and outgoing shortwave and longwave radiation at the land surface.",
+        compute! = (out, Y, p, t) ->
+            compute_net_radiation!(out, Y, p, t, land_model),
+    )
+
+    # Surface temperature
+    add_diagnostic_variable!(
+        short_name = "tsfc",
+        long_name = "Surface Temperature",
+        standard_name = "surface_temperature",
+        units = "K",
+        comments = "Temperature of the land surface.",
+        compute! = (out, Y, p, t) ->
+            compute_surface_temperature!(out, Y, p, t, land_model),
+    )
+
+    # Surface specific humidity
+    add_diagnostic_variable!(
+        short_name = "qsfc",
+        long_name = "Surface Specific Humidity",
+        standard_name = "surface_specific_humidity",
+        units = "",
+        comments = "Ratio of water vapor mass to total moist air parcel mass.",
+        compute! = (out, Y, p, t) ->
+            compute_surface_specific_humidity!(out, Y, p, t, land_model),
+    )
+
+    # Latent heat flux
+    add_diagnostic_variable!(
+        short_name = "lhf",
+        long_name = "Latent Heat Flux",
+        standard_name = "latent_heat_flux",
+        units = "W m^-2",
+        comments = "Exchange of energy at the land-atmosphere interface due to water evaporation or sublimation.",
+        compute! = (out, Y, p, t) ->
+            compute_latent_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    # Aerodynamic resistance
+    add_diagnostic_variable!(
+        short_name = "rae",
+        long_name = "Aerodynamic Resistance",
+        standard_name = "aerodynamic_resistance",
+        units = "m s^-1",
+        comments = "Effiency of turbulent transport controlling the land-atmosphere exchange of sensible and latent heat.",
+        compute! = (out, Y, p, t) ->
+            compute_aerodynamic_resistance!(out, Y, p, t, land_model),
+    )
+
+    # Sensible heat flux
+    add_diagnostic_variable!(
+        short_name = "shf",
+        long_name = "Sensible Heat Flux",
+        standard_name = "sensible_heat_flux",
+        units = "W m^-2",
+        comments = "Exchange of energy at the land-atmosphere interface due to temperature difference.",
+        compute! = (out, Y, p, t) ->
+            compute_sensible_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    # Vapor flux
+    add_diagnostic_variable!(
+        short_name = "vflux",
+        long_name = "Liquid water evaporation",
+        standard_name = "vapor_flux",
+        units = "m s^-1",
+        comments = "Flux of water from the land surface to the atmosphere. E.g., evaporation or sublimation.",
+        compute! = (out, Y, p, t) ->
+            compute_vapor_flux!(out, Y, p, t, land_model),
+    )
+
+    # Surface air density
+    add_diagnostic_variable!(
+        short_name = "rhosfc",
+        long_name = "Surface Air Density",
+        standard_name = "surface_air_density",
+        units = "kg m^−3",
+        comments = "Density of air at the land-atmosphere interface.",
+        compute! = (out, Y, p, t) ->
+            compute_surface_air_density!(out, Y, p, t, land_model),
+    )
+
+    # Stored in Y
+
+    # Soil temperature (3D) at depth
+    add_diagnostic_variable!(
+        short_name = "tsoil",
+        long_name = "Soil temperature",
+        standard_name = "soil_temperature",
+        units = "K",
+        comments = "Soil temperature at multiple soil depth.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_temperature!(out, Y, p, t, land_model),
+    )
+
+    # Surbsurface water storage
+    add_diagnostic_variable!(
+        short_name = "wsoil",
+        long_name = "subsurface Water Storage",
+        standard_name = "subsurface_water_storage",
+        units = "m",
+        comments = "Soil water content.",
+        compute! = (out, Y, p, t) ->
+            compute_subsurface_water_storage!(out, Y, p, t, land_model),
+    )
+
+    # Surface water content
+    add_diagnostic_variable!(
+        short_name = "wsfc",
+        long_name = "Surface Water Content",
+        standard_name = "surface_water_content",
+        units = "m",
+        comments = "Water at the soil surface.",
+        compute! = (out, Y, p, t) ->
+            compute_surface_water_content!(out, Y, p, t, land_model),
+    )
+
+    # Surface snow water content
+    add_diagnostic_variable!(
+        short_name = "ssfc",
+        long_name = "Snow Water Equivalent",
+        standard_name = "snow_water_equivalent",
+        units = "m",
+        comments = "Snow at the soil surface, expressed in water equivalent.",
+        compute! = (out, Y, p, t) ->
+            compute_snow_water_equivalent!(out, Y, p, t, land_model),
+    )
+
+    ###### SoilCanopyModel ######
+
+    add_diagnostic_variable!(
+        short_name = "slw",
+        long_name = "Soil Liquid Water",
+        standard_name = "soil_liquid_water",
+        units = "m^3 m^-3",
+        comments = "Soil moisture, the liquid water volume per soil volume. This liquid water is located in the soil pores.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_water_liquid!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "infil",
+        long_name = "Infiltration",
+        standard_name = "infiltration",
+        units = "m s^-1", # double check
+        comments = "The flux of liquid water into the soil.",
+        compute! = (out, Y, p, t) ->
+            compute_infiltration!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "scd",
+        long_name = "Soil CO2 Diffusivity",
+        standard_name = "soil_co2_diffusivity",
+        units = "", # need to add
+        comments = "The diffusivity of CO2 in the porous phase of the soil. Depends on soil texture, moisture, and temperature.",
+        compute! = (out, Y, p, t) ->
+            compute_soilco2_diffusivity!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "scms",
+        long_name = "Soil CO2 Microbial Source",
+        standard_name = "soil_co2_microbial_source",
+        units = "", # check
+        comments = "The production of CO2 by microbes in the soil. Vary by layers of soil depth.",
+        compute! = (out, Y, p, t) ->
+            compute_soilco2_source_microbe!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "gs",
+        long_name = "Stomatal Conductance",
+        standard_name = "stomatal_conductance",
+        units = "m s^-1",
+        comments = "The conductance of leaves. This depends on stomatal opening. It varies with factors such as soil moisture or atmospheric water demand.",
+        compute! = (out, Y, p, t) ->
+            compute_stomatal_conductance!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "mt",
+        long_name = "Medlyn Term",
+        standard_name = "medlyn_term",
+        units = "", # check
+        comments = "",
+        compute! = (out, Y, p, t) ->
+            compute_medlyn_term!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "trans",
+        long_name = "Canopy Transpiration",
+        standard_name = "canopy_transpiration",
+        units = "", # check
+        comments = "The water evaporated from the canopy due to leaf transpiration.",
+        compute! = (out, Y, p, t) ->
+            compute_canopy_transpiration!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!( # not actually computed, but read from input or atmosphere model, and stored in p...
+        short_name = "rain",
+        long_name = "Rainfall",
+        standard_name = "rainfall",
+        units = "m",
+        comments = "Precipitation of liquid water.",
+        compute! = (out, Y, p, t) ->
+            compute_rainfall!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "an",
+        long_name = "Leaf Net Photosynthesis",
+        standard_name = "leaf_net_photosynthesis",
+        units = "", # check
+        comments = "Net photosynthesis of a leaf.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthesis_net_leaf!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "gpp",
+        long_name = "Gross Primary Productivity",
+        standard_name = "gross_primary_productivity",
+        units = "",
+        comments = "Net photosynthesis of the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthesis_net_canopy!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "rd",
+        long_name = "Leaf Respiration",
+        standard_name = "leaf_dark_respiration",
+        units = "",
+        comments = "Leaf respiration, called dark respiration because usually measured in the abscence of radiation.",
+        compute! = (out, Y, p, t) ->
+            compute_respiration_leaf!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "vcmax25",
+        long_name = "Vcmax25",
+        standard_name = "vcmax25",
+        units = "",
+        comments = "The parameter vcmax at 25 degree celsius. Important for the Farquhar model of leaf photosynthesis.",
+        compute! = (out, Y, p, t) -> compute_vcmax25!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "par",
+        long_name = "Photosynthetically Active Radiation",
+        standard_name = "photosynthetically_active_radiation",
+        units = "",
+        comments = "The subset of total radiation that activates photosynthesis reaching the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthetically_active_radiation!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "apar",
+        long_name = "Absorbed Photosynthetically Active Radiation",
+        standard_name = "absorbed_photosynthetically_active_radiation",
+        units = "",
+        comments = "The amount of photosynthetically active radiation absorbed by the leaf. The rest if reflected or transmitted.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthetically_active_radiation_absorbed!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "rpar",
+        long_name = "Reflected Photosynthetically Active Radiation",
+        standard_name = "reflected_photosynthetically_active_radiation",
+        units = "",
+        comments = "The amount of photosynthetically active radiation reflected by leaves.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthetically_active_radiation_reflected!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "tpar",
+        long_name = "Transmitted Photosynthetically Active Radiation",
+        standard_name = "transmitted_photosynthetically_active_radiation",
+        units = "",
+        comments = "The amount of photosynthetically active radiation transmitted by leaves.",
+        compute! = (out, Y, p, t) ->
+            compute_photosynthetically_active_radiation_transmitted!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "nir",
+        long_name = "Near Infrared Radiation",
+        standard_name = "near_infrared_radiation",
+        units = "W m^-2",
+        comments = "The amount of near infrared radiation reaching the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_near_infrared_radiation!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "anir",
+        long_name = "Absorbed Near Infrared Radiation",
+        standard_name = "absorbed_near_infrared_radiation",
+        units = "W m^-2",
+        comments = "The amount of near infrared radiation reaching the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_near_infrared_radiation_absorbed!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "rnir",
+        long_name = "Reflected Near Infrared Radiation",
+        standard_name = "reflected_near_infrared_radiation",
+        units = "W m^-2",
+        comments = "The amount of near infrared radiation reaching the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_near_infrared_radiation_reflected!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "tnir",
+        long_name = "Transmitted Near Infrared Radiation",
+        standard_name = "transmitted_near_infrared_radiation",
+        units = "W m^-2",
+        comments = "The amount of near infrared radiation reaching the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_near_infrared_radiation_transmitted!(
+                out,
+                Y,
+                p,
+                t,
+                land_model,
+            ),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "swn",
+        long_name = "Net Shortwave Radiation",
+        standard_name = "net_shortwave_radiation",
+        units = "W m^-2",
+        comments = "The net (in minus out) radiation at the surface.",
+        compute! = (out, Y, p, t) ->
+            compute_radiation_shortwave_net!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "lwn",
+        long_name = "Net Longwave Radiation",
+        standard_name = "net_longwave_radiation",
+        units = "W m^-2",
+        comments = "The net (in minus out) radiation at the surface.",
+        compute! = (out, Y, p, t) ->
+            compute_radiation_longwave_net!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "ra",
+        long_name = "Autotrophic Respiration",
+        standard_name = "autotrophic_respiration",
+        units = "mol m^-2 s^-1",
+        comments = "Canopy autotrophic respiration, the sum of leaves, stems and roots respiration.",
+        compute! = (out, Y, p, t) ->
+            compute_autotrophic_respiration!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "soilco2",
+        long_name = "Soil CO2 concentration",
+        standard_name = "soil_co2",
+        units = "",
+        comments = "Concentration of CO2 in the air of soil pores.",
+        compute! = (out, Y, p, t) -> compute_soilco2!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "soilrn",
+        long_name = "Soil Net Radiation",
+        standard_name = "soil_net_radiation",
+        units = "W m^-2",
+        comments = "Net radiation at the soil surface.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_net_radiation!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "soillhf",
+        long_name = "Soil Latent Heat Flux",
+        standard_name = "soil_Latent_Heat_Flux",
+        units = "W m^-2",
+        comments = "Soil evaporation.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_latent_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "soilshf",
+        long_name = "Soil Sensible Heat Flux",
+        standard_name = "soil_sensible_Heat_Flux",
+        units = "W m^-2",
+        comments = "Soil sensible heat flux.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_sensible_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "soilrae",
+        long_name = "Soil Aerodynamic Resistance",
+        standard_name = "soil_aerodynamic_resistance",
+        units = "",
+        comments = "Soil aerodynamic resistance.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_aerodynamic_resistance!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "hr",
+        long_name = "Heterotrophic Respiration",
+        standard_name = "heterotrophic_respiration",
+        units = "mol m^-2 s^-1",
+        comments = "CO2 efflux at the soil surface due to microbial decomposition of soil organic matter.",
+        compute! = (out, Y, p, t) ->
+            compute_heterotrophic_respiration!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "shc",
+        long_name = "Soil Hydraulic Conductivity",
+        standard_name = "soil_hydraulic_conductivity",
+        units = "",
+        comments = "Soil hydraulic conductivity.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_hydraulic_conductivity!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "stc",
+        long_name = "Soil Thermal Conductivity",
+        standard_name = "soil_thermal_conductivity",
+        units = "",
+        comments = "Soil thermal conductivity.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_thermal_conductivity!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "swp",
+        long_name = "Soil Water Potential",
+        standard_name = "soil_water_potential",
+        units = "",
+        comments = "Soil water potential.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_water_potential!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "sza",
+        long_name = "Solar Zenith Angle",
+        standard_name = "solar_zenith_angle",
+        units = "",
+        comments = "Solar zenith angle.",
+        compute! = (out, Y, p, t) ->
+            compute_solar_zenith_angle!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "msf",
+        long_name = "Moisture Stress Factor",
+        standard_name = "moisture_stress_factor",
+        units = "",
+        comments = "Sensitivity of plants conductance to soil water content.",
+        compute! = (out, Y, p, t) ->
+            compute_moisture_stress_factor!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "cwp",
+        long_name = "Canopy Water Potential",
+        standard_name = "canopy_water_potential",
+        units = "",
+        comments = "The water potential of the canopy.",
+        compute! = (out, Y, p, t) ->
+            compute_canopy_water_potential!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "fa",
+        long_name = "Cross Section",
+        standard_name = "cross_section",
+        units = "",
+        comments = "The area of stem relative to ground area.", #??
+        compute! = (out, Y, p, t) ->
+            compute_cross_section!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "far",
+        long_name = "Root Cross Section",
+        standard_name = "Cross Section",
+        units = "",
+        comments = "The area of roots relative to ground area.", #??
+        compute! = (out, Y, p, t) ->
+            compute_cross_section_roots!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "ai",
+        long_name = "Area Index",
+        standard_name = "area_index",
+        units = "",
+        comments = "The area index.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_area_index!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "clhf",
+        long_name = "Canopy Latent Heat Flux",
+        standard_name = "canopy_latent_heat_flux",
+        units = "",
+        comments = "Canopy evaporation.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_canopy_latent_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "cshf",
+        long_name = "Canopy Sensible Heat Flux",
+        standard_name = "canopy_sensible_heat_flux",
+        units = "",
+        comments = "Canopy sensible heat flux.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_canopy_sensible_heat_flux!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "crae",
+        long_name = "Canopy Aerodynamic Resistance",
+        standard_name = "canopy_aerodynamic_resistance",
+        units = "",
+        comments = "Canopy aerodynamic_resistance.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_canopy_aerodynamic_resistance!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "ct",
+        long_name = "Canopy Temperature",
+        standard_name = "canopy_temperature",
+        units = "",
+        comments = "Canopy temperature.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_canopy_temperature!(out, Y, p, t, land_model),
+    )
+
+    add_diagnostic_variable!(
+        short_name = "si",
+        long_name = "Soil Ice",
+        standard_name = "soil_ice",
+        units = "m^3 m^-3",
+        comments = "soil ice.", #?? of steam, leaves, roots?
+        compute! = (out, Y, p, t) ->
+            compute_soil_ice!(out, Y, p, t, land_model),
+    )
+end
diff --git a/src/diagnostics/diagnostic.jl b/src/diagnostics/diagnostic.jl
new file mode 100644
index 0000000000..4894dd3f9c
--- /dev/null
+++ b/src/diagnostics/diagnostic.jl
@@ -0,0 +1,99 @@
+# ClimaLand diagnostics contains # a dictionary `ALL_DIAGNOSTICS` with all the
+# diagnostics we know how to compute, keyed over their short name. If you want
+# to add more diagnostics, look at the included files. You can add your own file
+# if you want to define several new diagnostics that are conceptually related.
+# The dictionary `ALL_DIAGNOSTICS` should be considered an implementation
+# detail, use the getters/setters.
+
+const ALL_DIAGNOSTICS = Dict{String, DiagnosticVariable}()
+
+"""
+
+    add_diagnostic_variable!(; short_name,
+                               long_name,
+                               standard_name,
+                               units,
+                               description,
+                               compute!)
+
+
+Add a new variable to the `ALL_DIAGNOSTICS` dictionary (this function mutates the state of
+`ClimaLand.ALL_DIAGNOSTICS`).
+
+If possible, please follow the naming scheme outline in
+https://airtable.com/appYNLuWqAgzLbhSq/shrKcLEdssxb8Yvcp/tblL7dJkC3vl5zQLb
+
+Keyword arguments
+=================
+
+- `short_name`: Name used to identify the variable in the output files and in the file
+                names. Short but descriptive. `ClimaLand` diagnostics are identified by the
+                short name. We follow the Coupled Model Intercomparison Project conventions.
+
+- `long_name`: Name used to identify the variable in the output files.
+
+- `standard_name`: Standard name, as in
+                   http://cfconventions.org/Data/cf-standard-names/71/build/cf-standard-name-table.html
+
+- `units`: Physical units of the variable.
+
+- `comments`: More verbose explanation of what the variable is, or comments related to how
+              it is defined or computed.
+
+- `compute!`: Function that compute the diagnostic variable from the state.
+                             It has to take two arguments: the `integrator`, and a
+                             pre-allocated area of memory where to write the result of the
+                             computation. If no pre-allocated area is available, a new
+                             one will be allocated. To avoid extra allocations, this
+                             function should perform the calculation in-place (i.e., using
+                             `.=`).
+
+"""
+function add_diagnostic_variable!(;
+    short_name,
+    long_name,
+    standard_name = "",
+    units,
+    comments = "",
+    compute!,
+)
+    haskey(ALL_DIAGNOSTICS, short_name) && @warn(
+        "overwriting diagnostic `$short_name` entry containing fields\n" *
+        "$(map(
+            field -> "$(getfield(ALL_DIAGNOSTICS[short_name], field))",
+            filter(field -> field != :compute!, fieldnames(DiagnosticVariable)),
+        ))"
+    )
+
+    ALL_DIAGNOSTICS[short_name] = DiagnosticVariable(;
+        short_name,
+        long_name,
+        standard_name,
+        units,
+        comments,
+        compute!,
+    )
+end
+
+"""
+
+    get_diagnostic_variable!(short_name)
+
+Return a `DiagnosticVariable` from its `short_name`, if it exists.
+"""
+function get_diagnostic_variable(short_name)
+    haskey(ALL_DIAGNOSTICS, short_name) ||
+        error("diagnostic $short_name does not exist")
+
+    return ALL_DIAGNOSTICS[short_name]
+end
+
+# Do you want to define more diagnostics? Add them here
+include("bucket_compute_methods.jl")
+include("soilcanopy_compute_methods.jl")
+
+# define_diagnostics.jl contains the list of all the diagnostics
+include("define_diagnostics.jl")
+
+# Default diagnostics and higher level interfaces
+include("default_diagnostics.jl")
diff --git a/src/diagnostics/soilcanopy_compute_methods.jl b/src/diagnostics/soilcanopy_compute_methods.jl
new file mode 100644
index 0000000000..80653a5d91
--- /dev/null
+++ b/src/diagnostics/soilcanopy_compute_methods.jl
@@ -0,0 +1,579 @@
+# stored in p
+
+function compute_soil_net_radiation!(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soil.R_n)
+    else
+        out .= p.soil.R_n
+    end
+end
+
+function compute_soil_latent_heat_flux!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soil.turbulent_flux.lhf) # is this different from canopy.energy.lhf?
+    else
+        out .= p.soil.turbulent_flux.lhf
+    end
+end
+
+function compute_soil_aerodynamic_resistance!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soil.turbulent_fluxes.r_ae)
+    else
+        out .= p.soil.turbulent_fluxes.r_ae
+    end
+end
+
+function compute_soil_sensible_heat_flux!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soil.turbulent_fluxes.shf)
+    else
+        out .= p.soil.turbulent_fluxes.shf
+    end
+end
+
+function compute_vapor_flux!(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soil.turbulent_fluxes.vapor_flux)
+    else
+        out .= p.soil.turbulent_fluxes.vapor_flux
+    end
+end
+
+function compute_soil_temperature!(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.T) # or is it p.soil.T?
+    else
+        out .= p.drivers.T
+    end
+end
+
+function compute_soil_water_liquid(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soil.θ_l) # or is it in Y.soil?
+    else
+        out .= p.soil.θ_l
+    end
+end
+
+function compute_infiltration(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soil.infiltration)
+    else
+        out .= p.soil.infiltration
+    end
+end
+
+function compute_soilco2_diffusivity(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soilco2.D) # NOTE: we will need a method to compute surface co2 efflux
+    else
+        out .= p.soilco2.D
+    end
+end
+
+function compute_soilco2_source_microbe(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soilco2.Sm)
+    else
+        out .= p.soilco2.Sm
+    end
+end
+
+function compute_stomatal_conductance(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.conductance.gs) # doublecheck: stomata, not canopy
+    else
+        out .= p.canopy.conductance.gs
+    end
+end
+
+function compute_medlyn_term(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.conductance.medlyn_term)
+    else
+        out .= p.canopy.conductance.medlyn_term
+    end
+end
+
+function compute_canopy_transpiration(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.transpiration) # doublecheck: canopy, not leaf
+    else
+        out .= p.canopy.transpiration
+    end
+end
+
+function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.P_liq) # I guess this is read and put it p. not computed. curious if we should handle this differently.
+    else
+        out .= p.drivers.P_liq
+    end
+end
+
+function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.P_snow) # following comment above, we could have a default getting only model output, and one also getting some inputs like drivers
+    else
+        out .= p.drivers.P_snow
+    end
+end
+
+function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.P) # not sure if precip or pressure
+    else
+        out .= p.drivers.P
+    end
+end
+
+function compute_wind_speed(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.u)
+    else
+        out .= p.drivers.u
+    end
+end
+
+function compute_specific_humidity(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.q) # check if this is correct. Also, check if Bucket has it and same name or not.
+    else
+        out .= p.drivers.q
+    end
+end
+
+function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.c_co2)
+    else
+        out .= p.drivers.c_co2
+    end
+end
+
+function compute_radiation_shortwave_down(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.drivers.SW_d)
+    else
+        out .= p.drivers.SW_d
+    end
+end
+
+function compute_radiation_longwave_down(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.drivers.LW_d)
+    else
+        out .= p.drivers.LW_d
+    end
+end
+
+function compute_photosynthesis_net_leaf(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.photosynthesis.An)
+    else
+        out .= p.canopy.photosynthesis.An
+    end
+end
+
+function compute_photosynthesis_net_canopy(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+) # could be gross primary productivity, but this is consistent with leaf
+    if isnothing(out)
+        return copy(p.canopy.photosynthesis.GPP)
+    else
+        out .= p.canopy.photosynthesis.GPP
+    end
+end
+
+function compute_respiration_leaf(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.photosynthesis.Rd)
+    else
+        out .= p.canopy.photosynthesis.Rd
+    end
+end
+
+function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.photosynthesis.vcmax25)
+    else
+        out .= p.canopy.photosynthesis.vcmax25
+    end
+end
+
+function compute_photosynthetically_active_radiation(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.par)
+    else
+        out .= p.canopy.radiative_transfer.par
+    end
+end
+
+function compute_photosynthetically_active_radiation_absorbed(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.apar)
+    else
+        out .= p.canopy.radive_transfer.apar
+    end
+end
+
+function compute_photosynthetically_active_radiation_reflected(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.rpar)
+    else
+        out .= p.canopy.radiative_transfer.rpar
+    end
+end
+
+function compute_photosynthetically_active_radiation_transmitted(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.tpar)
+    else
+        out .= p.canopy.radiative_transfer.tpar
+    end
+end
+
+function compute_near_infrared_radiation(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.nir)
+    else
+        out .= p.canopy.radiative_transfer.nir
+    end
+end
+
+function compute_near_infrared_radiation_absorbed(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.anir)
+    else
+        out .= p.canopy.radiative_transfer.anir
+    end
+end
+
+function compute_near_infrared_radiation_reflected(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.rnir)
+    else
+        out .= p.canopy.radiative_transfer.rnir
+    end
+end
+
+function compute_near_infrared_radiation_transmitted(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.tnir)
+    else
+        out .= p.canopy.radiative_transfer.tnir
+    end
+end
+
+function compute_radiation_shortwave_net(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.SW_n)
+    else
+        out .= p.canopy.radiative_transfer.SW_n
+    end
+end
+
+function compute_radiation_longwave_net(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.radiative_transfer.LW_n)
+    else
+        out .= p.canopy.radiative_transfer.LW_n
+    end
+end
+
+function compute_autotrophic_respiration(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.autotrophic_respiration.Ra)
+    else
+        out .= p.canopy.autotrophic_respiration.Ra
+    end
+end
+
+function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(Y.soilco2.C)
+    else
+        out .= Y.soilco2.C
+    end
+end
+
+function compute_heterotrophic_respiration(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soilco2.top_bc)
+    else
+        p.soilco2.top_bc
+    end
+end
+
+function compute_soil_hydraulic_conductivity(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soil.K)
+    else
+        p.soil.K
+    end
+end
+
+function compute_soil_water_potential(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.soil.ψ)
+    else
+        p.soil.ψ
+    end
+end
+
+function compute_soil_thermal_conductivity(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.soil.κ)
+    else
+        p.soil.κ
+    end
+end
+
+function compute_solar_zenith_angle(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.drivers.θs)
+    else
+        p.drivers.θs
+    end
+end
+
+function compute_moisture_stress_factor(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.hydraulics.β)
+    else
+        p.canopy.hydraulics.β
+    end
+end
+
+function compute_canopy_water_potential(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.hydraulics.ψ)
+    else
+        p.canopy.hydraulics.ψ
+    end
+end
+
+function compute_cross_section(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.hydraulics.fa)
+    else
+        p.canopy.hydraulics.fa
+    end
+end
+
+function compute_cross_section_roots(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.hydraulics.fa_roots)
+    else
+        p.canopy.hydraulics.fa_roots
+    end
+end
+
+function compute_area_index(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(p.canopy.hydraulics.area_index)
+    else
+        p.canopy.hydraulics.area_index
+    end
+end
+
+function compute_canopy_latent_heat_flux(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.energy.lhf)
+    else
+        p.canopy.canopy.lhf
+    end
+end
+
+function compute_canopy_sensible_heat_flux(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.energy.shf)
+    else
+        p.canopy.canopy.shf
+    end
+end
+
+function compute_canopy_aerodynamic_resistance(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
+    if isnothing(out)
+        return copy(p.canopy.energy.r_ae)
+    else
+        p.canopy.canopy.r_ae
+    end
+end
+
+function compute_canopy_temperature(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(Y.canopy.energy.T)
+    else
+        Y.canopy.energy.T
+    end
+end
+
+function compute_soil_ice(out, Y, p, t, land_model::SoilCanopyModel)
+    if isnothing(out)
+        return copy(Y.soil.θ_i)
+    else
+        Y.soil.θ_i
+    end
+end
diff --git a/src/diagnostics/standard_diagnostic_frequencies.jl b/src/diagnostics/standard_diagnostic_frequencies.jl
new file mode 100644
index 0000000000..abf80f4d47
--- /dev/null
+++ b/src/diagnostics/standard_diagnostic_frequencies.jl
@@ -0,0 +1,262 @@
+"""
+    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]
diff --git a/test/Project.toml b/test/Project.toml
index 9d8e99120a..03175b535e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,10 +3,13 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
 ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
+ClimaDiagnostics = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 ClimaLand = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
+ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
diff --git a/test/diagnostics/diagnostics_tests.jl b/test/diagnostics/diagnostics_tests.jl
new file mode 100644
index 0000000000..337db9320e
--- /dev/null
+++ b/test/diagnostics/diagnostics_tests.jl
@@ -0,0 +1,17 @@
+using Test
+using ClimaLand
+
+@test isdefined(ClimaLand.Diagnostics, :compute_albedo!)
+
+# Just to trigger the error
+out = Y = p = t = land_model = nothing
+
+@test_throws ErrorException("Cannot compute albedo with model = Nothing") ClimaLand.Diagnostics.compute_albedo!(
+    out,
+    Y,
+    p,
+    t,
+    land_model,
+)
+
+@test_throws ErrorException ClimaLand.Diagnostics.get_diagnostic_variable("Foo")