https://sealevel.info/Comment-on-Stallinga2023/comment_on_stallinga2023_supp_files/calc_est_co2_removal_rates_v04.pl

#!/usr/bin/perl

# Simulation of CO2 removal from the atmosphere following achievement of zero emissions,
# based on linear relation between CO2 level and removal rate (fossil only), per:
# https://sealevel.info/Global_Carbon_Budget_2023v1.1_with_removal_rate.xlsx and
# https://sealevel.info/Global_Carbon_Budget_2023v1.1_with_removal_rate_plot2.png

# By Dave Burton, www.sealevel.info.  Uncopyrighted.


$equilibrium_level = 285;  # 11.105300 / 0.039005411
$sink_rate = 0.01832;  # 0.039005411/2.12940
$starting_level = $co2level = 421.08;
$starting_year = $year = 2023;

# input is the current CO2 level, in ppmv
# output is the amount of CO2 removed in one year by natural sinks, assuming ENSO neutral conditions
sub removal_rate {
  local($co2level) = shift;
  local($removalrate) = 0;
  local($co2elevation) = $co2level - $equilibrium_level;
  if ($co2level <= $equilibrium_level) {
    $removalrate = 0;
  } else {
    $removalrate = $co2elevation * $sink_rate;
  }
  return $removalrate;
}

# DEBUG PRINT
# print "Removal rates as a function of CO2 level:\n";
# for ($i=270; $i<=500; $i++) {
#   $removalrate = &removal_rate($i);
#   printf("$i %2.4f\n", $removalrate);
# }
# print "\n";

# SIMULATE DECLINE IN CO2 LEVEL IF EMISSIONS SUDDENLY WENT TO ZERO

$efoldinglife = 1 / $sink_rate;
$ln2 = 0.69314718056;  # ln(2)
$e = 2.71828182846;
$halflife = $efoldinglife * $ln2;
$halflevel = ($equilibrium_level + $starting_level) / 2;

printf("Simulated CO2 level, starting at %6.2f ppmv in $year, and zero emissions,\n",  $co2level);
print "based on linear relation between CO2 level and removal rate (fossil only), per:\n";
print "https://sealevel.info/Global_Carbon_Budget_2023v1.1_with_removal_rate.xlsx and\n";
print "https://sealevel.info/Global_Carbon_Budget_2023v1.1_with_removal_rate_plot2.png\n\n";

printf("Equilibrium CO2 level is %6.2f, half-way down to that is %5.2f ppmv.\n", $equilibrium_level, $halflevel);
printf("Half-life is %5.2f years, e-folding lifetime (adjustment time) is %5.2f years.\n", $halflife, $efoldinglife);

$msg1_done = $msg2_done = $msg3_done = 0;

print "YEAR  CO2\n";
while ($co2level > 330) {
   printf("$year %6.1f\n", $co2level);

   if (($co2level <= $halflevel) && !$msg3_done) {
     printf( "** %5.1f ppmv: half of the %5.1f ppmv CO2 above equilibrium is gone **\n",
             $halflevel, ($starting_level - $equilibrium_level) );
     $msg3_done = 1;
   }
   $year += 1;
   $removalrate = &removal_rate( $co2level );
   $co2level -= $removalrate;

   if (($year - $starting_year) >= $halflife) {
     if (!$msg1_done) {
       printf( "** date %6.1f: One half-life (%2.1f years) reached **\n",
               ($starting_year + $halflife), $halflife);
       $msg1_done = 1;
     }
   }
   if (($year - $starting_year) >= $efoldinglife) {
     if (!$msg2_done) {
       printf( "** date %6.1f: e-folding life (%2.1f years) and level (%6.2f ppmv) reached **\n",
               ($starting_year + $efoldinglife), $efoldinglife,
               ((($starting_level - $equilibrium_level) / $e) + $equilibrium_level) );
       $msg2_done = 1;
     }
   }
}

exit(0);