diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 79933105275..adeeb2c0ded 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -87,6 +87,10 @@ steps:
         command: "julia --color=yes --project=.buildkite experiments/integrated/fluxnet/run_fluxnet.jl US-Ha1"
         artifact_paths: "experiments/integrated/fluxnet/US-Ha1/out/*png"
 
+      - label: "emr_test"
+        command: "julia --color=yes --project=.buildkite experiments/integrated/fluxnet/run_any_fluxnet.jl AU-Emr"
+        artifact_paths: "experiments/integrated/fluxnet/AU-Emr/out/*png"
+
       - label: "ozark_pft"
         command: "julia --color=yes --project=.buildkite experiments/integrated/fluxnet/ozark_pft.jl"
         artifact_paths: "experiments/integrated/fluxnet/US-MOz/pft/out/*png"
diff --git a/Artifacts.toml b/Artifacts.toml
index bd9505b77f6..c67fb66d8e5 100644
--- a/Artifacts.toml
+++ b/Artifacts.toml
@@ -50,17 +50,17 @@ git-tree-sha1 = "136f1db3ed969614fb589c5250545a3ed1e8aaab"
     url = "https://caltech.box.com/shared/static/of1admpndmikoumtgk5j3yvt92v71awk.gz"
 
 ["clm_data_0.9x1.25"]
-git-tree-sha1 = "283c62220fca8c4afe36a62348d3e9a159af2ee9"
+git-tree-sha1 = "38488be5dd476890deb6fbd2b26e3e1cee2d449e"
 
     [["clm_data_0.9x1.25".download]]
-    sha256 = "88d5899a729de800017a74bd8a7278582c2c55c0d8730a529bae384425d70e4e"
-    url = "https://caltech.box.com/shared/static/mxs3l1c1dppjwy81assk8w8ofn9d70pu.gz"
+    sha256 = "cb5ead78868243879261d26cbb42d665bfd98135ebbebc052487d4d13ebf3a98"
+    url = "https://caltech.box.com/shared/static/6ms2qzyiekorksewmfyie6mbonlq5ctp.gz"
 ["clm_data_0.125x0.125"]
-git-tree-sha1 = "6284ddefbc7937d9c1fb68fa731ff3f00b68e917"
+git-tree-sha1 = "b8816f0af2c6d182111a156c3b8eed2ed83efea7"
 
     [["clm_data_0.125x0.125".download]]
-    sha256 = "8d2a61baad53346515f4a66c4342f8aafce37ca9520020c17f99cc4f4aa98b15"
-    url = "https://caltech.box.com/shared/static/uteaw1l6fg7wmwhb3ey1c0jq4r8vg0ts.gz"
+    sha256 = "993b5b97b6ec22ec484db4e523b65a231d736ead8b5e0759ad3abd00050bb58c"
+    url = "https://caltech.box.com/shared/static/fx9138ysguyg4g65l4tvibbk16er0nqo.gz"
 [era5_land_forcing_data2008]
 git-tree-sha1 = "93d2e93f491e77cb8fba2a1b8b3946f38bde469e"
 [era5_land_forcing_data2008_lowres]
@@ -195,3 +195,6 @@ git-tree-sha1 = "c35ba0e899040cb8153226ab751f69100f475d39"
     [[mizoguchi_soil_freezing_data.download]]
     sha256 = "0027cc080ba45ba33dc790b176ec2854353ce7dce4eae4bef72963b0dd944e0b"
     url = "https://caltech.box.com/shared/static/tn1bnqjmegyetw5kzd2ixq5pbnb05s3u.gz"
+
+[fluxnet2015]
+git-tree-sha1 = "94d6eb8e3bc5fc361a43597ab4fb6941bdbe2850"
diff --git a/docs/src/tutorials/calibration/perfect_model_site_level_calibration.jl b/docs/src/tutorials/calibration/perfect_model_site_level_calibration.jl
index b8b30b0758e..a39452caed5 100644
--- a/docs/src/tutorials/calibration/perfect_model_site_level_calibration.jl
+++ b/docs/src/tutorials/calibration/perfect_model_site_level_calibration.jl
@@ -78,8 +78,10 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID)
 (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val))
 
 # Get maximum simulation start and end dates
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
 stop_date = DateTime(2010, 4, 1, 6, 30)  # Set the stop date manually
 Δt = 450.0  # seconds
 
diff --git a/docs/src/tutorials/integrated/changing_snowy_land_parameterizations.jl b/docs/src/tutorials/integrated/changing_snowy_land_parameterizations.jl
index 7187a200f57..4351addab0e 100644
--- a/docs/src/tutorials/integrated/changing_snowy_land_parameterizations.jl
+++ b/docs/src/tutorials/integrated/changing_snowy_land_parameterizations.jl
@@ -53,8 +53,11 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID);
 (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val));
 # Set a start and stop date of the simulation in UTC, as well as
 # a timestep in seconds
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset);
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset);
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 Δt = 450.0;
 
 # Setup the domain for the model. This corresponds to
diff --git a/docs/src/tutorials/integrated/fluxnet_data.jl b/docs/src/tutorials/integrated/fluxnet_data.jl
index 7b4324cc014..917b2109eb3 100644
--- a/docs/src/tutorials/integrated/fluxnet_data.jl
+++ b/docs/src/tutorials/integrated/fluxnet_data.jl
@@ -48,8 +48,11 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID)
 
 # It is also useful to know the bounds of the data,
 # in UTC, to use as the start and stop date of the simulation.
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 
 # Now we can construct the forcing objects. Under the hood, this
 # function finds the local path to the fluxtower data (and downloads it
diff --git a/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl b/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl
index 324137e740d..47654e6933c 100644
--- a/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl
+++ b/docs/src/tutorials/integrated/snowy_land_fluxnet_tutorial.jl
@@ -48,8 +48,11 @@ site_ID_val = FluxnetSimulations.replace_hyphen(site_ID);
 (; atmos_h) = FluxnetSimulations.get_fluxtower_height(FT, Val(site_ID_val));
 # Set a start and stop date of the simulation in UTC, as well as
 # a timestep in seconds
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset);
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 Δt = 450.0;
 
 # Setup the domain for the model. This corresponds to
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 498f4c787cc..1f80a65a858 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -103,8 +103,11 @@ prognostic_land_components = (:canopy, :soil, :soilco2)
 
 # Set up the timestepping information for the simulation
 dt = Float64(450) # 7.5 minutes
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 # Define the PFT land cover percentages for the Ozark site. Currently we only
 # use the dominant PFT, which for Ozark is deciduous broadleaf temperate trees.
 pft_pcts = [
diff --git a/experiments/integrated/fluxnet/ozark_pmodel.jl b/experiments/integrated/fluxnet/ozark_pmodel.jl
index 15bccef8a06..369a21b574e 100644
--- a/experiments/integrated/fluxnet/ozark_pmodel.jl
+++ b/experiments/integrated/fluxnet/ozark_pmodel.jl
@@ -95,8 +95,11 @@ dt = Float64(450) # 7.5 minutes
 
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
 (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet(
diff --git a/experiments/integrated/fluxnet/ozark_soilsnow.jl b/experiments/integrated/fluxnet/ozark_soilsnow.jl
index 697936cc589..724e044a5ad 100644
--- a/experiments/integrated/fluxnet/ozark_soilsnow.jl
+++ b/experiments/integrated/fluxnet/ozark_soilsnow.jl
@@ -97,8 +97,11 @@ dt = Float64(900)
 
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 
 # Height of sensor on flux tower
 atmos_h = FT(32)
diff --git a/experiments/integrated/fluxnet/run_any_fluxnet.jl b/experiments/integrated/fluxnet/run_any_fluxnet.jl
new file mode 100644
index 00000000000..7bb6f6e0cb8
--- /dev/null
+++ b/experiments/integrated/fluxnet/run_any_fluxnet.jl
@@ -0,0 +1,306 @@
+import SciMLBase
+import ClimaComms
+ClimaComms.@import_required_backends
+using ClimaCore
+import ClimaParams as CP
+using Dates
+using ClimaDiagnostics
+using ClimaUtilities
+
+using ClimaLand
+using ClimaLand.Domains: Column
+using ClimaLand.Snow
+using ClimaLand.Soil
+using ClimaLand.Soil.Biogeochemistry
+using ClimaLand.Canopy
+using ClimaLand.Canopy.PlantHydraulics
+import ClimaLand
+import ClimaLand.Parameters as LP
+using ClimaLand.Simulations: LandSimulation, solve!
+import ClimaLand.FluxnetSimulations as FluxnetSimulations
+using CairoMakie, ClimaAnalysis, GeoMakie, Printf, StatsBase
+import ClimaLand.LandSimVis as LandSimVis
+
+# Parse site_ID from command line
+if length(ARGS) < 1
+    error("Usage: julia --project run_any_fluxnet.jl <SITE_ID> (e.g. US-MOz)")
+end
+site_ID = ARGS[1]
+
+const FT = Float64
+toml_dict = LP.create_toml_dict(FT)
+climaland_dir = pkgdir(ClimaLand)
+prognostic_land_components = (:canopy, :snow, :soil, :soilco2)
+
+# Get the default domain info for a generic 10m soil column, 20 layers
+(; dz_tuple, nelements, zmin, zmax) = FluxnetSimulations.get_domain_info(FT)
+
+# Get location info from FLUXNET metadata
+(; time_offset, lat, long, atmos_h) = FluxnetSimulations.get_location(site_ID)
+
+# Construct the ClimaLand domain to run the simulation on
+land_domain = Column(;
+    zlim = (zmin, zmax),
+    nelements = nelements,
+    dz_tuple = dz_tuple,
+    longlat = (long, lat),
+)
+surface_domain = ClimaLand.Domains.obtain_surface_domain(land_domain)
+
+# Get PFT from CLM dominant PFT map
+pft = ClimaLand.Canopy.clm_dominant_pft(land_domain.space.surface)
+
+# Get parameters from global maps + PFT + metadata
+(;
+    soil_ν,
+    soil_K_sat,
+    soil_S_s,
+    soil_hydrology_cm,
+    θ_r,
+    ν_ss_quartz,
+    ν_ss_om,
+    ν_ss_gravel,
+    z_0m_soil,
+    z_0b_soil,
+    soil_ϵ,
+    soil_albedo,
+    Ω,
+    χl,
+    G_Function,
+    α_PAR_leaf,
+    λ_γ_PAR,
+    τ_PAR_leaf,
+    α_NIR_leaf,
+    τ_NIR_leaf,
+    ϵ_canopy,
+    ac_canopy,
+    g1,
+    Drel,
+    g0,
+    Vcmax25,
+    SAI,
+    f_root_to_shoot,
+    K_sat_plant,
+    ψ63,
+    Weibull_param,
+    a,
+    conductivity_model,
+    retention_model,
+    plant_ν,
+    plant_S_s,
+    rooting_depth,
+    n_stem,
+    n_leaf,
+    h_leaf,
+    h_stem,
+    h_canopy,
+) = FluxnetSimulations.get_parameters(FT, site_ID, land_domain, pft)
+
+# Set up the timestepping information for the simulation
+dt = Float64(450) # 7.5 minutes
+
+# Get data dates and forcing from the flux tower site
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
+(; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet(
+    site_ID,
+    lat,
+    long,
+    time_offset,
+    atmos_h,
+    start_date,
+    toml_dict,
+    FT,
+)
+
+# Set up the soil model with CLM albedo from global maps
+forcing = (; atmos, radiation)
+retention_parameters = (;
+    ν = soil_ν,
+    θ_r,
+    K_sat = soil_K_sat,
+    hydrology_cm = soil_hydrology_cm,
+)
+composition_parameters = (; ν_ss_om, ν_ss_quartz, ν_ss_gravel)
+
+soil = Soil.EnergyHydrology{FT}(
+    land_domain,
+    forcing,
+    toml_dict;
+    prognostic_land_components,
+    additional_sources = (ClimaLand.RootExtraction{FT}(),),
+    albedo = soil_albedo,
+    runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(),
+    retention_parameters,
+    composition_parameters,
+    S_s = soil_S_s,
+    z_0m = z_0m_soil,
+    z_0b = z_0b_soil,
+    emissivity = soil_ϵ,
+)
+
+# Soil microbes model
+co2_prognostic_soil = Soil.Biogeochemistry.PrognosticMet(soil.parameters)
+drivers = Soil.Biogeochemistry.SoilDrivers(co2_prognostic_soil, atmos)
+soilco2 =
+    Soil.Biogeochemistry.SoilCO2Model{FT}(land_domain, drivers, toml_dict)
+
+# Canopy model components
+# Radiative transfer
+radiation_parameters = (;
+    Ω,
+    G_Function,
+    α_PAR_leaf,
+    τ_PAR_leaf,
+    α_NIR_leaf,
+    τ_NIR_leaf,
+)
+radiative_transfer = Canopy.TwoStreamModel{FT}(
+    surface_domain,
+    toml_dict;
+    radiation_parameters,
+    ϵ_canopy,
+)
+
+# PModel conductance
+conductance = PModelConductance{FT}(toml_dict)
+
+# PModel photosynthesis - auto-detect C3/C4 from CLM maps
+surface_space = land_domain.space.surface
+clm_photo_params = ClimaLand.Canopy.clm_photosynthesis_parameters(surface_space)
+is_c3 = ClimaCore.Fields.field2array(clm_photo_params.is_c3)[1]
+photosynthesis = PModel{FT}(surface_domain, toml_dict; is_c3)
+
+# Soil moisture stress using soil retention parameters
+soil_moisture_stress = PiecewiseMoistureStressModel{FT}(
+    land_domain,
+    toml_dict;
+    soil_params = retention_parameters,
+)
+
+# Plant hydraulics
+LAI =
+    ClimaLand.Canopy.prescribed_lai_modis(surface_space, start_date, stop_date)
+maxLAI = FluxnetSimulations.get_maxLAI_at_site(start_date, lat, long)
+RAI = maxLAI * f_root_to_shoot
+hydraulics = Canopy.PlantHydraulicsModel{FT}(
+    surface_domain,
+    toml_dict;
+    n_stem,
+    n_leaf,
+    h_stem,
+    h_leaf,
+    ν = plant_ν,
+    S_s = plant_S_s,
+    conductivity_model,
+    retention_model,
+)
+height = h_stem + h_leaf
+biomass =
+    Canopy.PrescribedBiomassModel{FT}(; LAI, SAI, RAI, rooting_depth, height)
+
+# Energy model
+energy = Canopy.BigLeafEnergyModel{FT}(toml_dict; ac_canopy)
+
+ground = ClimaLand.PrognosticGroundConditions{FT}()
+canopy_forcing = (; atmos, radiation, ground)
+
+# Combine the components into a CanopyModel
+canopy = Canopy.CanopyModel{FT}(
+    surface_domain,
+    canopy_forcing,
+    LAI,
+    toml_dict;
+    prognostic_land_components,
+    radiative_transfer,
+    photosynthesis,
+    conductance,
+    soil_moisture_stress,
+    hydraulics,
+    energy,
+    biomass,
+)
+
+# Snow model
+snow = Snow.SnowModel(
+    FT,
+    surface_domain,
+    forcing,
+    toml_dict,
+    dt;
+    prognostic_land_components,
+)
+
+# Integrated land model
+land = LandModel{FT}(canopy, snow, soil, soilco2)
+set_ic! = FluxnetSimulations.make_set_fluxnet_initial_conditions(
+    site_ID,
+    start_date,
+    time_offset,
+    land,
+)
+
+# Diagnostics
+output_vars = [
+    "gpp",
+    "shf",
+    "lhf",
+    "swu",
+    "lwu",
+    "swc",
+    "swe",
+    "tsoil",
+    "sco2",
+    "so2",
+    "soc",
+    "scms",
+]
+diags = ClimaLand.default_diagnostics(
+    land,
+    start_date;
+    output_writer = ClimaDiagnostics.Writers.DictWriter(),
+    output_vars,
+    reduction_period = :halfhourly,
+)
+
+simulation = LandSimulation(
+    start_date,
+    stop_date,
+    dt,
+    land;
+    user_callbacks = (),
+    set_ic!,
+    updateat = Second(dt),
+    diagnostics = diags,
+)
+
+@time solve!(simulation)
+
+# Save comparison plots
+comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset)
+savedir = joinpath(
+    pkgdir(ClimaLand),
+    "experiments/integrated/fluxnet/$(site_ID)/pmodel_generic/out",
+)
+mkpath(savedir)
+LandSimVis.make_diurnal_timeseries(
+    land_domain,
+    diags,
+    start_date;
+    savedir,
+    short_names = ["gpp", "shf", "lhf", "swu", "lwu"],
+    spinup_date = start_date + Day(20),
+    comparison_data,
+)
+LandSimVis.make_timeseries(
+    land_domain,
+    diags,
+    start_date;
+    savedir,
+    short_names = ["swc", "tsoil", "swe"],
+    spinup_date = start_date + Day(20),
+    comparison_data,
+)
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index 0428fb5b344..743303e434b 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -101,8 +101,11 @@ dt = Float64(450) # 7.5 minutes
 
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
-(start_date, stop_date) =
-    FluxnetSimulations.get_data_dates(site_ID, time_offset)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
 (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet(
     site_ID,
     lat,
diff --git a/experiments/integrated/fluxnet/vaira_paper.jl b/experiments/integrated/fluxnet/vaira_paper.jl
index 3921abe300a..2e1acdc45f0 100644
--- a/experiments/integrated/fluxnet/vaira_paper.jl
+++ b/experiments/integrated/fluxnet/vaira_paper.jl
@@ -96,8 +96,11 @@ dt = Float64(450) # 7.5 minutes
 
 # This reads in the data from the flux tower site and creates
 # the atmospheric and radiative driver structs for the model
-(start_date, _) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
-stop_date = start_date + Year(1)
+(data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+# Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+# so TimeVaryingInputs have data before t=0 even if initial rows are missing
+start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+stop_date = min(start_date + Year(1), DateTime(2020, 12, 31, 23, 59, 59))
 (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet(
     site_ID,
     lat,
diff --git a/experiments/integrated/performance/conservation/ozark_conservation.jl b/experiments/integrated/performance/conservation/ozark_conservation.jl
index 6fae242099a..492ea73b69a 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation.jl
@@ -126,8 +126,11 @@ for float_type in (Float32, Float64)
     prognostic_land_components = (:canopy, :soil, :soilco2)
 
     # Get the atmospheric and radiation forcing data
-    (start_date, stop_date) =
-        FluxnetSimulations.get_data_dates(site_ID, time_offset)
+    (data_start, data_stop) = FluxnetSimulations.get_data_dates(site_ID, time_offset)
+    # Constrain to 2000-2020 (MODIS LAI availability) and skip first day
+    # so TimeVaryingInputs have data before t=0 even if initial rows are missing
+    start_date = max(data_start + Day(1), DateTime(2000, 1, 1))
+    stop_date = min(data_stop, DateTime(2020, 12, 31, 23, 59, 59))
     (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet(
         site_ID,
         lat,
diff --git a/ext/FluxnetSimulationsExt.jl b/ext/FluxnetSimulationsExt.jl
index 2ec50608469..daf69952c15 100644
--- a/ext/FluxnetSimulationsExt.jl
+++ b/ext/FluxnetSimulationsExt.jl
@@ -6,12 +6,15 @@ import ClimaUtilities.TimeManager: ITime, date
 using Thermodynamics
 using Dates
 using DelimitedFiles
-using DocStringExtensions
 using Insolation
 import ClimaLand.Parameters as LP
 using ClimaLand
 using ClimaLand.Canopy
 using ClimaLand.PlantHydraulics
+using ClimaLand.Soil
+using ClimaLand.Domains
+using ClimaCore
+using NCDatasets
 export prescribed_forcing_fluxnet,
     make_set_fluxnet_initial_conditions,
     get_comparison_data,
@@ -19,7 +22,11 @@ export prescribed_forcing_fluxnet,
     get_data_dt,
     replace_hyphen,
     get_parameters,
-    get_domain_info
+    get_domain_info,
+    get_location
+
+# Include metadata retrieval functions.
+include("fluxnet_simulations/get_fluxnet_metadata.jl")
 
 # Include site-specific configurations, as well as the default generic site.
 include("fluxnet_simulations/generic_site.jl")
diff --git a/ext/fluxnet_simulations/generic_site.jl b/ext/fluxnet_simulations/generic_site.jl
index e2d07c46fd0..c695b0a0115 100644
--- a/ext/fluxnet_simulations/generic_site.jl
+++ b/ext/fluxnet_simulations/generic_site.jl
@@ -2,7 +2,182 @@
 #        MODULE FUNCTIONS         #
 ###################################
 
-# TODO
+"""
+    get_domain_info(FT; dz_bottom = FT(1.5), dz_top = FT(0.1),
+        nelements = 20, zmin = FT(-10),  zmax = FT(0))
+
+Gets and returns primary domain information for a generic Fluxnet site, using
+default values corresponding to a 10m deep soil column with 20 layers, with
+a resolution of 10 cm at the top of the domain and 1.5m at the bottom of the domain.
+"""
+function FluxnetSimulations.get_domain_info(
+    FT;
+    dz_bottom = FT(1.5),
+    dz_top = FT(0.1),
+    nelements = 20,
+    zmin = FT(-10),
+    zmax = FT(0),
+)
+
+    dz_tuple = (dz_bottom, dz_top)
+
+    return (; dz_tuple, nelements, zmin, zmax)
+end
+
+"""
+    get_location(site_ID::String; kwargs)
+
+Returns geographical and time information for a generic Fluxnet site.
+The values are fetched from a list of Fluxnet sites metadata, and can
+be overwritten by passing the corresponding keyword arguments to this
+function.
+"""
+function FluxnetSimulations.get_location(
+    site_ID::String;
+    site_info = FluxnetSimulationsExt.get_site_info(site_ID),
+    lat = site_info[:lat],
+    long = site_info[:long],
+    time_offset = site_info[:time_offset],
+)
+    atmos_h = site_info[:atmospheric_sensor_height][1] # take the first height as default
+    return (; time_offset, lat, long, atmos_h)
+end
+
+"""
+    get_parameters(FT, site_ID, domain, pft::Pft; kwargs...)
+
+Gets parameters for a generic Fluxnet site based on PFT (plant functional type),
+supplemented with additional autofilled values from global maps.
+
+Data sources:
+
+Hydraulics parameters:
+    - Holtzman, Nataniel, et al. 2023, https://doi.org/10.1029/2023WR035481
+"""
+function FluxnetSimulations.get_parameters(
+    FT,
+    site_ID,
+    domain,
+    pft::Pft = ClimaLand.Canopy.clm_dominant_pft(domain.space.surface),
+    ;
+    soil_params_gupta = ClimaLand.Soil.soil_vangenuchten_parameters(
+        domain.space.subsurface,
+        FT,
+    ),
+    soil_params_grids = ClimaLand.Soil.soil_composition_parameters(
+        domain.space.subsurface,
+        FT,
+    ),
+    soil_ν = ClimaCore.Fields.field2array(soil_params_gupta[:ν])[1],
+    soil_K_sat = ClimaCore.Fields.field2array(soil_params_gupta[:K_sat])[1],
+    soil_S_s = FT(1e-2),
+    soil_hydrology_cm = ClimaCore.MatrixFields.column_field2array(
+        soil_params_gupta[:hydrology_cm],
+    )[1],
+    θ_r = ClimaCore.Fields.field2array(soil_params_gupta[:θ_r])[1],
+    ν_ss_quartz = ClimaCore.Fields.field2array(
+        soil_params_grids[:ν_ss_quartz],
+    )[1],
+    ν_ss_om = ClimaCore.Fields.field2array(soil_params_grids[:ν_ss_om])[1],
+    ν_ss_gravel = ClimaCore.Fields.field2array(
+        soil_params_grids[:ν_ss_gravel],
+    )[1],
+    z_0m_soil = FT(0.01),
+    z_0b_soil = FT(0.01),
+    soil_ϵ = FT(0.98),
+    soil_albedo_params = ClimaLand.Soil.clm_soil_albedo_parameters(
+        domain.space.surface,
+    ),
+    soil_albedo = ClimaLand.Soil.CLMTwoBandSoilAlbedo{FT}(;
+        NamedTuple{keys(soil_albedo_params)}((
+            ClimaCore.Fields.field2array(v)[1] for
+            v in values(soil_albedo_params)
+        ))...,
+    ),
+    Ω = pft.parameters.Ω,
+    χl = pft.parameters.χl,
+    G_Function = CLMGFunction(χl),
+    α_PAR_leaf = pft.parameters.α_PAR_leaf,
+    λ_γ_PAR = FT(5e-7),
+    α_NIR_leaf = pft.parameters.α_NIR_leaf,
+    τ_PAR_leaf = pft.parameters.τ_PAR_leaf,
+    τ_NIR_leaf = pft.parameters.τ_NIR_leaf,
+    ϵ_canopy = pft.parameters.ϵ_canopy,
+    ac_canopy = pft.parameters.ac_canopy,
+    g1 = pft.parameters.g1,
+    Drel = FT(1.6),
+    g0 = FT(0.0),
+    Vcmax25 = pft.parameters.Vcmax25,
+    pc = FT(-2.0e6),
+    sc = FT(5e-6),
+    SAI = FT(0),
+    f_root_to_shoot = pft.parameters.f_root_to_shoot,
+    K_sat_plant = pft.parameters.K_sat_plant,
+    ψ63 = pft.parameters.ψ63,
+    Weibull_param = FT(4),
+    a = FT(0.1 * 0.0098),
+    conductivity_model = PlantHydraulics.Weibull{FT}(
+        K_sat_plant,
+        ψ63,
+        Weibull_param,
+    ),
+    retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a),
+    plant_ν = pft.parameters.plant_ν,
+    plant_S_s = FT(1e-2 * 0.0098),
+    rooting_depth = pft.parameters.rooting_depth,
+    n_stem = Int64(0),
+    n_leaf = Int64(1),
+    h_leaf = ClimaCore.Fields.field2array(
+        ClimaLand.Canopy.clm_canopy_height(domain.space.surface),
+    )[1],
+    h_canopy = FluxnetSimulationsExt.get_canopy_height(site_ID),
+    h_stem = ((h_canopy - h_leaf)) > 0 ? h_canopy - h_leaf : FT(0.0),
+)
+    return (;
+        soil_ν,
+        soil_K_sat,
+        soil_S_s,
+        soil_hydrology_cm,
+        θ_r,
+        ν_ss_quartz,
+        ν_ss_om,
+        ν_ss_gravel,
+        z_0m_soil,
+        z_0b_soil,
+        soil_ϵ,
+        soil_albedo,
+        Ω,
+        χl,
+        G_Function,
+        α_PAR_leaf,
+        λ_γ_PAR,
+        τ_PAR_leaf,
+        α_NIR_leaf,
+        τ_NIR_leaf,
+        ϵ_canopy,
+        ac_canopy,
+        g1,
+        Drel,
+        g0,
+        Vcmax25,
+        SAI,
+        f_root_to_shoot,
+        K_sat_plant,
+        ψ63,
+        Weibull_param,
+        a,
+        conductivity_model,
+        retention_model,
+        plant_ν,
+        plant_S_s,
+        rooting_depth,
+        n_stem,
+        n_leaf,
+        h_leaf,
+        h_stem,
+        h_canopy,
+    )
+end
 
 ###################################
 #            UTILITIES            #
diff --git a/ext/fluxnet_simulations/get_fluxnet_metadata.jl b/ext/fluxnet_simulations/get_fluxnet_metadata.jl
new file mode 100644
index 00000000000..8b6c1981b3a
--- /dev/null
+++ b/ext/fluxnet_simulations/get_fluxnet_metadata.jl
@@ -0,0 +1,121 @@
+using DelimitedFiles
+using ClimaLand
+
+function parse_array_field(field)
+    if field == "" || field == "NaN"
+        return Float64[]
+    elseif field isa Float64
+        return isnan(field) ? Float64[] : [field]
+    else
+        return parse.(Float64, split(String(field), ";"))
+    end
+end
+
+
+"""
+    get_site_info(site_ID)
+
+This function retrieves the site information from an artifact for a given fluxnet site ID.
+
+Returns a Float64-valued namedtuple containing:
+- `lat`: Latitude of the site (deg)
+- `long`: Longitude of the site (deg)
+- `time_offset`: Time offset from UTC (hours). Note that the ClimaLand convention is for
+    UTC = local_time + time_offset.
+- `atmospheric_sensor_height` (Vector): Heights of atmospheric sensors (m)
+"""
+function get_site_info(site_ID; fluxnet2015_metadata_path = nothing)
+    if fluxnet2015_metadata_path === nothing
+        # Use the default path from ClimaLand.Artifacts
+        fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path()
+        fluxnet2015_metadata_path =
+            joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv")
+    end
+    raw = readdlm(fluxnet2015_metadata_path, ',', Any)
+    header = raw[1, :]
+    data = raw[2:end, :]
+
+    # Find the row matching site_ID
+    row_idx = findfirst(row -> row[1] == site_ID, eachrow(data))
+    if isnothing(row_idx)
+        error("Site ID $site_ID not found in metadata.")
+    end
+
+    site_metadata = data[row_idx, :]
+
+    varnames =
+        ["latitude", "longitude", "utc_offset", "atmospheric_sensor_heights"]
+    column_name_map = Dict(
+        varname => findfirst(header[:] .== varname) for varname in varnames
+    )
+
+    nan_if_empty(x) =
+        (x isa String && isempty(x)) ? NaN :
+        (x isa AbstractFloat && isnan(x)) ? NaN : x
+
+    for (varname, column) in column_name_map
+        field = site_metadata[column]
+        if field == "" || (field isa AbstractFloat && isnan(field))
+            @warn "Field $(varname) is missing for site $site_ID"
+        end
+    end
+
+    return (;
+        lat = nan_if_empty(site_metadata[column_name_map["latitude"]]),
+        long = nan_if_empty(site_metadata[column_name_map["longitude"]]),
+        time_offset = Int64(
+            -1 * nan_if_empty(site_metadata[column_name_map["utc_offset"]]),
+        ),
+        atmospheric_sensor_height = parse_array_field(
+            site_metadata[column_name_map["atmospheric_sensor_heights"]],
+        ),
+    )
+end
+
+
+"""
+    get_canopy_height(site_ID)
+
+This function retrieves a canopy height from an artifact for a given fluxnet site ID.
+
+Returns a Float64-valued namedtuple containing:
+- `canopy_height`: Canopy height of the site (m)
+"""
+function get_canopy_height(site_ID; fluxnet2015_metadata_path = nothing)
+    if fluxnet2015_metadata_path === nothing
+        # Use the default path from ClimaLand.Artifacts
+        fluxnet2015_data_path = ClimaLand.Artifacts.fluxnet2015_data_path()
+        fluxnet2015_metadata_path =
+            joinpath(fluxnet2015_data_path, "metadata_DD_clean.csv")
+    end
+
+    raw = readdlm(fluxnet2015_metadata_path, ',', Any)
+    header = raw[1, :]
+    data = raw[2:end, :]
+
+    site_info = get_site_info(site_ID; fluxnet2015_metadata_path)
+
+    # Find the row matching site_ID
+    row_idx = findfirst(row -> row[1] == site_ID, eachrow(data))
+    if isnothing(row_idx)
+        error("Site ID $site_ID not found in metadata.")
+    end
+
+    site_metadata = data[row_idx, :]
+
+    varname = "canopy_height"
+
+    col_idx = findfirst(header[:] .== varname)
+
+    nan_if_empty(x) =
+        (x isa String && isempty(x)) ? NaN :
+        (x isa AbstractFloat && isnan(x)) ? NaN : x
+
+    canopy_height = site_metadata[col_idx]
+    if canopy_height == "" ||
+       (canopy_height isa AbstractFloat && isnan(canopy_height))
+        @warn "Field $(varname) is missing for site $site_ID"
+    end
+
+    return nan_if_empty(canopy_height)
+end
diff --git a/src/Artifacts.jl b/src/Artifacts.jl
index 0d616b367fc..de6df00ee44 100644
--- a/src/Artifacts.jl
+++ b/src/Artifacts.jl
@@ -310,11 +310,48 @@ Citation: Siyan Ma, Liukang Xu, Joseph Verfaillie, Dennis Baldocchi (2023), Amer
 AmeriFlux CC-BY-4.0 License
 """
 function experiment_fluxnet_data_path(site_ID; context = nothing)
-    @assert site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")
+    try
+        full_fluxnet_path = @clima_artifact("fluxnet2015", context)
+        dirs = filter(
+            d -> isdir(joinpath(full_fluxnet_path, d)),
+            readdir(full_fluxnet_path),
+        )
+
+        match = findfirst(d -> occursin(site_ID, d), dirs)
+        @assert match !== nothing "No Fluxnet data found for site ID: $site_ID"
+
+        site_dir = dirs[match]
+        parts = split(site_dir, "FULLSET")
+        site_path_hh = join([parts[1], "FULLSET_HH", parts[2]], "") * ".csv"
+        site_path_hr = join([parts[1], "FULLSET_HR", parts[2]], "") * ".csv"
+
+        if isfile(joinpath(full_fluxnet_path, site_dir, site_path_hh))
+            data_path = joinpath(full_fluxnet_path, site_dir, site_path_hh)
+        elseif isfile(joinpath(full_fluxnet_path, site_dir, site_path_hr))
+            data_path = joinpath(full_fluxnet_path, site_dir, site_path_hr)
+        else
+            error(
+                "There exists a directory $site_dir for site ID $site_ID, but found no data files for half-hourly or hourly data.",
+            )
+        end
+
+        return data_path
+    catch
+        @info "Either the full fluxnet2015 dataset does not exist locally, or the site ID was not found.
+            Falling back to the fluxnet_sites artifact which contains only US-MOz, US-Var, US-NR1, and US-Ha1."
+
+        @assert site_ID ∈ ("US-MOz", "US-Var", "US-NR1", "US-Ha1")
+
+        folder_path = @clima_artifact("fluxnet_sites", context)
+        data_path = joinpath(folder_path, "$(site_ID).csv")
+        return data_path
+    end
+
+end
+
 
-    folder_path = @clima_artifact("fluxnet_sites", context)
-    data_path = joinpath(folder_path, "$(site_ID).csv")
-    return data_path
+function fluxnet2015_data_path(; context = nothing)
+    return @clima_artifact("fluxnet2015", context)
 end
 
 """
diff --git a/src/standalone/Vegetation/spatially_varying_parameters.jl b/src/standalone/Vegetation/spatially_varying_parameters.jl
index 56a685cc373..b75ccb71e14 100644
--- a/src/standalone/Vegetation/spatially_varying_parameters.jl
+++ b/src/standalone/Vegetation/spatially_varying_parameters.jl
@@ -354,3 +354,59 @@ function clm_canopy_height(
         return min.(canopy_height, max_height)
     end
 end
+
+"""
+    clm_dominant_pft(
+        surface_space;
+        regridder_type = :InterpolationsRegridder,
+        extrapolation_bc = (
+            Interpolations.Periodic(),
+            Interpolations.Flat(),
+            Interpolations.Flat(),
+        ),
+        interpolation_method = Interpolations.Constant(),
+        lowres = ClimaLand.Domains.use_lowres_clm(surface_space),
+    )
+
+Reads spatially varying dominant PFT, from a NetCDF file
+based on CLM data, and regrids it to the grid defined by the
+`surface_space` of the Clima simulation. Returns a Pft() object.
+
+The values correspond to the value of the dominant PFT at each point.
+
+The NetCDF files are stored in ClimaArtifacts and more detail on their origin
+is provided there. The keyword arguments `regridder_type` and `extrapolation_bc`
+affect the regridding by (1) changing how we interpolate to ClimaCore points which
+are not in the data, and (2) changing how extrapolate to points beyond the range of the
+data. The keyword argument lowres is a flag that determines if the 0.9x1.25 or 0.125x0.125
+resolution CLM data artifact is used. If the lowres flag is not provided, the clm artifact
+with the closest resolution to the surface_space is used.
+"""
+function clm_dominant_pft(
+    surface_space;
+    regridder_type = :InterpolationsRegridder,
+    extrapolation_bc = (
+        Interpolations.Periodic(),
+        Interpolations.Flat(),
+        Interpolations.Flat(),
+    ),
+    interpolation_method = Interpolations.Constant(),
+    lowres = ClimaLand.Domains.use_lowres_clm(surface_space),
+)
+    context = ClimaComms.context(surface_space)
+    clm_artifact_path = Artifacts.clm_data_folder_path(; context, lowres)
+
+    dominant_pft_idx = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "dominant_PFT_map.nc"),
+        "dominant_PFT",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+
+    dominant_pft_idx =
+        trunc(Int64, ClimaCore.Fields.field2array(dominant_pft_idx)[1])
+
+    dominant_pft = ClimaLand.default_pfts[dominant_pft_idx]
+    return dominant_pft
+end