diff --git a/experiments/ClimaEarth/cli_options.jl b/experiments/ClimaEarth/cli_options.jl
index adc8fef516..8336030d4f 100644
--- a/experiments/ClimaEarth/cli_options.jl
+++ b/experiments/ClimaEarth/cli_options.jl
@@ -116,6 +116,10 @@ function argparse_settings()
         help = "Type of temperature anomaly for bucket model. [`amip`, `aquaplanet` (default)]"
         arg_type = String
         default = "aquaplanet"
+        "--use_land_diagnostics"
+        help = "Boolean flag indicating whether to compute and output land model diagnostics [`true` (default), `false`]"
+        arg_type = Bool
+        default = true
     end
     return s
 end
diff --git a/experiments/ClimaEarth/components/land/climaland_bucket.jl b/experiments/ClimaEarth/components/land/climaland_bucket.jl
index 5cddcdb024..4ee84197b7 100644
--- a/experiments/ClimaEarth/components/land/climaland_bucket.jl
+++ b/experiments/ClimaEarth/components/land/climaland_bucket.jl
@@ -63,6 +63,7 @@ function bucket_init(
     t_start::Float64,
     energy_check::Bool,
     surface_elevation,
+    use_land_diagnostics::Bool,
 ) where {FT}
     if config != "sphere"
         println(
@@ -132,7 +133,7 @@ function bucket_init(
     @. Y.bucket.T = T_sfc_0 + temp_anomaly(coords.subsurface)
     # `surface_elevation` is a ClimaCore.Fields.Field(`half` level)
     orog_adjusted_T = CC.Fields.field_values(Y.bucket.T) .- lapse_rate .* CC.Fields.field_values(surface_elevation)
-    # Adjust T based on surface elevation (p.bucket.T_sfc is then set using the 
+    # Adjust T based on surface elevation (p.bucket.T_sfc is then set using the
     # set_initial_cache! function)
     parent(Y.bucket.T) .= parent(orog_adjusted_T)
 
@@ -150,11 +151,15 @@ function bucket_init(
     prob = SciMLBase.ODEProblem(bucket_ode_function, Y, tspan, p)
 
     # Add diagnostics
-    netcdf_writer = CD.Writers.NetCDFWriter(domain.space.subsurface, artifacts_dir)
-    scheduled_diagnostics = CL.default_diagnostics(model, date_ref, output_writer = netcdf_writer)
+    if use_land_diagnostics
+        netcdf_writer = CD.Writers.NetCDFWriter(domain.space.subsurface, artifacts_dir)
+        scheduled_diagnostics = CL.default_diagnostics(model, date_ref, output_writer = netcdf_writer)
 
-    diagnostic_handler = CD.DiagnosticsHandler(scheduled_diagnostics, Y, p, t_start; dt = dt)
-    diag_cb = CD.DiagnosticsCallback(diagnostic_handler)
+        diagnostic_handler = CD.DiagnosticsHandler(scheduled_diagnostics, Y, p, t_start; dt = dt)
+        diag_cb = CD.DiagnosticsCallback(diagnostic_handler)
+    else
+        diag_cb = nothing
+    end
 
     integrator = SciMLBase.init(
         prob,
diff --git a/experiments/ClimaEarth/run_amip.jl b/experiments/ClimaEarth/run_amip.jl
index c7ebf9583a..eaf9f09424 100644
--- a/experiments/ClimaEarth/run_amip.jl
+++ b/experiments/ClimaEarth/run_amip.jl
@@ -146,6 +146,7 @@ restart_t = Int(config_dict["restart_t"])
 evolving_ocean = config_dict["evolving_ocean"]
 dt_rad = config_dict["dt_rad"]
 use_coupler_diagnostics = config_dict["use_coupler_diagnostics"]
+use_land_diagnostics = config_dict["use_land_diagnostics"]
 
 #=
 ## Setup Communication Context
@@ -281,6 +282,7 @@ if mode_name == "amip"
         t_start = t_start,
         energy_check = energy_check,
         surface_elevation,
+        use_land_diagnostics,
     )
 
     ## ocean stub
@@ -368,6 +370,7 @@ elseif mode_name in ("slabplanet", "slabplanet_aqua", "slabplanet_terra")
         t_start = t_start,
         energy_check = energy_check,
         surface_elevation,
+        use_land_diagnostics,
     )
 
     ## ocean model
@@ -417,6 +420,7 @@ elseif mode_name == "slabplanet_eisenman"
         t_start = t_start,
         energy_check = energy_check,
         surface_elevation,
+        use_land_diagnostics,
     )
 
     ## ocean stub (here set to zero area coverage)
diff --git a/experiments/ClimaEarth/run_cloudy_slabplanet.jl b/experiments/ClimaEarth/run_cloudy_slabplanet.jl
index bb0ddd583f..3c2e44efda 100644
--- a/experiments/ClimaEarth/run_cloudy_slabplanet.jl
+++ b/experiments/ClimaEarth/run_cloudy_slabplanet.jl
@@ -235,6 +235,7 @@ land_sim = bucket_init(
     t_start = tspan[1],
     energy_check = false,
     surface_elevation,
+    use_land_diagnostics = true,
 )
 
 ocean_sim = ocean_init(