diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 11516dc..f506265 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -2222,11 +2222,12 @@ module evolution
 !
   subroutine evolve_ssp324()
 
-    use blocks     , only : block_data, list_data
-    use boundaries , only : boundary_fluxes
-    use coordinates, only : nb, ne
-    use equations  , only : errors, ibp, cmax
-    use sources    , only : update_sources
+    use blocks         , only : block_data, list_data
+    use boundaries     , only : boundary_fluxes
+    use coordinates    , only : nb, ne
+    use equations      , only : errors, ibp, cmax
+    use iso_fortran_env, only : error_unit
+    use sources        , only : update_sources
 
     implicit none
 
@@ -2238,7 +2239,9 @@ module evolution
 
     real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00,                   &
                                twothird = 2.0d+00 / 3.0d+00
-!
+
+    character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()'
+
 !-------------------------------------------------------------------------------
 !
 #ifdef PROFILE
@@ -2301,6 +2304,7 @@ module evolution
         end do
 
         call update_variables(tm, dtm, status)
+        if (status /= 0) go to 100
 
       end do ! i = 1, 3
 
@@ -2318,6 +2322,7 @@ module evolution
       end do
 
       call update_variables(tm, dtm, status)
+      if (status /= 0) go to 100
 
 != 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate
 !
@@ -2344,6 +2349,7 @@ module evolution
       end do
 
       call update_variables(tm, dtm, status)
+      if (status /= 0) go to 100
 
 != 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation
 !
@@ -2385,16 +2391,21 @@ module evolution
       fc      = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi)
       dte     = dt * fc
 
-      if (fc > fac .or. nrej >= mrej) then
+      100 continue
+
+      if ((fc > fac .or. nrej >= mrej) .and. status == 0) then
         test = .false.
 
         errs(3) = errs(2)
         errs(2) = errs(1)
 
       else
-        dt  = dte
-
-        nrej        = nrej        + 1 ! rejection count in the current step
+        if (status == 0) then
+          dt   = dte
+          nrej = nrej + 1 ! rejection count in the current step
+        else
+          dt     = 2.5d-01 * dt
+        end if
         nrejections = nrejections + 1
 
 ! since the solution was rejected, we have to revert the primitive variables
@@ -2409,6 +2420,11 @@ module evolution
         end do
 
         call update_variables(tm, dtm, status)
+        if (status /= 0) then
+          write(error_unit,"('[', a, ']: ', a)") trim(loc),                    &
+                           "Possible bug: the revert to the previous state " //&
+                           "found it unphysical."
+        end if
 
       end if