diff --git a/config/ci_configs/amip_component_dts.yml b/config/ci_configs/amip_component_dts.yml
index 1651e876fe..06c7f81849 100644
--- a/config/ci_configs/amip_component_dts.yml
+++ b/config/ci_configs/amip_component_dts.yml
@@ -8,6 +8,7 @@ dt_seaice: "37.5secs"
 dz_bottom: 30
 energy_check: false
 h_elem: 16
+land_model: "integrated"
 mode_name: "amip"
 moist: "equil"
 precip_model: "0M"
diff --git a/ext/ClimaCouplerMakieExt/debug_plots.jl b/ext/ClimaCouplerMakieExt/debug_plots.jl
index d1df6b74cb..0f32731ad3 100644
--- a/ext/ClimaCouplerMakieExt/debug_plots.jl
+++ b/ext/ClimaCouplerMakieExt/debug_plots.jl
@@ -122,17 +122,20 @@ Plot useful coupler fields (in `field_names`) and save plots to a directory.
 If `cs_fields_ref` is provided (e.g., using a copy of cs.fields from the initialization),
 plot the anomalies of the fields with respect to `cs_fields_ref`.
 """
-function Plotting.debug(cs_fields::CC.Fields.Field, dir, cs_fields_ref = nothing) # TODO should this dispatch on NamedTuple?
+function Plotting.debug(cs_fields::CC.Fields.Field, dir, cs_fields_ref = nothing)
     field_names = propertynames(cs_fields)
     fig = Makie.Figure(size = (1500, 800))
     min_square_len = ceil(Int, sqrt(length(field_names)))
+    has_nan = false
     for i in 1:min_square_len, j in 1:min_square_len
         field_index = (i - 1) * min_square_len + j
         if field_index <= length(field_names)
             field_name = field_names[field_index]
             field = getproperty(cs_fields, field_name)
 
-            title = string(field_name) * Plotting.print_extrema(field)
+            extrema_str, field_has_nan = Plotting.print_extrema(field)
+            has_nan = has_nan || field_has_nan
+            title = string(field_name) * extrema_str
             ax = Makie.Axis(fig[i, j * 2 - 1]; title)
             Plotting.debug_plot!(ax, fig, field, i, j)
         end
@@ -147,13 +150,20 @@ function Plotting.debug(cs_fields::CC.Fields.Field, dir, cs_fields_ref = nothing
                 field_name = field_names[field_index]
                 field = getproperty(cs_fields, field_name)
 
-                title = string(field_name) * Plotting.print_extrema(field)
+                extrema_str, field_has_nan = Plotting.print_extrema(field)
+                has_nan = has_nan || field_has_nan
+                title = string(field_name) * extrema_str
                 ax = Makie.Axis(fig[i, j * 2 - 1]; title)
                 Plotting.debug_plot!(ax, fig, field, i, j)
             end
         end
         Makie.save(joinpath(dir, "debug_coupler_anomalies.png"), fig)
     end
+
+    # Check for NaN errors after plots are saved
+    if has_nan
+        error("NaN values found in coupler fields extrema")
+    end
 end
 
 """
@@ -165,17 +175,25 @@ function Plotting.debug(sim::Interfacer.AbstractComponentSimulation, dir)
     field_names = Plotting.debug_plot_fields(sim)
     fig = Makie.Figure(size = (1500, 800))
     min_square_len = ceil(Int, sqrt(length(field_names)))
+    has_nan = false
     for i in 1:min_square_len, j in 1:min_square_len
         field_index = (i - 1) * min_square_len + j
         if field_index <= length(field_names)
             field_name = field_names[field_index]
             field = Interfacer.get_field(sim, Val(field_name))
-            title = string(field_name) * Plotting.print_extrema(field)
+            extrema_str, field_has_nan = Plotting.print_extrema(field)
+            has_nan = has_nan || field_has_nan
+            title = string(field_name) * extrema_str
             ax = Makie.Axis(fig[i, j * 2 - 1]; title)
             Plotting.debug_plot!(ax, fig, field, i, j)
         end
     end
     Makie.save(joinpath(dir, "debug_$(nameof(sim)).png"), fig)
+
+    # Check for NaN errors after plots are saved
+    if has_nan
+        error("NaN values found in field extrema of $(nameof(sim))")
+    end
 end
 
 """
@@ -226,15 +244,23 @@ Plotting.debug_plot!(ax, fig, field, i, j) = nothing # fallback method
 """
     Plotting.print_extrema(field::Union{CC.Fields.Field, Vector, StaticArrays.SVector, Number})
 
-Return the minimum and maximum values of a field as a string.
+Return a tuple `(string, has_nan::Bool)` where:
+- `string` is the minimum and maximum values of a field formatted as a string
+- `has_nan` is true if the extrema contain NaN values
 """
 function Plotting.print_extrema(
     field::Union{CC.Fields.Field, Vector, StaticArrays.SVector, Number},
 )
     ext_vals = extrema(field)
-    min = Printf.@sprintf("%.2E", ext_vals[1])
-    max = Printf.@sprintf("%.2E", ext_vals[2])
-    return " [$min, $max]"
+    min_val = ext_vals[1]
+    max_val = ext_vals[2]
+
+    # Check for NaN values
+    has_nan = isnan(min_val) || isnan(max_val)
+
+    min = Printf.@sprintf("%.2E", min_val)
+    max = Printf.@sprintf("%.2E", max_val)
+    return (" [$min, $max]", has_nan)
 end
 
 """
diff --git a/ext/ClimaCouplerOceananigansMakieExt.jl b/ext/ClimaCouplerOceananigansMakieExt.jl
index d2aaa3d4ff..5863d3ffe7 100644
--- a/ext/ClimaCouplerOceananigansMakieExt.jl
+++ b/ext/ClimaCouplerOceananigansMakieExt.jl
@@ -44,13 +44,20 @@ end
 """
     print_extrema(field)
 
-Return a string containing the minimum and maximum values of an Oceananigans
-field, formatted using scientific notation with 2 significant digits.
+Return a tuple `(string, has_nan::Bool)` where:
+- `string` is the minimum and maximum values of an Oceananigans field formatted as a string
+- `has_nan` is true if the extrema contain NaN values
 """
 function Plotting.print_extrema(field::OC.Field)
-    min = Printf.@sprintf("%.2E", minimum(field))
-    max = Printf.@sprintf("%.2E", maximum(field))
-    return " [$min, $max]"
+    min_val = minimum(field)
+    max_val = maximum(field)
+
+    # Check for NaN values
+    has_nan = isnan(min_val) || isnan(max_val)
+
+    min = Printf.@sprintf("%.2E", min_val)
+    max = Printf.@sprintf("%.2E", max_val)
+    return (" [$min, $max]", has_nan)
 end
 function Plotting.print_extrema(operation::OC.AbstractOperations.AbstractOperation)
     evaluated_field = OC.Field(operation)