diff --git a/sources/evolution.F90 b/sources/evolution.F90
index cc03629..ba99d4e 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -1601,7 +1601,7 @@ module evolution
     type(block_data), pointer :: pdata
 
     integer      :: status
-    real(kind=8) :: tm, dtm
+    real(kind=8) :: t, ds
 !
 !-------------------------------------------------------------------------------
 !
@@ -1609,76 +1609,94 @@ module evolution
     call start_timer(imu)
 #endif /* PROFILE */
 
+    status = 1
+
+    do while(status /= 0)
+
+      ds = 5.0d-01 * dt
+
 != 1st step: U(1) = U(n) + dt * F[U(n)]
 !
-    tm  = time + dt
-    dtm = dt
+      t = time + dt
 
-    pdata => list_data
-    do while (associated(pdata))
+      pdata => list_data
+      do while (associated(pdata))
 
-      call update_increment(pdata)
-      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
+        call update_increment(pdata)
+        call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
 
-      pdata => pdata%next
-    end do
+        pdata => pdata%next
+      end do
 
-    call boundary_fluxes()
+      call boundary_fluxes()
 
-    pdata => list_data
-    do while (associated(pdata))
+      pdata => list_data
+      do while (associated(pdata))
 
-      pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
+        pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
 
-      pdata%u => pdata%uu(:,:,:,:,2)
+        pdata%u => pdata%uu(:,:,:,:,2)
 
-      pdata => pdata%next
-    end do
+        pdata => pdata%next
+      end do
 
-    call update_variables(tm, dtm, status)
+      call update_variables(t, dt, status)
+
+      if (status /= 0) go to 100
 
 != 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)]
 !
-    tm  = time + dt
-    dtm = 0.5d+00 * dt
+      pdata => list_data
+      do while (associated(pdata))
 
-    pdata => list_data
-    do while (associated(pdata))
+        call update_increment(pdata)
+        call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
 
-      call update_increment(pdata)
-      call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
+        pdata => pdata%next
+      end do
+
+      call boundary_fluxes()
+
+      pdata => list_data
+      do while (associated(pdata))
+
+        pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1)                   &
+                                       + pdata%uu(:,:,:,:,2))                  &
+                                  + ds * pdata%du(:,:,:,:)
+
+        pdata => pdata%next
+      end do
+
+      if (ibp > 0) then
+        adecay(:) = exp(aglm(:) * cmax * dt)
+
+        pdata => list_data
+        do while (associated(pdata))
+
+          pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
+
+          pdata => pdata%next
+        end do
+      end if
+
+      call update_variables(t, ds, status)
+
+      100 continue
+
+      if (status /= 0) dt = 2.5d-01 * dt
 
-      pdata => pdata%next
     end do
 
-    call boundary_fluxes()
-
     pdata => list_data
     do while (associated(pdata))
 
-      pdata%uu(:,:,:,:,1) = 0.5d+00 * (pdata%uu(:,:,:,:,1)                     &
-                                     + pdata%uu(:,:,:,:,2)                     &
-                               +  dt * pdata%du(:,:,:,:))
+      pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
 
       pdata%u => pdata%uu(:,:,:,:,1)
 
       pdata => pdata%next
     end do
 
-    if (ibp > 0) then
-      adecay(:) = exp(aglm(:) * cmax * dt)
-
-      pdata => list_data
-      do while (associated(pdata))
-
-        pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
-
-        pdata => pdata%next
-      end do
-    end if
-
-    call update_variables(tm, dtm, status)
-
 #ifdef PROFILE
     call stop_timer(imu)
 #endif /* PROFILE */