diff --git a/examples/COPSE/COPSE_bergman2004_bergman2004.jl b/examples/COPSE/COPSE_bergman2004_bergman2004.jl
index ebf3c95..d3b632b 100644
--- a/examples/COPSE/COPSE_bergman2004_bergman2004.jl
+++ b/examples/COPSE/COPSE_bergman2004_bergman2004.jl
@@ -17,7 +17,7 @@ include("compare_output.jl")
 include("copse_bergman2004_bergman2004_expts.jl")
 
 
-# load archived model output 
+# load archived model output
 comparemodel = CompareOutput.copse_output_load("bergman2004","")
 # comparemodel = nothing
 
@@ -42,24 +42,26 @@ if use_TEMP_DAE
     println("integrate, DAE")
     # first run includes JIT time
     @time PALEOmodel.ODE.integrateDAE(
-        run, initial_state, modeldata, (-600e6, 0), 
+        run, initial_state, modeldata, (-600e6, 0),
         solvekwargs=(
+            dtmin=0.0,
             saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
         )
     )
-    
-    # sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?       
+
+    # sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
     # @time PALEOmodel.ODE.integrateDAEForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=IDA(linear_solver=:KLU))
 else
     println("integrate, ODE")
-    # first run includes JIT time    
+    # first run includes JIT time
     @time PALEOmodel.ODE.integrate(
-        run, initial_state, modeldata, (-600e6, 0), 
+        run, initial_state, modeldata, (-600e6, 0),
         solvekwargs=(
+            dtmin=0.0,
             saveat=1e6, # save output every 1e6 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/
         )
-    )  
-    # sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?       
+    )
+    # sparse jacobian with IDA(linear_solver=:KLU) fails at ~-350e6 yr ?
     # @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-650e6, 0), alg=CVODE_BDF(linear_solver=:KLU))
 end
 
diff --git a/examples/COPSE/COPSE_reloaded_reloaded.ipynb b/examples/COPSE/COPSE_reloaded_reloaded.ipynb
index 733868c..be04bbd 100644
--- a/examples/COPSE/COPSE_reloaded_reloaded.ipynb
+++ b/examples/COPSE/COPSE_reloaded_reloaded.ipynb
@@ -100,8 +100,8 @@
     "initial_deriv = similar(initial_state)\n",
     "\n",
     "PALEOmodel.SolverFunctions.ModelODE(modeldata)(\n",
-    "    initial_deriv, \n",
-    "    initial_state, \n",
+    "    initial_deriv,\n",
+    "    initial_state,\n",
     "    nothing,\n",
     "    0.0\n",
     ")\n",
@@ -139,12 +139,11 @@
     "\n",
     "# first run includes JIT time\n",
     "@time PALEOmodel.ODE.integrate(\n",
-    "    paleorun, initial_state, modeldata, (-1000e6, 0), \n",
+    "    paleorun, initial_state, modeldata, (-1000e6, 0),\n",
     "    solvekwargs=(\n",
     "        reltol=1e-4,\n",
     "    )\n",
-    "); \n",
-    "  "
+    ");\n"
    ]
   },
   {
@@ -303,7 +302,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# model-model plot \n",
+    "# model-model plot\n",
     "copse_reloaded_reloaded_plot([output_baseline_6C, output_baseline_3C]);\n"
    ]
   },
@@ -379,7 +378,7 @@
     "println(\"state and sms vars:\")\n",
     "println()\n",
     "\n",
-    "for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv)) \n",
+    "for (var, var_sms) in zip(PB.get_vars(modeldata.solver_view_all.stateexplicit), PB.get_vars(modeldata.solver_view_all.stateexplicit_deriv))\n",
     "    println(var.name, \"\\t\\t\", PB.get_var_type(var), \"\\t\\t\", PB.get_attribute(var, :units), \"\\t\\t\", PB.get_attribute(var, :description))\n",
     "    println(var_sms.name, \"\\t\\t\", PB.get_var_type(var_sms), \"\\t\\t\", PB.get_attribute(var_sms, :units), \"\\t\\t\", PB.get_attribute(var_sms, :description))\n",
     "end\n",
diff --git a/examples/COPSE/COPSE_reloaded_reloaded.jl b/examples/COPSE/COPSE_reloaded_reloaded.jl
index 41fdbe2..b019573 100644
--- a/examples/COPSE/COPSE_reloaded_reloaded.jl
+++ b/examples/COPSE/COPSE_reloaded_reloaded.jl
@@ -1,7 +1,7 @@
 using Logging
 import DataFrames
 
-using Plots 
+using Plots
 
 import PALEOboxes as PB
 import PALEOmodel
@@ -19,7 +19,7 @@ global_logger(ConsoleLogger(stderr,Logging.Info))
 
 @info "Start $(@__FILE__)"
 
-# load archived model output 
+# load archived model output
 # comparemodel=nothing
 
 # Baseline Phanerozoic configuration
@@ -43,13 +43,13 @@ tspan=(-1000e6, 0)
 # OOE oscillations with VCI and linear mocb
 # comparemodel = nothing
 # model = copse_reloaded_reloaded_expts(
-#     "reloaded", 
+#     "reloaded",
 #     [
 #         "VCI",
 #         "mocbProdLinear",
 #         ("set_par", "global", "force_LIPs", "co2releasefield", "CO2max"),
 #     ],
-# ) 
+# )
 # tspan=(-1000e6, 1e6)
 
 ######################################
@@ -69,18 +69,19 @@ println("integrate, no jacobian")
 run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())
 
 @time PALEOmodel.ODE.integrate(
-    run, initial_state, modeldata, tspan, 
+    run, initial_state, modeldata, tspan,
     solvekwargs=( # https://diffeq.sciml.ai/stable/basics/common_solver_opts/
+    dtmin=0.0,
         reltol=1e-4,
         # reltol=1e-5,
         # saveat=1e6,
         # dtmax=1e4,
     )
-) 
+)
 
 # println("integrate, jacobian from autodifferentiation")
 # @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-5,))
-                                   
+
 # PB.TestUtils.bench_model(run.model)
 
 ##############################
@@ -99,10 +100,10 @@ pager=PALEOmodel.PlotPager(9, (legend_background_color=nothing, ))
 copse_reloaded_reloaded_plot(run.output, pager=pager)
 
 #####################################
-# Check diff against comparison model 
+# Check diff against comparison model
 ########################################
 
-if !isnothing(comparemodel)   
+if !isnothing(comparemodel)
     diff = CompareOutput.compare_copse_output(run.output, comparemodel)
     firstpoint=50
     show(sort(DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean), :std), allrows=true)
@@ -110,7 +111,7 @@ if !isnothing(comparemodel)
     println()
 
     # Overlay plots for comparison
-    CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager) 
+    CompareOutput.copse_reloaded_comparecopse(run.output, comparemodel, include_Sr=true, pager=pager)
 end
 
 ####################################################################################################
diff --git a/examples/COPSE/runtests.jl b/examples/COPSE/runtests.jl
index 12e969d..0677011 100644
--- a/examples/COPSE/runtests.jl
+++ b/examples/COPSE/runtests.jl
@@ -13,7 +13,7 @@ include("compare_output.jl")
 
 @testset "COPSE examples" begin
 skipped_testsets = [
-    # "COPSE_reloaded_reloaded", 
+    # "COPSE_reloaded_reloaded",
     # "COPSE_reloaded_bergman2004",
     # "COPSE_bergman2004_bergman2004",
     # "COPSE_reloaded_reloaded.ipynb",
@@ -24,7 +24,7 @@ skipped_testsets = [
 
     include("copse_reloaded_reloaded_expts.jl")
 
-    # load archived model output 
+    # load archived model output
     comparemodel = CompareOutput.copse_output_load("reloaded", "reloaded_baseline")
 
     model = copse_reloaded_reloaded_expts("reloaded", ["baseline"])
@@ -33,24 +33,25 @@ skipped_testsets = [
     run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
 
     @time PALEOmodel.ODE.integrate(
-        run, initial_state, modeldata, (-1000e6, 0), 
+        run, initial_state, modeldata, (-1000e6, 0),
         solvekwargs=(
+            dtmin=0.0,
             reltol=1e-4,
         )
-    ) 
+    )
     # @time PALEOmodel.ODE.integrateForwardDiff(run, initial_state, modeldata, (-1000e6, 0), jac_ad_t_sparsity=0.0, solvekwargs=(reltol=1e-4,))
 
     # conservation checks
-    conschecks = [  
+    conschecks = [
         ("total_C",     :v,             1e-6 ),
         ("total_C",     :v_moldelta,    1e-6),
         ("total_S",     :v,             1e-6),
         ("total_S",     :v_moldelta,    1e-6),
         ("total_redox", nothing,        1e-6)
-    ]    
+    ]
     for (varname, propertyname, rtol) in conschecks
         startval, endval = PB.get_property(
-            PB.get_data(run.output, "global."*varname), 
+            PB.get_data(run.output, "global."*varname),
             propertyname=propertyname,
         )[[1, end]]
         println("check $varname $startval $endval $rtol")
@@ -62,7 +63,7 @@ skipped_testsets = [
     firstpoint=50
     diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
 
-    stdlimits = [   
+    stdlimits = [
         (:mccb_delta_diff,      1.5e-3),
         (:Sr_delta_diff,        4.0e-7),
         (:pCO2PAL_reldiff,      6.0e-4),
@@ -80,7 +81,7 @@ skipped_testsets = [
 
 end
 
-!("COPSE_reloaded_bergman2004" in skipped_testsets) && 
+!("COPSE_reloaded_bergman2004" in skipped_testsets) &&
 @testset "COPSE_reloaded_bergman2004" begin
     include("copse_reloaded_bergman2004_expts.jl")
 
@@ -93,21 +94,22 @@ end
     run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
 
     @time PALEOmodel.ODE.integrate(
-        run, initial_state, modeldata, (-1000e6, 0), 
+        run, initial_state, modeldata, (-1000e6, 0),
         solvekwargs=(
             reltol=1e-4,
+            dtmin=0.0,
             # saveat=1e6,
         ),
     )
 
     # conservation checks
-    conschecks = [  
+    conschecks = [
         ("total_C",     :v,             1e-6 ),
         ("total_C",     :v_moldelta,    1e-6),
         ("total_S",     :v,             1e-6),
         ("total_S",     :v_moldelta,    1e-6),
         ("total_redox", nothing,        1e-6)
-    ]    
+    ]
     for (varname, propertyname, rtol) in conschecks
         startval, endval = PB.get_property(
             PB.get_data(run.output, "global."*varname),
@@ -122,7 +124,7 @@ end
     firstpoint = 100
     diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
 
-    stdlimits = [   
+    stdlimits = [
         (:mccb_delta_diff,      4.0e-3),
         (:pCO2PAL_reldiff,      3.0e-4),
         (:pO2PAL_reldiff,       4.0e-4),
@@ -139,7 +141,7 @@ end
 
 end
 
-!("COPSE_bergman2004_bergman2004" in skipped_testsets) && 
+!("COPSE_bergman2004_bergman2004" in skipped_testsets) &&
 @testset "COPSE_bergman2004_bergman2004" begin
     include("copse_bergman2004_bergman2004_expts.jl")
 
@@ -150,22 +152,23 @@ end
     initial_state, modeldata = PALEOmodel.initialize!(model)
 
     run = PALEOmodel.Run(model=model, output=PALEOmodel.OutputWriters.OutputMemory())
-    
+
     @time PALEOmodel.ODE.integrate(
         run, initial_state, modeldata, (-600e6, 0),
         solvekwargs=(
+            dtmin=0.0,
             saveat=1e6,
         )
     )
 
     # conservation checks
-    conschecks = [  
+    conschecks = [
         ("total_C",     :v,             1e-6 ),
         ("total_C",     :v_moldelta,    1e-6),
         ("total_S",     :v,             1e-6),
         ("total_S",     :v_moldelta,    1e-6),
         ("total_redox", nothing,        1e-6)
-    ]    
+    ]
     for (varname, propertyname, rtol) in conschecks
         startval, endval = PB.get_property(
             PB.get_data(run.output, "global."*varname),
@@ -180,7 +183,7 @@ end
     firstpoint = 1
     diffsummary = DataFrames.describe(diff[firstpoint:end, :], :std, :min, :max, :mean)
 
-    stdlimits = [   
+    stdlimits = [
         (:mccb_delta_diff,      3.0e-2),
         (:pCO2PAL_reldiff,      9.0e-3),
         (:pO2PAL_reldiff,       6.0e-3),
@@ -197,7 +200,7 @@ end
 
 end
 
-!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) && 
+!("COPSE_reloaded_reloaded.ipynb" in skipped_testsets) &&
 @testset "COPSE_reloaded_reloaded.ipynb" begin
     # unicodeplots()
     gr() # headless environment will require ENV["GKSwstype"] = "100"