diff --git a/config/default.cfg b/config/default.cfg
index 5b05a14565ca4fd7d3baf19295d200fe9030651c..03b465d7494b80d82effc814bde9123683d18005 100755
--- a/config/default.cfg
+++ b/config/default.cfg
@@ -367,6 +367,7 @@ cfg$gms$cm_APscen              <- "SSP2"  # def <- SSP2
 cfg$gms$cm_Full_Integration    <- "off"   # def = off
 
 cfg$gms$c_CO2priceDependent_AdjCosts    <- "on"   # def = on
+cfg$gms$c_scaleEmiHistorical   <- "off"    # def = off
 
 # MAGICC configuration (AJS)
 # Calibrate year 2000 temperature to HADCRUT4 data (which is very close to AR5).
diff --git a/core/datainput.gms b/core/datainput.gms
index 9cbd9afdd004272b025a40e665f31d556b71ac1f..375ef20ad167dcb5ba843323e7cbb2a5487b71f9 100644
--- a/core/datainput.gms
+++ b/core/datainput.gms
@@ -290,6 +290,24 @@ table f_dataetaglob(tall,all_te)                      "global eta data"
 $include "./core/input/generisdata_varying_eta.prn"
 ;
 
+* Read in mac historical emissions to calibrate MAC reference emissions
+parameter p_histEmiMac(tall,all_regi,all_enty)    "historical emissions per MAC"
+/
+$ondelim
+$include "./core/input/p_histEmiMac.cs4r"
+$offdelim
+/
+;
+* Read in historical emissions per sector to calibrate MAC reference emissions
+parameter p_histEmiSector(tall,all_regi,all_enty,emi_sectors,sector_types)    "historical emissions per sector"
+/
+$ondelim
+$include "./core/input/p_histEmiSector.cs4r"
+$offdelim
+/
+;
+
+
 ***---------------------------------------------------------------------------
 *** Import and set regional data
 ***---------------------------------------------------------------------------
diff --git a/core/declarations.gms b/core/declarations.gms
index 775fe5372940044d006bf7ba99ea1b789ec1ce5b..c4ebdb439204428f1bcc7bf5a7adf4c090add4a8 100644
--- a/core/declarations.gms
+++ b/core/declarations.gms
@@ -227,6 +227,9 @@ pm_cumDeprecFactor_old(ttot,all_regi,all_in)         "investment depreciation wi
 pm_cumDeprecFactor_new(ttot,all_regi,all_in)         "investment depreciation within a period, applied to the investment of t"
 
 p_Mport2005correct(all_regi,all_enty)                "correction factor to match fossil supply and internal region energy demand in the initial year"
+
+p_histEmiMac(tall,all_regi,all_enty)                 "historical emissions per MAC; from Eurostat and CEDS, to correct CH4 and N2O reporting"
+p_histEmiSector(tall,all_regi,all_enty,emi_sectors,sector_types) "historical emissions per sector; from Eurostat and CEDS, to correct CH4 and N2O reporting"
 ;
 
 ***----------------------------------------------------------------------------------------
diff --git a/core/postsolve.gms b/core/postsolve.gms
index bf7f542cd343c11974c8559c41d8e8bc7e55e440..d7d13b8280702abec9493fb52fbd4b071a59fe80 100644
--- a/core/postsolve.gms
+++ b/core/postsolve.gms
@@ -211,6 +211,10 @@ display s_actualbudgetco2;
 	 );
 );
 
+*** ---------------------------------------------------------------------------------------------------------------
+*** ENGAGE peakBudg formulation that works with several CO2 price path realizations of module 45 ---------------------
+*** it results in a peak budget with zero net CO2 emissions afterwards
+*** ---------------------------------------------------------------------------------------------------------------
 if(cm_iterative_target_adj eq 7,
 *JeS/CB* Update tax levels/ multigasbudget values to reach the peak CO2 budget, but make sure CO2 emissions afterward are close to zero on the global level
  
@@ -237,20 +241,20 @@ display p_actualbudgetco2;
 *** use multiplicative for budgets higher than 1600 Gt; for lower budgets, use multiplicative adjustment only for first 3 iterations, 
 			if(ord(iteration) lt 3 or c_budgetCO2 > 1600,
 			    !! change in CO2 price through adjustment: new price - old price; needed for adjustment option 2
-				p_taxCO2eq_iterationdiff(t,regi)$(t.val le cm_peakBudgYr) = pm_taxCO2eq(t,regi) * min(max((s_actualbudgetco2/c_budgetCO2)** (25/(2 * iteration.val + 23)),0.5+iteration.val/208),2 - iteration.val/102)  - pm_taxCO2eq(t,regi);
+				p_taxCO2eq_iterationdiff(t,regi) = pm_taxCO2eq(t,regi) * min(max((s_actualbudgetco2/c_budgetCO2)** (25/(2 * iteration.val + 23)),0.5+iteration.val/208),2 - iteration.val/102)  - pm_taxCO2eq(t,regi);
 				pm_taxCO2eq(t,regi)$(t.val le cm_peakBudgYr) = pm_taxCO2eq(t,regi) + p_taxCO2eq_iterationdiff(t,regi) ;
 				p_taxCO2eq_until2150(t,regi) = p_taxCO2eq_until2150(t,regi) + p_taxCO2eq_iterationdiff(t,regi) ;
 *** then switch to triangle-approximation based on last two iteration data points			
 			else
 			    !! change in CO2 price through adjustment: new price - old price; the two instances of "pm_taxCO2eq" cancel out -> only the difference term
 				!! until cm_peakBudgYr: expolinear price trajectory
-				p_taxCO2eq_iterationdiff_tmp(t,regi)$(t.val le cm_peakBudgYr) = 
+				p_taxCO2eq_iterationdiff_tmp(t,regi) = 
 				                      max(p_taxCO2eq_iterationdiff(t,regi) * min(max((c_budgetCO2 - s_actualbudgetco2)/(s_actualbudgetco2 - s_actualbudgetco2_last),-2),2),-pm_taxCO2eq(t,regi)/2);
 				pm_taxCO2eq(t,regi)$(t.val le cm_peakBudgYr) = pm_taxCO2eq(t,regi) + 
 				                      max(p_taxCO2eq_iterationdiff(t,regi) * min(max((c_budgetCO2 - s_actualbudgetco2)/(s_actualbudgetco2 - s_actualbudgetco2_last),-2),2),-pm_taxCO2eq(t,regi)/2);
-			    p_taxCO2eq_until2150(t,regi)$(t.val le cm_peakBudgYr) = p_taxCO2eq_until2150(t,regi) + 
+			    p_taxCO2eq_until2150(t,regi) = p_taxCO2eq_until2150(t,regi) + 
 				                      max(p_taxCO2eq_iterationdiff(t,regi) * min(max((c_budgetCO2 - s_actualbudgetco2)/(s_actualbudgetco2 - s_actualbudgetco2_last),-2),2),-p_taxCO2eq_until2150(t,regi)/2);
-				p_taxCO2eq_iterationdiff(t,regi)$(t.val le cm_peakBudgYr) = p_taxCO2eq_iterationdiff_tmp(t,regi);
+				p_taxCO2eq_iterationdiff(t,regi) = p_taxCO2eq_iterationdiff_tmp(t,regi);
 				!! after cm_peakBudgYr: adjustment so that emissions become zero: increase/decrease tax in each time step after cm_peakBudgYr by percentage of that year's total CO2 emissions of 2015 emissions
 			);
       o_taxCO2eq_iterDiff_Itr(iteration,regi) = p_taxCO2eq_iterationdiff("2030",regi);
@@ -297,8 +301,7 @@ display p_actualbudgetco2;
 		    display "shift peakBudgYr right";
             o_peakBudgYr_Itr(iteration+1) =  pm_ttot_val(ttot + 1);  !! ttot+1 is the new peakBudgYr
 			loop(t$(t.val ge pm_ttot_val(ttot + 1)),
-              pm_taxCO2eq(t,regi) = p_taxCO2eq_until2150(ttot+1,regi) 
-			                        + (t.val - pm_ttot_val(ttot + 1)) * cm_taxCO2inc_after_peakBudgYr * sm_DptCO2_2_TDpGtC;  !! increase by cm_taxCO2inc_after_peakBudgYr per year 
+              pm_taxCO2eq(t,regi) = p_taxCO2eq_until2150(t,regi);
             );
 		  );
         
diff --git a/core/preloop.gms b/core/preloop.gms
index 1aff0111c6e0e2077e4f2bf093ac84292cf3e466..f1306937adb9a445e856771e70d0063126ab35b4 100644
--- a/core/preloop.gms
+++ b/core/preloop.gms
@@ -106,5 +106,30 @@ pm_vintage_in(regi,"1",te) = pm_vintage_in(regi,"1",te) * max((pm_histfegrowth(r
 pm_vintage_in(regi,"6",te) = pm_vintage_in(regi,"6",te) * max(((pm_histfegrowth(regi,entyFe)- 0.005 + 1/fm_dataglob("lifetime",te))/(1/fm_dataglob("lifetime",te)) + 1) * 0.75,0.2);
 );
 
+$ifthen setGlobal c_scaleEmiHistorical
+*re-scale MAgPie reference emissions to be inline with eurostat data (MagPie overestimates non-CO2 GHG emissions by a factor of 50% more)
+display p_macBaseMagpie;
+loop(enty$(sameas(enty,"ch4rice") OR sameas(enty,"ch4animals") OR sameas(enty,"ch4anmlwst")),
+  p_macBaseMagpie(ttot,regi,enty)$(p_histEmiSector("2005",regi,"ch4","agriculture","process") AND (ttot.val ge 2005)) =
+   p_macBaseMagpie(ttot,regi,enty) *
+    ( (p_histEmiSector("2005",regi,"ch4","agriculture","process")+p_histEmiSector("2005",regi,"ch4","lulucf","process")) !!no rescaling needed - REMIND-internal unit is Mt CH4
+      /
+      (sum(enty2$(sameas(enty2,"ch4rice") OR sameas(enty2,"ch4animals") OR sameas(enty2,"ch4anmlwst")), p_macBaseMagpie("2005",regi,enty2)) + p_macBaseExo("2005",regi,"ch4agwaste"))
+    )
+  ;
+);
+loop(enty$(sameas(enty,"n2ofertin") OR sameas(enty,"n2ofertcr") OR sameas(enty,"n2oanwstc") OR sameas(enty,"n2oanwstm") OR sameas(enty,"n2oanwstp")),
+  p_macBaseMagpie(ttot,regi,enty)$(p_histEmiSector("2005",regi,"n2o","agriculture","process") AND (ttot.val ge 2005)) =
+    p_macBaseMagpie(ttot,regi,enty) *
+    ( p_histEmiSector("2005",regi,"n2o","agriculture","process")/( 44 / 28) !! rescaling to Mt N (internal unit for N2O emissions)
+* eurostat uses 298 to convert N2O to CO2eq
+      /
+      (sum(enty2$(sameas(enty,"n2ofertin") OR sameas(enty2,"n2ofertcr") OR sameas(enty2,"n2oanwstc") OR sameas(enty2,"n2oanwstm") OR sameas(enty2,"n2oanwstp")), p_macBaseMagpie("2005",regi,enty2)) + p_macBaseExo("2005",regi,"n2oagwaste"))
+    )
+  ;
+);
+display p_macBaseMagpie;
+$endif
+
 
 *** EOF ./core/preloop.gms
diff --git a/core/presolve.gms b/core/presolve.gms
index c88b227b4d1d74339d42a56699e583b3e523c2bc..65e53a0a46d43105c3d32500df44c8c705b90fa5 100644
--- a/core/presolve.gms
+++ b/core/presolve.gms
@@ -326,7 +326,25 @@ p_macLevFree(ttot,regi,enty)$( ttot.val gt 2005 )
       p_macUse2005(regi,enty)
     )$( (ttot.val gt 2040) )
 ;
+
+$ifthen setGlobal c_scaleEmiHistorical
+
+**p_macLevFree(ttot,regi,emiMacMagpie(enty))=0;
+* Set minimum abatment levels based on historical emissions
+p_macLevFree("2010",regi,enty)$(p_histEmiSector("2005",regi,"ch4","agriculture","process") AND (sameas(enty,"ch4rice") OR sameas(enty,"ch4animals") OR sameas(enty,"ch4anmlwst"))) = max( 0, 1 - (p_histEmiSector("2010",regi,"ch4","agriculture","process")+p_histEmiSector("2010",regi,"ch4","lulucf","process"))/(p_histEmiSector("2005",regi,"ch4","agriculture","process")+p_histEmiSector("2005",regi,"ch4","lulucf","process")));
+p_macLevFree(ttot,regi,enty)$((ttot.val ge 2015) AND p_histEmiSector("2005",regi,"ch4","agriculture","process") AND (sameas(enty,"ch4rice") OR sameas(enty,"ch4animals") OR sameas(enty,"ch4anmlwst"))) = max( 0, 1 - (p_histEmiSector("2015",regi,"ch4","agriculture","process")+p_histEmiSector("2015",regi,"ch4","lulucf","process"))/(p_histEmiSector("2005",regi,"ch4","agriculture","process")+p_histEmiSector("2005",regi,"ch4","lulucf","process")) );
+p_macLevFree("2010",regi,enty)$(p_histEmiSector("2005",regi,"ch4","agriculture","process") AND (sameas(enty,"n2ofertin") OR sameas(enty,"n2ofertcr") OR sameas(enty,"n2oanwstc") OR sameas(enty,"n2oanwstm") OR sameas(enty,"n2oanwstp"))) = max( 0, 1 - (p_histEmiSector("2010",regi,"n2o","agriculture","process")+p_histEmiSector("2010",regi,"n2o","lulucf","process"))/(p_histEmiSector("2005",regi,"n2o","agriculture","process")+p_histEmiSector("2005",regi,"n2o","lulucf","process")) );
+p_macLevFree(ttot,regi,emiMacMagpie(enty))$((ttot.val ge 2015) AND p_histEmiSector("2005",regi,"n2o","agriculture","process") AND (sameas(enty,"n2ofertin") OR sameas(enty,"n2ofertcr") OR sameas(enty,"n2oanwstc") OR sameas(enty,"n2oanwstm") OR sameas(enty,"n2oanwstp"))) = max( 0, 1 - (p_histEmiSector("2015",regi,"n2o","agriculture","process")+p_histEmiSector("2015",regi,"n2o","lulucf","process"))/(p_histEmiSector("2005",regi,"n2o","agriculture","process")+p_histEmiSector("2005",regi,"n2o","lulucf","process")) );
+
+p_macLevFree("2010",regi,enty)$((p_histEmiMac("2010",regi,enty)) AND (sameas(enty,"ch4wstl") OR sameas(enty,"ch4wsts"))) = max( 0, 1 - (p_histEmiMac("2010",regi,enty)) /vm_macBase.l("2010",regi,enty) );
+p_macLevFree(ttot,regi,enty)$((ttot.val ge 2015) AND (p_histEmiMac("2015",regi,enty)) AND (sameas(enty,"ch4wstl") OR sameas(enty,"ch4wsts"))) = max( 0, 1 - (p_histEmiMac("2015",regi,enty))/vm_macBase.l("2015",regi,enty) );
+
+$else
+
 p_macLevFree(ttot,regi,emiMacMagpie(enty))=0;
+
+$endif
+
 pm_macAbatLev(ttot,regi,enty) = 0.0;
 pm_macAbatLev("2005",regi,enty) = p_macUse2005(regi,enty);
 pm_macAbatLev("2010",regi,enty) = p_macLevFree("2010",regi,enty);
diff --git a/core/sets.gms b/core/sets.gms
index b2a7c7652e02f57a054f1aab8aed80a99b669d3e..db162edeb5b376f0402413b91c80de4b35c889de 100755
--- a/core/sets.gms
+++ b/core/sets.gms
@@ -1707,6 +1707,28 @@ sectorExogEmi(all_sectorEmi) "sectors with exogenous emissions"
     extraction
     indprocess
 /
+emi_sectors  "comprehensive sector set used for more detailed emissions accounting (REMIND-EU) and for CH4 tier 1 scaling - potentially to be integrated with similar set all_exogEmi"
+/
+        power   "public electricity and heat production"
+        refining "petroleum refining"
+        solids  "manufacture of solid fuels and other energy industries"
+        extraction "fugitive emissions from fuel extraction"
+        build   "Commercial sector, institutional sector and households"
+        indst   "industry (including industrial processes)"
+        trans   "transportation"
+        agriculture "agriculture (plus forestry and fishing energy use)"
+        waste   "waste management"
+        cdr     "Transport, capture and storage of CO2"
+        lulucf  "Land use,  land use change,  and forestry (LULUCF)"
+        bunkers "International bunkers (maritime and aviation)"
+        other   "other sectors and multilateral operations"
+        indirect
+/
+sector_types "differentiation of energy and process emissions in each sector"
+/
+        energy "fuel combustion part (and emissions) of the sector activity"
+        process "process sepecific part (and emissions) of the sector activity"
+/
 ccsCo2(all_enty)    "only cco2 (???)"
 /
         cco2
diff --git a/main.gms b/main.gms
index 996e103dff1e53d2632fa9c5ec123fc7e30eff76..5e3c8604c1103dc919fe8f1b20a9167624861e24 100644
--- a/main.gms
+++ b/main.gms
@@ -378,6 +378,7 @@ cm_carbonprice_temperatureLimit       = 1.8;   !! def = 1.8
 cm_DiscRateScen = 0;!! def = 0
 cm_noReboundEffect = 0;
 $setGlobal cm_EsubGrowth  low  !! def = low
+$setGlobal c_scaleEmiHistorical  off  !! def = off
 
 
 $setGlobal c_regi_nucscen  all !! def = all