diff --git a/docs/src/diagnostics/developers_diagnostics.md b/docs/src/diagnostics/developers_diagnostics.md
index a1ed75c1f1..c639d5d7dd 100644
--- a/docs/src/diagnostics/developers_diagnostics.md
+++ b/docs/src/diagnostics/developers_diagnostics.md
@@ -17,7 +17,7 @@ Internally, this is done by using the [`ClimaDiagnostics.jl`](https://github.com
  `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)`.
+ - `default_diagnostics.jl` defines default diagnostics functions to use on a model simulation. For example, `default_diagnostics(land_model::BucketModel, 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.
@@ -66,7 +66,7 @@ For each model, we define a function `default_diagnostics` which will define wha
 on what schedule (for example, hourly average). For example,
 
 ```Julia
-function default_diagnostics(land_model::BucketModel, t_start; output_writer)
+function default_diagnostics(land_model::BucketModel{FT}; output_writer) where {FT}
 
     define_diagnostics!(land_model)
 
@@ -87,7 +87,7 @@ function default_diagnostics(land_model::BucketModel, t_start; output_writer)
     ]
 
     default_outputs =
-        hourly_averages(bucket_diagnostics...; output_writer, t_start)
+        hourly_averages(FT, bucket_diagnostics...; output_writer)
     return [default_outputs...]
 end
 ```
@@ -103,11 +103,10 @@ If `average_period = :hourly`, `default_outputs` calls `hourly_averages`, et cet
 We defined some functions of diagnostic schedule that may often be used in `standard_diagnostic_frequencies.jl`, for example
 
 ```Julia
-hourly_averages(short_names...; output_writer, t_start) = common_diagnostics(
-    60 * 60 * one(t_start),
+hourly_averages(FT, short_names...; output_writer) = common_diagnostics(
+    60 * 60 * one(FT),
     (+),
     output_writer,
-    t_start,
     short_names...;
     pre_output_hook! = average_pre_output_hook!,
 )
diff --git a/docs/src/diagnostics/users_diagnostics.md b/docs/src/diagnostics/users_diagnostics.md
index 143ca80922..ea54a10987 100644
--- a/docs/src/diagnostics/users_diagnostics.md
+++ b/docs/src/diagnostics/users_diagnostics.md
@@ -46,11 +46,11 @@ providing the space and output_dir defined in steps 1. and 2.
 Now that you defined your model and your writter, you can create a callback function to be called when solving your model. For example:
 
 ```Julia
-t0 = 0 # the starting time of your simulation
+t0 = 0 # the start date of your simulation
 
-reference_date = DateTime(2024) # reference_date is the DateTime of your starting time
+start_date = DateTime(2024) # start_date is the DateTime of your start date
 
-diags = ClimaLand.default_diagnostics(model, t0, reference_date; output_writer = nc_writer)
+diags = ClimaLand.default_diagnostics(model, start_date; output_writer = nc_writer)
 
 diagnostic_handler =
     ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt)
@@ -118,11 +118,10 @@ add_diagnostic_variable!(
 ### Define how to schedule your variables. For example, you want the seasonal maximum of your variables, where season is defined as 90 days.
 
 ```Julia
-seasonal_maxs(short_names...; output_writer, t_start) = common_diagnostics(
-    90 * 24 * 60 * 60 * one(t_start),
+seasonal_maxs(FT, short_names...; output_writer) = common_diagnostics(
+    90 * 24 * 60 * 60 * one(FT),
     max,
     output_writer,
-    t_start,
     short_names...,
 )
 ```
@@ -134,7 +133,7 @@ Now, you can call your schedule with your variables.
 ```Julia
 my_custom_diagnostics = ["lhf", "bor"]
 
-diags = seasonal_maxs(my_custom_diagnostics...; output_writer, t_start)
+diags = seasonal_maxs(FT, my_custom_diagnostics...; output_writer)
 ```
 
 ### Analyze your simulation output
diff --git a/docs/tutorials/shared_utilities/driver_tutorial.jl b/docs/tutorials/shared_utilities/driver_tutorial.jl
index 2f87c0557a..0708a325cc 100644
--- a/docs/tutorials/shared_utilities/driver_tutorial.jl
+++ b/docs/tutorials/shared_utilities/driver_tutorial.jl
@@ -24,12 +24,12 @@
 # shortwave and longwave flux (W/m^2). The radiative driver is also where
 # a function which computes the zenith angle for the site is stored.
 
-# Both drivers store the reference time for the data/simulation.
+# Both drivers store the start date for the data/simulation.
 # This is the DateTime object which corresponds to the time at which t=0
 # in the simulation. Additionally, for site-level runs, both drivers store the
 # forcing data as a spline function fit to the data which takes the time
 # `t` as an argument, where `t` is the simulation time measured in seconds since
-# the reference time. The reference time should be in UTC.
+# the start date. The start date should be in UTC.
 
 # Note: for coupled runs, corresponding types `CoupledAtmosphere`
 # and `CoupledRadiativeFluxes` exist. However, these are not defined
@@ -54,9 +54,9 @@ time_offset = 7;
 # Site latitude and longitude
 lat = 38.7441; # degree
 long = -92.2000; # degree
-# Compute the reference time in UTC, and convert local datetime
-# vector into a vector of seconds since the reference time
-ref_time = local_datetime[1] + Dates.Hour(time_offset);
+# Compute the start date in UTC, and convert local datetime
+# vector into a vector of seconds since the start date
+start_date = local_datetime[1] + Dates.Hour(time_offset);
 data_dt = 3600.0;
 seconds = 0:data_dt:((length(local_datetime) - 1) * data_dt);
 
@@ -81,16 +81,16 @@ earth_param_set = LP.LandParameters(Float64);
 insol_params = earth_param_set.insol_params # parameters of Earth's orbit required to compute the insolation
 function zenith_angle(
     t,
-    ref_time;
+    start_date;
     latitude = lat,
     longitude = long,
     insol_params = insol_params,
 )
-    current_datetime = ref_time + Dates.Second(round(t)) # Time in UTC
+    current_datetime = start_date + Dates.Second(round(t)) # Time in UTC
 
     d, δ, η_UTC = (Insolation.helper_instantaneous_zenith_angle(
         current_datetime,
-        ref_time,
+        start_date,
         insol_params,
     ))
 
@@ -110,7 +110,7 @@ radiation = ClimaLand.PrescribedRadiativeFluxes(
     Float64,
     SW_d,
     LW_d,
-    ref_time;
+    start_date;
     θs = zenith_angle,
 );
 
diff --git a/docs/tutorials/standalone/Bucket/bucket_tutorial.jl b/docs/tutorials/standalone/Bucket/bucket_tutorial.jl
index ab2fad4df1..c9dbb2b98a 100644
--- a/docs/tutorials/standalone/Bucket/bucket_tutorial.jl
+++ b/docs/tutorials/standalone/Bucket/bucket_tutorial.jl
@@ -156,7 +156,7 @@ using ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
 # We also want to plot the solution
 using Plots
 
-# And we need to use the DateTime type to store reference times
+# And we need to use the DateTime type to store start dates
 using Dates
 
 FT = Float32;
@@ -203,7 +203,7 @@ z_0b = FT(1e-3);
 ρc_soil = FT(2e6);
 # Snow melt timescale
 τc = FT(3600);
-# Simulation start time, end time, and timestep
+# Simulation start date, end time, and timestep
 t0 = 0.0;
 tf = 7 * 86400;
 Δt = 3600.0;
@@ -214,7 +214,7 @@ bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc);
 # The PrescribedAtmosphere and PrescribedRadiation need to take in a reference
 # time, the date of the start of the simulation. In this tutorial we will
 # consider this January 1, 2005.
-ref_time = DateTime(2005);
+start_date = DateTime(2005);
 
 # To drive the system in standalone mode,
 # the user must provide
@@ -247,7 +247,7 @@ bucket_atmos = PrescribedAtmosphere(
     TimeVaryingInput(u_atmos),
     TimeVaryingInput(q_atmos),
     TimeVaryingInput(P_atmos),
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set,
 );
@@ -262,7 +262,7 @@ bucket_rad = PrescribedRadiativeFluxes(
     FT,
     TimeVaryingInput(SW_d),
     TimeVaryingInput(LW_d),
-    ref_time,
+    start_date,
 );
 
 
diff --git a/docs/tutorials/standalone/Soil/evaporation.jl b/docs/tutorials/standalone/Soil/evaporation.jl
index ff0052d3fc..6f07ae3c29 100644
--- a/docs/tutorials/standalone/Soil/evaporation.jl
+++ b/docs/tutorials/standalone/Soil/evaporation.jl
@@ -39,14 +39,14 @@ thermo_params = LP.thermodynamic_parameters(earth_param_set);
 # Our assumption is that in the lab experiment there was no radiative heating
 # or cooling of the soil.
 
-ref_time = DateTime(2005) # required argument, but not used in this case
+start_date = DateTime(2005) # required argument, but not used in this case
 SW_d = (t) -> 0
 LW_d = (t) -> 301.15^4 * 5.67e-8
 radiation = PrescribedRadiativeFluxes(
     FT,
     TimeVaryingInput(SW_d),
     TimeVaryingInput(LW_d),
-    ref_time,
+    start_date,
 );
 # Set up atmospheric conditions that result in the potential evaporation
 # rate obsereved in the experiment.
@@ -74,7 +74,7 @@ atmos = PrescribedAtmosphere(
     TimeVaryingInput(u_atmos),
     TimeVaryingInput(q_atmos),
     TimeVaryingInput(P_atmos),
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set;
     gustiness = gustiness,
diff --git a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl
index 08bcc86011..476ae5a17c 100644
--- a/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl
+++ b/docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl
@@ -71,14 +71,14 @@ params = ClimaLand.Soil.EnergyHydrologyParameters(
     d_ds,
 );
 
-ref_time = DateTime(2005)
+start_date = DateTime(2005)
 SW_d = (t) -> 0
 LW_d = (t) -> 294.15^4 * 5.67e-8
 radiation = PrescribedRadiativeFluxes(
     FT,
     TimeVaryingInput(SW_d),
     TimeVaryingInput(LW_d),
-    ref_time,
+    start_date,
 )
 # Atmos
 T_air = FT(301.15)
@@ -104,7 +104,7 @@ atmos = PrescribedAtmosphere(
     TimeVaryingInput(u_atmos),
     TimeVaryingInput(q_atmos),
     TimeVaryingInput(P_atmos),
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set;
     gustiness = gustiness,
diff --git a/docs/tutorials/standalone/Soil/sublimation.jl b/docs/tutorials/standalone/Soil/sublimation.jl
index f90fa3820a..48b4a35ede 100644
--- a/docs/tutorials/standalone/Soil/sublimation.jl
+++ b/docs/tutorials/standalone/Soil/sublimation.jl
@@ -25,14 +25,14 @@ FT = Float64;
 earth_param_set = LP.LandParameters(FT)
 thermo_params = LP.thermodynamic_parameters(earth_param_set);
 
-ref_time = DateTime(2005)
+start_date = DateTime(2005)
 SW_d = (t) -> 0
 LW_d = (t) -> 270.0^4 * 5.67e-8
 radiation = PrescribedRadiativeFluxes(
     FT,
     TimeVaryingInput(SW_d),
     TimeVaryingInput(LW_d),
-    ref_time,
+    start_date,
 );
 # Set up atmospheric conditions that result in the potential evaporation
 # rate obsereved in the experiment.
@@ -60,7 +60,7 @@ atmos = PrescribedAtmosphere(
     TimeVaryingInput(u_atmos),
     TimeVaryingInput(q_atmos),
     TimeVaryingInput(P_atmos),
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set;
     gustiness = gustiness,
diff --git a/experiments/benchmarks/bucket.jl b/experiments/benchmarks/bucket.jl
index 2bb4f933ac..7ec653bd36 100644
--- a/experiments/benchmarks/bucket.jl
+++ b/experiments/benchmarks/bucket.jl
@@ -56,7 +56,7 @@ function setup_prob(t0, tf, Δt; nelements = (100, 10))
         npolynomial = 1,
         dz_tuple = FT.((1.0, 0.05)),
     )
-    ref_time = DateTime(2005)
+    start_date = DateTime(2005)
 
     # Initialize parameters
     σS_c = FT(0.2)
@@ -69,7 +69,7 @@ function setup_prob(t0, tf, Δt; nelements = (100, 10))
 
     surface_space = bucket_domain.space.surface
     # Construct albedo parameter object using temporal map
-    albedo = PrescribedSurfaceAlbedo{FT}(ref_time, t0, surface_space)
+    albedo = PrescribedSurfaceAlbedo{FT}(start_date, surface_space)
 
     bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc)
 
@@ -91,7 +91,7 @@ function setup_prob(t0, tf, Δt; nelements = (100, 10))
         TimeVaryingInput(u_atmos),
         TimeVaryingInput(q_atmos),
         TimeVaryingInput(P_atmos),
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
@@ -106,7 +106,7 @@ function setup_prob(t0, tf, Δt; nelements = (100, 10))
         FT,
         TimeVaryingInput(SW_d),
         TimeVaryingInput(LW_d),
-        ref_time,
+        start_date,
     )
 
     model = BucketModel(
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 4a27d246dd..6024e0af00 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -71,8 +71,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
 
-    ref_time = DateTime(2021)
-    t_start = 0.0
+    start_date = DateTime(2021)
 
     # Forcing data
     era5_artifact_path =
@@ -81,8 +80,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "rf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -91,8 +89,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -101,24 +98,21 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "ws",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "q",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sp",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
 
@@ -126,8 +120,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "t2m",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     h_atmos = FT(10)
@@ -139,7 +132,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         u_atmos,
         q_atmos,
         P_atmos,
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
@@ -149,8 +142,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "ssrd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
@@ -158,28 +150,27 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "strd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
 
     function zenith_angle(
         t,
-        ref_time;
+        start_date;
         latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
         longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
         insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
     ) where {FT}
         # This should be time in UTC
-        current_datetime = ref_time + Dates.Second(round(t))
+        current_datetime = start_date + Dates.Second(round(t))
 
         # Orbital Data uses Float64, so we need to convert to our sim FT
         d, δ, η_UTC =
             FT.(
                 Insolation.helper_instantaneous_zenith_angle(
                     current_datetime,
-                    ref_time,
+                    start_date,
                     insol_params,
                 )
             )
@@ -193,7 +184,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         ).:1
     end
     radiation =
-        PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle)
+        PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date; θs = zenith_angle)
 
     soil_params_artifact_path =
         ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
@@ -492,8 +483,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
         "lai",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (;
             preprocess_func = (data) -> data > 0.05 ? data : 0.0,
diff --git a/experiments/benchmarks/richards.jl b/experiments/benchmarks/richards.jl
index b778c225bf..3174ebb2f1 100644
--- a/experiments/benchmarks/richards.jl
+++ b/experiments/benchmarks/richards.jl
@@ -74,7 +74,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
 
-    ref_time = DateTime(2021)
+    start_date = DateTime(2021)
 
     # Read in f_max data and land sea mask
     infile_path = ClimaLand.Artifacts.topmodel_data_path()
@@ -195,15 +195,13 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     # Below, the preprocess_func argument is used to
     # 1. Convert precipitation to be negative (as it is downwards)
     # 2. Convert accumulations over an hour to a rate per second
-    ref_time = DateTime(2021)
-    t_start = 0.0
+    start_date = DateTime(2021)
     # Precipitation:
     precip = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "tp",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
diff --git a/experiments/integrated/fluxnet/met_drivers_FLUXNET.jl b/experiments/integrated/fluxnet/met_drivers_FLUXNET.jl
index fae2015ac2..4eec1e2392 100644
--- a/experiments/integrated/fluxnet/met_drivers_FLUXNET.jl
+++ b/experiments/integrated/fluxnet/met_drivers_FLUXNET.jl
@@ -124,7 +124,7 @@ LW_IN = TimeVaryingInput(seconds, FT.(drivers.LW_IN.values[:]); context)
 SW_IN = TimeVaryingInput(seconds, FT.(drivers.SW_IN.values[:]); context)
 snow_precip = TimeVaryingInput((t) -> FT(0))
 
-# Construct the drivers. For the reference time we will use the UTC time at the
+# Construct the drivers. For the start date we will use the UTC time at the
 # start of the simulation
 atmos = ClimaLand.PrescribedAtmosphere(
     precip,
@@ -141,20 +141,20 @@ atmos = ClimaLand.PrescribedAtmosphere(
 
 function zenith_angle(
     t,
-    ref_time;
+    start_date;
     latitude = lat,
     longitude = long,
     insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
 ) where {FT}
     # This should be time in UTC
-    current_datetime = ref_time + Dates.Second(round(t))
+    current_datetime = start_date + Dates.Second(round(t))
 
     # Orbital Data uses Float64, so we need to convert to our sim FT
     d, δ, η_UTC =
         FT.(
             Insolation.helper_instantaneous_zenith_angle(
                 current_datetime,
-                ref_time,
+                start_date,
                 insol_params,
             )
         )
diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
index 86baac684b..010c7c663d 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -56,8 +56,7 @@ domain = ClimaLand.Domains.SphericalShell(;
 surface_space = domain.space.surface
 subsurface_space = domain.space.subsurface
 
-ref_time = DateTime(2021);
-t_start = 0.0
+start_date = DateTime(2021);
 
 # Forcing data
 era5_artifact_path =
@@ -66,8 +65,7 @@ precip = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "rf",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
 )
@@ -76,8 +74,7 @@ snow_precip = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "sf",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
 )
@@ -86,24 +83,21 @@ u_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "ws",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 q_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "q",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 P_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "sp",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 
@@ -111,8 +105,7 @@ T_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "t2m",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 h_atmos = FT(10);
@@ -124,7 +117,7 @@ atmos = PrescribedAtmosphere(
     u_atmos,
     q_atmos,
     P_atmos,
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set,
 );
@@ -134,8 +127,7 @@ SW_d = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "ssrd",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
 )
@@ -143,28 +135,27 @@ LW_d = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "strd",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
 )
 
 function zenith_angle(
     t,
-    ref_time;
+    start_date;
     latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
     longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
     insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
 ) where {FT}
     # This should be time in UTC
-    current_datetime = ref_time + Dates.Second(round(t))
+    current_datetime = start_date + Dates.Second(round(t))
 
     # Orbital Data uses Float64, so we need to convert to our sim FT
     d, δ, η_UTC =
         FT.(
             Insolation.helper_instantaneous_zenith_angle(
                 current_datetime,
-                ref_time,
+                start_date,
                 insol_params,
             )
         )
@@ -172,7 +163,7 @@ function zenith_angle(
     Insolation.instantaneous_zenith_angle.(d, δ, η_UTC, longitude, latitude).:1
 end
 radiation =
-    PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle);
+    PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date; θs = zenith_angle);
 
 include(
     joinpath(
@@ -240,8 +231,7 @@ LAIfunction = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
     "lai",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (;
         preprocess_func = (data) -> data > 0.05 ? data : 0.0,
@@ -374,8 +364,7 @@ nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, output_dir)
 
 diags = ClimaLand.default_diagnostics(
     land,
-    t0,
-    ref_time;
+    start_date;
     output_writer = nc_writer,
     average_period = :hourly,
 )
diff --git a/experiments/integrated/performance/conservation/ozark_conservation.jl b/experiments/integrated/performance/conservation/ozark_conservation.jl
index 91fcc08cee..3c40fd18a4 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation.jl
@@ -78,7 +78,7 @@ for float_type in (Float32, Float64)
         cache_Tair =
             [parent(sv.saveval[k].drivers.T)[1] for k in 1:length(sv.t)]
         @assert mean(
-            abs.(radiation.θs.(sv.t, radiation.ref_time) .- cache_θs),
+            abs.(radiation.θs.(sv.t, radiation.start_date) .- cache_θs),
         ) < eps(FT)
         T_mutable = Vector{FT}(undef, 1)
         atmos_T = map(sv.t) do time
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 0236b9248f..22c41dfea1 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -105,9 +105,9 @@ land_domain = ClimaLand.Domains.SphericalShell(;
 );
 canopy_domain = ClimaLand.Domains.obtain_surface_domain(land_domain)
 sfc_cds = ClimaCore.Fields.coordinate_field(land_domain.space.surface)
-# First pick the reference time and start time of the simulation, since time varying input depends on that
-t0 = Float64(0)# start at reference time
-ref_time = DateTime("202001010000", "yyyymmddHHMM")
+# First pick the start date and time of the simulation, since time varying input depends on that
+t0 = Float64(0) # start at start date
+start_date = DateTime("202001010000", "yyyymmddHHMM")
 # Time varying input
 LAIfunction = TimeVaryingInput((t) -> 2.0)
 # Atmospheric and radiative forcing
@@ -120,7 +120,7 @@ LW_IN = TimeVaryingInput((t) -> 5.67e-8 * 298.0^4)
 SW_IN = TimeVaryingInput((t) -> 500.0)
 snow_precip = TimeVaryingInput((t) -> 0.0)
 atmos_h = FT(32)
-# Construct the drivers. For the reference time we will use the UTC time at the
+# Construct the drivers. For the start date we will use the UTC time at the
 # start of the simulation
 atmos = ClimaLand.PrescribedAtmosphere(
     precip,
@@ -129,25 +129,25 @@ atmos = ClimaLand.PrescribedAtmosphere(
     atmos_u,
     atmos_q,
     atmos_p,
-    ref_time,
+    start_date,
     atmos_h,
     earth_param_set;
 )
 function zenith_angle(
     t,
-    ref_time;
+    start_date;
     cd_field = sfc_cds,
     insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
 ) where {FT}
     # This should be time in UTC
-    current_datetime = ref_time + Dates.Second(round(t))
+    current_datetime = start_date + Dates.Second(round(t))
 
     # Orbital Data uses Float64, so we need to convert to our sim FT
     d, δ, η_UTC =
         FT.(
             Insolation.helper_instantaneous_zenith_angle(
                 current_datetime,
-                ref_time,
+                start_date,
                 insol_params,
             )
         )
@@ -165,7 +165,7 @@ radiation = ClimaLand.PrescribedRadiativeFluxes(
     FT,
     SW_IN,
     LW_IN,
-    ref_time,
+    start_date,
     θs = zenith_angle,
 )
 
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index 46163bf2ee..466d75200d 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -72,8 +72,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
 
-    ref_time = DateTime(2021)
-    t_start = t0
+    start_date = DateTime(2021)
     # Forcing data
     era5_artifact_path =
         ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context)    # Precipitation:
@@ -81,8 +80,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "rf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -91,8 +89,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -101,24 +98,21 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "ws",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "q",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sp",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
 
@@ -126,8 +120,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "t2m",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     h_atmos = FT(10)
@@ -139,7 +132,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         u_atmos,
         q_atmos,
         P_atmos,
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
@@ -149,8 +142,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "ssrd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
@@ -158,28 +150,27 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "strd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
 
     function zenith_angle(
         t,
-        ref_time;
+        start_date;
         latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
         longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
         insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
     ) where {FT}
         # This should be time in UTC
-        current_datetime = ref_time + Dates.Second(round(t))
+        current_datetime = start_date + Dates.Second(round(t))
 
         # Orbital Data uses Float64, so we need to convert to our sim FT
         d, δ, η_UTC =
             FT.(
                 Insolation.helper_instantaneous_zenith_angle(
                     current_datetime,
-                    ref_time,
+                    start_date,
                     insol_params,
                 )
             )
@@ -193,7 +184,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         ).:1
     end
     radiation =
-        PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle)
+        PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date; θs = zenith_angle)
 
     soil_params_artifact_path =
         ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
@@ -493,8 +484,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
         "lai",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (;
             preprocess_func = (data) -> data > 0.05 ? data : 0.0,
@@ -620,8 +610,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
 
     diags = ClimaLand.default_diagnostics(
         land,
-        t0,
-        ref_time;
+        start_date;
         output_writer = nc_writer,
         output_vars = :short,
     )
diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
index 66f22c3acd..c603ebfd30 100644
--- a/experiments/long_runs/soil.jl
+++ b/experiments/long_runs/soil.jl
@@ -70,8 +70,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
 
-    ref_time = DateTime(2021)
-    t_start = t0
+    start_date = DateTime(2021)
     # Forcing data
     era5_artifact_path =
         ClimaLand.Artifacts.era5_land_forcing_data2021_folder_path(; context)    # Precipitation:
@@ -79,8 +78,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "rf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -89,8 +87,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sf",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
     )
@@ -99,24 +96,21 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "ws",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "q",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sp",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
 
@@ -124,8 +118,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "t2m",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
     h_atmos = FT(10)
@@ -137,7 +130,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         u_atmos,
         q_atmos,
         P_atmos,
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
@@ -147,8 +140,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "ssrd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
@@ -156,28 +148,27 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "strd",
         surface_space;
-        reference_date = ref_time,
-        t_start,
+        reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
     )
 
     function zenith_angle(
         t,
-        ref_time;
+        start_date;
         latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
         longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
         insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
     ) where {FT}
         # This should be time in UTC
-        current_datetime = ref_time + Dates.Second(round(t))
+        current_datetime = start_date + Dates.Second(round(t))
 
         # Orbital Data uses Float64, so we need to convert to our sim FT
         d, δ, η_UTC =
             FT.(
                 Insolation.helper_instantaneous_zenith_angle(
                     current_datetime,
-                    ref_time,
+                    start_date,
                     insol_params,
                 )
             )
@@ -191,7 +182,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         ).:1
     end
     radiation =
-        PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = zenith_angle)
+        PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date; θs = zenith_angle)
 
     soil_params_artifact_path =
         ClimaLand.Artifacts.soil_params_artifact_folder_path(; context)
@@ -422,8 +413,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
 
     diags = ClimaLand.default_diagnostics(
         soil,
-        t0,
-        ref_time;
+        start_date;
         output_writer = nc_writer,
         average_period = :monthly,
     )
diff --git a/experiments/standalone/Bucket/bucket_era5.jl b/experiments/standalone/Bucket/bucket_era5.jl
index d76ab40de6..207b08b5ad 100644
--- a/experiments/standalone/Bucket/bucket_era5.jl
+++ b/experiments/standalone/Bucket/bucket_era5.jl
@@ -103,7 +103,7 @@ else
         dz_tuple = FT.((1.0, 0.05)),
     )
 end
-ref_time = DateTime(2021);
+start_date = DateTime(2021);
 
 # Set up parameters
 σS_c = FT(0.2);
@@ -122,7 +122,6 @@ tf = 14 * 86400;
 device_suffix =
     typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ?
     "cpu" : "gpu"
-t_start = t0
 surface_space = bucket_domain.space.surface
 α_snow = FT(0.8)
 albedo = PrescribedBaregroundAlbedo{FT}(α_snow, surface_space);
@@ -141,8 +140,7 @@ precip = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "rf",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
 )
@@ -151,8 +149,7 @@ snow_precip = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "sf",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
 )
@@ -161,24 +158,21 @@ u_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "ws",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 q_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
     "q",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 P_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "sp",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 
@@ -186,8 +180,7 @@ T_atmos = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "t2m",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
 )
 h_atmos = FT(10);
@@ -200,7 +193,7 @@ bucket_atmos = PrescribedAtmosphere(
     u_atmos,
     q_atmos,
     P_atmos,
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set,
 );
@@ -213,8 +206,7 @@ SW_d = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "ssrd",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
 )
@@ -222,13 +214,12 @@ LW_d = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "strd",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
 )
 
-bucket_rad = PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time);
+bucket_rad = PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date);
 
 model = BucketModel(
     parameters = bucket_parameters,
diff --git a/experiments/standalone/Bucket/global_bucket_function.jl b/experiments/standalone/Bucket/global_bucket_function.jl
index 71317c020d..aca5438400 100644
--- a/experiments/standalone/Bucket/global_bucket_function.jl
+++ b/experiments/standalone/Bucket/global_bucket_function.jl
@@ -97,7 +97,7 @@ tf = 7 * 86400;
 Δt = 3600.0;
 
 bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc);
-ref_time = DateTime(2005);
+start_date = DateTime(2005);
 
 # Precipitation:
 precip = (t) -> 0;
@@ -117,7 +117,7 @@ bucket_atmos = PrescribedAtmosphere(
     TimeVaryingInput(u_atmos),
     TimeVaryingInput(q_atmos),
     TimeVaryingInput(P_atmos),
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set,
 );
@@ -132,7 +132,7 @@ bucket_rad = PrescribedRadiativeFluxes(
     FT,
     TimeVaryingInput(SW_d),
     TimeVaryingInput(LW_d),
-    ref_time,
+    start_date,
 );
 
 model = BucketModel(
@@ -168,12 +168,8 @@ space = bucket_domain.space.subsurface
 
 nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(space, output_dir)
 
-diags = ClimaLand.default_diagnostics(
-    model,
-    t0,
-    ref_time;
-    output_writer = nc_writer,
-)
+diags =
+    ClimaLand.default_diagnostics(model, start_date; output_writer = nc_writer)
 
 diagnostic_handler =
     ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt)
diff --git a/experiments/standalone/Bucket/global_bucket_temporalmap.jl b/experiments/standalone/Bucket/global_bucket_temporalmap.jl
index d73d590ddb..09fbfd03f9 100644
--- a/experiments/standalone/Bucket/global_bucket_temporalmap.jl
+++ b/experiments/standalone/Bucket/global_bucket_temporalmap.jl
@@ -97,7 +97,7 @@ function setup_prob(t0, tf, Δt)
         npolynomial = 1,
         dz_tuple = FT.((1.0, 0.05)),
     )
-    ref_time = DateTime(2005)
+    start_date = DateTime(2005)
 
     # Initialize parameters
     σS_c = FT(0.2)
@@ -110,7 +110,7 @@ function setup_prob(t0, tf, Δt)
 
     surface_space = bucket_domain.space.surface
     # Construct albedo parameter object using temporal map
-    albedo = PrescribedSurfaceAlbedo{FT}(ref_time, t0, surface_space)
+    albedo = PrescribedSurfaceAlbedo{FT}(start_date, surface_space)
 
     bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc)
 
@@ -132,7 +132,7 @@ function setup_prob(t0, tf, Δt)
         TimeVaryingInput(u_atmos),
         TimeVaryingInput(q_atmos),
         TimeVaryingInput(P_atmos),
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
@@ -147,7 +147,7 @@ function setup_prob(t0, tf, Δt)
         FT,
         TimeVaryingInput(SW_d),
         TimeVaryingInput(LW_d),
-        ref_time,
+        start_date,
     )
 
 
diff --git a/experiments/standalone/Snow/process_snowmip.jl b/experiments/standalone/Snow/process_snowmip.jl
index f0c0435b0f..072132127a 100644
--- a/experiments/standalone/Snow/process_snowmip.jl
+++ b/experiments/standalone/Snow/process_snowmip.jl
@@ -51,23 +51,23 @@ snowf = -met_data["Snowf"][:][mask] ./ 1000.0 # convert to dS/dt
 SW_d = TimeVaryingInput(seconds, SWdown; context)
 LW_d = TimeVaryingInput(seconds, LWdown; context)
 
-ref_time = timestamp[mask][1]
+start_date = timestamp[mask][1]
 function zenith_angle(
     t,
-    ref_time;
+    start_date;
     latitude = FT(lat),
     longitude = FT(long),
     insol_params::Insolation.Parameters.InsolationParameters{FT} = param_set.insol_params,
 ) where {FT}
     # This should be time in UTC
-    current_datetime = ref_time + Dates.Second(round(t))
+    current_datetime = start_date + Dates.Second(round(t))
 
     # Orbital Data uses Float64, so we need to convert to our sim FT
     d, δ, η_UTC =
         FT.(
             Insolation.helper_instantaneous_zenith_angle(
                 current_datetime,
-                ref_time,
+                start_date,
                 insol_params,
             )
         )
@@ -87,7 +87,7 @@ radiation = ClimaLand.PrescribedRadiativeFluxes(
     FT,
     SW_d,
     LW_d,
-    ref_time;
+    start_date;
     θs = zenith_angle,
 )
 
@@ -109,7 +109,7 @@ atmos = PrescribedAtmosphere(
     u_atmos,
     q_atmos,
     P_atmos,
-    ref_time,
+    start_date,
     h_atmos,
     earth_param_set,
 )
diff --git a/experiments/standalone/Soil/richards_runoff.jl b/experiments/standalone/Soil/richards_runoff.jl
index 9f8ca9bf17..be56342490 100644
--- a/experiments/standalone/Soil/richards_runoff.jl
+++ b/experiments/standalone/Soil/richards_runoff.jl
@@ -157,15 +157,13 @@ era5_artifact_path =
 # Below, the preprocess_func argument is used to
 # 1. Convert precipitation to be negative (as it is downwards)
 # 2. Convert accumulations over an hour to a rate per second
-ref_time = DateTime(2021);
-t_start = 0.0
+start_date = DateTime(2021);
 # Precipitation:
 precip = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
     "tp",
     surface_space;
-    reference_date = ref_time,
-    t_start,
+    reference_date = start_date,
     regridder_type,
     file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
 )
diff --git a/lib/ClimaLandSimulations/experiments/ozark.jl b/lib/ClimaLandSimulations/experiments/ozark.jl
index a2e6a2679c..bd000d774c 100644
--- a/lib/ClimaLandSimulations/experiments/ozark.jl
+++ b/lib/ClimaLandSimulations/experiments/ozark.jl
@@ -7,7 +7,7 @@ sv_test, sol_test, Y_test, p_test = run_fluxnet(
     params = ozark_default_params(; hetero_resp = hetero_resp_ozark()),
 )
 
-# defaults, except start time
+# defaults, except start date
 sv, sol, Y, p = run_fluxnet(
     "US-MOz";
     setup = make_setup(;
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_drivers.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_drivers.jl
index 61e0a68b83..f68619597c 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_drivers.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_drivers.jl
@@ -123,7 +123,7 @@ function make_drivers(site_ID, setup, config, params, context)
     SW_IN = TimeVaryingInput(seconds, FT.(drivers.SW_IN.values[:]); context)
     snow_precip = TimeVaryingInput((t) -> FT(0))
 
-    # Construct the drivers. For the reference time we will use the UTC time at the
+    # Construct the drivers. For the start date we will use the UTC time at the
     # start of the simulation
     atmos = ClimaLand.PrescribedAtmosphere(
         precip,
@@ -140,20 +140,20 @@ function make_drivers(site_ID, setup, config, params, context)
 
     function zenith_angle(
         t,
-        ref_time;
+        start_date;
         latitude = config.lat,
         longitude = config.long,
         insol_params::Insolation.Parameters.InsolationParameters{FT} = earth_param_set.insol_params,
     ) where {FT}
         # This should be time in UTC
-        current_datetime = ref_time + Dates.Second(round(t))
+        current_datetime = start_date + Dates.Second(round(t))
 
         # Orbital Data uses Float64, so we need to convert to our sim FT
         d, δ, η_UTC =
             FT.(
                 Insolation.helper_instantaneous_zenith_angle(
                     current_datetime,
-                    ref_time,
+                    start_date,
                     insol_params,
                 )
             )
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
index 5f80beee14..bc9c562b63 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/fluxnet_utilities/make_setup.jl
@@ -10,7 +10,7 @@ setup site specific domain and time stepping:
 - n_leaf, the number of leaves
 - h_stem, the height of each stem (m)
 - h_leaf, the height of each leaf (m)
-- t0, the start time of the simulation (s)
+- t0, the start date of the simulation (s)
 - dt, the time step (s)
 - n, the number of time step between saving outputs
 """
diff --git a/lib/ClimaLandSimulations/test/runtests.jl b/lib/ClimaLandSimulations/test/runtests.jl
index 655115dee5..2021faf7eb 100644
--- a/lib/ClimaLandSimulations/test/runtests.jl
+++ b/lib/ClimaLandSimulations/test/runtests.jl
@@ -19,7 +19,7 @@ end
           @NamedTuple{t::Vector{Float64}, saveval::Vector{NamedTuple}}
 end
 
-@testset "Fluxnet single site, custom start time" begin
+@testset "Fluxnet single site, custom start date" begin
     sv, sol, Y, p = run_fluxnet(
         "US-MOz";
         setup = make_setup(;
diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
index 14a26b3c08..f9cf0c035d 100644
--- a/src/diagnostics/default_diagnostics.jl
+++ b/src/diagnostics/default_diagnostics.jl
@@ -14,7 +14,6 @@ export default_diagnostics
                                 period,
                                 reduction,
                                 output_writer,
-                                t_start,
                                 reference_date,
                                 short_names...;
                                 pre_output_hook! = nothing,
@@ -26,7 +25,6 @@ function common_diagnostics(
     period,
     reduction,
     output_writer,
-    t_start,
     reference_date,
     short_names...;
     pre_output_hook! = nothing,
@@ -35,8 +33,8 @@ function common_diagnostics(
         map(short_names) do short_name
             output_schedule_func =
                 period isa Period ?
-                EveryCalendarDtSchedule(period; t_start, reference_date) :
-                EveryDtSchedule(period; t_start)
+                EveryCalendarDtSchedule(period; reference_date) :
+                EveryDtSchedule(period)
             return ScheduledDiagnostic(
                 variable = get_diagnostic_variable(short_name),
                 compute_schedule_func = EveryStepSchedule(),
@@ -53,11 +51,10 @@ include("standard_diagnostic_frequencies.jl")
 
 # Bucket
 function default_diagnostics(
-    land_model::BucketModel,
-    t_start,
+    land_model::BucketModel{FT},
     reference_date;
     output_writer,
-)
+) where {FT}
 
     define_diagnostics!(land_model)
 
@@ -78,9 +75,9 @@ function default_diagnostics(
     ]
 
     default_outputs = hourly_averages(
+        FT,
         bucket_diagnostics...;
         output_writer,
-        t_start,
         reference_date,
     )
 
@@ -89,13 +86,12 @@ end
 
 # SoilCanopyModel
 function default_diagnostics(
-    land_model::SoilCanopyModel,
-    t_start,
+    land_model::SoilCanopyModel{FT},
     reference_date;
     output_writer,
     output_vars = :long,
     average_period = :daily,
-)
+) where {FT}
 
     define_diagnostics!(land_model)
 
@@ -163,23 +159,23 @@ function default_diagnostics(
 
     if average_period == :hourly
         default_outputs = hourly_averages(
+            FT,
             soilcanopy_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     elseif average_period == :daily
         default_outputs = daily_averages(
+            FT,
             soilcanopy_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     elseif average_period == :monthly
         default_outputs = monthly_averages(
+            FT,
             soilcanopy_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     end
@@ -190,12 +186,11 @@ end
 
 # SoilModel
 function default_diagnostics(
-    land_model::EnergyHydrology,
-    t_start,
+    land_model::EnergyHydrology{FT},
     reference_date;
     output_writer,
     average_period = :daily,
-)
+) where {FT}
 
     define_diagnostics!(land_model)
 
@@ -203,23 +198,23 @@ function default_diagnostics(
 
     if average_period == :hourly
         default_outputs = hourly_averages(
+            FT,
             soil_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     elseif average_period == :daily
         default_outputs = daily_averages(
+            FT,
             soil_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     elseif average_period == :monthly
         default_outputs = monthly_averages(
+            FT,
             soil_diagnostics...;
             output_writer,
-            t_start,
             reference_date,
         )
     end
diff --git a/src/diagnostics/land_compute_methods.jl b/src/diagnostics/land_compute_methods.jl
index 3532d949c8..41d9e33eb7 100644
--- a/src/diagnostics/land_compute_methods.jl
+++ b/src/diagnostics/land_compute_methods.jl
@@ -113,7 +113,7 @@ end
 
 ## Drivers Module ##
 
-@diagnostic_compute "soil_organic_carbon" SoilCanopyModel p.drivers.soc # need to fix this in src/shared_utilities/drivers 
+@diagnostic_compute "soil_organic_carbon" SoilCanopyModel p.drivers.soc # need to fix this in src/shared_utilities/drivers
 @diagnostic_compute "pressure" SoilCanopyModel p.drivers.P
 @diagnostic_compute "rainfall" SoilCanopyModel p.drivers.P_liq
 @diagnostic_compute "radiation_longwave_down" SoilCanopyModel p.drivers.LW_d
@@ -151,7 +151,7 @@ function compute_heterotrophic_respiration!(
     else
         out .= p.soilco2.top_bc .* FT(83.26)
     end
-end # Convert from kg C to mol CO2. 
+end # Convert from kg C to mol CO2.
 # To convert from kg C to mol CO2, we need to multiply by:
 # [3.664 kg CO2/ kg C] x [10^3 g CO2/ kg CO2] x [1 mol CO2/44.009 g CO2] = 83.26 mol CO2/kg C
 
diff --git a/src/diagnostics/standard_diagnostic_frequencies.jl b/src/diagnostics/standard_diagnostic_frequencies.jl
index b68b1ca02a..2e48a04e62 100644
--- a/src/diagnostics/standard_diagnostic_frequencies.jl
+++ b/src/diagnostics/standard_diagnostic_frequencies.jl
@@ -1,272 +1,262 @@
 """
-    monthly_maxs(short_names...; output_writer, t_start, reference_date)
+    monthly_maxs(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the monthly max for the given variables.
 """
-monthly_maxs(short_names...; output_writer, t_start, reference_date) =
-    common_diagnostics(Month(1), max, output_writer, t_start, short_names...)
+monthly_maxs(FT, short_names...; output_writer, reference_date) =
+    common_diagnostics(Month(1), max, output_writer, short_names...)
 
 """
-    monthly_max(short_names; output_writer, t_start, reference_date)
+    monthly_max(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the monthly max for the given variable.
 """
-monthly_max(short_names; output_writer, t_start, reference_date) =
-    monthly_maxs(short_names; output_writer, t_start, reference_date)[1]
+monthly_max(FT, short_names; output_writer, reference_date) =
+    monthly_maxs(FT, short_names; output_writer, reference_date)[1]
 
 """
-    monthly_mins(short_names...; output_writer, t_start, reference_date)
+    monthly_mins(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the monthly min for the given variables.
 """
-monthly_mins(short_names...; output_writer, t_start, reference_date) =
-    common_diagnostics(Month(1), min, output_writer, t_start, short_names...)
+monthly_mins(FT, short_names...; output_writer, reference_date) =
+    common_diagnostics(Month(1), min, output_writer, short_names...)
 
 """
-    monthly_min(short_names; output_writer, t_start, reference_date)
+    monthly_min(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the monthly min for the given variable.
 """
-monthly_min(short_names; output_writer, t_start, reference_date) =
-    monthly_mins(short_names; output_writer, t_start, reference_date)[1]
+monthly_min(FT, short_names; output_writer, reference_date) =
+    monthly_mins(FT, short_names; output_writer, reference_date)[1]
 
 """
-    monthly_averages(short_names...; output_writer, t_start, reference_date)
+    monthly_averages(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the monthly average for the given variables.
 """
 # An average is just a sum with a normalization before output
-monthly_averages(short_names...; output_writer, t_start, reference_date) =
+monthly_averages(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
         Month(1),
         (+),
         output_writer,
-        t_start,
         reference_date,
         short_names...;
         pre_output_hook! = average_pre_output_hook!,
     )
 
 """
-    monthly_average(short_names; output_writer, t_start, reference_date)
+    monthly_average(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that compute the monthly average for the given variable.
 """
 # An average is just a sum with a normalization before output
-monthly_average(short_names; output_writer, t_start, reference_date) =
-    monthly_averages(short_names; output_writer, t_start, reference_date)[1]
+monthly_average(FT, short_names; output_writer, reference_date) =
+    monthly_averages(FT, short_names; output_writer, reference_date)[1]
 
 """
-    tendaily_maxs(short_names...; output_writer, t_start, reference_date)
+    tendaily_maxs(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the max over ten days for the given variables.
 """
-tendaily_maxs(short_names...; output_writer, t_start, reference_date) =
+tendaily_maxs(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        10 * 24 * 60 * 60 * one(t_start),
+        10 * 24 * 60 * 60 * one(FT),
         max,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    tendaily_max(short_names; output_writer, t_start, reference_date)
+    tendaily_max(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the max over ten days for the given variable.
 """
-tendaily_max(short_names; output_writer, t_start, reference_date) =
-    tendaily_maxs(short_names; output_writer, t_start, reference_date)[1]
+tendaily_max(FT, short_names; output_writer, reference_date) =
+    tendaily_maxs(FT, short_names; output_writer, reference_date)[1]
 
 """
-    tendaily_mins(short_names...; output_writer, t_start, reference_date)
+    tendaily_mins(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the min over ten days for the given variables.
 """
-tendaily_mins(short_names...; output_writer, t_start, reference_date) =
+tendaily_mins(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        10 * 24 * 60 * 60 * one(t_start),
+        10 * 24 * 60 * 60 * one(FT),
         min,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    tendaily_min(short_names; output_writer, t_start, reference_date)
+    tendaily_min(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the min over ten days for the given variable.
 """
-tendaily_min(short_names; output_writer, t_start, reference_date) =
-    tendaily_mins(short_names; output_writer, t_start, reference_date)[1]
+tendaily_min(FT, short_names; output_writer, reference_date) =
+    tendaily_mins(FT, short_names; output_writer, reference_date)[1]
 
 """
-    tendaily_averages(short_names...; output_writer, t_start, reference_date)
+    tendaily_averages(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the average over ten days for the given variables.
 """
 # An average is just a sum with a normalization before output
-tendaily_averages(short_names...; output_writer, t_start, reference_date) =
+tendaily_averages(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        10 * 24 * 60 * 60 * one(t_start),
+        10 * 24 * 60 * 60 * one(FT),
         (+),
         output_writer,
-        t_start,
         reference_date,
         short_names...;
         pre_output_hook! = average_pre_output_hook!,
     )
 
 """
-    tendaily_average(short_names; output_writer, t_start, reference_date)
+    tendaily_average(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that compute the average over ten days for the given variable.
 """
 # An average is just a sum with a normalization before output
-tendaily_average(short_names; output_writer, t_start, reference_date) =
-    tendaily_averages(short_names; output_writer, t_start, reference_date)[1]
+tendaily_average(FT, short_names; output_writer, reference_date) =
+    tendaily_averages(FT, short_names; output_writer, reference_date)[1]
 
 """
-    daily_maxs(short_names...; output_writer, t_start, reference_date)
+    daily_maxs(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the daily max for the given variables.
 """
-daily_maxs(short_names...; output_writer, t_start, reference_date) =
+daily_maxs(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        24 * 60 * 60 * one(t_start),
+        24 * 60 * 60 * one(FT),
         max,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    daily_max(short_names; output_writer, t_start, reference_date)
+    daily_max(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the daily max for the given variable.
 """
-daily_max(short_names; output_writer, t_start, reference_date) =
-    daily_maxs(short_names; output_writer, t_start, reference_date)[1]
+daily_max(FT, short_names; output_writer, reference_date) =
+    daily_maxs(FT, short_names; output_writer, reference_date)[1]
 
 """
-    daily_mins(short_names...; output_writer, t_start, reference_date)
+    daily_mins(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the daily min for the given variables.
 """
-daily_mins(short_names...; output_writer, t_start, reference_date) =
+daily_mins(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        24 * 60 * 60 * one(t_start),
+        24 * 60 * 60 * one(FT),
         min,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    daily_min(short_names; output_writer, t_start, reference_date)
+    daily_min(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the daily min for the given variable.
 """
-daily_min(short_names; output_writer, t_start, reference_date) =
-    daily_mins(short_names; output_writer, t_start, reference_date)[1]
+daily_min(FT, short_names; output_writer, reference_date) =
+    daily_mins(FT, short_names; output_writer, reference_date)[1]
 
 """
-    daily_averages(short_names...; output_writer, t_start, reference_date)
+    daily_averages(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the daily average for the given variables.
 """
 # An average is just a sum with a normalization before output
-daily_averages(short_names...; output_writer, t_start, reference_date) =
+daily_averages(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        24 * 60 * 60 * one(t_start),
+        24 * 60 * 60 * one(FT),
         (+),
         output_writer,
-        t_start,
         reference_date,
         short_names...;
         pre_output_hook! = average_pre_output_hook!,
     )
 
 """
-    daily_average(short_names; output_writer, t_start, reference_date)
+    daily_average(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that compute the daily average for the given variable.
 """
 # An average is just a sum with a normalization before output
-daily_average(short_names; output_writer, t_start, reference_date) =
-    daily_averages(short_names; output_writer, t_start, reference_date)[1]
+daily_average(FT, short_names; output_writer, reference_date) =
+    daily_averages(FT, short_names; output_writer, reference_date)[1]
 
 """
-    hourly_maxs(short_names...; output_writer, t_start, reference_date)
+    hourly_maxs(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the hourly max for the given variables.
 """
-hourly_maxs(short_names...; output_writer, t_start, reference_date) =
+hourly_maxs(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        60 * 60 * one(t_start),
+        60 * 60 * one(FT),
         max,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    hourly_max(short_names; output_writer, t_start, reference_date)
+    hourly_max(FT, short_names; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the hourly max for the given variable.
 """
-hourly_max(short_names; output_writer, t_start, reference_date) =
-    hourly_maxs(short_names; output_writer, t_start, reference_date)[1]
+hourly_max(FT, short_names; output_writer, reference_date) =
+    hourly_maxs(FT, short_names; output_writer, reference_date)[1]
 
 """
-    hourly_mins(short_names...; output_writer, t_start, reference_date)
+    hourly_mins(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the hourly min for the given variables.
 """
-hourly_mins(short_names...; output_writer, t_start, reference_date) =
+hourly_mins(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        60 * 60 * one(t_start),
+        60 * 60 * one(FT),
         min,
         output_writer,
-        t_start,
         reference_date,
         short_names...,
     )
 
 """
-    hourly_mins(short_names...; output_writer, t_start, reference_date)
+    hourly_mins(FT, short_names...; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the hourly min for the given variable.
 """
-hourly_min(short_names; output_writer, t_start, reference_date) =
-    hourly_mins(short_names; output_writer, t_start, reference_date)[1]
+hourly_min(FT, short_names; output_writer, reference_date) =
+    hourly_mins(FT, short_names; output_writer, reference_date)[1]
 
 # An average is just a sum with a normalization before output
 """
-    hourly_averages(short_names...; output_writer, t_start, reference_date)
+    hourly_averages(FT, short_names...; output_writer,  reference_date)
 
 Return a list of `ScheduledDiagnostics` that compute the hourly average for the given variables.
 """
-hourly_averages(short_names...; output_writer, t_start, reference_date) =
+hourly_averages(FT, short_names...; output_writer, reference_date) =
     common_diagnostics(
-        60 * 60 * one(t_start),
+        60 * 60 * one(FT),
         (+),
         output_writer,
-        t_start,
         reference_date,
         short_names...;
         pre_output_hook! = average_pre_output_hook!,
     )
 
 """
-    hourly_average(short_names...; output_writer, t_start, reference_date)
+    hourly_average(FT, short_names...; output_writer,  reference_date)
 
 Return a `ScheduledDiagnostics` that computes the hourly average for the given variable.
 """
-hourly_average(short_names; output_writer, t_start, reference_date) =
-    hourly_averages(short_names; output_writer, t_start, reference_date)[1]
+hourly_average(FT, short_names; output_writer, reference_date) =
+    hourly_averages(FT, short_names; output_writer, reference_date)[1]
diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl
index 6ce71264e9..df784d4cc4 100644
--- a/src/shared_utilities/drivers.jl
+++ b/src/shared_utilities/drivers.jl
@@ -125,8 +125,8 @@ struct PrescribedAtmosphere{
     P::RA
     "CO2 concentration in atmosphere (mol/mol)"
     c_co2::CA
-    "Reference time - the datetime corresponding to t=0 for the simulation"
-    ref_time::DT
+    "Start date - the datetime corresponding to t=0 for the simulation"
+    start_date::DT
     "Reference height (m), relative to surface elevation"
     h::FT
     "Minimum wind speed (gustiness; m/s)"
@@ -140,14 +140,14 @@ struct PrescribedAtmosphere{
         u,
         q,
         P,
-        ref_time,
+        start_date,
         h::FT,
         earth_param_set;
         gustiness = FT(1),
         c_co2 = TimeVaryingInput((t) -> 4.2e-4),
     ) where {FT}
         thermo_params = LP.thermodynamic_parameters(earth_param_set)
-        args = (liquid_precip, snow_precip, T, u, q, P, c_co2, ref_time)
+        args = (liquid_precip, snow_precip, T, u, q, P, c_co2, start_date)
         return new{typeof(h), typeof.(args)..., typeof(thermo_params)}(
             args...,
             h,
@@ -412,12 +412,12 @@ struct PrescribedRadiativeFluxes{
     SW_d::SW
     "Downward longwave radiation function of time (W/m^2): positive indicates towards surface"
     LW_d::LW
-    "Reference time - the datetime corresponding to t=0 for the simulation"
-    ref_time::DT
+    "Start date - the datetime corresponding to t=0 for the simulation"
+    start_date::DT
     "Sun zenith angle, in radians"
     θs::T
-    function PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time; θs = nothing)
-        args = (SW_d, LW_d, ref_time, θs)
+    function PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date; θs = nothing)
+        args = (SW_d, LW_d, start_date, θs)
         return new{FT, typeof.(args)...}(args...)
     end
 end
@@ -894,7 +894,7 @@ function make_update_drivers(r::PrescribedRadiativeFluxes{FT}) where {FT}
         evaluate!(p.drivers.SW_d, r.SW_d, t)
         evaluate!(p.drivers.LW_d, r.LW_d, t)
         if !isnothing(r.θs)
-            p.drivers.θs .= FT.(r.θs(t, r.ref_time))
+            p.drivers.θs .= FT.(r.θs(t, r.start_date))
         else
             p.drivers.θs .= FT(0)
         end
diff --git a/src/standalone/Bucket/Bucket.jl b/src/standalone/Bucket/Bucket.jl
index fa37ae2174..ddb33e1330 100644
--- a/src/standalone/Bucket/Bucket.jl
+++ b/src/standalone/Bucket/Bucket.jl
@@ -156,8 +156,7 @@ end
 
 """
     PrescribedSurfaceAlbedo{FT}(
-        date_ref::Union{DateTime, DateTimeNoLeap},
-        t_start,
+        start_date::Union{DateTime, DateTimeNoLeap},
         Space::ClimaCore.Spaces.AbstractSpace;
         get_infile = ClimaLand.Artifacts.cesm2_albedo_dataset_path,
         varname = "sw_alb"
@@ -171,8 +170,7 @@ and download the data if it doesn't already exist on the machine.
 The input data file must have a time component.
 """
 function PrescribedSurfaceAlbedo{FT}(
-    date_ref::Union{DateTime, DateTimeNoLeap},
-    t_start,
+    start_date::Union{DateTime, DateTimeNoLeap},
     space::ClimaCore.Spaces.AbstractSpace;
     albedo_file_path = ClimaLand.Artifacts.cesm2_albedo_dataset_path(),
     varname = "sw_alb",
@@ -187,8 +185,7 @@ function PrescribedSurfaceAlbedo{FT}(
         albedo_file_path,
         varname,
         space;
-        reference_date = date_ref,
-        t_start,
+        reference_date = start_date,
         regridder_type,
     )
 
diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl
index 41dd15ea71..1b0db859f5 100644
--- a/src/standalone/Vegetation/Canopy.jl
+++ b/src/standalone/Vegetation/Canopy.jl
@@ -460,7 +460,7 @@ function ClimaLand.make_update_aux(
         compute_NIR!(inc_nir, RT, canopy.radiation, p, t)
         K = extinction_coeff.(G_Function, θs)
         DOY = Dates.dayofyear(
-            canopy.atmos.ref_time + Dates.Second(floor(Int64, t)),
+            canopy.atmos.start_date + Dates.Second(floor(Int64, t)),
         )
         @. frac_diff = diffuse_fraction(
             DOY,
diff --git a/test/standalone/Bucket/albedo_types.jl b/test/standalone/Bucket/albedo_types.jl
index 651558d6ba..a4d7841b52 100644
--- a/test/standalone/Bucket/albedo_types.jl
+++ b/test/standalone/Bucket/albedo_types.jl
@@ -131,22 +131,22 @@ end
     surface_coords = Fields.coordinate_field(space)
 
     infile_path = ClimaLand.Artifacts.cesm2_albedo_dataset_path()
-    date_ref_noleap = NCDataset(infile_path, "r") do ds
+    start_date_noleap = NCDataset(infile_path, "r") do ds
         ds["time"][1]
     end
     # Converting from NoLeap
-    date_ref = Dates.DateTime(
-        Dates.year(date_ref_noleap),
-        Dates.month(date_ref_noleap),
-        Dates.day(date_ref_noleap),
-        Dates.hour(date_ref_noleap),
-        Dates.minute(date_ref_noleap),
-        Dates.second(date_ref_noleap),
-        Dates.millisecond(date_ref_noleap),
+    start_date = Dates.DateTime(
+        Dates.year(start_date_noleap),
+        Dates.month(start_date_noleap),
+        Dates.day(start_date_noleap),
+        Dates.hour(start_date_noleap),
+        Dates.minute(start_date_noleap),
+        Dates.second(start_date_noleap),
+        Dates.millisecond(start_date_noleap),
     )
     t_start = Float64(0)
 
-    albedo = PrescribedSurfaceAlbedo{FT}(date_ref, t_start, space)
+    albedo = PrescribedSurfaceAlbedo{FT}(start_date, space)
 
     Y = (; bucket = (; W = Fields.zeros(space)))
     p = (; bucket = (; α_sfc = Fields.zeros(space)))
@@ -155,7 +155,7 @@ end
     varname = "sw_alb"
     file_dates = DataHandling.available_dates(albedo.albedo.data_handler)
 
-    new_date = date_ref + Second(t_start)
+    new_date = start_date + Second(t_start)
     t_curr = t_start
     for i in 1:5
         @assert new_date == file_dates[i]
@@ -205,14 +205,14 @@ end
             albedo = PrescribedBaregroundAlbedo{FT}(α_snow, surface_space)
 
             # Radiation
-            ref_time = DateTime(2005, 1, 15, 12)
+            start_date = DateTime(2005, 1, 15, 12)
             SW_d = (t) -> 0.0
             LW_d = (t) -> 5.67e-8 * 280.0^4.0
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             precip = (t) -> 0 # no precipitation
@@ -228,7 +228,7 @@ end
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
@@ -294,7 +294,7 @@ end
         Dates.millisecond(date),
     )
     file_dates = to_datetime.(file_dates_noleap)
-    date_ref = file_dates[1]
+    start_date_in_file = file_dates[1]
 
     bucket_domains = [
         Column(; zlim = FT.((-100.0, 0.0)), nelements = 10),
@@ -309,16 +309,17 @@ end
     for bucket_domain in bucket_domains
         space = bucket_domain.space.surface
         if bucket_domain isa SphericalShell
-            albedo_model = PrescribedSurfaceAlbedo{FT}(date_ref, t_start, space)
+            albedo_model =
+                PrescribedSurfaceAlbedo{FT}(start_date_in_file, space)
             # Radiation
-            ref_time = DateTime(2005, 1, 15, 12)
+            start_date = DateTime(2005, 1, 15, 12)
             SW_d = (t) -> 0
             LW_d = (t) -> 5.67e-8 * 280.0^4.0
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             precip = (t) -> 0 # no precipitation
@@ -327,7 +328,7 @@ end
             q_atmos = (t) -> 0.0 # no atmos water
             h_atmos = FT(1e-8)
             P_atmos = (t) -> 101325
-            ref_time = DateTime(2005, 1, 15, 12)
+            start_date = DateTime(2005, 1, 15, 12)
             bucket_atmos = PrescribedAtmosphere(
                 TimeVaryingInput(precip),
                 TimeVaryingInput(precip),
@@ -335,7 +336,7 @@ end
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
@@ -359,13 +360,13 @@ end
             set_initial_cache!(p, Y, FT(0.0))
             data_manual = DataHandling.regridded_snapshot(
                 albedo_model.albedo.data_handler,
-                date_ref,
+                start_date_in_file,
             )
 
             @test p.bucket.α_sfc == data_manual
 
             update_aux! = make_update_aux(model)
-            new_date = date_ref + Second(t_start)
+            new_date = start_date_in_file + Second(t_start)
             t_curr = t_start
             for i in 1:5
                 @assert new_date == file_dates[i]
@@ -387,8 +388,7 @@ end
             @test_throws "Using an albedo map requires a global run." PrescribedSurfaceAlbedo{
                 FT,
             }(
-                date_ref,
-                t_start,
+                start_date_in_file,
                 space,
             )
         end
diff --git a/test/standalone/Bucket/snow_bucket_tests.jl b/test/standalone/Bucket/snow_bucket_tests.jl
index f991a76dc8..e49a893c4c 100644
--- a/test/standalone/Bucket/snow_bucket_tests.jl
+++ b/test/standalone/Bucket/snow_bucket_tests.jl
@@ -59,14 +59,14 @@ for FT in (Float32, Float64)
 
     @testset "Small negative values FT = $FT " begin
         d = bucket_domains[1]
-        ref_time = DateTime(2005)
+        start_date = DateTime(2005)
         SW_d = (t) -> 20
         LW_d = (t) -> 20
         bucket_rad = PrescribedRadiativeFluxes(
             FT,
             TimeVaryingInput(SW_d),
             TimeVaryingInput(LW_d),
-            ref_time,
+            start_date,
         )
         # Atmos
         liquid_precip = (t) -> 0.0 # precipitation
@@ -84,7 +84,7 @@ for FT in (Float32, Float64)
             TimeVaryingInput(u_atmos),
             TimeVaryingInput(q_atmos),
             TimeVaryingInput(P_atmos),
-            ref_time,
+            start_date,
             h_atmos,
             earth_param_set,
         )
@@ -178,14 +178,14 @@ for FT in (Float32, Float64)
     for i in 1:3
         @testset "Conservation of water and energy I (snow present), FT = $FT" begin
             # Radiation
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             SW_d = (t) -> 20
             LW_d = (t) -> 20
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             liquid_precip = (t) -> -1e-8 # precipitation
@@ -203,7 +203,7 @@ for FT in (Float32, Float64)
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
@@ -292,14 +292,14 @@ for FT in (Float32, Float64)
     for i in 1:3
         @testset "Conservation of water and energy II (no snow to start), FT = $FT" begin
             # Radiation
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             SW_d = (t) -> 20
             LW_d = (t) -> 20
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             liquid_precip = (t) -> -1e-8 # precipitation
@@ -317,7 +317,7 @@ for FT in (Float32, Float64)
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
@@ -406,14 +406,14 @@ for FT in (Float32, Float64)
     @testset "Conservation of water and energy - nonuniform evaporation, FT = $FT" begin
         i = 3
         # Radiation
-        ref_time = DateTime(2005)
+        start_date = DateTime(2005)
         SW_d = (t) -> 20
         LW_d = (t) -> 20
         bucket_rad = PrescribedRadiativeFluxes(
             FT,
             TimeVaryingInput(SW_d),
             TimeVaryingInput(LW_d),
-            ref_time,
+            start_date,
         )
         # Atmos
         liquid_precip = (t) -> -1e-8 # precipitation
@@ -431,7 +431,7 @@ for FT in (Float32, Float64)
             TimeVaryingInput(u_atmos),
             TimeVaryingInput(q_atmos),
             TimeVaryingInput(P_atmos),
-            ref_time,
+            start_date,
             h_atmos,
             earth_param_set,
         )
diff --git a/test/standalone/Bucket/soil_bucket_tests.jl b/test/standalone/Bucket/soil_bucket_tests.jl
index 59c04a9866..99130dbbde 100644
--- a/test/standalone/Bucket/soil_bucket_tests.jl
+++ b/test/standalone/Bucket/soil_bucket_tests.jl
@@ -62,14 +62,14 @@ for FT in (Float32, Float64)
 
         @testset "Zero flux tendency, FT = $FT" begin
             # Radiation
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             SW_d = (t) -> 0
             LW_d = (t) -> 5.67e-8 * 280.0^4.0
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             precip = (t) -> 0 # no precipitation
@@ -85,7 +85,7 @@ for FT in (Float32, Float64)
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
@@ -162,14 +162,14 @@ for FT in (Float32, Float64)
 
         @testset "Energy + Moisture Conservation, FT = $FT" begin
             # Radiation
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             SW_d = (t) -> 10
             LW_d = (t) -> 300
             bucket_rad = PrescribedRadiativeFluxes(
                 FT,
                 TimeVaryingInput(SW_d),
                 TimeVaryingInput(LW_d),
-                ref_time,
+                start_date,
             )
             # Atmos
             precip = (t) -> -1e-6
@@ -186,7 +186,7 @@ for FT in (Float32, Float64)
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set,
             )
diff --git a/test/standalone/Snow/snow.jl b/test/standalone/Snow/snow.jl
index 8c486649d7..94dea91937 100644
--- a/test/standalone/Snow/snow.jl
+++ b/test/standalone/Snow/snow.jl
@@ -16,7 +16,7 @@ import ClimaLand.Parameters as LP
     FT = Float32
     earth_param_set = LP.LandParameters(FT)
 
-    ref_time = DateTime(2005)
+    start_date = DateTime(2005)
     param_set = LP.LandParameters(FT)
     Δt = FT(180.0)
     parameters = SnowParameters{FT}(Δt; earth_param_set = param_set)
@@ -24,7 +24,7 @@ import ClimaLand.Parameters as LP
     "Radiation"
     SW_d = TimeVaryingInput((t) -> eltype(t)(20.0))
     LW_d = TimeVaryingInput((t) -> eltype(t)(20.0))
-    rad = ClimaLand.PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time)
+    rad = ClimaLand.PrescribedRadiativeFluxes(FT, SW_d, LW_d, start_date)
     "Atmos"
     precip = TimeVaryingInput((t) -> eltype(t)(0)) # no precipitation
     T_atmos = TimeVaryingInput((t) -> eltype(t)(290.0))
@@ -39,7 +39,7 @@ import ClimaLand.Parameters as LP
         u_atmos,
         q_atmos,
         P_atmos,
-        ref_time,
+        start_date,
         h_atmos,
         earth_param_set,
     )
diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl
index 8f6bba7b78..47a5ae7d0b 100644
--- a/test/standalone/Soil/climate_drivers.jl
+++ b/test/standalone/Soil/climate_drivers.jl
@@ -48,14 +48,14 @@ for FT in (Float32, Float64)
         z_0m = FT(0.001)
         z_0b = z_0m
         # Radiation
-        ref_time = DateTime(2005)
+        start_date = DateTime(2005)
         SW_d = (t) -> 500
         LW_d = (t) -> 5.67e-8 * 280.0^4.0
         radiation = PrescribedRadiativeFluxes(
             FT,
             TimeVaryingInput(SW_d),
             TimeVaryingInput(LW_d),
-            ref_time,
+            start_date,
         )
         # Atmos
         precip = (t) -> 1e-8
@@ -72,7 +72,7 @@ for FT in (Float32, Float64)
             TimeVaryingInput(u_atmos),
             TimeVaryingInput(q_atmos),
             TimeVaryingInput(P_atmos),
-            ref_time,
+            start_date,
             h_atmos,
             earth_param_set,
         )
diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl
index 3a79dcf5d1..9f52df0fb2 100644
--- a/test/standalone/Vegetation/canopy_model.jl
+++ b/test/standalone/Vegetation/canopy_model.jl
@@ -50,17 +50,17 @@ import ClimaParams
 
             function zenith_angle(
                 t,
-                ref_time;
+                start_date;
                 latitude = lat,
                 longitude = long,
                 insol_params = earth_param_set.insol_params,
             )
-                current_datetime = ref_time + Dates.Second(round(t))
+                current_datetime = start_date + Dates.Second(round(t))
                 d, δ, η_UTC =
                     FT.(
                         Insolation.helper_instantaneous_zenith_angle(
                             current_datetime,
-                            ref_time,
+                            start_date,
                             insol_params,
                         )
                     )
@@ -95,7 +95,7 @@ import ClimaParams
             P_atmos = t -> 1e5 # Pa
             h_atmos = h_int # m
             c_atmos = (t) -> 4.11e-4 # mol/mol
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             atmos = PrescribedAtmosphere(
                 TimeVaryingInput(liquid_precip),
                 TimeVaryingInput(snow_precip),
@@ -103,7 +103,7 @@ import ClimaParams
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set;
                 c_co2 = TimeVaryingInput(c_atmos),
@@ -112,7 +112,7 @@ import ClimaParams
                 FT,
                 TimeVaryingInput(shortwave_radiation),
                 TimeVaryingInput(longwave_radiation),
-                ref_time;
+                start_date;
                 θs = zenith_angle,
             )
 
@@ -560,17 +560,17 @@ end
 
             function zenith_angle(
                 t,
-                ref_time;
+                start_date;
                 latitude = lat,
                 longitude = long,
                 insol_params = earth_param_set.insol_params,
             )
-                current_datetime = ref_time + Dates.Second(round(t))
+                current_datetime = start_date + Dates.Second(round(t))
                 d, δ, η_UTC =
                     FT.(
                         Insolation.helper_instantaneous_zenith_angle(
                             current_datetime,
-                            ref_time,
+                            start_date,
                             insol_params,
                         )
                     )
@@ -605,7 +605,7 @@ end
             P_atmos = t -> 1e5 # Pa
             h_atmos = h_int # m
             c_atmos = (t) -> 4.11e-4 # mol/mol
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             atmos = PrescribedAtmosphere(
                 TimeVaryingInput(liquid_precip),
                 TimeVaryingInput(snow_precip),
@@ -613,7 +613,7 @@ end
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set;
                 c_co2 = TimeVaryingInput(c_atmos),
@@ -622,7 +622,7 @@ end
                 FT,
                 TimeVaryingInput(shortwave_radiation),
                 TimeVaryingInput(longwave_radiation),
-                ref_time;
+                start_date;
                 θs = zenith_angle,
             )
 
@@ -836,17 +836,17 @@ end
 
             function zenith_angle(
                 t,
-                ref_time;
+                start_date;
                 latitude = lat,
                 longitude = long,
                 insol_params = earth_param_set.insol_params,
             )
-                current_datetime = ref_time + Dates.Second(round(t))
+                current_datetime = start_date + Dates.Second(round(t))
                 d, δ, η_UTC =
                     FT.(
                         Insolation.helper_instantaneous_zenith_angle(
                             current_datetime,
-                            ref_time,
+                            start_date,
                             insol_params,
                         )
                     )
@@ -881,7 +881,7 @@ end
             P_atmos = t -> 1e5 # Pa
             h_atmos = h_int # m
             c_atmos = (t) -> 4.11e-4 # mol/mol
-            ref_time = DateTime(2005)
+            start_date = DateTime(2005)
             atmos = PrescribedAtmosphere(
                 TimeVaryingInput(liquid_precip),
                 TimeVaryingInput(snow_precip),
@@ -889,7 +889,7 @@ end
                 TimeVaryingInput(u_atmos),
                 TimeVaryingInput(q_atmos),
                 TimeVaryingInput(P_atmos),
-                ref_time,
+                start_date,
                 h_atmos,
                 earth_param_set;
                 c_co2 = TimeVaryingInput(c_atmos),
@@ -898,7 +898,7 @@ end
                 FT,
                 TimeVaryingInput(shortwave_radiation),
                 TimeVaryingInput(longwave_radiation),
-                ref_time;
+                start_date;
                 θs = zenith_angle,
             )
 
diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl
index 2671e84234..6c3d380af5 100644
--- a/test/standalone/Vegetation/plant_hydraulics_test.jl
+++ b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -140,17 +140,17 @@ for FT in (Float32, Float64)
 
         function zenith_angle(
             t,
-            ref_time;
+            start_date;
             latitude = lat,
             longitude = long,
             insol_params = earth_param_set.insol_params,
         )
-            current_datetime = ref_time + Dates.Second(round(t))
+            current_datetime = start_date + Dates.Second(round(t))
             d, δ, η_UTC =
                 FT.(
                     Insolation.helper_instantaneous_zenith_angle(
                         current_datetime,
-                        ref_time,
+                        start_date,
                         insol_params,
                     )
                 )
@@ -185,7 +185,7 @@ for FT in (Float32, Float64)
         P_atmos = t -> 1e5 # Pa
         h_atmos = h_int # m
         c_atmos = (t) -> 4.11e-4 # mol/mol
-        ref_time = DateTime(2005)
+        start_date = DateTime(2005)
         atmos = PrescribedAtmosphere(
             TimeVaryingInput(liquid_precip),
             TimeVaryingInput(snow_precip),
@@ -193,7 +193,7 @@ for FT in (Float32, Float64)
             TimeVaryingInput(u_atmos),
             TimeVaryingInput(q_atmos),
             TimeVaryingInput(P_atmos),
-            ref_time,
+            start_date,
             h_atmos,
             earth_param_set;
             c_co2 = TimeVaryingInput(c_atmos),
@@ -202,7 +202,7 @@ for FT in (Float32, Float64)
             FT,
             TimeVaryingInput(shortwave_radiation),
             TimeVaryingInput(longwave_radiation),
-            ref_time;
+            start_date;
             θs = zenith_angle,
         )
         Δz = FT(1.0) # height of compartments
@@ -434,17 +434,17 @@ for FT in (Float32, Float64)
 
         function zenith_angle(
             t,
-            ref_time;
+            start_date;
             latitude = lat,
             longitude = long,
             insol_params = earth_param_set.insol_params,
         )
-            current_datetime = ref_time + Dates.Second(round(t))
+            current_datetime = start_date + Dates.Second(round(t))
             d, δ, η_UTC =
                 FT.(
                     Insolation.helper_instantaneous_zenith_angle(
                         current_datetime,
-                        ref_time,
+                        start_date,
                         insol_params,
                     )
                 )
@@ -479,7 +479,7 @@ for FT in (Float32, Float64)
         P_atmos = t -> 1e5 # Pa
         h_atmos = h_int # m
         c_atmos = (t) -> 4.11e-4 # mol/mol
-        ref_time = DateTime(2005)
+        start_date = DateTime(2005)
         atmos = PrescribedAtmosphere(
             TimeVaryingInput(liquid_precip),
             TimeVaryingInput(snow_precip),
@@ -487,7 +487,7 @@ for FT in (Float32, Float64)
             TimeVaryingInput(u_atmos),
             TimeVaryingInput(q_atmos),
             TimeVaryingInput(P_atmos),
-            ref_time,
+            start_date,
             h_atmos,
             earth_param_set;
             c_co2 = TimeVaryingInput(c_atmos),
@@ -496,7 +496,7 @@ for FT in (Float32, Float64)
             FT,
             TimeVaryingInput(shortwave_radiation),
             TimeVaryingInput(longwave_radiation),
-            ref_time;
+            start_date;
             θs = zenith_angle,
         )