diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 1739c0aed0465773e0e292d141a7056729bd21cf..d3739984ba9543f1c5154ee884e5384c152677eb 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -744,6 +744,7 @@ struct preferences_t {
     double artificial_orbital_angular_momentum_accretion_rate;
     double artificial_accretion_start_time;
     double artificial_accretion_end_time;
+    double artificial_mass_accretion_rate_by_stellar_type[NUMBER_OF_STELLAR_TYPES];
     double minimum_donor_menv_for_comenv;
     double magnetic_braking_factor;
     double magnetic_braking_gamma;
diff --git a/src/binary_star_functions/artificial_mass_accretion_rates.c b/src/binary_star_functions/artificial_mass_accretion_rates.c
index 64dee80dacb774e4c4cc972d8f78d2fecdf60603..61b87c40a96c7476c7c9f0ffb877a24dabb80ccd 100644
--- a/src/binary_star_functions/artificial_mass_accretion_rates.c
+++ b/src/binary_star_functions/artificial_mass_accretion_rates.c
@@ -12,8 +12,17 @@ void artificial_mass_accretion_rates(struct stardata_t * const stardata)
     {
         Foreach_star(star)
         {
+            /*
+             * Accrete at a fixed rate independent of stellar type
+             */
             star->derivative[DERIVATIVE_STELLAR_MASS_ARTIFICIAL] =
                 stardata->preferences->artificial_mass_accretion_rate[star->starnum];
+
+            /*
+             * Accrete at a fixed rate depending on stellar type
+             */
+            star->derivative[DERIVATIVE_STELLAR_MASS_ARTIFICIAL] =
+                stardata->preferences->artificial_mass_accretion_rate_by_stellar_type[star->stellar_type];
         }
     }
 }
diff --git a/src/evolution/evolution_difflog.c b/src/evolution/evolution_difflog.c
index c584dc6a233f9be6c56866f62f70b37d0ca2a035..801204cdee2d087d4baf17571b0e62b1dfbd7e10 100644
--- a/src/evolution/evolution_difflog.c
+++ b/src/evolution/evolution_difflog.c
@@ -775,7 +775,7 @@ static Boolean cf_diffstats(struct stardata_t * RESTRICT const stardata,
         }
 
         /* core collapse supernovae */
-        if(stardata->model.supernova)//now->SN_type[k] != SN_NONE)
+        if(stardata->model.supernova)
         {
             flags[SN_LABEL] = TRUE;
         }
diff --git a/src/logging/log_flags.def b/src/logging/log_flags.def
index e54491e8be3ad753879ffdb001287614620f0646..e3069df7e83cf91e11331514a1061ce52c847f4b 100644
--- a/src/logging/log_flags.def
+++ b/src/logging/log_flags.def
@@ -9,6 +9,7 @@
         X( PROBABILITY )                        \
         X( PERIASTRON_COLLISION )               \
         X( SN_AND_KICK )                        \
+        X( MCHAND )                             \
         X( OFF_MS )                             \
         X( STELLAR_TYPE_CHANGE )                \
         X( COMENV_POST )                        \
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 1a7aa9ae08e9e199da3519777482bf66e59404b2..e7932170ce710e34d579ea4d0c6b1315f6947072 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -1301,6 +1301,15 @@ BATCHMODE_ARGS
         Var(stardata->preferences->artificial_mass_accretion_rate),
         1.0
         CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+    ARG_SECTION_STARS,
+        "artificial_mass_accretion_rate_by_stellar_type%d",
+        "Constant mass accretion rate for stellar type <n>.",
+        ARG_FLOAT_SCANF ,
+        WTTS_USE_DEFAULT,
+        Var(stardata->preferences->artificial_mass_accretion_rate_by_stellar_type),
+        1.0
+        CMD_LINE_ARG_T_OPTIONAL_ARGS },
 {
     ARG_SECTION_STARS,
         "artificial_angular_momentum_accretion_rate%d",
diff --git a/src/stellar_structure/stellar_structure_AGB.c b/src/stellar_structure/stellar_structure_AGB.c
index 3867874e4a4bc8fe1d5b5dd4e873bd864e3fcf41..27f0508c012cd93bd767880359cad0d3b524337a 100644
--- a/src/stellar_structure/stellar_structure_AGB.c
+++ b/src/stellar_structure/stellar_structure_AGB.c
@@ -253,6 +253,15 @@ Stellar_type stellar_structure_AGB(struct star_t * const oldstar,
         (carbon_ignition != CARBON_IGNITION_NONE &&
          Boolean_(newstar->core_mass[CORE_CO] > mcmax));
 
+
+    if(supernova == TRUE &&
+       newstar->stellar_type == TPAGB)
+    {
+        Append_logstring(LOG_MCHAND,
+                         "0TPAGB core MCh reached");
+    }
+
+
     Dprint("supernova? %d\n",supernova);
 
     const Boolean envelope_lost =
diff --git a/src/stellar_structure/stellar_structure_WD.c b/src/stellar_structure/stellar_structure_WD.c
index 13c2b91b98626d189bcb4cf61713e8b31478e03a..2ae4f8331d14c65d7d36eb601a30e2befbf16b18 100644
--- a/src/stellar_structure/stellar_structure_WD.c
+++ b/src/stellar_structure/stellar_structure_WD.c
@@ -125,6 +125,9 @@ void stellar_structure_WD(struct star_t * const newstar,
                newstar->mass,
                MCh);
 
+        Append_logstring(LOG_MCHAND,
+                         "0WD MCh reached");
+
         struct star_t * news;
 
         if(newstar->stellar_type==ONeWD)
diff --git a/src/timestep/timestep_RLOF.c b/src/timestep/timestep_RLOF.c
index a433c60d27c43b1f075dea743fd2e06fce8a191a..6a58c1e5fd496f87a1eeecd34395f6dcc4e1c6a2 100644
--- a/src/timestep/timestep_RLOF.c
+++ b/src/timestep/timestep_RLOF.c
@@ -109,6 +109,7 @@ void timestep_RLOF(Timestep_prototype_args)
     if(Is_not_zero(stardata->model.RLOF_recommended_timestep))
     {
         Limit_timestep(*dt,
+                       stardata->preferences->timestep_multipliers[DT_LIMIT_RLOF]*
                        stardata->model.RLOF_recommended_timestep,
                        star,
                        DT_LIMIT_RLOF);
diff --git a/src/timestep/timestep_limits.c b/src/timestep/timestep_limits.c
index dee8310405e8525a9277f1f54c9357b8ed664fbf..368f40356d88d5eae9b5f02e5afe5fd27966934c 100644
--- a/src/timestep/timestep_limits.c
+++ b/src/timestep/timestep_limits.c
@@ -145,6 +145,45 @@ void timestep_limits(Timestep_prototype_args)
         }
     }
 
+
+    /*
+     * Mass gain: this usually requires that we are
+     * in a binary, but artificial accretion could
+     * also be set, so check it here.
+     */
+
+    {
+        const double accretion_rate = Mdot_gain(star);
+
+        if(Is_not_zero(accretion_rate))
+        {
+            if(Use_timestep(DT_LIMIT_MASS_GAIN))
+            {
+                const double dtlim = fsolver * 1e-6 *
+                    stardata->preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] *
+                    star->mass / accretion_rate;
+                Limit_timestep(*dt,
+                               dtlim,
+                               star,
+                               DT_LIMIT_MASS_GAIN);
+            }
+
+#ifdef SECOND_DERIVATIVES
+            if(Use_timestep(DT_LIMIT_MASS_GAIN2))
+            {
+                const double dtlim2 =
+                    stardata->preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN2] *
+                    SECOND_DERIVATIVE_F * Second_derivative(DT_LIMIT_MASS_GAIN) * 1e-6;
+                Limit_timestep(*dt,
+                               dtlim2,
+                               star,
+                               DT_LIMIT_MASS_GAIN2);
+            }
+#endif//SECOND_DERIVATIVES
+        }
+    }
+
+
     if(System_is_binary)
     {
         /*
@@ -250,34 +289,6 @@ void timestep_limits(Timestep_prototype_args)
 
         timestep_RLOF(Timestep_call_args);
 
-        {
-            const double accretion_rate = Mdot_gain(star);
-
-            if(Is_not_zero(accretion_rate))
-            {
-                if(Use_timestep(DT_LIMIT_MASS_GAIN))
-                {
-                    const double dtlim = fsolver * 0.1 * star->mass / accretion_rate;
-                    Limit_timestep(*dt,
-                                   dtlim,
-                                   star,
-                                   DT_LIMIT_MASS_GAIN);
-                }
-
-#ifdef SECOND_DERIVATIVES
-                if(Use_timestep(DT_LIMIT_MASS_GAIN2))
-                {
-                    const double dtlim2 =  SECOND_DERIVATIVE_F * Second_derivative(DT_LIMIT_MASS_GAIN) * 1e-6;
-                    Limit_timestep(*dt,
-                                   dtlim2,
-                                   star,
-                                   DT_LIMIT_MASS_GAIN2);
-                }
-#endif//SECOND_DERIVATIVES
-            }
-        }
-
-
         if(stardata->preferences->individual_novae == TRUE &&
            WHITE_DWARF(star->stellar_type) &&
            (star->nova != NOVA_STATE_NONE ||
diff --git a/src/timestep/timestep_set_default_multipliers.c b/src/timestep/timestep_set_default_multipliers.c
index a00fc77c1bb8c5f300c00f3c87ea173c5dad98bd..ad7b8f4fbd1167e7f7692e6c0ec8bbcc07fc5ff4 100644
--- a/src/timestep/timestep_set_default_multipliers.c
+++ b/src/timestep/timestep_set_default_multipliers.c
@@ -62,7 +62,7 @@ void timestep_set_default_multipliers(struct preferences_t * preferences)
     preferences->timestep_multipliers[DT_LIMIT_NOVAE] = 0.1; /* 39 */
     preferences->timestep_multipliers[DT_LIMIT_ARTIFICIAL_ACCRETION] = 1.0; /* 40 */
     preferences->timestep_multipliers[DT_LIMIT_SN] = 1.0; /* 41 */
-    preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 1.0; /* 42 */
+    preferences->timestep_multipliers[DT_LIMIT_MASS_GAIN] = 0.01; /* 42 */
     preferences->timestep_multipliers[DT_LIMIT_MASS_LOSS] = 1.0; /* 43 */
 #ifdef EXPONENTIAL_TIDES
     /* exponential tides are more stable, so can have a longer multiplier */
diff --git a/tbse b/tbse
index 219fdd075d548915c8dbc8a708f6ead06a1e8ef6..e1c57ba2399541bec0696a2a89aabc51aa5bc0a9 100755
--- a/tbse
+++ b/tbse
@@ -382,10 +382,17 @@ ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE1=0.0
 ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE2=0.0
 ARTIFICIAL_ORBITAL_ANGULAR_MOMENTUM_ACCRETION_RATE=0.0
 # Start and end times for the accretion, ignored if set to
-# ARTIFICIAL_ACCRETION_IGNORE = -1
+# ARTIFICIAL_ACCRETION_IGNORE == -1
 ARTIFICIAL_ACCRETION_START_TIME=ARTIFICIAL_ACCRETION_IGNORE
 ARTIFICIAL_ACCRETION_END_TIME=ARTIFICIAL_ACCRETION_IGNORE
 
+# You can also accrete by stellar type, e.g.
+# to accrete on COWDs (type 11) at a rate 1e-10 Msun/yr
+# do this. This still respects the start and end time
+# set above.
+#ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE="--artificial_mass_accretion_rate_by_stellar_type11 1e-10"
+ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE=""
+
 # Mass loss at the helium flash
 # If >=0, this is an amount in Msun
 # If <0, this is an algorithm 
@@ -2519,6 +2526,7 @@ $WARMUP \
 --individual_novae $INDIVIDUAL_NOVAE \
 --artificial_mass_accretion_rate1 $ARTIFICIAL_MASS_ACCRETION_RATE1 \
 --artificial_mass_accretion_rate2 $ARTIFICIAL_MASS_ACCRETION_RATE2 \
+$ARTIFICIAL_MASS_ACCRETION_BY_STELLAR_TYPE \
 --artificial_angular_momentum_accretion_rate1 $ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE1 \
 --artificial_angular_momentum_accretion_rate2 $ARTIFICIAL_ANGULAR_MOMENTUM_ACCRETION_RATE2 \
 --artificial_accretion_start_time $ARTIFICIAL_ACCRETION_START_TIME \
@@ -2980,7 +2988,7 @@ else
     ############################################################
     
     eval "echo $ARGS_PREFIX $ARGS $ARGS_POSTFIX | xargs $EXECUTABLE"
-
+    
     STATUS=$?
     if [ "$VERBOSE" -eq "1" ];
     then