diff --git a/fair/forward.py b/fair/forward.py
index 360ffaeb..b67d66fa 100644
--- a/fair/forward.py
+++ b/fair/forward.py
@@ -299,9 +299,41 @@ def fair_scm(
         R_i = np.zeros(carbon_boxes_shape)
 
     if restart_in:
-        R_i[0]=restart_in[0]
-        T_j[0]=restart_in[1]
-        C_acc[0] = restart_in[2]
+        R_minus1 = restart_in[0]
+        T_j_minus1 = restart_in[1]
+        C_acc_minus1 = restart_in[2]
+        E_minus1 = restart_in[3]
+        C_minus1 = np.sum(R_minus1,axis=-1) + C_0[0]
+        # Calculate the parametrised iIRF and check if it is over the
+        # maximum allowed value
+        iirf[0] = rc * C_acc_minus1 + rt*np.sum(T_j_minus1) + r0
+        if iirf[0] >= iirf_max:
+            iirf[0] = iirf_max
+            
+        # Linearly interpolate a solution for alpha
+        time_scale_sf = (
+          root(iirf_interp_funct,0.16,args=(a,tau,iirf[0])))['x']
+
+        # Multiply default timescales by scale factor
+        tau_new = tau * time_scale_sf
+    
+        R_i[0,:] = R_minus1*np.exp(-1.0/tau_new) + a*emissions[0] / ppm_gtc
+        # Sum the boxes to get the total concentration anomaly
+        C[0,0] = np.sum(R_i[0,:],axis=-1) + C_0[0]
+        # Calculate the additional carbon uptake
+        C_acc[0] =  C_acc_minus1 + 0.5*(emissions[0] + E_minus1) - (
+          C[0,0] - C_minus1)*ppm_gtc
+          
+        if np.isscalar(other_rf):
+            F[0,0] = co2_log(C[0,0], C_pi[0], F2x) + other_rf
+        else:
+            F[0,0] = co2_log(C[0,0], C_pi[0], F2x) + other_rf[0]
+
+        F[0,0] = F[0,0] * scale[0]
+
+        T_j[0,:] = forc_to_temp(T_j_minus1, q[0,:], d, F[0,:])
+        T[0]=np.sum(T_j[0,:],axis=-1)
+
     else:
         # Initialise the carbon pools to be correct for first timestep in
         # numerical method
@@ -559,16 +591,14 @@ def fair_scm(
                 T_j[t,:] = forc_to_temp(T_j[t-1,:], q[t,:], d, F[t,:])
                 T[t]=np.sum(T_j[t,:],axis=-1)
 
-#    if emissions_driven:
-#        # add delta CO2 concentrations to initial value
-#        C[:,0] = C[:,0] + C_0[0]
-
     if not useMultigas:
         C = np.squeeze(C)
         F = np.squeeze(F)
-        
+        E_minus1 = emissions[-1]
+    else:
+        E_minus1 = np.sum(emissions[-1,1:3])
     if restart_out:
-        restart_out_val=(R_i[-1],T_j[-1],C_acc[-1])
+        restart_out_val=(R_i[-1],T_j[-1],C_acc[-1],E_minus1)
         return C, F, T, restart_out_val
     else:
         return C, F, T
diff --git a/tests/integration/integration_test.py b/tests/integration/integration_test.py
index 15965749..5e155784 100644
--- a/tests/integration/integration_test.py
+++ b/tests/integration/integration_test.py
@@ -213,47 +213,3 @@ def test_co2only_scale_error():
         emissions = rcp85.Emissions.co2,
         useMultigas=False,
         scale=1.0)
-
-
-def test_restart_co2_co2():
-    # intended behaviour: the last year of run 1 = the first year of run 2
-    C1, F1, T1, restart = fair.forward.fair_scm(
-        emissions   = rcp45.Emissions.co2,
-        useMultigas = False,
-        restart_out = True
-    )
-
-    C2, F2, T2 = fair.forward.fair_scm(
-        emissions   = np.zeros(10),
-        useMultigas = False,
-        restart_in  = restart
-    )
-
-    print C1[-1], C2
-    print F1[-1], F2
-    print T1[-1], T2
-    assert C1[-1] == C2[0]
-    assert F1[-1] == F2[0]
-    assert T1[-1] == T2[0]
-
-
-def test_restart_all_co2():
-    # intended behaviour: the last year of run 1 = the first year of run 2
-    C1, F1, T1, restart = fair.forward.fair_scm(
-        emissions   = rcp45.Emissions.emissions,
-        useMultigas = True,
-        restart_out = True
-    )
-
-    C2, F2, T2 = fair.forward.fair_scm(
-        emissions   = np.zeros(10),
-        useMultigas = False,
-        restart_in  = restart
-    )
-
-    print C1[-1,0], C2
-    print F1[-1,0], F2
-    print T1[-1], T2
-    assert C1[-1,0] == C2[0]
-    assert F1[-1,0] == F2[0]
-    assert T1[-1] == T2[0]
diff --git a/tests/unit/unit_test.py b/tests/unit/unit_test.py
index a93690d4..2fc89ec2 100644
--- a/tests/unit/unit_test.py
+++ b/tests/unit/unit_test.py
@@ -86,7 +86,7 @@ def test_steady():
         steady.emissions()
 
 
-def test_restart_continuous():
+def test_restart_co2_continuous():
     # Tests to check that a CO2-only run with a restart produces the same
     # results as a CO2-only run without a restart.
     
@@ -107,5 +107,6 @@ def test_restart_continuous():
         restart_in  = restart
         )
         
-    print C.shape, C1.shape, C2.shape
     assert np.all(C == np.concatenate((C1, C2)))
+    assert np.all(F == np.concatenate((F1, F2)))
+    assert np.all(T == np.concatenate((T1, T2)))