diff --git a/modules/15_climate/box/equations.gms b/modules/15_climate/box/equations.gms
index 80cb4efe7dd5b6844bae645d2d4cffd19c2b5112..697035b64d75572b5a2460d1f91e5fa01fa33dd2 100644
--- a/modules/15_climate/box/equations.gms
+++ b/modules/15_climate/box/equations.gms
@@ -8,129 +8,192 @@
 ta2ttot10(ta10, ttot)$((ttot.val lt ta10.val) AND (pm_ttot_val(ttot+1) gt ta10.val)AND (ttot.val ge 2005)) = Yes;
 display ttot2ta10, ta2ttot10;
 
-
-*** taken from box-model.inc ---------------------------
-
-***     carbon cycle
-
-* AL changed towards an IRF Model with 2 or 3 time scales
+***---------------------------------------------------------------------------
+*'  Carbon Cycle. CO2 concentration is calculated using and Impulse-Response-Function Model with 3 time scales
+***---------------------------------------------------------------------------
 q15_cc(ta10)$(ord(ta10) ge 1)..
-                   v15_conc(ta10,'CO2') =e= s15_c0 + ((s15_c2000-s15_c0)/(s15_ca0*s15_cq0+s15_ca1*s15_cq1+s15_ca2*s15_cq2+s15_ca3*s15_cq3))*
-                                                                 (s15_ca0*s15_cq0 + s15_ca1*s15_cq1*exp(-(ord(ta10)-1)/s15_ctau1) + s15_ca2*s15_cq2*exp(-(ord(ta10)-1)/s15_ctau2) + s15_ca3*s15_cq3*exp(-(ord(ta10)-1)/s15_ctau3))+
-                                                                                                 s15_cconvi*sum((tx)$(ord(tx)lt ord(ta10)),v15_emi(tx,'CO2')*p15_epsilon(ta10-ord(tx)));
-
-***     N20 concentration (from ACC2)
-
-q15_concN2OQ(ta10) $(Ord(ta10) GT 1)     .. (v15_conc(ta10,'N2O')-v15_conc(ta10-1,'N2O'))/s15_DELTAT =E= 0.5*
-                                 ((1/s15_CNVN2O*(v15_emi(ta10,'N2O')+ s15_NATN2O)-(1/s15_TAUN2O*(v15_conc(ta10,'N2O')/s15_CONN2O2000R)**(-s15_SENTAUN2O)*v15_conc(ta10,'N2O')))+
-                                 (1/s15_CNVN2O*(v15_emi(ta10-1,'N2O')+ s15_NATN2O)-(1/s15_TAUN2O*(v15_conc(ta10-1,'N2O')/s15_CONN2O2000R)**(-s15_SENTAUN2O)*v15_conc(ta10-1,'N2O'))));
-*** end of stuff taken from box-model.inc ---------------------------
-
-*** taken from box-model.inc ---------------------------
-q15_concCH4Q(ta10) $(Ord(ta10) GT 1)     .. (v15_conc(ta10,'CH4')-v15_conc(ta10-1,'CH4'))/s15_DELTAT =E= 0.5*
-                                 ((1/s15_CNVCH4*(v15_emi(ta10,'CH4')+ s15_NATCH4)-
-                                 (p15_conroh(ta10)/s15_TAUCH4OH+1/s15_TAUCH4SS)*v15_conc(ta10,'CH4'))+
-                                 (1/s15_CNVCH4*(v15_emi(ta10-1,'CH4')+ s15_NATCH4)-
-                                 (p15_conroh(ta10-1)/s15_TAUCH4OH+1/s15_TAUCH4SS)*v15_conc(ta10-1,'CH4')));
-
-***    CH4 radiative forcing (from ACC2)
-
-q15_forcCH4Q(ta10)                    .. v15_forcComp(ta10,'CH4') =E= s15_RHOCH4*(SQRT(v15_conc(ta10,'CH4'))-SQRT(s15_CONCH4PRE))
-                                 -s15_OVERLFAC1*LOG(1+s15_OVERLFAC2*(v15_conc(ta10,'CH4')*s15_CONN2OPRE)**s15_OVERLEXP1+
-                                 s15_OVERLFAC3*v15_conc(ta10,'CH4')*(v15_conc(ta10,'CH4')*s15_CONN2OPRE)**s15_OVERLEXP2)
-                                 +s15_OVERLFAC1*LOG(1+s15_OVERLFAC2*(s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP1+
-                                 s15_OVERLFAC3*s15_CONCH4PRE*(s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP2);
-***    N2O radiative forcing (from ACC2)
-q15_forcN2OQ(ta10)                    .. v15_forcComp(ta10,'N2O') =E= s15_RHON2O*(SQRT(v15_conc(ta10,'N2O'))-SQRT(s15_CONN2OPRE))
-                                 -s15_OVERLFAC1*LOG(1+s15_OVERLFAC2*(s15_CONCH4PRE*v15_conc(ta10,'N2O'))**s15_OVERLEXP1+
-                                 s15_OVERLFAC3*s15_CONCH4PRE*(s15_CONCH4PRE*v15_conc(ta10,'N2O'))**s15_OVERLEXP2)
-                                 +s15_OVERLFAC1*LOG(1+s15_OVERLFAC2*(s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP1+
-                                 s15_OVERLFAC3*s15_CONCH4PRE*(s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP2);
-
-*AL* changed towards the ACC2 approach to use the foghg
-***     so2 radiative forcing module
-
-    q15_forcso2(ta10)..
-                    v15_forcComp(ta10,'SO2') =e= s15_dso1990 * v15_emi(ta10,'SO2') / s15_so1990 +
-                             s15_iso1990 * log(1 + v15_emi(ta10,'SO2')/s15_enatso2) /
-                                       log(1 + s15_so1990/s15_enatso2) ;
+    v15_conc(ta10,'CO2')
+	=e=
+	s15_c0 + ((s15_c2000-s15_c0)/(s15_ca0*s15_cq0+s15_ca1*s15_cq1+s15_ca2*s15_cq2+s15_ca3*s15_cq3))
+	         * (s15_ca0*s15_cq0 + s15_ca1*s15_cq1*exp(-(ord(ta10)-1)/s15_ctau1) + s15_ca2*s15_cq2*exp(-(ord(ta10)-1)/s15_ctau2) + s15_ca3*s15_cq3*exp(-(ord(ta10)-1)/s15_ctau3))
+		   + s15_cconvi * sum((tx)$(ord(tx)lt ord(ta10)),v15_emi(tx,'CO2')*p15_epsilon(ta10-ord(tx)));
+
+***---------------------------------------------------------------------------
+*'  CH4 concentration
+***---------------------------------------------------------------------------
+q15_concCH4Q(ta10)$(Ord(ta10) GT 1)..
+    (v15_conc(ta10,'CH4')-v15_conc(ta10-1,'CH4'))/s15_DELTAT
+	=e=
+	0.5 * (
+	      (1/s15_CNVCH4 * (v15_emi(ta10,'CH4')+ s15_NATCH4) - (p15_conroh(ta10) / s15_TAUCH4OH + 1 / s15_TAUCH4SS) * v15_conc(ta10,'CH4'))
+		  + (1/s15_CNVCH4*(v15_emi(ta10-1,'CH4') + s15_NATCH4) - (p15_conroh(ta10-1) / s15_TAUCH4OH + 1 / s15_TAUCH4SS) * v15_conc(ta10-1,'CH4'))
+		  );
+
+***---------------------------------------------------------------------------
+*'  N20 concentration (from ACC2)
+***---------------------------------------------------------------------------
+q15_concN2OQ(ta10)$(Ord(ta10) GT 1)..
+    (v15_conc(ta10,'N2O')-v15_conc(ta10-1,'N2O'))/s15_DELTAT
+	=e=
+	0.5 * (
+	(1/s15_CNVN2O * (v15_emi(ta10,'N2O')+ s15_NATN2O) - (1/s15_TAUN2O*(v15_conc(ta10,'N2O')/s15_CONN2O2000R)**(-s15_SENTAUN2O)*v15_conc(ta10,'N2O')))
+	+ (1/s15_CNVN2O * (v15_emi(ta10-1,'N2O') + s15_NATN2O) - (1/s15_TAUN2O*(v15_conc(ta10-1,'N2O')/s15_CONN2O2000R)**(-s15_SENTAUN2O)*v15_conc(ta10-1,'N2O')))
+	);
+
+***---------------------------------------------------------------------------
+*'    CO2 radiative forcing
+***---------------------------------------------------------------------------
+q15_forcco2(ta10)..
+    v15_forcComp(ta10,'CO2')
+	=e=
+	s15_fcodb * log(v15_conc(ta10,'CO2')/s15_c0) / log(2) ;
+
+***---------------------------------------------------------------------------
+*'    CH4 radiative forcing (from ACC2)
+***---------------------------------------------------------------------------
+q15_forcCH4Q(ta10)..
+    v15_forcComp(ta10,'CH4')
+	=e=
+	s15_RHOCH4 * (SQRT(v15_conc(ta10,'CH4')) - SQRT(s15_CONCH4PRE))
+    - s15_OVERLFAC1 * LOG(1 + s15_OVERLFAC2 *(v15_conc(ta10,'CH4') * s15_CONN2OPRE)**s15_OVERLEXP1
+	                        + s15_OVERLFAC3 * v15_conc(ta10,'CH4') * (v15_conc(ta10,'CH4')*s15_CONN2OPRE)**s15_OVERLEXP2
+							)
+    + s15_OVERLFAC1 * LOG(1 + s15_OVERLFAC2 * (s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP1
+	                        + s15_OVERLFAC3 * s15_CONCH4PRE * (s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP2
+							);
+							
+***---------------------------------------------------------------------------
+*'    N2O radiative forcing (from ACC2)
+***---------------------------------------------------------------------------
+q15_forcN2OQ(ta10)..
+    v15_forcComp(ta10,'N2O')
+	=e=
+	s15_RHON2O * (SQRT(v15_conc(ta10,'N2O'))-SQRT(s15_CONN2OPRE))
+    - s15_OVERLFAC1 * LOG(1 + s15_OVERLFAC2 * (s15_CONCH4PRE*v15_conc(ta10,'N2O'))**s15_OVERLEXP1
+	                      + s15_OVERLFAC3 * s15_CONCH4PRE * (s15_CONCH4PRE*v15_conc(ta10,'N2O'))**s15_OVERLEXP2
+						  )
+    + s15_OVERLFAC1 * LOG(1 + s15_OVERLFAC2 * (s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP1
+	                      + s15_OVERLFAC3 * s15_CONCH4PRE * (s15_CONCH4PRE*s15_CONN2OPRE)**s15_OVERLEXP2
+						  );
+
+***---------------------------------------------------------------------------
+*'    SO2 radiative forcing
+***---------------------------------------------------------------------------
+q15_forcso2(ta10)..
+    v15_forcComp(ta10,'SO2')
+	=e=
+	s15_dso1990 * v15_emi(ta10,'SO2') / s15_so1990 
+	+ s15_iso1990 * log(1 + v15_emi(ta10,'SO2')/s15_enatso2) / log(1 + s15_so1990/s15_enatso2) ;
 									   
-*** BC and OC from fossil fuels forcing
-
-	q15_forcbc(ta10)..
-					v15_forcComp(ta10,'BC') =e= (v15_emi(ta10,'BC') / s15_bc2005) * p15_oghgf_ffbc('2005');
+***---------------------------------------------------------------------------
+*'    BC and OC from fossil fuels radiative forcing; scales linear with emissions.
+***---------------------------------------------------------------------------
+q15_forcbc(ta10)..
+	v15_forcComp(ta10,'BC')
+	=e=
+	(v15_emi(ta10,'BC') / s15_bc2005) * p15_oghgf_ffbc('2005');
 					
-	q15_forcoc(ta10)..
-					v15_forcComp(ta10,'OC') =e= (v15_emi(ta10,'OC') / s15_oc2005) * p15_oghgf_ffoc('2005');					
-
-
-
-***     co2 radiative forcing module
-
-  q15_forcco2(ta10)..
-                   v15_forcComp(ta10,'CO2') =e= s15_fcodb * log(v15_conc(ta10,'CO2')/s15_c0) / log(2) ;
-
-***     total radiative forcing (foghg given exogenously)
-
-  q15_forctotal(ta10)..
-                   v15_forcComp(ta10,'TTL') =e= v15_forcComp(ta10,'CO2') +
-                                v15_forcComp(ta10,'SO2') +
-                                v15_forcComp(ta10,'CH4') +
-                                v15_forcComp(ta10,'N2O') +
-                                v15_forcComp(ta10,'BC') + 
-                                v15_forcComp(ta10,'OC') + 
-                                v15_forcComp(ta10,'oghg_nokyo') +
-                                v15_forcComp(ta10,'oghg_kyo');
+q15_forcoc(ta10)..
+	v15_forcComp(ta10,'OC')
+	=e=
+	(v15_emi(ta10,'OC') / s15_oc2005) * p15_oghgf_ffoc('2005');					
+
+***---------------------------------------------------------------------------
+*'    Total radiative forcing (foghg given exogenously)
+***---------------------------------------------------------------------------
+q15_forctotal(ta10)..
+    v15_forcComp(ta10,'TTL')
+	=e=
+	hv15_forcComp(ta10,'CO2') 
+	+ v15_forcComp(ta10,'SO2') 
+	+ v15_forcComp(ta10,'CH4') 
+	+ v15_forcComp(ta10,'N2O') 
+	+ v15_forcComp(ta10,'BC') 
+	+ v15_forcComp(ta10,'OC') 
+	+ v15_forcComp(ta10,'oghg_nokyo') 
+	+ v15_forcComp(ta10,'oghg_kyo');
 								
-*GL: Total forcing of Kyoto gases - required for formulation of policy targets
-  q15_forc_kyo(ta10)..
-                   v15_forcKyo(ta10) =e= v15_forcComp(ta10,'CO2') +
-                                v15_forcComp(ta10,'CH4') +
-                                v15_forcComp(ta10,'N2O') +
-                                 v15_forcComp(ta10,'oghg_kyo');
-
-*JeS* rcp forcing
-	q15_forc_rcp(ta10)..
-					v15_forcRcp(ta10) =e= v15_forcComp(ta10,'CO2') +
-                                v15_forcComp(ta10,'CH4') +
-                                v15_forcComp(ta10,'N2O') +
-								v15_forcComp(ta10,'oghg_kyo') + 
-								v15_forcComp(ta10,'SO2') + 
-								v15_forcComp(ta10,'BC') + 
-								v15_forcComp(ta10,'OC') +
-								v15_forcComp(ta10,'oghg_nokyo_rcp');
+***---------------------------------------------------------------------------
+*'    Total radiative forcing of Kyoto gases
+***---------------------------------------------------------------------------
+q15_forc_kyo(ta10)..
+    v15_forcKyo(ta10) 
+	=e=
+	v15_forcComp(ta10,'CO2') 
+	+ v15_forcComp(ta10,'CH4') 
+	+ v15_forcComp(ta10,'N2O') 
+	+ v15_forcComp(ta10,'oghg_kyo');
+
+***---------------------------------------------------------------------------
+*'    RCP forcing
+***---------------------------------------------------------------------------
+q15_forc_rcp(ta10)..
+	v15_forcRcp(ta10)
+	=e=
+	v15_forcComp(ta10,'CO2') 
+	+ v15_forcComp(ta10,'CH4') 
+	+ v15_forcComp(ta10,'N2O') 
+	+ v15_forcComp(ta10,'oghg_kyo') 
+	+ v15_forcComp(ta10,'SO2') 
+	+ v15_forcComp(ta10,'BC') 
+	+ v15_forcComp(ta10,'OC') 
+	+ v15_forcComp(ta10,'oghg_nokyo_rcp');
 																
-***     new temperature equations
-
-q15_clisys1(ta10) $(Ord(ta10) > 1) .. (v15_tempFast(ta10)-v15_tempFast(ta10-1))/s15_deltat_box =E= 0.5/s15_RPCTT1*
-                                         ((s15_RPCTA1*v15_forcComp(ta10,'TTL')/3.7-v15_tempFast(ta10))+
-                                          (s15_RPCTA1*v15_forcComp(ta10-1,'TTL')/3.7-v15_tempFast(ta10-1)));
-
-q15_clisys2(ta10) $(Ord(ta10) > 1) .. (v15_tempSlow(ta10)-v15_tempSlow(ta10-1))/s15_deltat_box =E= 0.5/s15_RPCTT2*
-                                         ((s15_RPCTA2*v15_forcComp(ta10,'TTL')/3.7-v15_tempSlow(ta10))+
-                                          (s15_RPCTA2*v15_forcComp(ta10-1,'TTL')/3.7-v15_tempSlow(ta10-1)));
-
-q15_clisys(ta10) .. v15_temp(ta10) =E= v15_tempFast(ta10) + v15_tempSlow(ta10);
-
-q15_clisys01(ta10) $(ord(ta10)=1)..
-                    v15_tempFast(ta10) =e= s15_RPCTA1 * s15_temp2000/s15_tsens;
-
-q15_clisys02(ta10) $(ord(ta10)=1)..
-                    v15_tempSlow(ta10) =e= s15_RPCTA2 * s15_temp2000/s15_tsens;
-
-*JeS* forcing overshoot for damage function
+***---------------------------------------------------------------------------
+*'    Temperature equations, consisting of a fast and a slow response function
+***---------------------------------------------------------------------------
+q15_clisys01(ta10)$(ord(ta10)=1)..
+    v15_tempFast(ta10)
+	=e=
+	s15_RPCTA1 * s15_temp2000 / s15_tsens;
+
+q15_clisys02(ta10)$(ord(ta10)=1)..
+    v15_tempSlow(ta10)
+	=e=
+	s15_RPCTA2 * s15_temp2000 / s15_tsens;
+
+q15_clisys1(ta10)$(Ord(ta10) > 1)..
+    (v15_tempFast(ta10) - v15_tempFast(ta10-1)) / s15_deltat_box
+    =e=
+    0.5 / s15_RPCTT1 * ( (s15_RPCTA1 * v15_forcComp(ta10,'TTL') / 3.7 - v15_tempFast(ta10))
+	                   + (s15_RPCTA1 * v15_forcComp(ta10-1,'TTL') / 3.7 - v15_tempFast(ta10-1))
+					   );
+
+q15_clisys2(ta10)$(Ord(ta10) > 1)..
+    (v15_tempSlow(ta10) - v15_tempSlow(ta10-1)) / s15_deltat_box
+	=e=
+	0.5 / s15_RPCTT2 * ( (s15_RPCTA2 * v15_forcComp(ta10,'TTL') / 3.7 - v15_tempSlow(ta10))
+	                   + (s15_RPCTA2 * v15_forcComp(ta10-1,'TTL') / 3.7 - v15_tempSlow(ta10-1))
+					   );
+
+q15_clisys(ta10)..
+    v15_temp(ta10)
+	=e=
+	v15_tempFast(ta10) + v15_tempSlow(ta10);
+
+***---------------------------------------------------------------------------
+*'    Forcing overshoot for damage function
+***---------------------------------------------------------------------------
 q15_forc_os(t)..
-        vm_forcOs(t) =E= v15_forcKyo(t)-s15_gr_forc_kyo+v15_slackForc(t);
-*** end of stuff taken from box-model.inc ---------------------------
-
-*-----------------------------------------------------------------------------------------
-*** ------------------------link to core--------------------------------------------------
-*-----------------------------------------------------------------------------------------
-
-q15_linkEMI(ttot2ta10(ttot, ta10),emis2climate10(enty,FOB10)).. vm_emiAllGlob(ttot,enty) =e= v15_emi(ta10,FOB10);
+    vm_forcOs(t)
+	=e=
+	v15_forcKyo(t) - s15_gr_forc_kyo + v15_slackForc(t);
+
+***-----------------------------------------------------------------------------------------
+*'     link to core
+***-----------------------------------------------------------------------------------------
+q15_linkEMI(ttot2ta10(ttot, ta10),emis2climate10(enty,FOB10))..
+    vm_emiAllGlob(ttot,enty)
+	=e=
+	v15_emi(ta10,FOB10);
 $IF %cm_so2_out_of_opt% == "on" q15_linkEMI_aer(ttot2ta10(ttot, ta10),emiaer2climate10(emiaer,FOB10)).. p15_so2emi(ttot,emiaer) =e= v15_emi(ta10,FOB10);
 
-*---interpolation for linking (annual resolution of climate modules)
-q15_interEMI(ta2ttot10(ta10, ttot),FOBEMI(FOB10)).. v15_emi(ta10,FOB10) =e= (1-p15_interpol(ta10)) * v15_emi(ttot,FOB10) + p15_interpol(ta10) *  v15_emi(ttot+1,FOB10);
+***-----------------------------------------------------------------------------------------
+*'     interpolation for linking (annual resolution of climate modules)
+***-----------------------------------------------------------------------------------------
+q15_interEMI(ta2ttot10(ta10, ttot),FOBEMI(FOB10))..
+    v15_emi(ta10,FOB10)
+	=e=
+	(1-p15_interpol(ta10)) * v15_emi(ttot,FOB10) + p15_interpol(ta10) *  v15_emi(ttot+1,FOB10);
 
 *** EOF ./modules/15_climate/box/equations.gms
diff --git a/modules/15_climate/box/realization.gms b/modules/15_climate/box/realization.gms
index 06cd376cc7592ebf22b103d09447a4ee5e6d7cfc..5ed71b30cd095a8a7012d1db3b9e020174705743 100644
--- a/modules/15_climate/box/realization.gms
+++ b/modules/15_climate/box/realization.gms
@@ -6,6 +6,9 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/15_climate/box.gms
 
+*' @description 
+*' In this realization, concentration, forcing, and temperature values are calculated using a simple model that can be used within the optimization routine.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "sets" $include "./modules/15_climate/box/sets.gms"
 $Ifi "%phase%" == "declarations" $include "./modules/15_climate/box/declarations.gms"
diff --git a/modules/15_climate/magicc/postsolve.gms b/modules/15_climate/magicc/postsolve.gms
index 46885aae531bd6804600cb93e779f6ff440c6a8c..77539f6e0c08ea369c9e1ab8ff3fbc989007bf5f 100644
--- a/modules/15_climate/magicc/postsolve.gms
+++ b/modules/15_climate/magicc/postsolve.gms
@@ -6,45 +6,50 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/15_climate/magicc/postsolve.gms
 
-*cb 20140218 HERE goes the code for the iterative adjustment of the emission budget for SSP runs
-*** emission budgets are adjusted, such that a predefined forcing target in 2100 is met
-*** the actual 2100 forcing after each iteration is calculated by a magicc run started from GAMS
+***---------------------------------------------------------------------------
+*' HERE goes the code for the iterative adjustment of the emission budget for SSP runs
+*' emission budgets are adjusted, such that a predefined forcing target in 2100 is met
+*' the actual 2100 forcing after each iteration is calculated by a magicc run started from GAMS
+***---------------------------------------------------------------------------
 
 *** Generate MAGICC scenario file
 $include "./core/magicc.gms";
-
-* execute MAGICC (this is cheap enough, ~2s)
+*** execute MAGICC (this is cheap enough, ~2s)
 Execute "Rscript run_magicc.R";
-* read in results
+*** read in results
 Execute "Rscript read_DAT_TOTAL_ANTHRO_RF.R";
 Execute_Loadpoint 'p15_forc_magicc'  p15_forc_magicc;
 Execute "Rscript read_DAT_SURFACE_TEMP.R";
 Execute_Loadpoint 'p15_magicc_temp' pm_globalMeanTemperature = pm_globalMeanTemperature;
-* MAGICC only reports unitl 2300:
+*** MAGICC only reports unitl 2300:
 pm_globalMeanTemperature(tall)$(tall.val gt 2300) = 0;
 
-
+***---------------------------------------------------------------------------
+*' calibrate temperature (GMT anomaly) to match HADCRUT4 in 2000. 
+*' This ensures that different MAGICC configurations start at the same observed temperature.
+***---------------------------------------------------------------------------
 $ifthen.cm_magicc_calibrateTemperature2000 %cm_magicc_calibrateTemperature2000% == "HADCRUT4"
-*AJS* calibrate temperature (GMT anomaly) to match HADCRUT4 in 2000. This ensures that different MAGICC configurations start at the same observed temperature.
-*OLD* The HADCRUT4 offset from 1861-1880 (the SR1.5 reference period) to 1985-2015 is  0.67 degree Celsius (median).
-*UPDATE* Use the 2010 offset relative to 2006-2015 reference period (AR6SR Table 2.2, footnote 1; email GL from 19122018)
+
+***---------------------------------------------------------------------------
+*' Calibrate temperature such that anomaly in 2006-2015 reference period is 0.97 (SR1.5 Table 2.2, footnote 1)
+***---------------------------------------------------------------------------
 s15_tempOffset2010 = sum(tall$(tall.val gt 2005 and tall.val le 2015),pm_globalMeanTemperature(tall))/10; 
 display s15_tempOffset2010;
 pm_globalMeanTemperature(tall) = pm_globalMeanTemperature(tall) - s15_tempOffset2010 + 0.97;
 display pm_globalMeanTemperature;
 
-*temperature convergence indicator
+*** temperature convergence indicator
 p15_gmt_conv = 100*smax(t,abs(pm_globalMeanTemperature(t)/max(p15_gmt0(t),1e-8) -1));
 display p15_gmt_conv;
-*save temp from last iteration
+*** save temp from last iteration
 p15_gmt0(tall) = pm_globalMeanTemperature(tall);
 
 $endif.cm_magicc_calibrateTemperature2000
 
-* offset from HADCRUT4 to zero temperature in 1900, instead of the default 1870 (20 year averages each).
+*** offset from HADCRUT4 to zero temperature in 1900, instead of the default 1870 (20 year averages each).
 pm_globalMeanTemperatureZeroed1900(tall)  = pm_globalMeanTemperature(tall) + 0.092; 
 
-* derive temperature impulse response (TIRF) from MAGICC pulse scenarios
+*** derive temperature impulse response (TIRF) from MAGICC pulse scenarios
 $ifthen.cm_magicc_tirf "%cm_magicc_temperatureImpulseResponse%" == "on"
 * the TIRF does not change much with emissions profile (see, e.g., figure in Schultes et al. 2017); 
 * thus only compute TIRF after each of the first 10 iterations, then only every fifth iteration. 
@@ -56,16 +61,22 @@ if( ((iteration.val le 10) or ( mod(iteration.val,5 ) eq 0)) ,
 *NOTE the MAGICC results (*.OUT files) are from  the last pulse experiment now, so take care if reading them in after this point.
 $endif.cm_magicc_tirf
 
-*** Iterative adjustment based on 
+***---------------------------------------------------------------------------
+*' Iterative adjustment of budgets or carbon taxes to meet forcing target 
+***---------------------------------------------------------------------------
 if (cm_iterative_target_adj eq 2, !! otherwise adjustment happens in core/postsolve.gms 
   
-*** Iterative adjustment for budget runs  
+***---------------------------------------------------------------------------
+*' Iterative adjustment for budget runs: scale current budget with the ratio of target forcing s15_gr_forc_os to current forcing p15_forc_magicc.
+*' The offset is only there to increase the speed of convergence, the values have no physical meaning.
+*' For low stabilization targets (rcp2.0, rcp2.6, rcp3.7) the target is the 2100 forcing target (s15_rcpCluster eq 1),
+*' for lower targets the forcing target is valid during the full century (s15_rcpCluster eq 0).
+***---------------------------------------------------------------------------
   if ((cm_emiscen eq 6),
    
       display sm_budgetCO2eqGlob, s15_gr_forc_os, p15_forc_magicc;
    
       if (s15_rcpCluster eq 1,
-display ' Liebe Lavinia ';
         sm_budgetCO2eqGlob 
         = 
           sm_budgetCO2eqGlob 
@@ -98,7 +109,12 @@ display ' Liebe Lavinia ';
      display sm_budgetCO2eqGlob;
      );
 
-*** Iterative adjustments for tax runs
+***---------------------------------------------------------------------------
+*' Iterative adjustment for carbon tax runs: scale current tax pathway with the ratio of target forcing s15_gr_forc_os to current forcing p15_forc_magicc.
+*' The offset is only there to increase the speed of convergence, the values have no physical meaning.
+*' For low stabilization targets (rcp2.0, rcp2.6, rcp3.7) the target is the 2100 forcing target (s15_rcpCluster eq 1),
+*' for lower targets the forcing target is valid during the full century (s15_rcpCluster eq 0).
+***---------------------------------------------------------------------------
   if (cm_emiscen eq 9, 
   
     display pm_taxCO2eq, s15_gr_forc_os, p15_forc_magicc;
diff --git a/modules/15_climate/magicc/realization.gms b/modules/15_climate/magicc/realization.gms
index 56e3656c2ab1cffafad212d71ac2d5a559c98427..06246441dae14f84e21256bdf60963210f4e7cb5 100644
--- a/modules/15_climate/magicc/realization.gms
+++ b/modules/15_climate/magicc/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/15_climate/magicc.gms
 
+*' @description 
+*' In this realization, concentration, forcing, and temperature values are calculated using a MAGICC6.4. 
+*' MAGICC is run in between iterations and can be used to adapt carbon tax pathways and budgets to meet a give climate target.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "sets" $include "./modules/15_climate/magicc/sets.gms"
 $Ifi "%phase%" == "declarations" $include "./modules/15_climate/magicc/declarations.gms"
diff --git a/modules/15_climate/module.gms b/modules/15_climate/module.gms
index 573400ef283de0406afe554ffe0b379c383da5b7..5c09159d12ec5de2a221a1a316f4d07b5d6bfee1 100644
--- a/modules/15_climate/module.gms
+++ b/modules/15_climate/module.gms
@@ -6,6 +6,12 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/15_climate/15_climate.gms
 
+*' @title climate
+*'
+*' @description  The 15_climate module calculates the resulting climate variables using either MAGICC6.4 or a stylized box model that can be used within the optimization routine.
+*'
+*' @authors Jessica Strefler, Michaja Pehl, Christoph Bertram
+
 *###################### R SECTION START (MODULETYPES) ##########################
 $Ifi "%climate%" == "box" $include "./modules/15_climate/box/realization.gms"
 $Ifi "%climate%" == "magicc" $include "./modules/15_climate/magicc/realization.gms"
diff --git a/modules/21_tax/on/equations.gms b/modules/21_tax/on/equations.gms
index 8a3cf98e27895543fe3e3121fc126d09a937c5b3..5baf7adc3c6de62c15f952addf8e5b322036da18 100644
--- a/modules/21_tax/on/equations.gms
+++ b/modules/21_tax/on/equations.gms
@@ -63,7 +63,7 @@ v21_taxrevCO2luc(t,regi) =g= ( pm_taxCO2eq(t,regi)  + pm_taxCO2eqSCC(t,regi) + p
                            - p21_taxrevCO2LUC0(t,regi);
 
 ***---------------------------------------------------------------------------
-*'  Calculation of CCS tax: tax rate (defined as fraction(or multiplier) of O&M costs) times sequestration
+*'  Calculation of CCS tax: tax rate (defined as fraction(or multiplier) of O&M costs) times amount of CO2 sequestration
 *'  Documentation of overall tax approach is above at q21_taxrev.
 ***---------------------------------------------------------------------------
 q21_taxrevCCS(t,regi)$(t.val ge max(2010,cm_startyear))..
@@ -74,7 +74,7 @@ v21_taxrevCCS(t,regi)
 	- p21_taxrevCCS0(t,regi);
 
 ***---------------------------------------------------------------------------
-*'  Calculation of net negative emissions tax: tax rate (defined as fraction of carbon price) times negative emissions
+*'  Calculation of net-negative emissions tax: tax rate (defined as fraction of carbon price) times net-negative emissions
 *'  Documentation of overall tax approach is above at q21_taxrev.
 ***---------------------------------------------------------------------------
 q21_taxrevNetNegEmi(t,regi)$(t.val ge max(2010,cm_startyear))..
@@ -82,7 +82,7 @@ v21_taxrevNetNegEmi(t,regi) =g=  cm_frac_NetNegEmi * pm_taxCO2eq(t,regi) * v21_e
                                  - p21_taxrevNetNegEmi0(t,regi);
 
 ***---------------------------------------------------------------------------
-*'  Auxiliary calculation of negative emissions: 
+*'  Auxiliary calculation of net-negative emissions: 
 *'  v21_emiAllco2neg and v21_emiAllco2neg_slack are defined as positive variables
 *'  so as long as vm_emiAll is positive, v21_emiAllco2neg_slack adjusts so that sum is zero
 *'  if vm_emiAll is negative, in order to minimize tax v21_emiAllco2neg_slack becomes zero
diff --git a/modules/33_CDR/DAC/equations.gms b/modules/33_CDR/DAC/equations.gms
index df0bcc2f40c4b7e5df666849fb063748b7e934c3..b87ab3e58aebc311b87ebb732991b8b7343805a4 100644
--- a/modules/33_CDR/DAC/equations.gms
+++ b/modules/33_CDR/DAC/equations.gms
@@ -5,7 +5,11 @@
 *** |  REMIND License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/DAC/equations.gms
-*** new default: gas for heat production is with CCS; assume 90% capture rate.
+
+***---------------------------------------------------------------------------
+*'  Calculation of (negative) CO2 emissions from direct air capture. The first part of the equation describes emissions captured from the ambient air, 
+*'  the second part calculates the CO2 captured from the gas used for heat production assuming 90% capture rate.
+***---------------------------------------------------------------------------
 q33_capconst_dac(t,regi)..
 	v33_emiDAC(t,regi)
 	=e=
@@ -13,26 +17,37 @@ q33_capconst_dac(t,regi)..
 	-  (1 / pm_eta_conv(t,regi,"gash2c")) * fm_dataemiglob("pegas","seh2","gash2c","cco2") * vm_otherFEdemand(t,regi,"fegas")
 	;
 
+***---------------------------------------------------------------------------
+*'  Sum of all CDR emissions other than BECCS and afforestation, which are calculated in the core.
+***---------------------------------------------------------------------------
 q33_emicdrregi(t,regi)..
 	vm_emiCdr(t,regi,"co2")
 	=e=
 	v33_emiDAC(t,regi);
 	
-*** fehes is only district heat. Use gas or h2 to generate heat for DAC plants. 
-*** vm_otherFEdemand(t,regi,"fegas") is calculated as the total energy demand for heat minus what is already covered by h2, i.e. vm_otherFEdemand(t,regi,"feh2s") and vice versa.
+***---------------------------------------------------------------------------
+*'  Calculation of energy demand of direct air capture. Heat can be provided by gas or H2; 
+*'  vm_otherFEdemand(t,regi,"fegas") is calculated as the total energy demand for heat minus what is already covered by h2, i.e. vm_otherFEdemand(t,regi,"feh2s") and vice versa.
+***---------------------------------------------------------------------------
 q33_otherFEdemand(t,regi,entyFe)..
     vm_otherFEdemand(t,regi,entyFe)
     =e=
     - vm_emiCdr(t,regi,"co2") * sm_EJ_2_TWa * p33_dac_fedem(entyFe)
     - vm_otherFEdemand(t,regi,"feh2s")$(sameas(entyFe,"fegas")) - vm_otherFEdemand(t,regi,"fegas")$(sameas(entyFe,"feh2s"))
     ;
-	
+
+***---------------------------------------------------------------------------
+*'  Preparation of captured emissions to enter the CCS chain.
+***---------------------------------------------------------------------------	
 q33_ccsbal(t,regi,ccs2te(ccsCo2(enty),enty2,te))..
 	sum(teCCS2rlf(te,rlf), vm_ccs_cdr(t,regi,enty,enty2,te,rlf))
 	=e=
 	-vm_emiCdr(t,regi,"co2")
 	;
 
+***---------------------------------------------------------------------------
+*'  Limit the amount of H2 from biomass to the demand without DAC.
+***---------------------------------------------------------------------------
 q33_H2bio_lim(t,regi,te)..	         
 	vm_prodSE(t,regi,"pebiolc","seh2",te)$pe2se("pebiolc","seh2",te)
 	=l=
diff --git a/modules/33_CDR/DAC/realization.gms b/modules/33_CDR/DAC/realization.gms
index 97af2c6b6cd3ad083f5d6054313a4e07fff83ac9..a6a4f870d6857f1a3f01e4baf09f9131c6cc60da 100644
--- a/modules/33_CDR/DAC/realization.gms
+++ b/modules/33_CDR/DAC/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/DAC.gms
 
+*' @description 
+*' In this realization, direct air capture can be used to remove CO2 from the atmosphere in addition to BECCS and afforestation. Based on Broehm et al. we assume an energy demand of 
+*' 2 GJ/tCO2 electricity and 10 GJ/tCO2 heat which can be met via gas or H2. If gas is used, the resulting CO2 is captured with a capture rate of 90%.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "sets" $include "./modules/33_CDR/DAC/sets.gms"
 $Ifi "%phase%" == "declarations" $include "./modules/33_CDR/DAC/declarations.gms"
diff --git a/modules/33_CDR/DAC/sets.gms b/modules/33_CDR/DAC/sets.gms
index 03262a27dd4f6e7f68947e90496da99b7965cd0a..a6fb8f62aa8953bc6096eff729704bba86d457d5 100644
--- a/modules/33_CDR/DAC/sets.gms
+++ b/modules/33_CDR/DAC/sets.gms
@@ -7,7 +7,7 @@
 *** SOF ./modules/33_CDR/DAC/sets.gms
 sets
 
-te_dyn33(all_te)  "???"
+te_dyn33(all_te)  "all technologies"
 /
 		dac		"direct air capture"
 /
diff --git a/modules/33_CDR/all/equations.gms b/modules/33_CDR/all/equations.gms
index df9c9bf19c337d3a3fa455e0e1863a1a4bdb8fb2..efefaa8c05cd18aad6230b78f9827f5bf2543df7 100644
--- a/modules/33_CDR/all/equations.gms
+++ b/modules/33_CDR/all/equations.gms
@@ -6,13 +6,19 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/all/equations.gms
 
+***---------------------------------------------------------------------------
+*'  Calculation of the amount of ground rock spread in timestep t.
+***---------------------------------------------------------------------------
 q33_capconst_grindrock(t,regi)..
 	sum(rlf2,sum(rlf, v33_grindrock_onfield(t,regi,rlf,rlf2)))
 	=l=
 	sum(teNoTransform2rlf_dyn33(te,rlf2), vm_capFac(t,regi,"rockgrind") * vm_cap(t,regi,"rockgrind",rlf2))
 	;
 	
-* JeS: timestep in 2060 not yet quite right!  
+***---------------------------------------------------------------------------
+*'  Calculation of the total amount of ground rock on the fields in timestep t. The first part of the equation describes the decay of the rocks added until that time,
+*'  the rest describes the newly added rocks.
+***---------------------------------------------------------------------------
 q33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)$(ttot.val ge max(2010, cm_startyear))..
 	v33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)
 	=e=
@@ -21,6 +27,9 @@ q33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)$(ttot.val ge max(2010, cm_startyea
 	v33_grindrock_onfield(ttot,regi,rlf,rlf2) * (sum(tall $ ((tall.val le ttot.val) $ (tall.val gt (ttot.val-pm_ts(ttot)/2))),exp(-p33_co2_rem_rate(rlf) * (ttot.val-tall.val))))
 ;  
 
+***---------------------------------------------------------------------------
+*'  Calculation of (negative) CO2 emissions from enhanced weathering. 
+***---------------------------------------------------------------------------
 q33_emiEW(t,regi)..
 	v33_emiEW(t,regi)
 	=e=
@@ -29,20 +38,34 @@ q33_emiEW(t,regi)..
 		)
 	;	
 
-*** new default: gas for heat production is with CCS; assume 90% capture rate.	
+***---------------------------------------------------------------------------
+*'  Calculation of (negative) CO2 emissions from direct air capture. The first part of the equation describes emissions captured from the ambient air, 
+*'  the second part calculates the CO2 captured from the gas used for heat production assuming 90% capture rate.
+***---------------------------------------------------------------------------
 q33_capconst_dac(t,regi)..
 	v33_emiDAC(t,regi)
 	=e=
 	-sum(teNoTransform2rlf_dyn33(te,rlf2), vm_capFac(t,regi,"dac") * vm_cap(t,regi,"dac",rlf2))
 	-  (1 / pm_eta_conv(t,regi,"gash2c")) * fm_dataemiglob("pegas","seh2","gash2c","cco2") * vm_otherFEdemand(t,regi,"fegas")	
 	;
-	
+
+***---------------------------------------------------------------------------
+*'  Sum of all CDR emissions other than BECCS and afforestation, which are calculated in the core.
+***---------------------------------------------------------------------------	
 q33_emicdrregi(t,regi)..
 	vm_emiCdr(t,regi,"co2")
 	=e=
 	v33_emiEW(t,regi) + v33_emiDAC(t,regi)
 	;
-	
+
+***---------------------------------------------------------------------------
+*'  Calculation of energy demand of DAC and EW. 
+*'  The first part of the equation describes the electricity demand for grinding, 
+*'  the second part the diesel demand for transportation and spreading on crop fields.
+*'  The third part is DAC electricity demand and the last part is DAC heat demand.
+*'  Heat for DAC can be provided by gas or H2; vm_otherFEdemand(t,regi,"fegas") is calculated as the total 
+*'  DAC energy demand for heat minus what is already covered by h2, i.e. vm_otherFEdemand(t,regi,"feh2s") and vice versa.
+***---------------------------------------------------------------------------	
 q33_otherFEdemand(t,regi,entyFe)..
 	vm_otherFEdemand(t,regi,entyFe)
 	=e=
@@ -52,12 +75,18 @@ q33_otherFEdemand(t,regi,entyFe)..
    - vm_otherFEdemand(t,regi,"feh2s")$(sameas(entyFe,"fegas")) - vm_otherFEdemand(t,regi,"fegas")$(sameas(entyFe,"feh2s"))
 	;	
 	
+***---------------------------------------------------------------------------
+*'  Limit the amount of H2 from biomass to the demand without DAC.
+***---------------------------------------------------------------------------	
 q33_H2bio_lim(t,regi,te)..	         
 	vm_prodSE(t,regi,"pebiolc","seh2",te)$pe2se("pebiolc","seh2",te)
 	=l=
     vm_prodFe(t,regi,"seh2","feh2s","tdh2s") - vm_otherFEdemand(t,regi,"feh2s")
 	;
 
+***---------------------------------------------------------------------------
+*'  O&M costs of EW, consisting of fix costs for mining, grinding and spreading, and transportation costs.
+***---------------------------------------------------------------------------	
 q33_omcosts(t,regi)..
 	vm_omcosts_cdr(t,regi)
 	=e=
@@ -68,18 +97,27 @@ q33_omcosts(t,regi)..
 		)
 	)
 	;
-	
+
+***---------------------------------------------------------------------------
+*'  Limit total amount of ground rock on the fields to regional maximum potentials.
+***---------------------------------------------------------------------------		
 q33_potential(t,regi,rlf)..	
 	sum(rlf2,v33_grindrock_onfield_tot(t,regi,rlf,rlf2))
 	=l=
 	f33_maxProdGradeRegiWeathering(regi,rlf);	
-	
+
+***---------------------------------------------------------------------------
+*'  Preparation of captured emissions to enter the CCS chain.
+***---------------------------------------------------------------------------		
 q33_ccsbal(t,regi,ccs2te(ccsCo2(enty),enty2,te))..
 	sum(teCCS2rlf(te,rlf), vm_ccs_cdr(t,regi,enty,enty2,te,rlf))
 	=e=
 	-v33_emiDAC(t,regi)
 	;	
-  
+
+***---------------------------------------------------------------------------
+*'  An annual limit for the maximum amount of rocks spred [Gt] can be set via cm_LimRock, e.g. due to sustainability concerns.
+***---------------------------------------------------------------------------  
 q33_LimEmiEW(t,regi)..
              sum(rlf,
                   sum(rlf2,v33_grindrock_onfield(t,regi,rlf,rlf2))
diff --git a/modules/33_CDR/all/realization.gms b/modules/33_CDR/all/realization.gms
index c4566f79a100451cedd0a22cafd071de29595966..db4241e6d5e64bb08fd83ba38bebf43093de4cdb 100644
--- a/modules/33_CDR/all/realization.gms
+++ b/modules/33_CDR/all/realization.gms
@@ -6,6 +6,11 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/all.gms
 
+*' @description 
+*' In this realization, direct air capture and enhanced weathering can be used to remove CO2 from the atmosphere in addition to BECCS and afforestation. Based on Broehm et al. we assume an energy demand of 
+*' 2 GJ/tCO2 electricity and 10 GJ/tCO2 heat for DAC which can be met via gas or H2. If gas is used, the resulting CO2 is captured with a capture rate of 90%.
+*' For EW, electricty is needed to grind the rocks and diesel is needed for transportation and spreading on crop fields.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "sets" $include "./modules/33_CDR/all/sets.gms"
 $Ifi "%phase%" == "declarations" $include "./modules/33_CDR/all/declarations.gms"
diff --git a/modules/33_CDR/all/sets.gms b/modules/33_CDR/all/sets.gms
index c88fc53f870810a973902643f3955f2dc20f2929..5d071d50cd21b47fde038e758638ad96f97ab6b5 100644
--- a/modules/33_CDR/all/sets.gms
+++ b/modules/33_CDR/all/sets.gms
@@ -7,7 +7,7 @@
 *** SOF ./modules/33_CDR/all/sets.gms
 sets
 
-te_dyn33(all_te)  "???"
+te_dyn33(all_te)  "all technologies"
 /
 		rockgrind		"grinding rock for enhanced weathering"
 		dac		"direct air capture"
diff --git a/modules/33_CDR/module.gms b/modules/33_CDR/module.gms
index 1709d94eff45472097e34fc36fdcad74efaa18c9..73a4a52accf2510e0e30b454bbd57080c3e73f32 100644
--- a/modules/33_CDR/module.gms
+++ b/modules/33_CDR/module.gms
@@ -6,6 +6,12 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/33_CDR.gms
 
+*' @title CDR
+*'
+*' @description  The 33_CDR module calculates CO2 removed from the atmosphere by options other than BECCS or afforestation, which are calculated in the core.
+*'
+*' @authors Jessica Strefler
+
 *###################### R SECTION START (MODULETYPES) ##########################
 $Ifi "%CDR%" == "DAC" $include "./modules/33_CDR/DAC/realization.gms"
 $Ifi "%CDR%" == "all" $include "./modules/33_CDR/all/realization.gms"
diff --git a/modules/33_CDR/off/realization.gms b/modules/33_CDR/off/realization.gms
index dc409de831fca38b129851ac1e7e98d620fc0ad2..988df6d91a79e53c260cc954e7376f5721874ab0 100644
--- a/modules/33_CDR/off/realization.gms
+++ b/modules/33_CDR/off/realization.gms
@@ -6,6 +6,9 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/off.gms
 
+*' @description 
+*' In this realization, no additional CDR option other than BECCS and afforestation is available.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/33_CDR/off/declarations.gms"
 $Ifi "%phase%" == "bounds" $include "./modules/33_CDR/off/bounds.gms"
diff --git a/modules/33_CDR/weathering/equations.gms b/modules/33_CDR/weathering/equations.gms
index ec21906251d2839fdedd1f8119efe49014e73733..4da7f9e13fae964cae09cb29594881325eb35f11 100644
--- a/modules/33_CDR/weathering/equations.gms
+++ b/modules/33_CDR/weathering/equations.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/weathering/equations.gms
 
+***---------------------------------------------------------------------------
+*'  Calculation of the energy demand of enhanced weathering. The first part of the equation describes the electricity demand for grinding, 
+*'  the second part the diesel demand for transportation and spreading on crop fields.
+***---------------------------------------------------------------------------
 q33_otherFEdemand(t,regi,entyFe)$(sameas(entyFe,"feels") OR sameas(entyFe,"fedie"))..
 	vm_otherFEdemand(t,regi,entyFe)$(sameas(entyFe,"feels") OR sameas(entyFe,"fedie"))
 	=e=
@@ -13,13 +17,19 @@ q33_otherFEdemand(t,regi,entyFe)$(sameas(entyFe,"feels") OR sameas(entyFe,"fedie
    + sum(rlf$(rlf.val le 2), s33_rockfield_fedem$(sameas(entyFe,"fedie")) * sm_EJ_2_TWa * sum(rlf2,v33_grindrock_onfield(t,regi,rlf,rlf2)))
 	;
 
+***---------------------------------------------------------------------------
+*'  Calculation of the amount of ground rock spread in timestep t.
+***---------------------------------------------------------------------------
 q33_capconst_grindrock(t,regi)..
 	sum(rlf2,sum(rlf$(rlf.val le 2), v33_grindrock_onfield(t,regi,rlf,rlf2)))
 	=l=
 	sum(teNoTransform2rlf_dyn33(te,rlf2), vm_capFac(t,regi,"rockgrind") * vm_cap(t,regi,"rockgrind",rlf2))
 	;
 	
-* JeS: timestep in 2060 not yet quite right!  
+***---------------------------------------------------------------------------
+*'  Calculation of the total amount of ground rock on the fields in timestep t. The first part of the equation describes the decay of the rocks added until that time,
+*'  the rest describes the newly added rocks.
+***---------------------------------------------------------------------------
 q33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)$((ttot.val ge max(2010, cm_startyear))$(rlf.val le 2))..
 	v33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)$(rlf.val le 2)
 	=e=
@@ -28,6 +38,9 @@ q33_grindrock_onfield_tot(ttot,regi,rlf,rlf2)$((ttot.val ge max(2010, cm_startye
 	v33_grindrock_onfield(ttot,regi,rlf,rlf2)$(rlf.val le 2) * (sum(tall $ ((tall.val le ttot.val) $ (tall.val gt (ttot.val-pm_ts(ttot)/2))),exp(-p33_co2_rem_rate(rlf)$(rlf.val le 2) * (ttot.val-tall.val))))
 ;  
 
+***---------------------------------------------------------------------------
+*'  Calculation of (negative) CO2 emissions from enhanced weathering. 
+***---------------------------------------------------------------------------
 q33_emiEW(t,regi)..
 	v33_emiEW(t,regi)
 	=e=
@@ -36,12 +49,18 @@ q33_emiEW(t,regi)..
 	)
 	;
 
+***---------------------------------------------------------------------------
+*'  Sum of all CDR emissions other than BECCS and afforestation, which are calculated in the core.
+***---------------------------------------------------------------------------
 q33_emicdrregi(t,regi)..
 	vm_emiCdr(t,regi,"co2")
 	=e=
 	v33_emiEW(t,regi)
 	;
-	
+
+***---------------------------------------------------------------------------
+*'  O&M costs of EW, consisting of fix costs for mining, grinding and spreading, and transportation costs.
+***---------------------------------------------------------------------------	
 q33_omcosts_onfield(t,regi)..
 	vm_omcosts_cdr(t,regi)
 	=e=
@@ -52,12 +71,18 @@ q33_omcosts_onfield(t,regi)..
 		)
 	)
 	;
-	
+
+***---------------------------------------------------------------------------
+*'  Limit total amount of ground rock on the fields to regional maximum potentials.
+***---------------------------------------------------------------------------	
 q33_potential(t,regi,rlf)$(rlf.val le 2)..	
 	sum(rlf2,v33_grindrock_onfield_tot(t,regi,rlf,rlf2)$(rlf.val le 2))
 	=l=
 	f33_maxProdGradeRegiWeathering(regi,rlf)$(rlf.val le 2);
 
+***---------------------------------------------------------------------------
+*'  An annual limit for the maximum amount of rocks spred [Gt] can be set via cm_LimRock, e.g. due to sustainability concerns.
+***---------------------------------------------------------------------------
 q33_LimEmiEW(t,regi)..
              sum(rlf,
                   sum(rlf2,v33_grindrock_onfield(t,regi,rlf,rlf2))
diff --git a/modules/33_CDR/weathering/realization.gms b/modules/33_CDR/weathering/realization.gms
index 7a11cdd92b6900cf866c69d5898dba715139c1cd..de5a4affb8d6efb0f285a6f43f31b0f447a8a8ad 100644
--- a/modules/33_CDR/weathering/realization.gms
+++ b/modules/33_CDR/weathering/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/33_CDR/weathering.gms
 
+*' @description 
+*' In this realization, enhanced weathering of rocks can be used to remove CO2 from the atmosphere in addition to BECCS and afforestation. 
+*' Electricty is needed to grind the rocks and diesel is needed for transportation and spreading on crop fields.
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "sets" $include "./modules/33_CDR/weathering/sets.gms"
 $Ifi "%phase%" == "declarations" $include "./modules/33_CDR/weathering/declarations.gms"
diff --git a/modules/33_CDR/weathering/sets.gms b/modules/33_CDR/weathering/sets.gms
index 159a043529ac5cbd11cd1079ccf345e0f28f14f2..acbb8acc93cbbf1dde4b397d8b9400122ddf35e7 100644
--- a/modules/33_CDR/weathering/sets.gms
+++ b/modules/33_CDR/weathering/sets.gms
@@ -7,7 +7,7 @@
 *** SOF ./modules/33_CDR/weathering/sets.gms
 sets
 
-te_dyn33(all_te)   "???"
+te_dyn33(all_te)   "all technologies"
 /
 		rockgrind		"grinding rock for enhanced weathering"
 /
diff --git a/modules/39_CCU/module.gms b/modules/39_CCU/module.gms
index fd7983dc793e80384f3b1ed0a99aaa0031d0f4e2..476e61236ebda7c80119b41a84206517d5bf3496 100644
--- a/modules/39_CCU/module.gms
+++ b/modules/39_CCU/module.gms
@@ -6,6 +6,12 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/39_CCU/39_CCU.gms
 
+*' @title CCU
+*'
+*' @description  The 39_CCU module calculates emissions from synthetic gas and liquids.
+*'
+*' @authors Laura Popin, Jessica Strefler
+
 *####################### R SECTION START (MODULETYPES) ##############################
 $Ifi "%CCU%" == "off" $include "./modules/39_CCU/off/realization.gms"
 $Ifi "%CCU%" == "on" $include "./modules/39_CCU/on/realization.gms"
diff --git a/modules/39_CCU/on/equations.gms b/modules/39_CCU/on/equations.gms
index 54c7c5cb9994b0af00d8de090c0b25b2df05310f..b0e053592076096d6ab5166806d2bf6739d2fc2e 100644
--- a/modules/39_CCU/on/equations.gms
+++ b/modules/39_CCU/on/equations.gms
@@ -8,9 +8,8 @@
 
 
 ***---------------------------------------------------------------------------
-*** LP
-*** Managing the C/H ratio in CCU-Technologies
-*** amount of C temporary used in CCU-products in relation to the amount of hydrogen necessary [GtC/y]
+*' Managing the C/H ratio in CCU-Technologies
+*' amount of C temporary used in CCU-products in relation to the amount of hydrogen necessary [GtC/y]
 ***---------------------------------------------------------------------------
 
 q39_emiCCU(t,regi) .. 
diff --git a/modules/41_emicapregi/AbilityToPay/realization.gms b/modules/41_emicapregi/AbilityToPay/realization.gms
index 1906f7b840b50d01313cde244ad5f7763cee234f..c6fa2500b1ceec7ed1cdcd9632ff175f4c2c2708 100644
--- a/modules/41_emicapregi/AbilityToPay/realization.gms
+++ b/modules/41_emicapregi/AbilityToPay/realization.gms
@@ -6,6 +6,9 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/AbilityToPay.gms
 
+*' @description
+*' Emission caps/permits are allocated according to the ability to pay principle  
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/AbilityToPay/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/AbilityToPay/datainput.gms"
diff --git a/modules/41_emicapregi/CandC/datainput.gms b/modules/41_emicapregi/CandC/datainput.gms
index 9bd775e59518eb2f0595b1ea57c74bd4a3056178..abc52a38d5d139ff3fffbe5526328232d386a21b 100644
--- a/modules/41_emicapregi/CandC/datainput.gms
+++ b/modules/41_emicapregi/CandC/datainput.gms
@@ -15,14 +15,6 @@ Execute_Loadpoint "input_ref" p41_co2eq = vm_co2eq.l;
 p41_shEmi2005(regi) = p41_co2eq("2005",regi) / sum(regi2, p41_co2eq("2005",regi2) );
 display p41_shEmi2005;
 
-*LB* safe formulation of c_polscen = 11
-$ontext
-***  contraction & convergence (reference year 2020)  +++++++++
-        p41_lambda(tall) $(tall.val<2050) = (tall.val-2020) / 30;
-        p41_lambda(tall) $(tall.val>2049) = 1;
-$offtext
-
-*gl* calculate share of global emissions according to different burden sharing rules
      pm_shPerm(t,regi) =  p41_lambda(t) * pm_pop(t,regi) / sum(regi2,pm_pop(t,regi2))
          + (1 - p41_lambda(t)) * p41_shEmi2005(regi) / sum(regi2, p41_shEmi2005(regi2));
 *** EOF ./modules/41_emicapregi/CandC/datainput.gms
diff --git a/modules/41_emicapregi/CandC/realization.gms b/modules/41_emicapregi/CandC/realization.gms
index f98e426555d0cb992618eee9137528fea7b3d6e6..8b42ff92735607950e46528bd1752d7fff75c70e 100644
--- a/modules/41_emicapregi/CandC/realization.gms
+++ b/modules/41_emicapregi/CandC/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/CandC.gms
 
+*' @description
+*' Emission caps/permits are allocated according to the contraction and convergence 
+*' rule (transition towards equal per capita allocation; with reference years 2005 and 2050)
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/CandC/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/CandC/datainput.gms"
diff --git a/modules/41_emicapregi/GDPint/datainput.gms b/modules/41_emicapregi/GDPint/datainput.gms
index 7a500bdc85f35dfebc419c6b9c28cd8be998cacb..5fee4af62d16c6cd8f62bad7668e3776b40a9526 100644
--- a/modules/41_emicapregi/GDPint/datainput.gms
+++ b/modules/41_emicapregi/GDPint/datainput.gms
@@ -5,7 +5,8 @@
 *** |  REMIND License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/GDPint/datainput.gms
-*gl* calculate share of global emissions according to different burden sharing rules
+
+*** calculate share of global emissions 
      pm_shPerm(t,regi)  =  pm_gdp_gdx(t,regi) / sum(regi2, pm_gdp_gdx(t,regi2));
 
 
diff --git a/modules/41_emicapregi/GDPint/realization.gms b/modules/41_emicapregi/GDPint/realization.gms
index e13d093f5e782fb0429467bb347490258513d56f..3218f37a8b1f56598e8def246bd3a8475fbc8a02 100644
--- a/modules/41_emicapregi/GDPint/realization.gms
+++ b/modules/41_emicapregi/GDPint/realization.gms
@@ -6,6 +6,9 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/GDPint.gms
 
+*' @description
+*' Emission caps/permits are allocated according to GDP intensity
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/GDPint/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/GDPint/datainput.gms"
diff --git a/modules/41_emicapregi/POPint/datainput.gms b/modules/41_emicapregi/POPint/datainput.gms
index b0177b42307dc5defc1937a148e01b13abb1f3b1..8dbde4f1cd5b46d9041b381e66b2e5d489b4f703 100644
--- a/modules/41_emicapregi/POPint/datainput.gms
+++ b/modules/41_emicapregi/POPint/datainput.gms
@@ -5,8 +5,8 @@
 *** |  REMIND License Exception, version 1.0 (see LICENSE file).
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/POPint/datainput.gms
-*gl* calculate share of global emissions according to different burden sharing rules
-*ML* cumulated POP shares
+
+*** calculate share of global emissions 
      pm_shPerm(t,regi)  =  sum(t2, pm_pop(t2,regi))/sum((t2,regi2),pm_pop(t2,regi2));
 
 *** EOF ./modules/41_emicapregi/POPint/datainput.gms
diff --git a/modules/41_emicapregi/POPint/realization.gms b/modules/41_emicapregi/POPint/realization.gms
index aee1a59e26fc63c705a10931bb738ac9e8d0e16f..24451006f65bd2b27345db6e9ce0caf5c2706f34 100644
--- a/modules/41_emicapregi/POPint/realization.gms
+++ b/modules/41_emicapregi/POPint/realization.gms
@@ -6,6 +6,11 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/POPint.gms
 
+*' @description
+*' Emission caps/permits are allocated according to each region's share on
+*' cumulated population 
+*' 
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/POPint/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/POPint/datainput.gms"
diff --git a/modules/41_emicapregi/PerCapitaConvergence/datainput.gms b/modules/41_emicapregi/PerCapitaConvergence/datainput.gms
index fa0e00ec1ca45100206b033525a1550d39728cfd..abc037eee1448b4811d0163a1479806f78999ff8 100644
--- a/modules/41_emicapregi/PerCapitaConvergence/datainput.gms
+++ b/modules/41_emicapregi/PerCapitaConvergence/datainput.gms
@@ -6,7 +6,7 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/PerCapitaConvergence/datainput.gms
 
-***   contraction & convergence (reference year 2005)  +++++++++
+***   contraction & convergence (reference year 2020)  +++++++++
 p41_convergenceyear = 2050;
 p41_lambda(tall) $(tall.val < p41_convergenceyear) = (tall.val - 2020) / (p41_convergenceyear-2020);
 p41_lambda(tall) $(tall.val>2049) = 1;
@@ -21,14 +21,7 @@ display p41_shEmi2020;
 *** calculate global emissions pathway in cost-optimal scenario
 pm_emicapglob(ttot) = sum(regi, p41_co2eq(ttot,regi));
 
-*LB* safe formulation of c_polscen = 11
-$ontext
-***  contraction & convergence (reference year 2020)  +++++++++
-        p41_lambda(tall) $(tall.val<2050) = (tall.val-2020) / 30;
-        p41_lambda(tall) $(tall.val>2049) = 1;
-$offtext
-
-*gl* calculate share of global emissions according to different burden sharing rules
+*** calculate share of global emissions 
      pm_shPerm(t,regi) =  p41_lambda(t) * pm_pop(t,regi) / sum(regi2,pm_pop(t,regi2))
          + (1 - p41_lambda(t)) * p41_shEmi2020(regi) / sum(regi2, p41_shEmi2020(regi2));
 
diff --git a/modules/41_emicapregi/PerCapitaConvergence/realization.gms b/modules/41_emicapregi/PerCapitaConvergence/realization.gms
index 63e098727c4fb9231a8b29a3a185c10649428c5b..c87dddf9fa8d4b552512e863b3a59b591eca0e07 100644
--- a/modules/41_emicapregi/PerCapitaConvergence/realization.gms
+++ b/modules/41_emicapregi/PerCapitaConvergence/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/PerCapitaConvergence.gms
 
+*' @description
+*' Emission caps/permits are allocated according to the contraction and convergence 
+*' rule (transition towards equal per capita allocation; with reference years 2020 and 2050)
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/PerCapitaConvergence/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/PerCapitaConvergence/datainput.gms"
diff --git a/modules/41_emicapregi/exog/realization.gms b/modules/41_emicapregi/exog/realization.gms
index 4cefd2b425207a9e704a566bd940c87006f16f72..1e16d3e8d8154bb05f46aa9c90411194b9e9ef40 100644
--- a/modules/41_emicapregi/exog/realization.gms
+++ b/modules/41_emicapregi/exog/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/exog.gms
 
+*' @description
+*' Emission caps/permits are allocated from an exogenous emission path that have 
+*' to be provided "manually"
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "declarations" $include "./modules/41_emicapregi/exog/declarations.gms"
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/exog/datainput.gms"
diff --git a/modules/41_emicapregi/module.gms b/modules/41_emicapregi/module.gms
index bb79ebc4754089cd506f724100fc76d880ff38e4..69196c2f937d4e5c6c9b179feefaabfba14a2b7a 100644
--- a/modules/41_emicapregi/module.gms
+++ b/modules/41_emicapregi/module.gms
@@ -6,6 +6,20 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/41_emicapregi.gms
 
+*' @title Regional Emission Caps 
+*'
+*' @description
+*' This module computes reginal emission caps both in absolute terms and as share of global emissions.
+*' In a setting with emissions trading these caps represent allocated permits and permit shares, respectively.
+*' The allocation of caps and permits is based on different burden sharing rules.
+
+*' @limitations
+*' Permit allocation and emissions trading yield less robust results under Nash (decentralized optimization)
+*' compared to Negishi (Social planner optimization).
+
+*' @authors Marian Leimbach, Christoph Bertram
+
+
 *###################### R SECTION START (MODULETYPES) ##########################
 $Ifi "%emicapregi%" == "AbilityToPay" $include "./modules/41_emicapregi/AbilityToPay/realization.gms"
 $Ifi "%emicapregi%" == "CandC" $include "./modules/41_emicapregi/CandC/realization.gms"
diff --git a/modules/41_emicapregi/none/realization.gms b/modules/41_emicapregi/none/realization.gms
index 519c95e74a85e9958bf4e32abc84e2905110047f..af18d75d711c2811db271c588f5d0a4298daa47a 100644
--- a/modules/41_emicapregi/none/realization.gms
+++ b/modules/41_emicapregi/none/realization.gms
@@ -6,6 +6,10 @@
 *** |  Contact: remind@pik-potsdam.de
 *** SOF ./modules/41_emicapregi/none.gms
 
+*' @description
+*' No allocation of regional emission caps/permits - applies to tax scenarios and 
+*' no-tax scenarios without permit trading  
+
 *####################### R SECTION START (PHASES) ##############################
 $Ifi "%phase%" == "datainput" $include "./modules/41_emicapregi/none/datainput.gms"
 $Ifi "%phase%" == "bounds" $include "./modules/41_emicapregi/none/bounds.gms"
diff --git a/scripts/output/comparison/compareScenarios.R b/scripts/output/comparison/compareScenarios.R
index fd1d73f284c33b62f5a0472d6e014124ae277471..b94620870fad0274945ed2df4c64a510934d6921 100644
--- a/scripts/output/comparison/compareScenarios.R
+++ b/scripts/output/comparison/compareScenarios.R
@@ -57,7 +57,7 @@ start_comp <- function(outputdirs,shortTerm,outfilename) {
     outfilename    <- jobname
     tmp.env <- new.env()
     script <- "scripts/run_submit/run_compareScenarios.R"
-    tmp.erorr <- try(sys.source(script,envir=tmp.env))
+    tmp.error <- try(sys.source(script,envir=tmp.env))
     if(!is.null(tmp.error)) warning("Script ",script," was stopped by an error and not executed properly!")
     rm(tmp.env)
   }