diff --git a/docs/src/diagnostics/users_diagnostics.md b/docs/src/diagnostics/users_diagnostics.md
index 0c9728b579..198569052d 100644
--- a/docs/src/diagnostics/users_diagnostics.md
+++ b/docs/src/diagnostics/users_diagnostics.md
@@ -62,16 +62,34 @@ Your diagnostics have now been written in netcdf files in your output folder.
 
 When defining a custom diagnostic, follow these steps:
  - Define how to compute your diagnostic variable from your model state and cache.
- For example, let's say you want the bowen ratio (ratio between sensible heat and latent heat) in the Bucket model.
- ```
- function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel)
+For example, let's say you want the Bowen ratio (ratio between sensible heat and latent heat) in the Bucket model.
+```
+function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel)
     if isnothing(out)
         return copy(p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf)
     else
         out .= p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf
     end
 end
- ```
+```
+It is good practice to add error messages to inform other users that a diagnostic variable makes sense only with a
+specific `land_model`. This can be accomplished by prepending the `@with_error` macro at the function declaration,
+as in
+```
+import ClimaLand.Diagnostics: @witherror
+
+@with_error function compute_bowen_ratio!(out, Y, p, t, land_model::BucketModel)
+    if isnothing(out)
+        return copy(p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf)
+    else
+        out .= p.bucket.turbulent_fluxes.shf / p.bucket.turbulent_fluxes.lhf
+    end
+end
+```
+So, when someone tries outputting the Bowen ratio with a different model (e.g., `SnowModel`), `ClimaLand` will produce the following message:
+```
+Cannot compute albedo with model = SnowModel
+```
  - Add that diagnostic variable to your list of variables
  ```
  add_diagnostic_variable!(
diff --git a/src/diagnostics/bucket_compute_methods.jl b/src/diagnostics/bucket_compute_methods.jl
index 9ab65958f6..b63ff5f405 100644
--- a/src/diagnostics/bucket_compute_methods.jl
+++ b/src/diagnostics/bucket_compute_methods.jl
@@ -1,6 +1,6 @@
 # stored in p
 
-function compute_albedo!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_albedo!(out, Y, p, t, land_model::BucketModel)
     if isnothing(out)
         return copy(p.bucket.α_sfc)
     else
@@ -8,7 +8,13 @@ function compute_albedo!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_net_radiation!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_net_radiation!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.R_n)
     else
@@ -16,7 +22,13 @@ function compute_net_radiation!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_surface_temperature!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_surface_temperature!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.T_sfc)
     else
@@ -24,7 +36,7 @@ function compute_surface_temperature!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_surface_specific_humidity!(
+@with_error function compute_surface_specific_humidity!(
     out,
     Y,
     p,
@@ -38,7 +50,13 @@ function compute_surface_specific_humidity!(
     end
 end
 
-function compute_latent_heat_flux!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_latent_heat_flux!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.turbulent_fluxes.lhf)
     else
@@ -46,7 +64,13 @@ function compute_latent_heat_flux!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_aerodynamic_resistance!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_aerodynamic_resistance!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.turbulent_fluxes.r_ae)
     else
@@ -54,7 +78,13 @@ function compute_aerodynamic_resistance!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_sensible_heat_flux!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_sensible_heat_flux!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.turbulent_fluxes.shf)
     else
@@ -62,7 +92,7 @@ function compute_sensible_heat_flux!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_vapor_flux!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_vapor_flux!(out, Y, p, t, land_model::BucketModel)
     if isnothing(out)
         return copy(p.bucket.turbulent_fluxes.vapor_flux)
     else
@@ -70,7 +100,13 @@ function compute_vapor_flux!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_surface_air_density!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_surface_air_density!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(p.bucket.ρ_sfc)
     else
@@ -80,7 +116,13 @@ end
 
 # stored in Y
 
-function compute_soil_temperature!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_soil_temperature!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(Y.bucket.T)
     else
@@ -88,7 +130,7 @@ function compute_soil_temperature!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_subsurface_water_storage!(
+@with_error function compute_subsurface_water_storage!(
     out,
     Y,
     p,
@@ -102,7 +144,13 @@ function compute_subsurface_water_storage!(
     end
 end
 
-function compute_surface_water_content!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_surface_water_content!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(Y.bucket.Ws)
     else
@@ -110,7 +158,13 @@ function compute_surface_water_content!(out, Y, p, t, land_model::BucketModel)
     end
 end
 
-function compute_snow_water_equivalent!(out, Y, p, t, land_model::BucketModel)
+@with_error function compute_snow_water_equivalent!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::BucketModel,
+)
     if isnothing(out)
         return copy(Y.bucket.σS)
     else
diff --git a/src/diagnostics/define_diagnostics.jl b/src/diagnostics/define_diagnostics.jl
index b4f24380c2..e833784d48 100644
--- a/src/diagnostics/define_diagnostics.jl
+++ b/src/diagnostics/define_diagnostics.jl
@@ -1,28 +1,3 @@
-# General helper functions for undefined diagnostics for a particular model
-error_diagnostic_variable(variable, land_model::T) where {T} =
-    error("Cannot compute $variable with model = $T")
-
-# generate_error_functions is helper macro that generates the error message
-# when the user tries calling something that is incompatible with the model
-macro generate_error_functions(variable_names...)
-    functions = Expr[]
-    for variable in variable_names
-        function_name_sym = Symbol("compute_", variable, "!")
-        body = esc(quote
-            function $function_name_sym(_, _, _, _, land_model)
-                error_diagnostic_variable($variable, land_model)
-            end
-        end)
-        push!(functions, body)
-    end
-    return quote
-        $(functions...)
-    end
-end
-
-# TODO: Automatically generate this list from the names of the diagnostics
-@generate_error_functions "soil_net_radiation" "soil_latent_heat_flux" "soil_aerodynamic_resistance" "soil_sensible_heat_flux" "vapor_flux" "soil_temperature" "soil_water_liquid" "infiltration" "soilco2_diffusivity" "soilco2_source_microbe" "stomatal_conductance" "medlyn_term" "canopy_transpiration" "rainfall" "snowfall" "pressure" "wind_speed" "specific_humidity" "air_co2" "radiation_shortwave_down" "radiation_longwave_down" "photosynthesis_net_leaf" "photosynthesis_net_canopy" "respiration_leaf" "vcmax25" "photosynthetically_active_radiation" "photosynthetically_active_radiation_absorbed" "photosynthetically_active_radiation_reflected" "photosynthetically_active_radiation_transmitted" "near_infrared_radiation" "near_infrared_radiation_absorbed" "near_infrared_radiation_reflected" "near_infrared_radiation_transmitted" "radiation_shortwave_net" "radiation_longwave_net" "autotrophic_respiration" "soilco2" "heterotrophic_respiration" "soil_hydraulic_conductivity" "soil_water_potential" "soil_thermal_conductivity" "solar_zenith_angle" "moisture_stress_factor" "canopy_water_potential" "cross_section" "cross_section_roots" "area_index" "canopy_latent_heat_flux" "canopy_sensible_heat_flux" "canopy_aerodynamic_resistance" "canopy_temperature" "soil_ice"
-
 """
     define_diagnostics!(land_model)
 
diff --git a/src/diagnostics/diagnostic.jl b/src/diagnostics/diagnostic.jl
index 4894dd3f9c..c99383ff28 100644
--- a/src/diagnostics/diagnostic.jl
+++ b/src/diagnostics/diagnostic.jl
@@ -88,6 +88,44 @@ function get_diagnostic_variable(short_name)
     return ALL_DIAGNOSTICS[short_name]
 end
 
+# General helper functions for undefined diagnostics for a particular model
+error_diagnostic_variable(variable, land_model::T) where {T} =
+    error("Cannot compute $variable with model = $T")
+
+# with_error is a helper macro that generates the error message
+# when the user tries calling something that is incompatible with the model.
+# It should be called when defining compute functions
+macro with_error(compute_function_expr)
+    compute_function_expr.head == :function ||
+        error("Cannot parse this function, head is not a :function")
+
+    # Two firsts:
+    # 1st: extract the function signature
+    # 2nd: extract the name
+    function_name = first(first(compute_function_expr.args).args)
+    function_name isa Symbol || error("Cannot parse this function!")
+
+    # qualified_name ensures that this macro can be used outside of this module while
+    # still defining compute functions in this module
+    qualified_name = GlobalRef(Diagnostics, function_name)
+    method_already_exists = hasmethod(qualified_name, (Any, Any, Any, Any, Any))
+    if method_already_exists
+        return nothing
+    else
+        # Assuming the convention that functions are called "compute_variable!",
+        # otherwise the error might look a little less informative
+        variable_name =
+            replace(string(function_name), "compute_" => "", "!" => "")
+        return esc(
+            quote
+                function $qualified_name(_, _, _, _, land_model)
+                    error_diagnostic_variable($variable_name, land_model)
+                end
+            end,
+        )
+    end
+end
+
 # Do you want to define more diagnostics? Add them here
 include("bucket_compute_methods.jl")
 include("soilcanopy_compute_methods.jl")
diff --git a/src/diagnostics/soilcanopy_compute_methods.jl b/src/diagnostics/soilcanopy_compute_methods.jl
index eef836b786..a02cfe71bf 100644
--- a/src/diagnostics/soilcanopy_compute_methods.jl
+++ b/src/diagnostics/soilcanopy_compute_methods.jl
@@ -1,6 +1,12 @@
 # stored in p
 
-function compute_soil_net_radiation!(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_soil_net_radiation!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.soil.R_n)
     else
@@ -8,7 +14,7 @@ function compute_soil_net_radiation!(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soil_latent_heat_flux!(
+@with_error function compute_soil_latent_heat_flux!(
     out,
     Y,
     p,
@@ -22,7 +28,7 @@ function compute_soil_latent_heat_flux!(
     end
 end
 
-function compute_soil_aerodynamic_resistance!(
+@with_error function compute_soil_aerodynamic_resistance!(
     out,
     Y,
     p,
@@ -36,7 +42,7 @@ function compute_soil_aerodynamic_resistance!(
     end
 end
 
-function compute_soil_sensible_heat_flux!(
+@with_error function compute_soil_sensible_heat_flux!(
     out,
     Y,
     p,
@@ -50,7 +56,13 @@ function compute_soil_sensible_heat_flux!(
     end
 end
 
-function compute_vapor_flux!(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_vapor_flux!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.soil.turbulent_fluxes.vapor_flux)
     else
@@ -58,7 +70,7 @@ function compute_vapor_flux!(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soil_temperature!(
+@with_error function compute_soil_temperature!(
     out,
     Y,
     p,
@@ -72,7 +84,7 @@ function compute_soil_temperature!(
     end
 end
 
-function compute_soil_water_liquid!(
+@with_error function compute_soil_water_liquid!(
     out,
     Y,
     p,
@@ -86,7 +98,7 @@ function compute_soil_water_liquid!(
     end
 end
 
-function compute_infiltration(
+@with_error function compute_infiltration(
     out,
     Y,
     p,
@@ -100,7 +112,13 @@ function compute_infiltration(
     end
 end
 
-function compute_soilco2_diffusivity(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_soilco2_diffusivity(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.soilco2.D) # NOTE: we will need a method to compute surface co2 efflux
     else
@@ -108,7 +126,7 @@ function compute_soilco2_diffusivity(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soilco2_source_microbe(
+@with_error function compute_soilco2_source_microbe(
     out,
     Y,
     p,
@@ -122,7 +140,13 @@ function compute_soilco2_source_microbe(
     end
 end
 
-function compute_stomatal_conductance(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_stomatal_conductance(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.conductance.gs) # doublecheck: stomata, not canopy
     else
@@ -130,7 +154,13 @@ function compute_stomatal_conductance(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_medlyn_term(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_medlyn_term(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.conductance.medlyn_term)
     else
@@ -138,7 +168,13 @@ function compute_medlyn_term(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_canopy_transpiration(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_canopy_transpiration(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.transpiration) # doublecheck: canopy, not leaf
     else
@@ -146,7 +182,7 @@ function compute_canopy_transpiration(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P_liq) # I guess this is read and put it p. not computed. curious if we should handle this differently.
     else
@@ -154,7 +190,7 @@ function compute_rainfall(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P_snow) # following comment above, we could have a default getting only model output, and one also getting some inputs like drivers
     else
@@ -162,7 +198,7 @@ function compute_snowfall(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.P) # not sure if precip or pressure
     else
@@ -170,7 +206,13 @@ function compute_pressure(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_wind_speed(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_wind_speed(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.drivers.u)
     else
@@ -178,7 +220,13 @@ function compute_wind_speed(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_specific_humidity(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_specific_humidity(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.drivers.q) # check if this is correct. Also, check if Bucket has it and same name or not.
     else
@@ -186,7 +234,7 @@ function compute_specific_humidity(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.drivers.c_co2)
     else
@@ -194,7 +242,7 @@ function compute_air_co2(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_radiation_shortwave_down(
+@with_error function compute_radiation_shortwave_down(
     out,
     Y,
     p,
@@ -208,7 +256,7 @@ function compute_radiation_shortwave_down(
     end
 end
 
-function compute_radiation_longwave_down(
+@with_error function compute_radiation_longwave_down(
     out,
     Y,
     p,
@@ -222,7 +270,7 @@ function compute_radiation_longwave_down(
     end
 end
 
-function compute_photosynthesis_net_leaf(
+@with_error function compute_photosynthesis_net_leaf(
     out,
     Y,
     p,
@@ -236,7 +284,7 @@ function compute_photosynthesis_net_leaf(
     end
 end
 
-function compute_photosynthesis_net_canopy!(
+@with_error function compute_photosynthesis_net_canopy!(
     out,
     Y,
     p,
@@ -250,7 +298,13 @@ function compute_photosynthesis_net_canopy!(
     end
 end
 
-function compute_respiration_leaf(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_respiration_leaf(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.photosynthesis.Rd)
     else
@@ -258,7 +312,7 @@ function compute_respiration_leaf(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(p.canopy.photosynthesis.vcmax25)
     else
@@ -266,7 +320,7 @@ function compute_vcmax25(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_photosynthetically_active_radiation(
+@with_error function compute_photosynthetically_active_radiation(
     out,
     Y,
     p,
@@ -280,7 +334,7 @@ function compute_photosynthetically_active_radiation(
     end
 end
 
-function compute_photosynthetically_active_radiation_absorbed(
+@with_error function compute_photosynthetically_active_radiation_absorbed(
     out,
     Y,
     p,
@@ -294,7 +348,7 @@ function compute_photosynthetically_active_radiation_absorbed(
     end
 end
 
-function compute_photosynthetically_active_radiation_reflected(
+@with_error function compute_photosynthetically_active_radiation_reflected(
     out,
     Y,
     p,
@@ -308,7 +362,7 @@ function compute_photosynthetically_active_radiation_reflected(
     end
 end
 
-function compute_photosynthetically_active_radiation_transmitted(
+@with_error function compute_photosynthetically_active_radiation_transmitted(
     out,
     Y,
     p,
@@ -322,7 +376,7 @@ function compute_photosynthetically_active_radiation_transmitted(
     end
 end
 
-function compute_near_infrared_radiation(
+@with_error function compute_near_infrared_radiation(
     out,
     Y,
     p,
@@ -336,7 +390,7 @@ function compute_near_infrared_radiation(
     end
 end
 
-function compute_near_infrared_radiation_absorbed(
+@with_error function compute_near_infrared_radiation_absorbed(
     out,
     Y,
     p,
@@ -350,7 +404,7 @@ function compute_near_infrared_radiation_absorbed(
     end
 end
 
-function compute_near_infrared_radiation_reflected(
+@with_error function compute_near_infrared_radiation_reflected(
     out,
     Y,
     p,
@@ -364,7 +418,7 @@ function compute_near_infrared_radiation_reflected(
     end
 end
 
-function compute_near_infrared_radiation_transmitted(
+@with_error function compute_near_infrared_radiation_transmitted(
     out,
     Y,
     p,
@@ -378,7 +432,7 @@ function compute_near_infrared_radiation_transmitted(
     end
 end
 
-function compute_radiation_shortwave_net(
+@with_error function compute_radiation_shortwave_net(
     out,
     Y,
     p,
@@ -392,7 +446,7 @@ function compute_radiation_shortwave_net(
     end
 end
 
-function compute_radiation_longwave_net(
+@with_error function compute_radiation_longwave_net(
     out,
     Y,
     p,
@@ -406,7 +460,7 @@ function compute_radiation_longwave_net(
     end
 end
 
-function compute_autotrophic_respiration(
+@with_error function compute_autotrophic_respiration(
     out,
     Y,
     p,
@@ -420,7 +474,7 @@ function compute_autotrophic_respiration(
     end
 end
 
-function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
     if isnothing(out)
         return copy(Y.soilco2.C)
     else
@@ -428,7 +482,7 @@ function compute_soilco2(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_heterotrophic_respiration(
+@with_error function compute_heterotrophic_respiration(
     out,
     Y,
     p,
@@ -442,7 +496,7 @@ function compute_heterotrophic_respiration(
     end
 end
 
-function compute_soil_hydraulic_conductivity(
+@with_error function compute_soil_hydraulic_conductivity(
     out,
     Y,
     p,
@@ -456,7 +510,13 @@ function compute_soil_hydraulic_conductivity(
     end
 end
 
-function compute_soil_water_potential(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_soil_water_potential(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.soil.ψ)
     else
@@ -464,7 +524,7 @@ function compute_soil_water_potential(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soil_thermal_conductivity(
+@with_error function compute_soil_thermal_conductivity(
     out,
     Y,
     p,
@@ -478,7 +538,13 @@ function compute_soil_thermal_conductivity(
     end
 end
 
-function compute_solar_zenith_angle(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_solar_zenith_angle(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.drivers.θs)
     else
@@ -486,7 +552,7 @@ function compute_solar_zenith_angle(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_moisture_stress_factor(
+@with_error function compute_moisture_stress_factor(
     out,
     Y,
     p,
@@ -500,7 +566,7 @@ function compute_moisture_stress_factor(
     end
 end
 
-function compute_canopy_water_potential(
+@with_error function compute_canopy_water_potential(
     out,
     Y,
     p,
@@ -514,7 +580,13 @@ function compute_canopy_water_potential(
     end
 end
 
-function compute_cross_section(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_cross_section(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.hydraulics.fa)
     else
@@ -522,7 +594,13 @@ function compute_cross_section(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_cross_section_roots(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_cross_section_roots(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.hydraulics.fa_roots)
     else
@@ -530,7 +608,13 @@ function compute_cross_section_roots(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_area_index!(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_area_index!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(p.canopy.hydraulics.area_index.leaf)
     else
@@ -538,7 +622,7 @@ function compute_area_index!(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_canopy_latent_heat_flux(
+@with_error function compute_canopy_latent_heat_flux(
     out,
     Y,
     p,
@@ -552,7 +636,7 @@ function compute_canopy_latent_heat_flux(
     end
 end
 
-function compute_canopy_sensible_heat_flux(
+@with_error function compute_canopy_sensible_heat_flux(
     out,
     Y,
     p,
@@ -566,7 +650,7 @@ function compute_canopy_sensible_heat_flux(
     end
 end
 
-function compute_canopy_aerodynamic_resistance(
+@with_error function compute_canopy_aerodynamic_resistance(
     out,
     Y,
     p,
@@ -580,7 +664,13 @@ function compute_canopy_aerodynamic_resistance(
     end
 end
 
-function compute_canopy_temperature!(out, Y, p, t, land_model::SoilCanopyModel)
+@with_error function compute_canopy_temperature!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel,
+)
     if isnothing(out)
         return copy(Y.canopy.energy.T)
     else
@@ -588,7 +678,7 @@ function compute_canopy_temperature!(out, Y, p, t, land_model::SoilCanopyModel)
     end
 end
 
-function compute_soil_ice!(
+@with_error function compute_soil_ice!(
     out,
     Y,
     p,
diff --git a/test/diagnostics/diagnostics_tests.jl b/test/diagnostics/diagnostics_tests.jl
index 337db9320e..29ff4acdb6 100644
--- a/test/diagnostics/diagnostics_tests.jl
+++ b/test/diagnostics/diagnostics_tests.jl
@@ -3,6 +3,15 @@ using ClimaLand
 
 @test isdefined(ClimaLand.Diagnostics, :compute_albedo!)
 
+@test !hasmethod(
+    ClimaLand.Diagnostics.compute_albedo!,
+    (Any, Any, Any, Any, Any),
+)
+
+# Define some diagnostics for a DummyModel
+struct DummyModel end
+ClimaLand.Diagnostics.define_diagnostics!(DummyModel())
+
 # Just to trigger the error
 out = Y = p = t = land_model = nothing