diff --git a/docs/make.jl b/docs/make.jl
index d4f25fb2..0d31334c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,8 +1,8 @@
 using Documenter
 
 makedocs(
-	sitename = "MimiPAGE2009.jl",
-	pages = [
+	sitename="MimiPAGE2009.jl",
+	pages=[
 		"Home" => "index.md",
 		"Getting started" => "gettingstarted.md",
 		"Model Structure" => "model-structure.md",
@@ -11,5 +11,5 @@ makedocs(
 )
 
 deploydocs(
-    repo = "github.com/anthofflab/MimiPAGE2009.jl.git"
+    repo="github.com/anthofflab/MimiPAGE2009.jl.git"
 )
diff --git a/src/MimiPAGE2009.jl b/src/MimiPAGE2009.jl
index 13640bc9..47e2439c 100644
--- a/src/MimiPAGE2009.jl
+++ b/src/MimiPAGE2009.jl
@@ -45,7 +45,7 @@ include("components/EquityWeighting.jl")
 
 function buildpage(m::Model, policy::String="policy-a")
 
-    #add all the components
+    # add all the components
     add_comp!(m, co2emissions)
     add_comp!(m, co2cycle)
     add_comp!(m, co2forcing)
@@ -63,11 +63,11 @@ function buildpage(m::Model, policy::String="policy-a")
     add_comp!(m, ClimateTemperature)
     add_comp!(m, SeaLevelRise)
 
-    #Socio-Economics
+    # Socio-Economics
     addpopulation(m)
     add_comp!(m, GDP)
 
-    #Abatement Costs 
+    # Abatement Costs 
 
     addabatementcostparameters(m, :CO2, policy)
     addabatementcostparameters(m, :CH4, policy)
@@ -91,7 +91,7 @@ function buildpage(m::Model, policy::String="policy-a")
     addabatementcosts(m, :Lin, policy)
     add_comp!(m, TotalAbatementCosts)
 
-    #Adaptation Costs
+    # Adaptation Costs
     addadaptationcosts_sealevel(m)
     addadaptationcosts_economic(m)
     addadaptationcosts_noneconomic(m)
@@ -108,10 +108,10 @@ function buildpage(m::Model, policy::String="policy-a")
     # Total costs component 
     add_comp!(m, TotalCosts)
 
-    #Equity weighting and Total Costs
+    # Equity weighting and Total Costs
     add_comp!(m, EquityWeighting)
 
-    #connect parameters together
+    # connect parameters together
     connect_param!(m, :co2cycle => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
     connect_param!(m, :co2cycle => :rt_g0_baseglobaltemp, :ClimateTemperature => :rt_g0_baseglobaltemp)
     connect_param!(m, :co2cycle => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)
@@ -202,7 +202,7 @@ function buildpage(m::Model, policy::String="policy-a")
     connect_param!(m, :NonMarketDamages => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
     connect_param!(m, :NonMarketDamages => :rgdp_per_cap_MarketRemainGDP, :MarketDamages => :rgdp_per_cap_MarketRemainGDP)
     connect_param!(m, :NonMarketDamages => :rcons_per_cap_MarketRemainConsumption, :MarketDamages => :rcons_per_cap_MarketRemainConsumption)
-    connect_param!(m, :NonMarketDamages =>:atl_adjustedtolerableleveloftemprise, :AdaptiveCostsNonEconomic =>:atl_adjustedtolerablelevel, ignoreunits=true)
+    connect_param!(m, :NonMarketDamages => :atl_adjustedtolerableleveloftemprise, :AdaptiveCostsNonEconomic => :atl_adjustedtolerablelevel, ignoreunits=true)
     connect_param!(m, :NonMarketDamages => :imp_actualreduction, :AdaptiveCostsNonEconomic => :imp_adaptedimpacts)
     connect_param!(m, :NonMarketDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)
 
diff --git a/src/climate_model.jl b/src/climate_model.jl
index 9e347448..672a955a 100644
--- a/src/climate_model.jl
+++ b/src/climate_model.jl
@@ -19,7 +19,7 @@ m = Model()
 set_dimension!(m, :time, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
 set_dimension!(m, :region, ["EU", "USA", "OECD","USSR","China","SEAsia","Africa","LatAmerica"])
 
-#add all the components
+# add all the components
 add_comp!(m,co2emissions)
 add_comp!(m, co2cycle)
 add_comp!(m, co2forcing)
@@ -36,7 +36,7 @@ add_comp!(m, SulphateForcing)
 add_comp!(m, TotalForcing)
 add_comp!(m, ClimateTemperature)
 
-#connect parameters together
+# connect parameters together
 
 set_param!(m, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.])
 set_param!(m, :y_year_0, 2008.)
diff --git a/src/components/AbatementCostParameters.jl b/src/components/AbatementCostParameters.jl
index dc100a69..5e168cbd 100644
--- a/src/components/AbatementCostParameters.jl
+++ b/src/components/AbatementCostParameters.jl
@@ -2,26 +2,26 @@
 
     region = Index()
     y_year = Parameter(index=[time], unit="year")
-    y_year_0 = Parameter(unit="year", default = 2008.)
+    y_year_0 = Parameter(unit="year", default=2008.)
 
-    #gas inputs
+    # gas inputs
     emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = Parameter(unit="%")
     q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = Parameter(unit="% of BAU emissions")
     c0init_MostNegativeCostCutbackinBaseYear = Parameter(unit="\$/ton")
     qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = Parameter(unit="% of BAU emissions")
     cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = Parameter(unit="\$/ton")
-    ies_InitialExperienceStockofCutbacks = Parameter(unit= "Million ton")
-    e0_baselineemissions = Parameter(index=[region], unit= "Mtonne/year")
+    ies_InitialExperienceStockofCutbacks = Parameter(unit="Million ton")
+    e0_baselineemissions = Parameter(index=[region], unit="Mtonne/year")
 
-    #regional inputs
-    emitf_uncertaintyinBAUemissfactor = Parameter(index=[region], unit= "none")
+    # regional inputs
+    emitf_uncertaintyinBAUemissfactor = Parameter(index=[region], unit="none")
     q0f_negativecostpercentagefactor = Parameter(index=[region], unit="none")
     cmaxf_maxcostfactor = Parameter(index=[region], unit="none")
 
-    bau_businessasusualemissions = Parameter(index=[time, region], unit = "%")
+    bau_businessasusualemissions = Parameter(index=[time, region], unit="%")
     yagg = Parameter(index=[time], unit="year") # from equity weighting
 
-    #inputs with single, uncertain values
+    # inputs with single, uncertain values
     q0propmult_cutbacksatnegativecostinfinalyear = Parameter(unit="none", default=.733333333333333334)
     qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = Parameter(unit="none", default=1.2666666666666666)
     c0mult_mostnegativecostinfinalyear = Parameter(unit="none", default=.8333333333333334)
@@ -33,34 +33,34 @@
     equity_prop_equityweightsproportion = Parameter(unit="none", default=1.)
 
     # Inputs from other components
-    cbe_absoluteemissionreductions = Parameter(index=[time, region], unit="Mtonne") #TODO: default here?
+    cbe_absoluteemissionreductions = Parameter(index=[time, region], unit="Mtonne") # TODO: default here?
     
-    #Variables
-    emit_UncertaintyinBAUEmissFactor = Variable(index=[region], unit = "%")
-    q0propinit_CutbacksinNegativeCostinBaseYear = Variable(index=[region], unit= "% of BAU emissions")
-    cmaxinit_MaxCutbackCostinBaseYear = Variable(index=[region], unit = "\$/ton")
-    zc_zerocostemissions = Variable(index=[time, region], unit= "%")
+    # Variables
+    emit_UncertaintyinBAUEmissFactor = Variable(index=[region], unit="%")
+    q0propinit_CutbacksinNegativeCostinBaseYear = Variable(index=[region], unit="% of BAU emissions")
+    cmaxinit_MaxCutbackCostinBaseYear = Variable(index=[region], unit="\$/ton")
+    zc_zerocostemissions = Variable(index=[time, region], unit="%")
     cumcbe_cumulativereductionssincebaseyear = Variable(index=[time, region], unit="Mtonne")
     cumcbe_g_totalreductions = Variable(index=[time], unit="Mtonne")
-    learnfac_learning= Variable(index=[time, region], unit= "none")
+    learnfac_learning = Variable(index=[time, region], unit="none")
     auto = Variable(unit="% per year")
-    autofac = Variable(index=[time], unit= "% per year")
-    c0g = Variable(unit= "% per year")
-    c0 = Variable(index=[time], unit= "\$/ton")
-    qmaxminusq0propg = Variable(unit= "% per year")
-    qmaxminusq0prop = Variable(unit = "% of BAU emissions")
-    q0propg = Variable(unit = "% per year")
+    autofac = Variable(index=[time], unit="% per year")
+    c0g = Variable(unit="% per year")
+    c0 = Variable(index=[time], unit="\$/ton")
+    qmaxminusq0propg = Variable(unit="% per year")
+    qmaxminusq0prop = Variable(unit="% of BAU emissions")
+    q0propg = Variable(unit="% per year")
     q0prop = Variable(index=[time, region], unit="% of BAU emissions")
-    q0_absolutecutbacksatnegativecost = Variable(index=[time, region], unit= "Mtonne")
+    q0_absolutecutbacksatnegativecost = Variable(index=[time, region], unit="Mtonne")
     qmax_maxreferencereductions = Variable(index=[time, region], unit="Mtonne")
-    cmax = Variable(index=[time,region], unit = "\$/tonne")
-    blo = Variable(index=[time, region], unit = "per Mtonne")
-    alo = Variable(index=[time, region], unit = "\$/tonne")
-    bhi = Variable(index=[time, region], unit = "per Mtonne")
-    ahi = Variable(index=[time, region], unit = "\$/tonne")
-    mc_marginalcost = Variable(index=[time, region], unit = "\$/tonne")
-    tcq0 = Variable(index=[time, region], unit = "\$million")
-    tc_totalcost = Variable(index=[time, region], unit= "\$million")
+    cmax = Variable(index=[time,region], unit="\$/tonne")
+    blo = Variable(index=[time, region], unit="per Mtonne")
+    alo = Variable(index=[time, region], unit="\$/tonne")
+    bhi = Variable(index=[time, region], unit="per Mtonne")
+    ahi = Variable(index=[time, region], unit="\$/tonne")
+    mc_marginalcost = Variable(index=[time, region], unit="\$/tonne")
+    tcq0 = Variable(index=[time, region], unit="\$million")
+    tc_totalcost = Variable(index=[time, region], unit="\$million")
 
     function run_timestep(p, v, d, t)
 
@@ -72,43 +72,43 @@
             v.cmaxinit_MaxCutbackCostinBaseYear[r] = p.cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear *
                 p.cmaxf_maxcostfactor[r]
 
-            v.zc_zerocostemissions[t,r] = (1+v.emit_UncertaintyinBAUEmissFactor[r]/100 * (p.y_year[t]-p.y_year_0)/(p.y_year[end]-p.y_year_0)) * p.bau_businessasusualemissions[t,r]
+            v.zc_zerocostemissions[t,r] = (1 + v.emit_UncertaintyinBAUEmissFactor[r] / 100 * (p.y_year[t] - p.y_year_0) / (p.y_year[end] - p.y_year_0)) * p.bau_businessasusualemissions[t,r]
 
             if is_first(t)
                 v.cumcbe_cumulativereductionssincebaseyear[t,r] = 0.
             else
-                v.cumcbe_cumulativereductionssincebaseyear[t,r] = v.cumcbe_cumulativereductionssincebaseyear[t-1, r] + p.cbe_absoluteemissionreductions[t-1, r] * p.yagg[t-1]
+                v.cumcbe_cumulativereductionssincebaseyear[t,r] = v.cumcbe_cumulativereductionssincebaseyear[t - 1, r] + p.cbe_absoluteemissionreductions[t - 1, r] * p.yagg[t - 1]
             end
         end
         
         v.cumcbe_g_totalreductions[t] = sum(v.cumcbe_cumulativereductionssincebaseyear[t,:])
 
-        v.auto = (1-p.automult_autonomoustechchange^(1/(p.y_year[end]-p.y_year_0)))*100
-        v.autofac[t] = (1-v.auto/100)^(p.y_year[t] - p.y_year_0)
+        v.auto = (1 - p.automult_autonomoustechchange^(1 / (p.y_year[end] - p.y_year_0))) * 100
+        v.autofac[t] = (1 - v.auto / 100)^(p.y_year[t] - p.y_year_0)
 
-        v.c0g = (p.c0mult_mostnegativecostinfinalyear^(1/(p.y_year[end]-p.y_year_0))-1)*100
-        v.c0[t] = p.c0init_MostNegativeCostCutbackinBaseYear* (1+v.c0g/100)^(p.y_year[t]-p.y_year_0)
+        v.c0g = (p.c0mult_mostnegativecostinfinalyear^(1 / (p.y_year[end] - p.y_year_0)) - 1) * 100
+        v.c0[t] = p.c0init_MostNegativeCostCutbackinBaseYear * (1 + v.c0g / 100)^(p.y_year[t] - p.y_year_0)
 
-        v.qmaxminusq0propg = (p.qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear ^(1/(p.y_year[end]-p.y_year_0))- 1)* 100
-        v.qmaxminusq0prop = p.qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear * (1+ v.qmaxminusq0propg/100)^(p.y_year[t]-p.y_year_0)
+        v.qmaxminusq0propg = (p.qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear^(1 / (p.y_year[end] - p.y_year_0)) - 1) * 100
+        v.qmaxminusq0prop = p.qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear * (1 + v.qmaxminusq0propg / 100)^(p.y_year[t] - p.y_year_0)
 
-        v.q0propg = (p.q0propmult_cutbacksatnegativecostinfinalyear^(1/(p.y_year[end]-p.y_year_0))-1)*100
+        v.q0propg = (p.q0propmult_cutbacksatnegativecostinfinalyear^(1 / (p.y_year[end] - p.y_year_0)) - 1) * 100
 
         for r in d.region
-            v.learnfac_learning[t,r] = ((p.cross_experiencecrossoverratio *v.cumcbe_g_totalreductions[t]+ (1-p.cross_experiencecrossoverratio)*v.cumcbe_cumulativereductionssincebaseyear[t,r] + p.ies_InitialExperienceStockofCutbacks)/ p.ies_InitialExperienceStockofCutbacks)^ -(log(1/(1-p.learn_learningrate))/log(2))
+            v.learnfac_learning[t,r] = ((p.cross_experiencecrossoverratio * v.cumcbe_g_totalreductions[t] + (1 - p.cross_experiencecrossoverratio) * v.cumcbe_cumulativereductionssincebaseyear[t,r] + p.ies_InitialExperienceStockofCutbacks) / p.ies_InitialExperienceStockofCutbacks)^-(log(1 / (1 - p.learn_learningrate)) / log(2))
 
-            v.q0prop[t,r] = v.q0propinit_CutbacksinNegativeCostinBaseYear[r]* (1+v.q0propg/100)^(p.y_year[t]-p.y_year_0)
+            v.q0prop[t,r] = v.q0propinit_CutbacksinNegativeCostinBaseYear[r] * (1 + v.q0propg / 100)^(p.y_year[t] - p.y_year_0)
 
-            v.q0_absolutecutbacksatnegativecost[t,r]= (v.q0prop[t,r]/100)* (v.zc_zerocostemissions[t,r]/100) * p.e0_baselineemissions[r]
+            v.q0_absolutecutbacksatnegativecost[t,r] = (v.q0prop[t,r] / 100) * (v.zc_zerocostemissions[t,r] / 100) * p.e0_baselineemissions[r]
 
-            v.qmax_maxreferencereductions[t,r] = (v.qmaxminusq0prop/100) * (v.zc_zerocostemissions[t,r]/100)* p.e0_baselineemissions[r] + v.q0_absolutecutbacksatnegativecost[t,r]
+            v.qmax_maxreferencereductions[t,r] = (v.qmaxminusq0prop / 100) * (v.zc_zerocostemissions[t,r] / 100) * p.e0_baselineemissions[r] + v.q0_absolutecutbacksatnegativecost[t,r]
 
-            v.cmax[t,r] = v.cmaxinit_MaxCutbackCostinBaseYear[r] * v.learnfac_learning[t,r]* v.autofac[t]
+            v.cmax[t,r] = v.cmaxinit_MaxCutbackCostinBaseYear[r] * v.learnfac_learning[t,r] * v.autofac[t]
 
-            v.blo[t,r] = -2*log((1+p.curve_below_curvatureofMACcurvebelowzerocost)/(1-p.curve_below_curvatureofMACcurvebelowzerocost))/ v.q0_absolutecutbacksatnegativecost[t,r]
-            v.alo[t,r] = v.c0[t]/(exp(-v.blo[t,r]*v.q0_absolutecutbacksatnegativecost[t,r])-1)
-            v.bhi[t,r] = 2*log((1+p.curve_above_curvatureofMACcurveabovezerocost)/(1-p.curve_above_curvatureofMACcurveabovezerocost))/ (v.qmax_maxreferencereductions[t,r] - v.q0_absolutecutbacksatnegativecost[t,r])
-            v.ahi[t,r] = v.cmax[t,r]/ (exp(v.bhi[t,r]*(v.qmax_maxreferencereductions[t,r]-v.q0_absolutecutbacksatnegativecost[t,r]))-1)
+            v.blo[t,r] = -2 * log((1 + p.curve_below_curvatureofMACcurvebelowzerocost) / (1 - p.curve_below_curvatureofMACcurvebelowzerocost)) / v.q0_absolutecutbacksatnegativecost[t,r]
+            v.alo[t,r] = v.c0[t] / (exp(-v.blo[t,r] * v.q0_absolutecutbacksatnegativecost[t,r]) - 1)
+            v.bhi[t,r] = 2 * log((1 + p.curve_above_curvatureofMACcurveabovezerocost) / (1 - p.curve_above_curvatureofMACcurveabovezerocost)) / (v.qmax_maxreferencereductions[t,r] - v.q0_absolutecutbacksatnegativecost[t,r])
+            v.ahi[t,r] = v.cmax[t,r] / (exp(v.bhi[t,r] * (v.qmax_maxreferencereductions[t,r] - v.q0_absolutecutbacksatnegativecost[t,r])) - 1)
         end
     end
 end
@@ -151,8 +151,8 @@ function addabatementcostparameters(model::Model, class::Symbol, policy::String=
         setdistinctparameter(model, componentname, :qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear, 70.)
         setdistinctparameter(model, componentname, :cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear, 333.333333333333)
         setdistinctparameter(model, componentname, :ies_InitialExperienceStockofCutbacks, 2000.)
-        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model,"data/e0_baselineLGemissions.csv"))
-        setdistinctparameter(model, componentname, :bau_businessasusualemissions, readpagedata(model,"data/bau_linemissions.csv"))
+        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model, "data/e0_baselineLGemissions.csv"))
+        setdistinctparameter(model, componentname, :bau_businessasusualemissions, readpagedata(model, "data/bau_linemissions.csv"))
     else
         error("Unknown class of abatement costs.")
     end
diff --git a/src/components/AbatementCosts.jl b/src/components/AbatementCosts.jl
index 87d5ed44..49d4c18c 100644
--- a/src/components/AbatementCosts.jl
+++ b/src/components/AbatementCosts.jl
@@ -2,48 +2,48 @@
     region = Index()
 
    # From the AbatementCostParameters
-    zc_zerocostemissions = Parameter(index=[time, region], unit= "%")
-    q0_absolutecutbacksatnegativecost = Parameter(index=[time, region], unit= "Mtonne")
-    blo = Parameter(index=[time, region], unit = "per Mtonne")
-    alo = Parameter(index=[time, region], unit = "\$/tonne")
-    bhi = Parameter(index=[time, region], unit = "per Mtonne")
-    ahi = Parameter(index=[time, region], unit = "\$/tonne")
+    zc_zerocostemissions = Parameter(index=[time, region], unit="%")
+    q0_absolutecutbacksatnegativecost = Parameter(index=[time, region], unit="Mtonne")
+    blo = Parameter(index=[time, region], unit="per Mtonne")
+    alo = Parameter(index=[time, region], unit="\$/tonne")
+    bhi = Parameter(index=[time, region], unit="per Mtonne")
+    ahi = Parameter(index=[time, region], unit="\$/tonne")
 
     # Driver of abatement costs
-    e0_baselineemissions = Parameter(index=[region], unit= "Mtonne/year")
-    er_emissionsgrowth = Parameter(index=[time, region], unit= "%")
+    e0_baselineemissions = Parameter(index=[region], unit="Mtonne/year")
+    er_emissionsgrowth = Parameter(index=[time, region], unit="%")
 
     # Intermediate outputs
-    cb_reductionsfromzerocostemissions = Variable(index=[time, region], unit= "%")
+    cb_reductionsfromzerocostemissions = Variable(index=[time, region], unit="%")
     cbe_absoluteemissionreductions = Variable(index=[time, region], unit="Mtonne") # Goes to AbatementCostParameters
 
     # Main costs results
-    mc_marginalcost = Variable(index=[time, region], unit = "\$/tonne")
-    tcq0 = Variable(index=[time, region], unit = "\$million")
-    tc_totalcost = Variable(index=[time, region], unit= "\$million")
+    mc_marginalcost = Variable(index=[time, region], unit="\$/tonne")
+    tcq0 = Variable(index=[time, region], unit="\$million")
+    tc_totalcost = Variable(index=[time, region], unit="\$million")
 
     function run_timestep(p, v, d, t)
 
         for r in d.region
             v.cb_reductionsfromzerocostemissions[t,r] = max(p.zc_zerocostemissions[t,r] - p.er_emissionsgrowth[t,r], 0)
-            v.cbe_absoluteemissionreductions[t,r] = v.cb_reductionsfromzerocostemissions[t,r]* p.e0_baselineemissions[r]/100
+            v.cbe_absoluteemissionreductions[t,r] = v.cb_reductionsfromzerocostemissions[t,r] * p.e0_baselineemissions[r] / 100
 
-            if v.cbe_absoluteemissionreductions[t,r]< p.q0_absolutecutbacksatnegativecost[t,r]
-                v.mc_marginalcost[t,r] = p.alo[t,r]* (exp(p.blo[t,r]*(v.cbe_absoluteemissionreductions[t,r]- p.q0_absolutecutbacksatnegativecost[t,r]))-1)
+            if v.cbe_absoluteemissionreductions[t,r] < p.q0_absolutecutbacksatnegativecost[t,r]
+                v.mc_marginalcost[t,r] = p.alo[t,r] * (exp(p.blo[t,r] * (v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r])) - 1)
             else
-                v.mc_marginalcost[t,r] = p.ahi[t,r]*(exp(p.bhi[t,r]*(v.cbe_absoluteemissionreductions[t,r]- p.q0_absolutecutbacksatnegativecost[t,r]))-1)
+                v.mc_marginalcost[t,r] = p.ahi[t,r] * (exp(p.bhi[t,r] * (v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r])) - 1)
             end
 
             if p.q0_absolutecutbacksatnegativecost[t,r] == 0.
                 v.tcq0[t,r] = 0.
             else
-                v.tcq0[t,r] = (p.alo[t,r]/p.blo[t,r])*(1-exp(-p.blo[t,r]* p.q0_absolutecutbacksatnegativecost[t,r]))- p.alo[t,r]*p.q0_absolutecutbacksatnegativecost[t,r]
+                v.tcq0[t,r] = (p.alo[t,r] / p.blo[t,r]) * (1 - exp(-p.blo[t,r] * p.q0_absolutecutbacksatnegativecost[t,r])) - p.alo[t,r] * p.q0_absolutecutbacksatnegativecost[t,r]
             end
 
-            if v.cbe_absoluteemissionreductions[t,r]<p.q0_absolutecutbacksatnegativecost[t,r]
-                v.tc_totalcost[t,r] = (p.alo[t,r]/p.blo[t,r])*(exp(p.blo[t,r]*(v.cbe_absoluteemissionreductions[t,r]- p.q0_absolutecutbacksatnegativecost[t,r]))- exp(-p.blo[t,r]*p.q0_absolutecutbacksatnegativecost[t,r])) - p.alo[t,r]*v.cbe_absoluteemissionreductions[t,r]
+            if v.cbe_absoluteemissionreductions[t,r] < p.q0_absolutecutbacksatnegativecost[t,r]
+                v.tc_totalcost[t,r] = (p.alo[t,r] / p.blo[t,r]) * (exp(p.blo[t,r] * (v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r])) - exp(-p.blo[t,r] * p.q0_absolutecutbacksatnegativecost[t,r])) - p.alo[t,r] * v.cbe_absoluteemissionreductions[t,r]
             else
-                v.tc_totalcost[t,r] = (p.ahi[t,r]/p.bhi[t,r])* (exp(p.bhi[t,r]*(v.cbe_absoluteemissionreductions[t,r]-p.q0_absolutecutbacksatnegativecost[t,r]))-1) - p.ahi[t,r]*(v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r]) + v.tcq0[t,r]
+                v.tc_totalcost[t,r] = (p.ahi[t,r] / p.bhi[t,r]) * (exp(p.bhi[t,r] * (v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r])) - 1) - p.ahi[t,r] * (v.cbe_absoluteemissionreductions[t,r] - p.q0_absolutecutbacksatnegativecost[t,r]) + v.tcq0[t,r]
             end
         end 
     end       
@@ -75,11 +75,11 @@ function addabatementcosts(model::Model, class::Symbol, policy::String="policy-a
             setdistinctparameter(model, componentname, :er_emissionsgrowth, readpagedata(model, "data/$policy/er_N2Oemissionsgrowth.csv"))
         end
     elseif class == :Lin
-        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model,"data/e0_baselineLGemissions.csv"))
+        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model, "data/e0_baselineLGemissions.csv"))
         if policy == "policy-a"
-            setdistinctparameter(model, componentname, :er_emissionsgrowth, readpagedata(model,"data/er_LGemissionsgrowth.csv"))
+            setdistinctparameter(model, componentname, :er_emissionsgrowth, readpagedata(model, "data/er_LGemissionsgrowth.csv"))
         else
-            setdistinctparameter(model, componentname, :er_emissionsgrowth, readpagedata(model,"data/$policy/er_LGemissionsgrowth.csv"))
+            setdistinctparameter(model, componentname, :er_emissionsgrowth, readpagedata(model, "data/$policy/er_LGemissionsgrowth.csv"))
         end
     else
         error("Unknown class of abatement costs.")
diff --git a/src/components/AdaptationCosts.jl b/src/components/AdaptationCosts.jl
index e499049c..f2d75263 100644
--- a/src/components/AdaptationCosts.jl
+++ b/src/components/AdaptationCosts.jl
@@ -10,13 +10,13 @@
 
     automult_autonomouschange = Parameter(unit="none")
     impmax_maximumadaptivecapacity = Parameter(index=[region], unit="driver")
-    #tolerability parameters
+    # tolerability parameters
     plateau_increaseintolerableplateaufromadaptation = Parameter(index=[region], unit="driver")
     pstart_startdateofadaptpolicy = Parameter(index=[region], unit="year")
     pyears_yearstilfulleffect = Parameter(index=[region], unit="year")
-    impred_eventualpercentreduction = Parameter(index=[region], unit= "%")
-    istart_startdate = Parameter(index=[region], unit = "year")
-    iyears_yearstilfulleffect = Parameter(index=[region], unit= "year")
+    impred_eventualpercentreduction = Parameter(index=[region], unit="%")
+    istart_startdate = Parameter(index=[region], unit="year")
+    iyears_yearstilfulleffect = Parameter(index=[region], unit="year")
 
     cp_costplateau_eu = Parameter(unit="%GDP/driver")
     ci_costimpact_eu = Parameter(unit="%GDP/%driver")
@@ -34,26 +34,26 @@
     function run_timestep(p, v, d, tt)
 
         # Hope (2009), p. 21, equation -5
-        auto_autonomouschangepercent = (1 - p.automult_autonomouschange^(1/(p.y_year[end] - p.y_year_0)))*100 # % per year
-        v.autofac_autonomouschangefraction[tt] = (1 - auto_autonomouschangepercent/100)^(p.y_year[tt] - p.y_year_0) # Varies by year
+        auto_autonomouschangepercent = (1 - p.automult_autonomouschange^(1 / (p.y_year[end] - p.y_year_0))) * 100 # % per year
+        v.autofac_autonomouschangefraction[tt] = (1 - auto_autonomouschangepercent / 100)^(p.y_year[tt] - p.y_year_0) # Varies by year
 
         for rr in d.region
-            #calculate adjusted tolerable level and max impact based on adaptation policy
+            # calculate adjusted tolerable level and max impact based on adaptation policy
             if (p.y_year[tt] - p.pstart_startdateofadaptpolicy[rr]) < 0
-                v.atl_adjustedtolerablelevel[tt,rr]= 0
-            elseif ((p.y_year[tt]-p.pstart_startdateofadaptpolicy[rr])/p.pyears_yearstilfulleffect[rr])<1.
-                v.atl_adjustedtolerablelevel[tt,rr]=
-                    ((p.y_year[tt]-p.pstart_startdateofadaptpolicy[rr])/p.pyears_yearstilfulleffect[rr]) *
+                v.atl_adjustedtolerablelevel[tt,rr] = 0
+            elseif ((p.y_year[tt] - p.pstart_startdateofadaptpolicy[rr]) / p.pyears_yearstilfulleffect[rr]) < 1.
+                v.atl_adjustedtolerablelevel[tt,rr] =
+                    ((p.y_year[tt] - p.pstart_startdateofadaptpolicy[rr]) / p.pyears_yearstilfulleffect[rr]) *
                     p.plateau_increaseintolerableplateaufromadaptation[rr]
             else
                 v.atl_adjustedtolerablelevel[tt,rr] = p.plateau_increaseintolerableplateaufromadaptation[rr]
             end
 
-            if (p.y_year[tt]- p.istart_startdate[rr]) < 0
+            if (p.y_year[tt] - p.istart_startdate[rr]) < 0
                 v.imp_adaptedimpacts[tt,rr] = 0
-            elseif ((p.y_year[tt]-p.istart_startdate[rr])/p.iyears_yearstilfulleffect[rr]) < 1
+            elseif ((p.y_year[tt] - p.istart_startdate[rr]) / p.iyears_yearstilfulleffect[rr]) < 1
                 v.imp_adaptedimpacts[tt,rr] =
-                    (p.y_year[tt]-p.istart_startdate[rr])/p.iyears_yearstilfulleffect[rr]*
+                    (p.y_year[tt] - p.istart_startdate[rr]) / p.iyears_yearstilfulleffect[rr] *
                     p.impred_eventualpercentreduction[rr]
             else
                 v.imp_adaptedimpacts[tt,rr] = p.impred_eventualpercentreduction[rr]
diff --git a/src/components/CH4cycle.jl b/src/components/CH4cycle.jl
index 55bd6a20..42c005dd 100644
--- a/src/components/CH4cycle.jl
+++ b/src/components/CH4cycle.jl
@@ -1,54 +1,54 @@
 @defcomp ch4cycle begin
-    e_globalCH4emissions=Parameter(index=[time],unit="Mtonne/year")
-    e_0globalCH4emissions=Parameter(unit="Mtonne/year", default=363.00000000000006)
-    c_CH4concentration=Variable(index=[time],unit="ppbv")
-    pic_preindustconcCH4=Parameter(unit="ppbv", default=700.)
-    exc_excessconcCH4=Variable(unit="ppbv")
-    c0_CH4concbaseyr=Parameter(unit="ppbv", default=1860.)
-    re_remainCH4=Variable(index=[time],unit="ppbv")
-    re_remainCH4base=Variable(unit="ppbv")
-    nte_natCH4emissions=Variable(index=[time],unit="Mtonne/year")
-    air_CH4fractioninatm=Parameter(unit="%", default=100.)
-    tea_CH4emissionstoatm=Variable(index=[time],unit="Mtonne/year")
-    teay_CH4emissionstoatm=Variable(index=[time],unit="Mtonne/t")
-    y_year=Parameter(index=[time],unit="year")
-    y_year_0=Parameter(unit="year")
-    res_CH4atmlifetime=Parameter(unit="year", default=10.5)
-    den_CH4density=Parameter(unit="Mtonne/ppbv", default=2.78)
-    stim_CH4emissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
-    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=0.9258270139190647)
-    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")
+    e_globalCH4emissions = Parameter(index=[time], unit="Mtonne/year")
+    e_0globalCH4emissions = Parameter(unit="Mtonne/year", default=363.00000000000006)
+    c_CH4concentration = Variable(index=[time], unit="ppbv")
+    pic_preindustconcCH4 = Parameter(unit="ppbv", default=700.)
+    exc_excessconcCH4 = Variable(unit="ppbv")
+    c0_CH4concbaseyr = Parameter(unit="ppbv", default=1860.)
+    re_remainCH4 = Variable(index=[time], unit="ppbv")
+    re_remainCH4base = Variable(unit="ppbv")
+    nte_natCH4emissions = Variable(index=[time], unit="Mtonne/year")
+    air_CH4fractioninatm = Parameter(unit="%", default=100.)
+    tea_CH4emissionstoatm = Variable(index=[time], unit="Mtonne/year")
+    teay_CH4emissionstoatm = Variable(index=[time], unit="Mtonne/t")
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
+    res_CH4atmlifetime = Parameter(unit="year", default=10.5)
+    den_CH4density = Parameter(unit="Mtonne/ppbv", default=2.78)
+    stim_CH4emissionfeedback = Parameter(unit="Mtonne/degreeC", default=0.)
+    rtl_g0_baselandtemp = Parameter(unit="degreeC", default=0.9258270139190647)
+    rtl_g_landtemperature = Parameter(index=[time], unit="degreeC")
 
     function run_timestep(p, v, d, t)
 
         if is_first(t)
-            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
-            nte_0=p.stim_CH4emissionfeedback*p.rtl_g0_baselandtemp
-            v.nte_natCH4emissions[t]=p.stim_CH4emissionfeedback*p.rtl_g0_baselandtemp
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_CH4emissionstoatm[t]=(p.e_globalCH4emissions[t]+v.nte_natCH4emissions[t])*p.air_CH4fractioninatm/100
-            tea_0=(p.e_0globalCH4emissions+nte_0)*p.air_CH4fractioninatm/100
-            v.teay_CH4emissionstoatm[t]=(v.tea_CH4emissionstoatm[t]+tea_0)*(p.y_year[t]-p.y_year_0)/2
-            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
-            v.exc_excessconcCH4=p.c0_CH4concbaseyr-p.pic_preindustconcCH4
-            #Eq. 2 from Hope (2006) - base-year remaining emissions
-            v.re_remainCH4base=v.exc_excessconcCH4*p.den_CH4density
-            v.re_remainCH4[t]=v.re_remainCH4base*exp(-(p.y_year[t]-p.y_year_0)/p.res_CH4atmlifetime)+
-                v.teay_CH4emissionstoatm[t]*p.res_CH4atmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_CH4atmlifetime))/(p.y_year[t]-p.y_year_0)
+            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
+            nte_0 = p.stim_CH4emissionfeedback * p.rtl_g0_baselandtemp
+            v.nte_natCH4emissions[t] = p.stim_CH4emissionfeedback * p.rtl_g0_baselandtemp
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_CH4emissionstoatm[t] = (p.e_globalCH4emissions[t] + v.nte_natCH4emissions[t]) * p.air_CH4fractioninatm / 100
+            tea_0 = (p.e_0globalCH4emissions + nte_0) * p.air_CH4fractioninatm / 100
+            v.teay_CH4emissionstoatm[t] = (v.tea_CH4emissionstoatm[t] + tea_0) * (p.y_year[t] - p.y_year_0) / 2
+            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
+            v.exc_excessconcCH4 = p.c0_CH4concbaseyr - p.pic_preindustconcCH4
+            # Eq. 2 from Hope (2006) - base-year remaining emissions
+            v.re_remainCH4base = v.exc_excessconcCH4 * p.den_CH4density
+            v.re_remainCH4[t] = v.re_remainCH4base * exp(-(p.y_year[t] - p.y_year_0) / p.res_CH4atmlifetime) +
+                v.teay_CH4emissionstoatm[t] * p.res_CH4atmlifetime * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.res_CH4atmlifetime)) / (p.y_year[t] - p.y_year_0)
         else
-            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
-            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
-            v.nte_natCH4emissions[t]=p.stim_CH4emissionfeedback*p.rtl_g_landtemperature[t-1] #askChrisHope
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_CH4emissionstoatm[t]=(p.e_globalCH4emissions[t]+v.nte_natCH4emissions[t])*p.air_CH4fractioninatm/100
-            #eq.7 from Hope (2006) - average emissions to atm over time period
-            v.teay_CH4emissionstoatm[t]=(v.tea_CH4emissionstoatm[t]+v.tea_CH4emissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2
-            #eq.10 from Hope (2006) - remaining emissions in atmosphere
-            v.re_remainCH4[t]=v.re_remainCH4[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CH4atmlifetime)+
-                v.teay_CH4emissionstoatm[t]*p.res_CH4atmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CH4atmlifetime))/(p.y_year[t]-p.y_year[t-1])
+            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
+            # Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
+            v.nte_natCH4emissions[t] = p.stim_CH4emissionfeedback * p.rtl_g_landtemperature[t - 1] # askChrisHope
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_CH4emissionstoatm[t] = (p.e_globalCH4emissions[t] + v.nte_natCH4emissions[t]) * p.air_CH4fractioninatm / 100
+            # eq.7 from Hope (2006) - average emissions to atm over time period
+            v.teay_CH4emissionstoatm[t] = (v.tea_CH4emissionstoatm[t] + v.tea_CH4emissionstoatm[t - 1]) * (p.y_year[t] - p.y_year[t - 1]) / 2
+            # eq.10 from Hope (2006) - remaining emissions in atmosphere
+            v.re_remainCH4[t] = v.re_remainCH4[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_CH4atmlifetime) +
+                v.teay_CH4emissionstoatm[t] * p.res_CH4atmlifetime * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_CH4atmlifetime)) / (p.y_year[t] - p.y_year[t - 1])
         end
 
-        #eq.11 from Hope(2006) - CH4 concentration
-        v.c_CH4concentration[t]=p.pic_preindustconcCH4+v.exc_excessconcCH4*v.re_remainCH4[t]/v.re_remainCH4base
+        # eq.11 from Hope(2006) - CH4 concentration
+        v.c_CH4concentration[t] = p.pic_preindustconcCH4 + v.exc_excessconcCH4 * v.re_remainCH4[t] / v.re_remainCH4base
     end
 end
diff --git a/src/components/CH4emissions.jl b/src/components/CH4emissions.jl
index d0c44354..10095ab2 100644
--- a/src/components/CH4emissions.jl
+++ b/src/components/CH4emissions.jl
@@ -1,19 +1,19 @@
 @defcomp ch4emissions begin
-    region=Index()
+    region = Index()
 
-    e_globalCH4emissions=Variable(index=[time],unit="Mtonne/year")
-    e0_baselineCH4emissions=Parameter(index=[region],unit="Mtonne/year")
-    e_regionalCH4emissions=Variable(index=[time,region],unit="Mtonne/year")
-    er_CH4emissionsgrowth=Parameter(index=[time,region],unit="%")
+    e_globalCH4emissions = Variable(index=[time], unit="Mtonne/year")
+    e0_baselineCH4emissions = Parameter(index=[region], unit="Mtonne/year")
+    e_regionalCH4emissions = Variable(index=[time,region], unit="Mtonne/year")
+    er_CH4emissionsgrowth = Parameter(index=[time,region], unit="%")
 
     function run_timestep(p, v, d, t)
 
-        #eq.4 in Hope (2006) - regional CH4 emissions as % change from baseline
+        # eq.4 in Hope (2006) - regional CH4 emissions as % change from baseline
         for r in d.region
-            v.e_regionalCH4emissions[t,r]=p.er_CH4emissionsgrowth[t,r]*p.e0_baselineCH4emissions[r]/100
+            v.e_regionalCH4emissions[t,r] = p.er_CH4emissionsgrowth[t,r] * p.e0_baselineCH4emissions[r] / 100
         end
 
-        #eq. 5 in Hope (2006) - global CH4 emissions are sum of regional emissions
-        v.e_globalCH4emissions[t]=sum(v.e_regionalCH4emissions[t,:])
+        # eq. 5 in Hope (2006) - global CH4 emissions are sum of regional emissions
+        v.e_globalCH4emissions[t] = sum(v.e_regionalCH4emissions[t,:])
     end
 end
diff --git a/src/components/CH4forcing.jl b/src/components/CH4forcing.jl
index b047fb91..58a3c279 100644
--- a/src/components/CH4forcing.jl
+++ b/src/components/CH4forcing.jl
@@ -1,23 +1,23 @@
 @defcomp ch4forcing begin
-    c_N2Oconcentration=Parameter(index=[time],unit="ppbv")
-    c_CH4concentration=Parameter(index=[time],unit="ppbv")
-    f0_CH4baseforcing=Parameter(unit="W/m2", default=0.550)
-    fslope_CH4forcingslope=Parameter(unit="W/m2", default=0.036)
-    c0_baseN2Oconc=Parameter(unit="ppbv", default=322.)
-    c0_baseCH4conc=Parameter(unit="ppbv", default=1860.)
-    f_CH4forcing=Variable(index=[time],unit="W/m2")
-    over_baseoverlap=Variable(unit="W/m2")
-    over=Variable(index=[time],unit="W/m2")
+    c_N2Oconcentration = Parameter(index=[time], unit="ppbv")
+    c_CH4concentration = Parameter(index=[time], unit="ppbv")
+    f0_CH4baseforcing = Parameter(unit="W/m2", default=0.550)
+    fslope_CH4forcingslope = Parameter(unit="W/m2", default=0.036)
+    c0_baseN2Oconc = Parameter(unit="ppbv", default=322.)
+    c0_baseCH4conc = Parameter(unit="ppbv", default=1860.)
+    f_CH4forcing = Variable(index=[time], unit="W/m2")
+    over_baseoverlap = Variable(unit="W/m2")
+    over = Variable(index=[time], unit="W/m2")
 
     function run_timestep(p, v, d, t)
 
-        #from p.16 in Hope 2009
+        # from p.16 in Hope 2009
         if is_first(t)
-            #calculate baseline forcing overlap in first time period
-            v.over_baseoverlap=-0.47*log(1+2.0e-5*(p.c0_baseN2Oconc*p.c0_baseCH4conc)^0.75+5.3e-15*p.c0_baseCH4conc*(p.c0_baseCH4conc*p.c0_baseN2Oconc)^1.52)
+            # calculate baseline forcing overlap in first time period
+            v.over_baseoverlap = -0.47 * log(1 + 2.0e-5 * (p.c0_baseN2Oconc * p.c0_baseCH4conc)^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c0_baseN2Oconc)^1.52)
         end
 
-        v.over[t]=-0.47*log(1+2.0e-5*(p.c_CH4concentration[t]*p.c0_baseN2Oconc)^0.75+5.3e-15*p.c_CH4concentration[t]*(p.c0_baseN2Oconc*p.c_CH4concentration[t])^1.52)
-        v.f_CH4forcing[t]=p.f0_CH4baseforcing+p.fslope_CH4forcingslope*(sqrt(p.c_CH4concentration[t])-sqrt(p.c0_baseCH4conc))+v.over[t]-v.over_baseoverlap
+        v.over[t] = -0.47 * log(1 + 2.0e-5 * (p.c_CH4concentration[t] * p.c0_baseN2Oconc)^0.75 + 5.3e-15 * p.c_CH4concentration[t] * (p.c0_baseN2Oconc * p.c_CH4concentration[t])^1.52)
+        v.f_CH4forcing[t] = p.f0_CH4baseforcing + p.fslope_CH4forcingslope * (sqrt(p.c_CH4concentration[t]) - sqrt(p.c0_baseCH4conc)) + v.over[t] - v.over_baseoverlap
     end
 end
diff --git a/src/components/CO2cycle.jl b/src/components/CO2cycle.jl
index 57b94b79..099fd0bc 100644
--- a/src/components/CO2cycle.jl
+++ b/src/components/CO2cycle.jl
@@ -1,77 +1,77 @@
 @defcomp co2cycle begin
-    e_globalCO2emissions=Parameter(index=[time],unit="Mtonne/year")
-    e0_globalCO2emissions=Parameter(unit="Mtonne/year", default=38191.0315797948)
-    c_CO2concentration=Variable(index=[time],unit="ppbv")
-    pic_preindustconcCO2=Parameter(unit="ppbv", default=278000.)
-    exc_excessconcCO2=Variable(unit="ppbv")
-    c0_CO2concbaseyr=Parameter(unit="ppbv", default=395000.)
-    re_remainCO2=Variable(index=[time],unit="Mtonne")
-    re_remainCO2base=Variable(unit="Mtonne")
-    renoccf_remainCO2wocc=Variable(index=[time],unit="Mtonne")
-    air_CO2fractioninatm=Parameter(unit="%", default=62.00)
-    stay_fractionCO2emissionsinatm=Parameter(default=0.3)
-    tea_CO2emissionstoatm=Variable(index=[time],unit="Mtonne/year")
-    teay_CO2emissionstoatm=Variable(index=[time],unit="Mtonne/t")
-    ccf_CO2feedback=Parameter(unit="%/degreeC", default=9.66666666666667)
-    ccfmax_maxCO2feedback=Parameter(unit="%", default=53.3333333333333)
-    cea_cumCO2emissionsatm=Variable(index=[time],unit="Mtonne")
-    ce_0_basecumCO2emissions=Parameter(unit="Mtonne", default=2050000.)
-    y_year=Parameter(index=[time],unit="year")
-    y_year_0=Parameter(unit="year")
-    res_CO2atmlifetime=Parameter(unit="year", default=73.3333333333333)
-    den_CO2density=Parameter(unit="Mtonne/ppbv", default=7.8)
-    rt_g0_baseglobaltemp=Parameter(unit="degreeC", default=0.735309967925382)
-    rt_g_globaltemperature=Parameter(index=[time],unit="degreeC")
+    e_globalCO2emissions = Parameter(index=[time], unit="Mtonne/year")
+    e0_globalCO2emissions = Parameter(unit="Mtonne/year", default=38191.0315797948)
+    c_CO2concentration = Variable(index=[time], unit="ppbv")
+    pic_preindustconcCO2 = Parameter(unit="ppbv", default=278000.)
+    exc_excessconcCO2 = Variable(unit="ppbv")
+    c0_CO2concbaseyr = Parameter(unit="ppbv", default=395000.)
+    re_remainCO2 = Variable(index=[time], unit="Mtonne")
+    re_remainCO2base = Variable(unit="Mtonne")
+    renoccf_remainCO2wocc = Variable(index=[time], unit="Mtonne")
+    air_CO2fractioninatm = Parameter(unit="%", default=62.00)
+    stay_fractionCO2emissionsinatm = Parameter(default=0.3)
+    tea_CO2emissionstoatm = Variable(index=[time], unit="Mtonne/year")
+    teay_CO2emissionstoatm = Variable(index=[time], unit="Mtonne/t")
+    ccf_CO2feedback = Parameter(unit="%/degreeC", default=9.66666666666667)
+    ccfmax_maxCO2feedback = Parameter(unit="%", default=53.3333333333333)
+    cea_cumCO2emissionsatm = Variable(index=[time], unit="Mtonne")
+    ce_0_basecumCO2emissions = Parameter(unit="Mtonne", default=2050000.)
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
+    res_CO2atmlifetime = Parameter(unit="year", default=73.3333333333333)
+    den_CO2density = Parameter(unit="Mtonne/ppbv", default=7.8)
+    rt_g0_baseglobaltemp = Parameter(unit="degreeC", default=0.735309967925382)
+    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
 
     function run_timestep(p, v, d, t)
 
         if is_first(t)
-            #CO2 emissions gain calculated based on PAGE 2009
-            gain=p.ccf_CO2feedback*p.rt_g0_baseglobaltemp
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            tea0=p.e0_globalCO2emissions*p.air_CO2fractioninatm/100
-            v.tea_CO2emissionstoatm[t]=(p.e_globalCO2emissions[t])*p.air_CO2fractioninatm/100
-            v.teay_CO2emissionstoatm[t]=(v.tea_CO2emissionstoatm[t]+tea0)*(p.y_year[t]-p.y_year_0)/2
-            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
-            v.exc_excessconcCO2=p.c0_CO2concbaseyr-p.pic_preindustconcCO2
-            #Eq. 2 from Hope (2006) - base-year remaining emissions
-            v.re_remainCO2base=v.exc_excessconcCO2*p.den_CO2density
-            #PAGE 2009 initial remaining emissions without CO2 feedback
-            renoccf0_remainCO2wocc=v.re_remainCO2base/(1+gain/100)
-            #eq. 8 from Hope (2006) - baseline cumulative emissions to atmosphere
-            ceabase=p.ce_0_basecumCO2emissions*p.air_CO2fractioninatm/100
-            #eq.9 from Hope(2006) - cumulative emissions in atmosphere
-            v.cea_cumCO2emissionsatm[t]=ceabase+v.teay_CO2emissionstoatm[t]
-            #eq.11 from Hope (2006) - anthropogenic remaining emissions
-            v.renoccf_remainCO2wocc[t]=p.stay_fractionCO2emissionsinatm*ceabase*
-                (1-exp(-(p.y_year[t]-p.y_year_0)/
-                p.res_CO2atmlifetime))+renoccf0_remainCO2wocc*
-                exp(-(p.y_year[t]-p.y_year_0)/p.res_CO2atmlifetime)+
-                v.teay_CO2emissionstoatm[t]*exp(-(p.y_year[t]-p.y_year_0)/
-                (2*p.res_CO2atmlifetime))
-            #Hope 2009 - remaining emissions with CO2 feedback
-            v.re_remainCO2[t]=v.renoccf_remainCO2wocc[t]*(1+gain/100)
+            # CO2 emissions gain calculated based on PAGE 2009
+            gain = p.ccf_CO2feedback * p.rt_g0_baseglobaltemp
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            tea0 = p.e0_globalCO2emissions * p.air_CO2fractioninatm / 100
+            v.tea_CO2emissionstoatm[t] = (p.e_globalCO2emissions[t]) * p.air_CO2fractioninatm / 100
+            v.teay_CO2emissionstoatm[t] = (v.tea_CO2emissionstoatm[t] + tea0) * (p.y_year[t] - p.y_year_0) / 2
+            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
+            v.exc_excessconcCO2 = p.c0_CO2concbaseyr - p.pic_preindustconcCO2
+            # Eq. 2 from Hope (2006) - base-year remaining emissions
+            v.re_remainCO2base = v.exc_excessconcCO2 * p.den_CO2density
+            # PAGE 2009 initial remaining emissions without CO2 feedback
+            renoccf0_remainCO2wocc = v.re_remainCO2base / (1 + gain / 100)
+            # eq. 8 from Hope (2006) - baseline cumulative emissions to atmosphere
+            ceabase = p.ce_0_basecumCO2emissions * p.air_CO2fractioninatm / 100
+            # eq.9 from Hope(2006) - cumulative emissions in atmosphere
+            v.cea_cumCO2emissionsatm[t] = ceabase + v.teay_CO2emissionstoatm[t]
+            # eq.11 from Hope (2006) - anthropogenic remaining emissions
+            v.renoccf_remainCO2wocc[t] = p.stay_fractionCO2emissionsinatm * ceabase *
+                (1 - exp(-(p.y_year[t] - p.y_year_0) /
+                p.res_CO2atmlifetime)) + renoccf0_remainCO2wocc *
+                exp(-(p.y_year[t] - p.y_year_0) / p.res_CO2atmlifetime) +
+                v.teay_CO2emissionstoatm[t] * exp(-(p.y_year[t] - p.y_year_0) /
+                (2 * p.res_CO2atmlifetime))
+            # Hope 2009 - remaining emissions with CO2 feedback
+            v.re_remainCO2[t] = v.renoccf_remainCO2wocc[t] * (1 + gain / 100)
         else
-            #CO2 emissions gain calculated based on PAGE 2009
-            gain=min(p.ccf_CO2feedback*p.rt_g_globaltemperature[t-1],p.ccfmax_maxCO2feedback)
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_CO2emissionstoatm[t]=(p.e_globalCO2emissions[t])*p.air_CO2fractioninatm/100
-            #eq.7 from Hope (2006) - total emissions over time period
-            v.teay_CO2emissionstoatm[t]=(v.tea_CO2emissionstoatm[t]+v.tea_CO2emissionstoatm[t-1])*
-                (p.y_year[t]-p.y_year[t-1])/2
-            #eq.9 from Hope(2006) - cumulative emissions in atmosphere
-            v.cea_cumCO2emissionsatm[t]=v.cea_cumCO2emissionsatm[t-1]+v.teay_CO2emissionstoatm[t]
-            #eq.11 from Hope (2006) - anthropogenic remaining emissions
-            v.renoccf_remainCO2wocc[t]=p.stay_fractionCO2emissionsinatm*v.cea_cumCO2emissionsatm[t-1]*
-                (1-exp(-(p.y_year[t]-p.y_year[t-1])/
-                p.res_CO2atmlifetime))+v.renoccf_remainCO2wocc[t-1]*
-                exp(-(p.y_year[t]-p.y_year[t-1])/p.res_CO2atmlifetime)+
-                v.teay_CO2emissionstoatm[t]*exp(-(p.y_year[t]-p.y_year[t-1])/
-                (2*p.res_CO2atmlifetime))
-            #Hope 2009 - remaining emissions with CO2 feedback
-            v.re_remainCO2[t]=v.renoccf_remainCO2wocc[t]*(1+gain/100)
+            # CO2 emissions gain calculated based on PAGE 2009
+            gain = min(p.ccf_CO2feedback * p.rt_g_globaltemperature[t - 1], p.ccfmax_maxCO2feedback)
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_CO2emissionstoatm[t] = (p.e_globalCO2emissions[t]) * p.air_CO2fractioninatm / 100
+            # eq.7 from Hope (2006) - total emissions over time period
+            v.teay_CO2emissionstoatm[t] = (v.tea_CO2emissionstoatm[t] + v.tea_CO2emissionstoatm[t - 1]) *
+                (p.y_year[t] - p.y_year[t - 1]) / 2
+            # eq.9 from Hope(2006) - cumulative emissions in atmosphere
+            v.cea_cumCO2emissionsatm[t] = v.cea_cumCO2emissionsatm[t - 1] + v.teay_CO2emissionstoatm[t]
+            # eq.11 from Hope (2006) - anthropogenic remaining emissions
+            v.renoccf_remainCO2wocc[t] = p.stay_fractionCO2emissionsinatm * v.cea_cumCO2emissionsatm[t - 1] *
+                (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) /
+                p.res_CO2atmlifetime)) + v.renoccf_remainCO2wocc[t - 1] *
+                exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_CO2atmlifetime) +
+                v.teay_CO2emissionstoatm[t] * exp(-(p.y_year[t] - p.y_year[t - 1]) /
+                (2 * p.res_CO2atmlifetime))
+            # Hope 2009 - remaining emissions with CO2 feedback
+            v.re_remainCO2[t] = v.renoccf_remainCO2wocc[t] * (1 + gain / 100)
         end
-        #eq.11 from Hope(2006) - CO2 concentration
-        v.c_CO2concentration[t]=p.pic_preindustconcCO2+v.exc_excessconcCO2 * v.re_remainCO2[t]/v.re_remainCO2base
+        # eq.11 from Hope(2006) - CO2 concentration
+        v.c_CO2concentration[t] = p.pic_preindustconcCO2 + v.exc_excessconcCO2 * v.re_remainCO2[t] / v.re_remainCO2base
     end
 end
diff --git a/src/components/CO2emissions.jl b/src/components/CO2emissions.jl
index 45e6e4e5..002f9552 100644
--- a/src/components/CO2emissions.jl
+++ b/src/components/CO2emissions.jl
@@ -1,18 +1,18 @@
 @defcomp co2emissions begin
-    region=Index()
+    region = Index()
 
-    e_globalCO2emissions=Variable(index=[time],unit="Mtonne/year")
-    e0_baselineCO2emissions=Parameter(index=[region],unit="Mtonne/year")
-    e_regionalCO2emissions=Variable(index=[time,region],unit="Mtonne/year")
-    er_CO2emissionsgrowth=Parameter(index=[time,region],unit="%")
+    e_globalCO2emissions = Variable(index=[time], unit="Mtonne/year")
+    e0_baselineCO2emissions = Parameter(index=[region], unit="Mtonne/year")
+    e_regionalCO2emissions = Variable(index=[time,region], unit="Mtonne/year")
+    er_CO2emissionsgrowth = Parameter(index=[time,region], unit="%")
     
     function run_timestep(p, v, d, t)
 
-        #eq.4 in Hope (2006) - regional CO2 emissions as % change from baseline
+        # eq.4 in Hope (2006) - regional CO2 emissions as % change from baseline
         for r in d.region
-            v.e_regionalCO2emissions[t,r]=p.er_CO2emissionsgrowth[t,r]*p.e0_baselineCO2emissions[r]/100
+            v.e_regionalCO2emissions[t,r] = p.er_CO2emissionsgrowth[t,r] * p.e0_baselineCO2emissions[r] / 100
         end
-        #eq. 5 in Hope (2006) - global CO2 emissions are sum of regional emissions
-        v.e_globalCO2emissions[t]=sum(v.e_regionalCO2emissions[t,:])
+        # eq. 5 in Hope (2006) - global CO2 emissions are sum of regional emissions
+        v.e_globalCO2emissions[t] = sum(v.e_regionalCO2emissions[t,:])
     end
 end
\ No newline at end of file
diff --git a/src/components/CO2forcing.jl b/src/components/CO2forcing.jl
index 3d4480a6..6f3f8789 100644
--- a/src/components/CO2forcing.jl
+++ b/src/components/CO2forcing.jl
@@ -1,13 +1,13 @@
 @defcomp co2forcing begin
-    c_CO2concentration=Parameter(index=[time],unit="ppbv")
-    f0_CO2baseforcing=Parameter(unit="W/m2", default=1.735)
-    fslope_CO2forcingslope=Parameter(unit="W/m2", default=5.5)
-    c0_baseCO2conc=Parameter(unit="ppbv", default=395000.)
-    f_CO2forcing=Variable(index=[time],unit="W/m2")
+    c_CO2concentration = Parameter(index=[time], unit="ppbv")
+    f0_CO2baseforcing = Parameter(unit="W/m2", default=1.735)
+    fslope_CO2forcingslope = Parameter(unit="W/m2", default=5.5)
+    c0_baseCO2conc = Parameter(unit="ppbv", default=395000.)
+    f_CO2forcing = Variable(index=[time], unit="W/m2")
 
     function run_timestep(p, v, d, t)
 
-        #eq.13 in Hope 2006
-        v.f_CO2forcing[t]=p.f0_CO2baseforcing+p.fslope_CO2forcingslope*log(p.c_CO2concentration[t]/p.c0_baseCO2conc)
+        # eq.13 in Hope 2006
+        v.f_CO2forcing[t] = p.f0_CO2baseforcing + p.fslope_CO2forcingslope * log(p.c_CO2concentration[t] / p.c0_baseCO2conc)
     end
 end
diff --git a/src/components/ClimateTemperature.jl b/src/components/ClimateTemperature.jl
index 18e50af1..4665b698 100644
--- a/src/components/ClimateTemperature.jl
+++ b/src/components/ClimateTemperature.jl
@@ -34,11 +34,11 @@
     rtl_g_landtemperature = Variable(index=[time], unit="degreeC")
     rto_g_oceantemperature = Variable(index=[time], unit="degreeC")
     rt_g_globaltemperature = Variable(index=[time], unit="degreeC")
-    rt_g0_baseglobaltemp=Variable(unit="degreeC") #needed for feedback in CO2 cycle component
-    rtl_g0_baselandtemp=Variable(unit="degreeC") #needed for feedback in CH4 and N2O cycles
+    rt_g0_baseglobaltemp = Variable(unit="degreeC") # needed for feedback in CO2 cycle component
+    rtl_g0_baselandtemp = Variable(unit="degreeC") # needed for feedback in CH4 and N2O cycles
 
     function init(p, v, d)
-        #calculate global baseline temperature from initial regional temperatures
+        # calculate global baseline temperature from initial regional temperatures
 
         ocean_prop_ortion = 1. - sum(p.area) / 510000000.
 
@@ -46,7 +46,7 @@
         v.rtl_g0_baselandtemp = sum(p.rtl_0_realizedtemperature' .* p.area') / sum(p.area)
 
         # initial ocean and global temperatures
-        rto_g0_baseoceantemp = v.rtl_g0_baselandtemp/ p.rlo_ratiolandocean
+        rto_g0_baseoceantemp = v.rtl_g0_baselandtemp / p.rlo_ratiolandocean
         v.rt_g0_baseglobaltemp = ocean_prop_ortion * rto_g0_baseoceantemp + (1. - ocean_prop_ortion) * v.rtl_g0_baselandtemp
     end
 
@@ -79,7 +79,7 @@
             end
         else
             for rr in d.region
-                v.rt_realizedtemperature[tt, rr] = v.rt_realizedtemperature[tt-1, rr] + (1 - exp(-(p.y_year[tt] - p.y_year[tt-1]) / p.frt_warminghalflife)) * (v.et_equilibriumtemperature[tt, rr] - v.rt_realizedtemperature[tt-1, rr])
+                v.rt_realizedtemperature[tt, rr] = v.rt_realizedtemperature[tt - 1, rr] + (1 - exp(-(p.y_year[tt] - p.y_year[tt - 1]) / p.frt_warminghalflife)) * (v.et_equilibriumtemperature[tt, rr] - v.rt_realizedtemperature[tt - 1, rr])
             end
         end
 
diff --git a/src/components/Discontinuity.jl b/src/components/Discontinuity.jl
index 3a860c0f..8d849b0b 100644
--- a/src/components/Discontinuity.jl
+++ b/src/components/Discontinuity.jl
@@ -1,77 +1,77 @@
 @defcomp Discontinuity begin
 
-  region = Index()
-  y_year=Parameter(index=[time], unit="year")
-  y_year_0 = Parameter(unit="year")
+    region = Index()
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
 
-  rand_discontinuity = Parameter(unit="unitless", default=.5)
+    rand_discontinuity = Parameter(unit="unitless", default=.5)
 
-  irefeqdis_eqdiscimpact=Variable(index=[region], unit="%")
-  wincf_weightsfactor=Parameter(index=[region], unit="unitless", default=[1, 0.8, 0.8, 0.4, 0.8, 0.8, 0.6, 0.6])
-  wdis_gdplostdisc=Parameter(unit="%", default=15.)
+    irefeqdis_eqdiscimpact = Variable(index=[region], unit="%")
+    wincf_weightsfactor = Parameter(index=[region], unit="unitless", default=[1, 0.8, 0.8, 0.4, 0.8, 0.8, 0.6, 0.6])
+    wdis_gdplostdisc = Parameter(unit="%", default=15.)
 
-  igdpeqdis_eqdiscimpact=Variable(index=[time,region], unit="%")
-  rgdp_per_cap_NonMarketRemainGDP=Parameter(index=[time,region], unit="\$/person")
-  GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
-  ipow_incomeexponent=Parameter(unit="unitless", default=-0.13333333333333333)
+    igdpeqdis_eqdiscimpact = Variable(index=[time,region], unit="%")
+    rgdp_per_cap_NonMarketRemainGDP = Parameter(index=[time,region], unit="\$/person")
+    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
+    ipow_incomeexponent = Parameter(unit="unitless", default=-0.13333333333333333)
 
-  igdp_realizeddiscimpact=Variable(index=[time,region], unit="%")
-  occurdis_occurrencedummy=Variable(index=[time], unit="unitless")
-  expfdis_discdecay=Variable(index=[time], unit="unitless")
+    igdp_realizeddiscimpact = Variable(index=[time,region], unit="%")
+    occurdis_occurrencedummy = Variable(index=[time], unit="unitless")
+    expfdis_discdecay = Variable(index=[time], unit="unitless")
 
-  distau_discontinuityexponent=Parameter(unit="unitless", default=90.)
+    distau_discontinuityexponent = Parameter(unit="unitless", default=90.)
 
-  idis_lossfromdisc=Variable(index=[time], unit="degreeC")
-  tdis_tolerabilitydisc=Parameter(unit="degreeC", default=3.)
-  rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
-  pdis_probability=Parameter(unit="%/degreeC", default=20.)
+    idis_lossfromdisc = Variable(index=[time], unit="degreeC")
+    tdis_tolerabilitydisc = Parameter(unit="degreeC", default=3.)
+    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
+    pdis_probability = Parameter(unit="%/degreeC", default=20.)
 
-  isatg_saturationmodification = Parameter(unit="unitless")
-  isat_satdiscimpact=Variable(index=[time,region], unit="%")
+    isatg_saturationmodification = Parameter(unit="unitless")
+    isat_satdiscimpact = Variable(index=[time,region], unit="%")
 
-  isat_per_cap_DiscImpactperCapinclSaturation=Variable(index=[time,region], unit="\$/person")
-  rcons_per_cap_DiscRemainConsumption=Variable(index=[time, region], unit = "\$/person")
-  rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit = "\$/person")   
+    isat_per_cap_DiscImpactperCapinclSaturation = Variable(index=[time,region], unit="\$/person")
+    rcons_per_cap_DiscRemainConsumption = Variable(index=[time, region], unit="\$/person")
+    rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit="\$/person")   
 
     function run_timestep(p, v, d, t)
 
         v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - p.tdis_tolerabilitydisc)
 
         if is_first(t)
-            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
+            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                 v.occurdis_occurrencedummy[t] = 1
             else
                 v.occurdis_occurrencedummy[t] = 0
             end
-            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year_0)/p.distau_discontinuityexponent)
+            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year_0) / p.distau_discontinuityexponent)
         else
-            if v.idis_lossfromdisc[t]*(p.pdis_probability/100) > p.rand_discontinuity
+            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                 v.occurdis_occurrencedummy[t] = 1
-            elseif v.occurdis_occurrencedummy[t-1] == 1
+            elseif v.occurdis_occurrencedummy[t - 1] == 1
                 v.occurdis_occurrencedummy[t] = 1
             else
                 v.occurdis_occurrencedummy[t] = 0
             end
-            v.expfdis_discdecay[t]=exp(-(p.y_year[t] - p.y_year[t-1])/p.distau_discontinuityexponent)
+            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year[t - 1]) / p.distau_discontinuityexponent)
         end
 
         for r in d.region
-            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor[r]*p.wdis_gdplostdisc
+            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor[r] * p.wdis_gdplostdisc
 
-            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent
+            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent
 
             if is_first(t)
-                v.igdp_realizeddiscimpact[t,r]=v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*v.igdpeqdis_eqdiscimpact[t,r]
+                v.igdp_realizeddiscimpact[t,r] = v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * v.igdpeqdis_eqdiscimpact[t,r]
             else
-                v.igdp_realizeddiscimpact[t,r]=v.igdp_realizeddiscimpact[t-1,r]+v.occurdis_occurrencedummy[t]*(1-v.expfdis_discdecay[t])*(v.igdpeqdis_eqdiscimpact[t,r]-v.igdp_realizeddiscimpact[t-1,r])
+                v.igdp_realizeddiscimpact[t,r] = v.igdp_realizeddiscimpact[t - 1,r] + v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * (v.igdpeqdis_eqdiscimpact[t,r] - v.igdp_realizeddiscimpact[t - 1,r])
             end
 
             if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                 v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
             else
-                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100-p.isatg_saturationmodification)*((v.igdp_realizeddiscimpact[t,r]-p.isatg_saturationmodification)/((100-p.isatg_saturationmodification)+(v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
+                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100 - p.isatg_saturationmodification) * ((v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification) / ((100 - p.isatg_saturationmodification) + (v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
             end
-            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r]/100)*p.rgdp_per_cap_NonMarketRemainGDP[t,r]
+            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r] / 100) * p.rgdp_per_cap_NonMarketRemainGDP[t,r]
             v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]
         end
     end
diff --git a/src/components/EquityWeighting.jl b/src/components/EquityWeighting.jl
index 01c4b606..2d849526 100644
--- a/src/components/EquityWeighting.jl
+++ b/src/components/EquityWeighting.jl
@@ -8,7 +8,7 @@
     # Impacts across all gases
     pop_population = Parameter(index=[time, region], unit="million person")
 
-    #Total and Per-Capita Abatement and Adaptation Costs
+    # Total and Per-Capita Abatement and Adaptation Costs
     tct_percap_totalcosts_total = Parameter(index=[time, region], unit="\$/person")
     act_adaptationcosts_total = Parameter(index=[time, region], unit="\$million")
     act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")
@@ -109,7 +109,7 @@
             if is_first(tt)
                 v.yp_yearsperiod[TimestepIndex(1)] = p.y_year[TimestepIndex(1)] - p.y_year_0
             else
-                v.yp_yearsperiod[tt] = p.y_year[tt] - p.y_year[tt-1]
+                v.yp_yearsperiod[tt] = p.y_year[tt] - p.y_year[tt - 1]
             end
 
             if is_first(tt)
diff --git a/src/components/GDP.jl b/src/components/GDP.jl
index c9e7371d..5a4fffdc 100644
--- a/src/components/GDP.jl
+++ b/src/components/GDP.jl
@@ -13,21 +13,21 @@
     # Parameters
     y_year_0          = Parameter(unit="year")
     y_year            = Parameter(index=[time], unit="year")
-    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year") #From p.32 of Hope 2009
-    gdp_0             = Parameter(index=[region], unit="\$M") #GDP in y_year_0
-    save_savingsrate  = Parameter(unit="%", default=15.00) #pp33 PAGE09 documentation, "savings rate".
+    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year") # From p.32 of Hope 2009
+    gdp_0             = Parameter(index=[region], unit="\$M") # GDP in y_year_0
+    save_savingsrate  = Parameter(unit="%", default=15.00) # pp33 PAGE09 documentation, "savings rate".
     pop0_initpopulation = Parameter(index=[region], unit="million person")
-    pop_population    = Parameter(index=[time,region],unit="million person")
+    pop_population    = Parameter(index=[time,region], unit="million person")
 
     # Saturation, used in impacts
-    isat0_initialimpactfxnsaturation = Parameter(unit="unitless", default=33.333333333333336) #pp34 PAGE09 documentation
+    isat0_initialimpactfxnsaturation = Parameter(unit="unitless", default=33.333333333333336) # pp34 PAGE09 documentation
     isatg_impactfxnsaturation = Variable(unit="unitless")
 
     function init(p, v, d)
 
-        v.isatg_impactfxnsaturation = p.isat0_initialimpactfxnsaturation * (1 - p.save_savingsrate/100)
+        v.isatg_impactfxnsaturation = p.isat0_initialimpactfxnsaturation * (1 - p.save_savingsrate / 100)
         for rr in d.region
-            v.cons_percap_consumption_0[rr] = (p.gdp_0[rr] / p.pop0_initpopulation[rr])*(1 - p.save_savingsrate / 100)
+            v.cons_percap_consumption_0[rr] = (p.gdp_0[rr] / p.pop0_initpopulation[rr]) * (1 - p.save_savingsrate / 100)
         end
     end
 
@@ -37,23 +37,23 @@
         if is_first(t)
             ylo_periodstart = p.y_year_0
         else
-            ylo_periodstart = (p.y_year[t] + p.y_year[t-1]) / 2
+            ylo_periodstart = (p.y_year[t] + p.y_year[t - 1]) / 2
         end
 
         if t.t == length(p.y_year)
             yhi_periodend = p.y_year[t]
         else
-            yhi_periodend = (p.y_year[t] + p.y_year[t+1]) / 2
+            yhi_periodend = (p.y_year[t] + p.y_year[t + 1]) / 2
         end
 
-        v.yagg_periodspan[t] = yhi_periodend- ylo_periodstart
+        v.yagg_periodspan[t] = yhi_periodend - ylo_periodstart
 
         for r in d.region
-            #eq.28 in Hope 2002
+            # eq.28 in Hope 2002
             if is_first(t)
-                v.gdp[t, r] = p.gdp_0[r] * (1 + (p.grw_gdpgrowthrate[t,r]/100))^(p.y_year[t] - p.y_year_0)
+                v.gdp[t, r] = p.gdp_0[r] * (1 + (p.grw_gdpgrowthrate[t,r] / 100))^(p.y_year[t] - p.y_year_0)
             else
-                v.gdp[t, r] = v.gdp[t-1, r] * (1 + (p.grw_gdpgrowthrate[t,r]/100))^(p.y_year[t] - p.y_year[t-1])
+                v.gdp[t, r] = v.gdp[t - 1, r] * (1 + (p.grw_gdpgrowthrate[t,r] / 100))^(p.y_year[t] - p.y_year[t - 1])
             end
             v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
             v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]
diff --git a/src/components/LGcycle.jl b/src/components/LGcycle.jl
index e9a79144..f5f8d4c5 100644
--- a/src/components/LGcycle.jl
+++ b/src/components/LGcycle.jl
@@ -1,54 +1,54 @@
 @defcomp LGcycle begin
-    e_globalLGemissions=Parameter(index=[time],unit="Mtonne/year")
-    e_0globalLGemissions=Parameter(unit="Mtonne/year", default=557.2112715473608)
-    c_LGconcentration=Variable(index=[time],unit="ppbv")
-    pic_preindustconcLG=Parameter(unit="ppbv", default=0.)
-    exc_excessconcLG=Variable(unit="ppbv")
-    c0_LGconcbaseyr=Parameter(unit="ppbv", default=0.11)
-    re_remainLG=Variable(index=[time],unit="Mtonne")
-    nte_natLGemissions=Variable(index=[time],unit="Mtonne/year")
-    air_LGfractioninatm=Parameter(unit="%", default=100.)
-    tea_LGemissionstoatm=Variable(index=[time],unit="Mtonne/year")
-    teay_LGemissionstoatm=Variable(index=[time],unit="Mtonne/t")
-    y_year=Parameter(index=[time],unit="year")
-    y_year_0=Parameter(unit="year")
-    res_LGatmlifetime=Parameter(unit="year", default=1000.)
-    den_LGdensity=Parameter(unit="Mtonne/ppbv", default=100000.)
-    stim_LGemissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
-    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=0.9258270139190647)
-    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")
-    re_remainLGbase=Variable(unit="Mtonne")
+    e_globalLGemissions = Parameter(index=[time], unit="Mtonne/year")
+    e_0globalLGemissions = Parameter(unit="Mtonne/year", default=557.2112715473608)
+    c_LGconcentration = Variable(index=[time], unit="ppbv")
+    pic_preindustconcLG = Parameter(unit="ppbv", default=0.)
+    exc_excessconcLG = Variable(unit="ppbv")
+    c0_LGconcbaseyr = Parameter(unit="ppbv", default=0.11)
+    re_remainLG = Variable(index=[time], unit="Mtonne")
+    nte_natLGemissions = Variable(index=[time], unit="Mtonne/year")
+    air_LGfractioninatm = Parameter(unit="%", default=100.)
+    tea_LGemissionstoatm = Variable(index=[time], unit="Mtonne/year")
+    teay_LGemissionstoatm = Variable(index=[time], unit="Mtonne/t")
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
+    res_LGatmlifetime = Parameter(unit="year", default=1000.)
+    den_LGdensity = Parameter(unit="Mtonne/ppbv", default=100000.)
+    stim_LGemissionfeedback = Parameter(unit="Mtonne/degreeC", default=0.)
+    rtl_g0_baselandtemp = Parameter(unit="degreeC", default=0.9258270139190647)
+    rtl_g_landtemperature = Parameter(index=[time], unit="degreeC")
+    re_remainLGbase = Variable(unit="Mtonne")
 
     function run_timestep(p, v, d, t)
         if is_first(t)
-            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
-            nte0=p.stim_LGemissionfeedback*p.rtl_g0_baselandtemp
-            v.nte_natLGemissions[t]=p.stim_LGemissionfeedback*p.rtl_g0_baselandtemp
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            tea0=(p.e_0globalLGemissions+nte0)*p.air_LGfractioninatm/100
-            v.tea_LGemissionstoatm[t]=(p.e_globalLGemissions[t]+v.nte_natLGemissions[t])*p.air_LGfractioninatm/100
-            v.teay_LGemissionstoatm[t]=(v.tea_LGemissionstoatm[t]+tea0)*(p.y_year[t]-p.y_year_0)/2
-            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
-            v.exc_excessconcLG=p.c0_LGconcbaseyr-p.pic_preindustconcLG
-            #Eq. 2 from Hope (2006) - base-year remaining emissions
-            v.re_remainLGbase=v.exc_excessconcLG*p.den_LGdensity
-            v.re_remainLG[t]=v.re_remainLGbase*exp(-(p.y_year[t]-p.y_year_0)/p.res_LGatmlifetime)+
-                v.teay_LGemissionstoatm[t]*p.res_LGatmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_LGatmlifetime))/(p.y_year[t]-p.y_year_0)
+            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
+            nte0 = p.stim_LGemissionfeedback * p.rtl_g0_baselandtemp
+            v.nte_natLGemissions[t] = p.stim_LGemissionfeedback * p.rtl_g0_baselandtemp
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            tea0 = (p.e_0globalLGemissions + nte0) * p.air_LGfractioninatm / 100
+            v.tea_LGemissionstoatm[t] = (p.e_globalLGemissions[t] + v.nte_natLGemissions[t]) * p.air_LGfractioninatm / 100
+            v.teay_LGemissionstoatm[t] = (v.tea_LGemissionstoatm[t] + tea0) * (p.y_year[t] - p.y_year_0) / 2
+            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
+            v.exc_excessconcLG = p.c0_LGconcbaseyr - p.pic_preindustconcLG
+            # Eq. 2 from Hope (2006) - base-year remaining emissions
+            v.re_remainLGbase = v.exc_excessconcLG * p.den_LGdensity
+            v.re_remainLG[t] = v.re_remainLGbase * exp(-(p.y_year[t] - p.y_year_0) / p.res_LGatmlifetime) +
+                v.teay_LGemissionstoatm[t] * p.res_LGatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.res_LGatmlifetime)) / (p.y_year[t] - p.y_year_0)
         else
-            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
-            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
-            v.nte_natLGemissions[t]=p.stim_LGemissionfeedback*p.rtl_g_landtemperature[t-1]
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_LGemissionstoatm[t]=(p.e_globalLGemissions[t]+v.nte_natLGemissions[t])*p.air_LGfractioninatm/100
-            #eq.7 from Hope (2006) - average emissions to atm over time period
-            v.teay_LGemissionstoatm[t]=(v.tea_LGemissionstoatm[t]+v.tea_LGemissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2
-            #eq.10 from Hope (2006) - remaining emissions in atmosphere
-            v.re_remainLG[t]=v.re_remainLG[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_LGatmlifetime)+
-                v.teay_LGemissionstoatm[t]*p.res_LGatmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_LGatmlifetime))/(p.y_year[t]-p.y_year[t-1])
+            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
+            # Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
+            v.nte_natLGemissions[t] = p.stim_LGemissionfeedback * p.rtl_g_landtemperature[t - 1]
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_LGemissionstoatm[t] = (p.e_globalLGemissions[t] + v.nte_natLGemissions[t]) * p.air_LGfractioninatm / 100
+            # eq.7 from Hope (2006) - average emissions to atm over time period
+            v.teay_LGemissionstoatm[t] = (v.tea_LGemissionstoatm[t] + v.tea_LGemissionstoatm[t - 1]) * (p.y_year[t] - p.y_year[t - 1]) / 2
+            # eq.10 from Hope (2006) - remaining emissions in atmosphere
+            v.re_remainLG[t] = v.re_remainLG[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_LGatmlifetime) +
+                v.teay_LGemissionstoatm[t] * p.res_LGatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_LGatmlifetime)) / (p.y_year[t] - p.y_year[t - 1])
         end
 
-    #eq.11 from Hope(2006) - LG concentration
-        v.c_LGconcentration[t]=p.pic_preindustconcLG+v.exc_excessconcLG*v.re_remainLG[t]/v.re_remainLGbase
+    # eq.11 from Hope(2006) - LG concentration
+        v.c_LGconcentration[t] = p.pic_preindustconcLG + v.exc_excessconcLG * v.re_remainLG[t] / v.re_remainLGbase
 
     end
 end
diff --git a/src/components/LGemissions.jl b/src/components/LGemissions.jl
index 998a1e22..c6009793 100644
--- a/src/components/LGemissions.jl
+++ b/src/components/LGemissions.jl
@@ -1,24 +1,24 @@
 # SF6 defined as "linear gas" or gas 4 in PAGE 2009; equations are the same as for gas 3 (CH4) in PAGE2002
 @defcomp LGemissions begin
-    region=Index()
+    region = Index()
 
     # global emissions
-    e_globalLGemissions=Variable(index=[time],unit="Mtonne/year")
+    e_globalLGemissions = Variable(index=[time], unit="Mtonne/year")
     # baseline emissions
-    e0_baselineLGemissions=Parameter(index=[region],unit="Mtonne/year")
+    e0_baselineLGemissions = Parameter(index=[region], unit="Mtonne/year")
     # regional emissions
-    e_regionalLGemissions=Variable(index=[time,region],unit="Mtonne/year")
+    e_regionalLGemissions = Variable(index=[time,region], unit="Mtonne/year")
     # growth rate by region
-    er_LGemissionsgrowth=Parameter(index=[time,region],unit="%")
+    er_LGemissionsgrowth = Parameter(index=[time,region], unit="%")
 
     function run_timestep(p, v, d, t)
 
-        #eq.4 in Hope (2006) - regional LG emissions as % change from baseline
+        # eq.4 in Hope (2006) - regional LG emissions as % change from baseline
         for r in d.region
-            v.e_regionalLGemissions[t,r]=p.er_LGemissionsgrowth[t,r]*p.e0_baselineLGemissions[r]/100
+            v.e_regionalLGemissions[t,r] = p.er_LGemissionsgrowth[t,r] * p.e0_baselineLGemissions[r] / 100
         end
 
-        #eq. 5 in Hope (2006) - global LG emissions are sum of regional emissions
-        v.e_globalLGemissions[t]=sum(v.e_regionalLGemissions[t,:])
+        # eq. 5 in Hope (2006) - global LG emissions are sum of regional emissions
+        v.e_globalLGemissions[t] = sum(v.e_regionalLGemissions[t,:])
     end
 end
diff --git a/src/components/LGforcing.jl b/src/components/LGforcing.jl
index db113396..16d99103 100644
--- a/src/components/LGforcing.jl
+++ b/src/components/LGforcing.jl
@@ -1,15 +1,15 @@
 @defcomp LGforcing begin
 
-    c_LGconcentration=Parameter(index=[time],unit="ppbv")
-    c0_LGconcbaseyr=Parameter(unit="ppbv", default=0.11)
-    f_LGforcing=Variable(index=[time],unit="W/m2")
-    f0_LGforcingbase=Parameter(unit="W/m2", default=0.022)
-    fslope_LGforcingslope=Parameter(unit="W/m2", default=0.2)
+    c_LGconcentration = Parameter(index=[time], unit="ppbv")
+    c0_LGconcbaseyr = Parameter(unit="ppbv", default=0.11)
+    f_LGforcing = Variable(index=[time], unit="W/m2")
+    f0_LGforcingbase = Parameter(unit="W/m2", default=0.022)
+    fslope_LGforcingslope = Parameter(unit="W/m2", default=0.2)
 
     function run_timestep(p, v, d, t)
 
-        #eq.13 in Hope 2006
-        v.f_LGforcing[t]=p.f0_LGforcingbase+p.fslope_LGforcingslope*(p.c_LGconcentration[t]-p.c0_LGconcbaseyr)
+        # eq.13 in Hope 2006
+        v.f_LGforcing[t] = p.f0_LGforcingbase + p.fslope_LGforcingslope * (p.c_LGconcentration[t] - p.c0_LGconcbaseyr)
 
     end
 end
diff --git a/src/components/MarketDamages.jl b/src/components/MarketDamages.jl
index dd35534b..b3664705 100644
--- a/src/components/MarketDamages.jl
+++ b/src/components/MarketDamages.jl
@@ -2,76 +2,76 @@
     region = Index()
     y_year = Parameter(index=[time], unit="year")
 
-    #incoming parameters from Climate
+    # incoming parameters from Climate
     rtl_realizedtemperature = Parameter(index=[time, region], unit="degreeC")
 
-    #tolerability variables
+    # tolerability variables
     atl_adjustedtolerableleveloftemprise = Parameter(index=[time,region], unit="degreeC")
-    imp_actualreduction = Parameter(index=[time, region], unit= "%")
+    imp_actualreduction = Parameter(index=[time, region], unit="%")
     i_regionalimpact = Variable(index=[time, region], unit="degreeC")
 
-    #impact Parameters
-    rcons_per_cap_SLRRemainConsumption = Parameter(index=[time, region], unit = "\$/person")
-    rgdp_per_cap_SLRRemainGDP = Parameter(index=[time, region], unit = "\$/person")
+    # impact Parameters
+    rcons_per_cap_SLRRemainConsumption = Parameter(index=[time, region], unit="\$/person")
+    rgdp_per_cap_SLRRemainGDP = Parameter(index=[time, region], unit="\$/person")
 
-    save_savingsrate = Parameter(unit= "%", default=15.)
-    wincf_weightsfactor =Parameter(index=[region], unit="unitless")
-    W_MarketImpactsatCalibrationTemp =Parameter(unit="%GDP", default=0.5)
-    ipow_MarketIncomeFxnExponent =Parameter(default=-0.13333333333333333)
-    iben_MarketInitialBenefit=Parameter(default=.1333333333333)
+    save_savingsrate = Parameter(unit="%", default=15.)
+    wincf_weightsfactor = Parameter(index=[region], unit="unitless")
+    W_MarketImpactsatCalibrationTemp = Parameter(unit="%GDP", default=0.5)
+    ipow_MarketIncomeFxnExponent = Parameter(default=-0.13333333333333333)
+    iben_MarketInitialBenefit = Parameter(default=.1333333333333)
     tcal_CalibrationTemp = Parameter(unit="degreeC", default=3.)
     GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
 
-    #impact variables
+    # impact variables
     isatg_impactfxnsaturation = Parameter(unit="unitless")
-    rcons_per_cap_MarketRemainConsumption = Variable(index=[time, region], unit = "\$/person")
-    rgdp_per_cap_MarketRemainGDP = Variable(index=[time, region], unit = "\$/person")
-    iref_ImpactatReferenceGDPperCap=Variable(index=[time, region])
-    igdp_ImpactatActualGDPperCap=Variable(index=[time, region])
-    impmax_maxtempriseforadaptpolicyM = Parameter(index=[region], unit= "degreeC")
+    rcons_per_cap_MarketRemainConsumption = Variable(index=[time, region], unit="\$/person")
+    rgdp_per_cap_MarketRemainGDP = Variable(index=[time, region], unit="\$/person")
+    iref_ImpactatReferenceGDPperCap = Variable(index=[time, region])
+    igdp_ImpactatActualGDPperCap = Variable(index=[time, region])
+    impmax_maxtempriseforadaptpolicyM = Parameter(index=[region], unit="degreeC")
 
-    isat_ImpactinclSaturationandAdaptation= Variable(index=[time,region])
+    isat_ImpactinclSaturationandAdaptation = Variable(index=[time,region])
     isat_per_cap_ImpactperCapinclSaturationandAdaptation = Variable(index=[time,region], unit="\$/person")
-    pow_MarketImpactExponent=Parameter(unit="", default=2.16666666666665)
+    pow_MarketImpactExponent = Parameter(unit="", default=2.16666666666665)
 
     function run_timestep(p, v, d, t)
 
         for r in d.region
-            #calculate tolerability
-            if (p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r]) < 0
+            # calculate tolerability
+            if (p.rtl_realizedtemperature[t,r] - p.atl_adjustedtolerableleveloftemprise[t,r]) < 0
                 v.i_regionalimpact[t,r] = 0
             else
-                v.i_regionalimpact[t,r] = p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r]
+                v.i_regionalimpact[t,r] = p.rtl_realizedtemperature[t,r] - p.atl_adjustedtolerableleveloftemprise[t,r]
             end
 
-            v.iref_ImpactatReferenceGDPperCap[t,r]= p.wincf_weightsfactor[r]*((p.W_MarketImpactsatCalibrationTemp + p.iben_MarketInitialBenefit * p.tcal_CalibrationTemp)*
-                (v.i_regionalimpact[t,r]/p.tcal_CalibrationTemp)^p.pow_MarketImpactExponent - v.i_regionalimpact[t,r] * p.iben_MarketInitialBenefit)
+            v.iref_ImpactatReferenceGDPperCap[t,r] = p.wincf_weightsfactor[r] * ((p.W_MarketImpactsatCalibrationTemp + p.iben_MarketInitialBenefit * p.tcal_CalibrationTemp) *
+                (v.i_regionalimpact[t,r] / p.tcal_CalibrationTemp)^p.pow_MarketImpactExponent - v.i_regionalimpact[t,r] * p.iben_MarketInitialBenefit)
 
-            v.igdp_ImpactatActualGDPperCap[t,r]= v.iref_ImpactatReferenceGDPperCap[t,r]*
-                (p.rgdp_per_cap_SLRRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_MarketIncomeFxnExponent
+            v.igdp_ImpactatActualGDPperCap[t,r] = v.iref_ImpactatReferenceGDPperCap[t,r] *
+                (p.rgdp_per_cap_SLRRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_MarketIncomeFxnExponent
 
             if v.igdp_ImpactatActualGDPperCap[t,r] < p.isatg_impactfxnsaturation
                 v.isat_ImpactinclSaturationandAdaptation[t,r] = v.igdp_ImpactatActualGDPperCap[t,r]
             else
-                v.isat_ImpactinclSaturationandAdaptation[t,r] = p.isatg_impactfxnsaturation+
-                    ((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)*
-                    ((v.igdp_ImpactatActualGDPperCap[t,r]-p.isatg_impactfxnsaturation)/
-                    (((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)+
-                    (v.igdp_ImpactatActualGDPperCap[t,r]-
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = p.isatg_impactfxnsaturation +
+                    ((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) *
+                    ((v.igdp_ImpactatActualGDPperCap[t,r] - p.isatg_impactfxnsaturation) /
+                    (((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) +
+                    (v.igdp_ImpactatActualGDPperCap[t,r] -
                     p.isatg_impactfxnsaturation)))
-                end
+            end
 
             if v.i_regionalimpact[t,r] < p.impmax_maxtempriseforadaptpolicyM[r]
-                v.isat_ImpactinclSaturationandAdaptation[t,r]=v.isat_ImpactinclSaturationandAdaptation[t,r]*(1-p.imp_actualreduction[t,r]/100)
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] * (1 - p.imp_actualreduction[t,r] / 100)
             else
                 v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] *
-                    (1-(p.imp_actualreduction[t,r]/100)* p.impmax_maxtempriseforadaptpolicyM[r] /
+                    (1 - (p.imp_actualreduction[t,r] / 100) * p.impmax_maxtempriseforadaptpolicyM[r] /
                     v.i_regionalimpact[t,r])
             end
 
-            v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptation[t,r]/100)*p.rgdp_per_cap_SLRRemainGDP[t,r]
+            v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptation[t,r] / 100) * p.rgdp_per_cap_SLRRemainGDP[t,r]
             v.rcons_per_cap_MarketRemainConsumption[t,r] = p.rcons_per_cap_SLRRemainConsumption[t,r] - v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r]
-            v.rgdp_per_cap_MarketRemainGDP[t,r] = v.rcons_per_cap_MarketRemainConsumption[t,r]/(1-p.save_savingsrate/100)
+            v.rgdp_per_cap_MarketRemainGDP[t,r] = v.rcons_per_cap_MarketRemainConsumption[t,r] / (1 - p.save_savingsrate / 100)
         end
 
     end
diff --git a/src/components/N2Ocycle.jl b/src/components/N2Ocycle.jl
index b56cf79b..5aa4b0cc 100644
--- a/src/components/N2Ocycle.jl
+++ b/src/components/N2Ocycle.jl
@@ -1,56 +1,56 @@
 @defcomp n2ocycle begin
-    e_globalN2Oemissions=Parameter(index=[time],unit="Mtonne/year")
-    e_0globalN2Oemissions=Parameter(unit="Mtonne/year", default=11.046520000000001)
-    c_N2Oconcentration=Variable(index=[time],unit="ppbv")
-    pic_preindustconcN2O=Parameter(unit="ppbv", default=270.)
-    exc_excessconcN2O=Variable(unit="ppbv")
-    c0_N2Oconcbaseyr=Parameter(unit="ppbv", default=322.)
-    re_remainN2O=Variable(index=[time],unit="Mtonne")
-    re_remainN2Obase=Variable(unit="Mtonne")
-    nte_natN2Oemissions=Variable(index=[time],unit="Mtonne/year")
-    air_N2Ofractioninatm=Parameter(unit="%", default=100.)
-    tea_N2Oemissionstoatm=Variable(index=[time],unit="Mtonne/year")
-    teay_N2Oemissionstoatm=Variable(index=[time],unit="Mtonne/t")
-    y_year=Parameter(index=[time],unit="year")
-    y_year_0=Parameter(unit="year")
-    res_N2Oatmlifetime=Parameter(unit="year", default=114.)
-    den_N2Odensity=Parameter(unit="Mtonne/ppbv", default=7.8)
-    stim_N2Oemissionfeedback=Parameter(unit="Mtonne/degreeC", default=0.)
-    rtl_g0_baselandtemp=Parameter(unit="degreeC", default=0.9258270139190647)
-    rtl_g_landtemperature=Parameter(index=[time],unit="degreeC")
+    e_globalN2Oemissions = Parameter(index=[time], unit="Mtonne/year")
+    e_0globalN2Oemissions = Parameter(unit="Mtonne/year", default=11.046520000000001)
+    c_N2Oconcentration = Variable(index=[time], unit="ppbv")
+    pic_preindustconcN2O = Parameter(unit="ppbv", default=270.)
+    exc_excessconcN2O = Variable(unit="ppbv")
+    c0_N2Oconcbaseyr = Parameter(unit="ppbv", default=322.)
+    re_remainN2O = Variable(index=[time], unit="Mtonne")
+    re_remainN2Obase = Variable(unit="Mtonne")
+    nte_natN2Oemissions = Variable(index=[time], unit="Mtonne/year")
+    air_N2Ofractioninatm = Parameter(unit="%", default=100.)
+    tea_N2Oemissionstoatm = Variable(index=[time], unit="Mtonne/year")
+    teay_N2Oemissionstoatm = Variable(index=[time], unit="Mtonne/t")
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
+    res_N2Oatmlifetime = Parameter(unit="year", default=114.)
+    den_N2Odensity = Parameter(unit="Mtonne/ppbv", default=7.8)
+    stim_N2Oemissionfeedback = Parameter(unit="Mtonne/degreeC", default=0.)
+    rtl_g0_baselandtemp = Parameter(unit="degreeC", default=0.9258270139190647)
+    rtl_g_landtemperature = Parameter(index=[time], unit="degreeC")
 
     function run_timestep(p, v, d, t)
-        #note that Hope (2009) states that Equations 1-12 for methane also apply to N2O
+        # note that Hope (2009) states that Equations 1-12 for methane also apply to N2O
 
         if is_first(t)
-            #eq.3 from Hope (2006) - natural emissions feedback, using global temperatures calculated in ClimateTemperature component
-            nte_0=p.stim_N2Oemissionfeedback*p.rtl_g0_baselandtemp
-            v.nte_natN2Oemissions[t]=p.stim_N2Oemissionfeedback*p.rtl_g0_baselandtemp
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_N2Oemissionstoatm[t]=(p.e_globalN2Oemissions[t]+v.nte_natN2Oemissions[t])*p.air_N2Ofractioninatm/100
-            tea_0=(p.e_0globalN2Oemissions+nte_0)*p.air_N2Ofractioninatm/100
-            v.teay_N2Oemissionstoatm[t]=(v.tea_N2Oemissionstoatm[t]+tea_0)*(p.y_year[t]-p.y_year_0)/2
-            #adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
-            v.exc_excessconcN2O=p.c0_N2Oconcbaseyr-p.pic_preindustconcN2O
-            #Eq. 2 from Hope (2006) - base-year remaining emissions
-            v.re_remainN2Obase=v.exc_excessconcN2O*p.den_N2Odensity
-            v.re_remainN2O[t]=v.re_remainN2Obase*exp(-(p.y_year[t]-p.y_year_0)/p.res_N2Oatmlifetime)+
-                v.teay_N2Oemissionstoatm[t]*p.res_N2Oatmlifetime*(1-exp(-(p.y_year[t]-p.y_year_0)/p.res_N2Oatmlifetime))/(p.y_year[t]-p.y_year_0)
+            # eq.3 from Hope (2006) - natural emissions feedback, using global temperatures calculated in ClimateTemperature component
+            nte_0 = p.stim_N2Oemissionfeedback * p.rtl_g0_baselandtemp
+            v.nte_natN2Oemissions[t] = p.stim_N2Oemissionfeedback * p.rtl_g0_baselandtemp
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_N2Oemissionstoatm[t] = (p.e_globalN2Oemissions[t] + v.nte_natN2Oemissions[t]) * p.air_N2Ofractioninatm / 100
+            tea_0 = (p.e_0globalN2Oemissions + nte_0) * p.air_N2Ofractioninatm / 100
+            v.teay_N2Oemissionstoatm[t] = (v.tea_N2Oemissionstoatm[t] + tea_0) * (p.y_year[t] - p.y_year_0) / 2
+            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
+            v.exc_excessconcN2O = p.c0_N2Oconcbaseyr - p.pic_preindustconcN2O
+            # Eq. 2 from Hope (2006) - base-year remaining emissions
+            v.re_remainN2Obase = v.exc_excessconcN2O * p.den_N2Odensity
+            v.re_remainN2O[t] = v.re_remainN2Obase * exp(-(p.y_year[t] - p.y_year_0) / p.res_N2Oatmlifetime) +
+                v.teay_N2Oemissionstoatm[t] * p.res_N2Oatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.res_N2Oatmlifetime)) / (p.y_year[t] - p.y_year_0)
         else
-            #eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
-            #Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
-            v.nte_natN2Oemissions[t]=p.stim_N2Oemissionfeedback*p.rtl_g_landtemperature[t-1]
-            #eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
-            v.tea_N2Oemissionstoatm[t]=(p.e_globalN2Oemissions[t]+v.nte_natN2Oemissions[t])*p.air_N2Ofractioninatm/100
-            #eq.7 from Hope (2006) - average emissions to atm over time period
-            v.teay_N2Oemissionstoatm[t]=(v.tea_N2Oemissionstoatm[t]+v.tea_N2Oemissionstoatm[t-1])*(p.y_year[t]-p.y_year[t-1])/2
-            #eq.10 from Hope (2006) - remaining emissions in atmosphere
-            v.re_remainN2O[t]=v.re_remainN2O[t-1]*exp(-(p.y_year[t]-p.y_year[t-1])/p.res_N2Oatmlifetime)+
-                v.teay_N2Oemissionstoatm[t]*p.res_N2Oatmlifetime*(1-exp(-(p.y_year[t]-p.y_year[t-1])/p.res_N2Oatmlifetime))/(p.y_year[t]-p.y_year[t-1])
+            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
+            # Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
+            v.nte_natN2Oemissions[t] = p.stim_N2Oemissionfeedback * p.rtl_g_landtemperature[t - 1]
+            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
+            v.tea_N2Oemissionstoatm[t] = (p.e_globalN2Oemissions[t] + v.nte_natN2Oemissions[t]) * p.air_N2Ofractioninatm / 100
+            # eq.7 from Hope (2006) - average emissions to atm over time period
+            v.teay_N2Oemissionstoatm[t] = (v.tea_N2Oemissionstoatm[t] + v.tea_N2Oemissionstoatm[t - 1]) * (p.y_year[t] - p.y_year[t - 1]) / 2
+            # eq.10 from Hope (2006) - remaining emissions in atmosphere
+            v.re_remainN2O[t] = v.re_remainN2O[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_N2Oatmlifetime) +
+                v.teay_N2Oemissionstoatm[t] * p.res_N2Oatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_N2Oatmlifetime)) / (p.y_year[t] - p.y_year[t - 1])
         end
 
-        #eq.11 from Hope(2006) - N2O concentration
-        v.c_N2Oconcentration[t]=p.pic_preindustconcN2O+v.exc_excessconcN2O*v.re_remainN2O[t]/v.re_remainN2Obase
+        # eq.11 from Hope(2006) - N2O concentration
+        v.c_N2Oconcentration[t] = p.pic_preindustconcN2O + v.exc_excessconcN2O * v.re_remainN2O[t] / v.re_remainN2Obase
 
     end
 end
diff --git a/src/components/N2Oemissions.jl b/src/components/N2Oemissions.jl
index 3b2fab8a..6526876c 100644
--- a/src/components/N2Oemissions.jl
+++ b/src/components/N2Oemissions.jl
@@ -1,20 +1,20 @@
 @defcomp n2oemissions begin
-    region=Index()
+    region = Index()
 
-    e_globalN2Oemissions=Variable(index=[time],unit="Mtonne/year")
-    e0_baselineN2Oemissions=Parameter(index=[region],unit="Mtonne/year")
-    e_regionalN2Oemissions=Variable(index=[time,region],unit="Mtonne/year")
-    er_N2Oemissionsgrowth=Parameter(index=[time,region],unit="%")
+    e_globalN2Oemissions = Variable(index=[time], unit="Mtonne/year")
+    e0_baselineN2Oemissions = Parameter(index=[region], unit="Mtonne/year")
+    e_regionalN2Oemissions = Variable(index=[time,region], unit="Mtonne/year")
+    er_N2Oemissionsgrowth = Parameter(index=[time,region], unit="%")
 
     function run_timestep(p, v, d, t)
-        #note that Hope (2009) states that Equations 1-12 for methane also apply to N2O
+        # note that Hope (2009) states that Equations 1-12 for methane also apply to N2O
 
-        #eq.4 in Hope (2006) - regional N2O emissions as % change from baseline
+        # eq.4 in Hope (2006) - regional N2O emissions as % change from baseline
         for r in d.region
-            v.e_regionalN2Oemissions[t,r]=p.er_N2Oemissionsgrowth[t,r]*p.e0_baselineN2Oemissions[r]/100
+            v.e_regionalN2Oemissions[t,r] = p.er_N2Oemissionsgrowth[t,r] * p.e0_baselineN2Oemissions[r] / 100
         end
 
-        #eq. 5 in Hope (2006) - global N2O emissions are sum of regional emissions
-        v.e_globalN2Oemissions[t]=sum(v.e_regionalN2Oemissions[t,:])
+        # eq. 5 in Hope (2006) - global N2O emissions are sum of regional emissions
+        v.e_globalN2Oemissions[t] = sum(v.e_regionalN2Oemissions[t,:])
     end
 end
diff --git a/src/components/N2Oforcing.jl b/src/components/N2Oforcing.jl
index 5daf67ae..6ff50e0b 100644
--- a/src/components/N2Oforcing.jl
+++ b/src/components/N2Oforcing.jl
@@ -1,24 +1,24 @@
 
 
 @defcomp n2oforcing begin
-    c_N2Oconcentration=Parameter(index=[time],unit="ppbv")
-    c_CH4concentration=Parameter(index=[time],unit="ppbv")
-    f0_N2Obaseforcing=Parameter(unit="W/m2", default=0.180)
-    fslope_N2Oforcingslope=Parameter(unit="W/m2", default=0.12)
-    c0_baseN2Oconc=Parameter(unit="ppbv", default=322.)
-    c0_baseCH4conc=Parameter(unit="ppbv", default=1860.)
-    f_N2Oforcing=Variable(index=[time],unit="W/m2")
-    over_baseoverlap=Variable(unit="W/m2")
+    c_N2Oconcentration = Parameter(index=[time], unit="ppbv")
+    c_CH4concentration = Parameter(index=[time], unit="ppbv")
+    f0_N2Obaseforcing = Parameter(unit="W/m2", default=0.180)
+    fslope_N2Oforcingslope = Parameter(unit="W/m2", default=0.12)
+    c0_baseN2Oconc = Parameter(unit="ppbv", default=322.)
+    c0_baseCH4conc = Parameter(unit="ppbv", default=1860.)
+    f_N2Oforcing = Variable(index=[time], unit="W/m2")
+    over_baseoverlap = Variable(unit="W/m2")
     
     function run_timestep(p, v, d, t)
 
-        #from p.16 in Hope 2009
+        # from p.16 in Hope 2009
         if is_first(t)
-            #calculate baseline forcing overlap in first time period
-            v.over_baseoverlap=-0.47*log(1+2.0e-5*(p.c0_baseN2Oconc*p.c0_baseCH4conc)^0.75+5.3e-15*p.c0_baseCH4conc*(p.c0_baseCH4conc*p.c0_baseN2Oconc)^1.52)
+            # calculate baseline forcing overlap in first time period
+            v.over_baseoverlap = -0.47 * log(1 + 2.0e-5 * (p.c0_baseN2Oconc * p.c0_baseCH4conc)^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c0_baseN2Oconc)^1.52)
         end
 
-        over=-0.47*log(1+2.0e-5*(p.c0_baseCH4conc*p.c_N2Oconcentration[t])^0.75+5.3e-15*p.c0_baseCH4conc*(p.c0_baseCH4conc*p.c_N2Oconcentration[t])^1.52)
-        v.f_N2Oforcing[t]=p.f0_N2Obaseforcing+p.fslope_N2Oforcingslope*(sqrt(p.c_N2Oconcentration[t])-sqrt(p.c0_baseN2Oconc))+over-v.over_baseoverlap
+        over = -0.47 * log(1 + 2.0e-5 * (p.c0_baseCH4conc * p.c_N2Oconcentration[t])^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c_N2Oconcentration[t])^1.52)
+        v.f_N2Oforcing[t] = p.f0_N2Obaseforcing + p.fslope_N2Oforcingslope * (sqrt(p.c_N2Oconcentration[t]) - sqrt(p.c0_baseN2Oconc)) + over - v.over_baseoverlap
     end
 end
diff --git a/src/components/NonMarketDamages.jl b/src/components/NonMarketDamages.jl
index 6f00624e..eed553fd 100644
--- a/src/components/NonMarketDamages.jl
+++ b/src/components/NonMarketDamages.jl
@@ -5,78 +5,78 @@
 
     y_year = Parameter(index=[time], unit="year")
 
-    #incoming parameters from Climate
+    # incoming parameters from Climate
     rtl_realizedtemperature = Parameter(index=[time, region], unit="degreeC")
 
-    #tolerability parameters
-    impmax_maxtempriseforadaptpolicyNM = Parameter(index=[region], unit= "degreeC")
+    # tolerability parameters
+    impmax_maxtempriseforadaptpolicyNM = Parameter(index=[region], unit="degreeC")
     atl_adjustedtolerableleveloftemprise = Parameter(index=[time,region], unit="degreeC")
-    imp_actualreduction = Parameter(index=[time, region], unit= "%")
+    imp_actualreduction = Parameter(index=[time, region], unit="%")
 
-    #tolerability variables
+    # tolerability variables
     i_regionalimpact = Variable(index=[time, region], unit="degreeC")
 
-    #impact Parameters
-    rcons_per_cap_MarketRemainConsumption = Parameter(index=[time, region], unit = "\$/person")
-    rgdp_per_cap_MarketRemainGDP = Parameter(index=[time, region], unit = "\$/person")
+    # impact Parameters
+    rcons_per_cap_MarketRemainConsumption = Parameter(index=[time, region], unit="\$/person")
+    rgdp_per_cap_MarketRemainGDP = Parameter(index=[time, region], unit="\$/person")
 
-    save_savingsrate = Parameter(unit= "%", default=15.)
-    wincf_weightsfactor =Parameter(index=[region], unit="unitless")
-    w_NonImpactsatCalibrationTemp =Parameter(unit="%GDP", default=0.5333333333333333)
-    ipow_NonMarketIncomeFxnExponent =Parameter(unit="unitless", default=0.)
-    iben_NonMarketInitialBenefit=Parameter(unit="%GDP/degreeC", default=0.08333333333333333)
+    save_savingsrate = Parameter(unit="%", default=15.)
+    wincf_weightsfactor = Parameter(index=[region], unit="unitless")
+    w_NonImpactsatCalibrationTemp = Parameter(unit="%GDP", default=0.5333333333333333)
+    ipow_NonMarketIncomeFxnExponent = Parameter(unit="unitless", default=0.)
+    iben_NonMarketInitialBenefit = Parameter(unit="%GDP/degreeC", default=0.08333333333333333)
     tcal_CalibrationTemp = Parameter(unit="degreeC", default=3.)
     GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
     pow_NonMarketExponent = Parameter(unit="", default=2.1666666666666665)
 
-    #impact variables
+    # impact variables
     isatg_impactfxnsaturation = Parameter(unit="unitless")
-    rcons_per_cap_NonMarketRemainConsumption = Variable(index=[time, region], unit = "\$/person")
-    rgdp_per_cap_NonMarketRemainGDP = Variable(index=[time, region], unit = "\$/person")
-    iref_ImpactatReferenceGDPperCap=Variable(index=[time, region], unit="%")
-    igdp_ImpactatActualGDPperCap=Variable(index=[time, region], unit="%")
-    isat_ImpactinclSaturationandAdaptation= Variable(index=[time,region], unit="\$")
+    rcons_per_cap_NonMarketRemainConsumption = Variable(index=[time, region], unit="\$/person")
+    rgdp_per_cap_NonMarketRemainGDP = Variable(index=[time, region], unit="\$/person")
+    iref_ImpactatReferenceGDPperCap = Variable(index=[time, region], unit="%")
+    igdp_ImpactatActualGDPperCap = Variable(index=[time, region], unit="%")
+    isat_ImpactinclSaturationandAdaptation = Variable(index=[time,region], unit="\$")
     isat_per_cap_ImpactperCapinclSaturationandAdaptation = Variable(index=[time,region], unit="\$/person")
 
     function run_timestep(p, v, d, t)
 
         for r in d.region
 
-            if p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r] < 0
+            if p.rtl_realizedtemperature[t,r] - p.atl_adjustedtolerableleveloftemprise[t,r] < 0
                 v.i_regionalimpact[t,r] = 0
             else
-                v.i_regionalimpact[t,r] = p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r]
+                v.i_regionalimpact[t,r] = p.rtl_realizedtemperature[t,r] - p.atl_adjustedtolerableleveloftemprise[t,r]
             end
 
-            v.iref_ImpactatReferenceGDPperCap[t,r]= p.wincf_weightsfactor[r]*
-                ((p.w_NonImpactsatCalibrationTemp + p.iben_NonMarketInitialBenefit *p.tcal_CalibrationTemp)*
-                    (v.i_regionalimpact[t,r]/p.tcal_CalibrationTemp)^p.pow_NonMarketExponent - v.i_regionalimpact[t,r] * p.iben_NonMarketInitialBenefit)
-
-            v.igdp_ImpactatActualGDPperCap[t,r]= v.iref_ImpactatReferenceGDPperCap[t,r]*
-                (p.rgdp_per_cap_MarketRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_NonMarketIncomeFxnExponent
-
-                if v.igdp_ImpactatActualGDPperCap[t,r] < p.isatg_impactfxnsaturation
-                    v.isat_ImpactinclSaturationandAdaptation[t,r] = v.igdp_ImpactatActualGDPperCap[t,r]
-                else
-                    v.isat_ImpactinclSaturationandAdaptation[t,r] = p.isatg_impactfxnsaturation+
-                        ((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)*
-                        ((v.igdp_ImpactatActualGDPperCap[t,r]-p.isatg_impactfxnsaturation)/
-                        (((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)+
-                        (v.igdp_ImpactatActualGDPperCap[t,r]-
+            v.iref_ImpactatReferenceGDPperCap[t,r] = p.wincf_weightsfactor[r] *
+                ((p.w_NonImpactsatCalibrationTemp + p.iben_NonMarketInitialBenefit * p.tcal_CalibrationTemp) *
+                    (v.i_regionalimpact[t,r] / p.tcal_CalibrationTemp)^p.pow_NonMarketExponent - v.i_regionalimpact[t,r] * p.iben_NonMarketInitialBenefit)
+
+            v.igdp_ImpactatActualGDPperCap[t,r] = v.iref_ImpactatReferenceGDPperCap[t,r] *
+                (p.rgdp_per_cap_MarketRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_NonMarketIncomeFxnExponent
+
+            if v.igdp_ImpactatActualGDPperCap[t,r] < p.isatg_impactfxnsaturation
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.igdp_ImpactatActualGDPperCap[t,r]
+            else
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = p.isatg_impactfxnsaturation +
+                        ((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) *
+                        ((v.igdp_ImpactatActualGDPperCap[t,r] - p.isatg_impactfxnsaturation) /
+                        (((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) +
+                        (v.igdp_ImpactatActualGDPperCap[t,r] -
                         p.isatg_impactfxnsaturation)))
-                    end
+            end
 
-                if v.i_regionalimpact[t,r] < p.impmax_maxtempriseforadaptpolicyNM[r]
-                    v.isat_ImpactinclSaturationandAdaptation[t,r]=v.isat_ImpactinclSaturationandAdaptation[t,r]*(1-p.imp_actualreduction[t,r]/100)
-                else
-                    v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] *
-                        (1-(p.imp_actualreduction[t,r]/100)* p.impmax_maxtempriseforadaptpolicyNM[r] /
+            if v.i_regionalimpact[t,r] < p.impmax_maxtempriseforadaptpolicyNM[r]
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] * (1 - p.imp_actualreduction[t,r] / 100)
+            else
+                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] *
+                        (1 - (p.imp_actualreduction[t,r] / 100) * p.impmax_maxtempriseforadaptpolicyNM[r] /
                         v.i_regionalimpact[t,r])
-                end
+            end
 
-            v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptation[t,r]/100)*p.rgdp_per_cap_MarketRemainGDP[t,r]
+            v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptation[t,r] / 100) * p.rgdp_per_cap_MarketRemainGDP[t,r]
             v.rcons_per_cap_NonMarketRemainConsumption[t,r] = p.rcons_per_cap_MarketRemainConsumption[t,r] - v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r]
-            v.rgdp_per_cap_NonMarketRemainGDP[t,r] = v.rcons_per_cap_NonMarketRemainConsumption[t,r]/(1-p.save_savingsrate/100)
+            v.rgdp_per_cap_NonMarketRemainGDP[t,r] = v.rcons_per_cap_NonMarketRemainConsumption[t,r] / (1 - p.save_savingsrate / 100)
         end
     end
 end
diff --git a/src/components/Population.jl b/src/components/Population.jl
index cb4aa03b..5d8f6011 100644
--- a/src/components/Population.jl
+++ b/src/components/Population.jl
@@ -17,16 +17,16 @@
         for rr in d.region
             # Eq.28 in Hope 2002 (defined for GDP, but also applies to population)
             if is_first(tt)
-                v.pop_population[tt, rr] = p.pop0_initpopulation[rr] * (1 + p.popgrw_populationgrowth[tt, rr]/100)^(p.y_year[tt] - p.y_year_0)
+                v.pop_population[tt, rr] = p.pop0_initpopulation[rr] * (1 + p.popgrw_populationgrowth[tt, rr] / 100)^(p.y_year[tt] - p.y_year_0)
             else
-                v.pop_population[tt, rr] = v.pop_population[tt-1, rr] * (1 + p.popgrw_populationgrowth[tt, rr]/100)^(p.y_year[tt] - p.y_year[tt-1])
+                v.pop_population[tt, rr] = v.pop_population[tt - 1, rr] * (1 + p.popgrw_populationgrowth[tt, rr] / 100)^(p.y_year[tt] - p.y_year[tt - 1])
             end
         end
     end
 end
 
 # Still need this function in order to set the parameters than depend on 
-# readpagedata, which takes model as an input. These cannot be set using 
+    # readpagedata, which takes model as an input. These cannot be set using 
 # the default keyword arg for now.
 function addpopulation(model::Model)
     populationcomp = add_comp!(model, Population)
diff --git a/src/components/SLRDamages.jl b/src/components/SLRDamages.jl
index 929b0230..f44b412a 100644
--- a/src/components/SLRDamages.jl
+++ b/src/components/SLRDamages.jl
@@ -6,80 +6,80 @@
     y_year = Parameter(index=[time], unit="year")
     y_year_0 = Parameter(unit="year")
 
-    #incoming parameters from SeaLevelRise
+    # incoming parameters from SeaLevelRise
     s_sealevel = Parameter(index=[time], unit="m")
-    #incoming parameters to calculate consumption per capita after Costs
+    # incoming parameters to calculate consumption per capita after Costs
     cons_percap_consumption = Parameter(index=[time, region], unit="\$/person")
-    tct_per_cap_totalcostspercap = Parameter(index=[time,region], unit= "\$/person")
+    tct_per_cap_totalcostspercap = Parameter(index=[time,region], unit="\$/person")
     act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")
 
-    #component parameters
-    impmax_maxSLRforadaptpolicySLR = Parameter(index=[region], unit= "m") 
+    # component parameters
+    impmax_maxSLRforadaptpolicySLR = Parameter(index=[region], unit="m") 
 
-    save_savingsrate = Parameter(unit= "%", default=15.00) #pp33 PAGE09 documentation, "savings rate".
-    wincf_weightsfactor =Parameter(index=[region], unit="unitless")
-    W_SatCalibrationSLR =Parameter(default=1.0) #pp33 PAGE09 documentation, "Sea level impact at calibration sea level rise"
-    ipow_SLRIncomeFxnExponent =Parameter(default=-0.30)
-    pow_SLRImpactFxnExponent=Parameter(default=0.7333333333333334)
-    iben_SLRInitialBenefit=Parameter(default=0.00)
+    save_savingsrate = Parameter(unit="%", default=15.00) # pp33 PAGE09 documentation, "savings rate".
+    wincf_weightsfactor = Parameter(index=[region], unit="unitless")
+    W_SatCalibrationSLR = Parameter(default=1.0) # pp33 PAGE09 documentation, "Sea level impact at calibration sea level rise"
+    ipow_SLRIncomeFxnExponent = Parameter(default=-0.30)
+    pow_SLRImpactFxnExponent = Parameter(default=0.7333333333333334)
+    iben_SLRInitialBenefit = Parameter(default=0.00)
     scal_calibrationSLR = Parameter(default=0.5)
     GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
 
-    #component variables
-    cons_percap_aftercosts = Variable(index=[time, region], unit = "\$/person")
-    gdp_percap_aftercosts = Variable(index=[time, region], unit = "\$/person")
+    # component variables
+    cons_percap_aftercosts = Variable(index=[time, region], unit="\$/person")
+    gdp_percap_aftercosts = Variable(index=[time, region], unit="\$/person")
 
     atl_adjustedtolerablelevelofsealevelrise = Parameter(index=[time,region], unit="m") # meter
     imp_actualreductionSLR = Parameter(index=[time, region], unit="%")
     i_regionalimpactSLR = Variable(index=[time, region], unit="m")
 
     iref_ImpactatReferenceGDPperCapSLR = Variable(index=[time, region])
-    igdp_ImpactatActualGDPperCapSLR=Variable(index=[time, region])
+    igdp_ImpactatActualGDPperCapSLR = Variable(index=[time, region])
 
     isatg_impactfxnsaturation = Parameter(unit="unitless")
-    isat_ImpactinclSaturationandAdaptationSLR= Variable(index=[time,region])
+    isat_ImpactinclSaturationandAdaptationSLR = Variable(index=[time,region])
     isat_per_cap_SLRImpactperCapinclSaturationandAdaptation = Variable(index=[time, region], unit="\$/person")
 
-    rcons_per_cap_SLRRemainConsumption = Variable(index=[time, region], unit = "\$/person") #include?
-    rgdp_per_cap_SLRRemainGDP = Variable(index=[time, region], unit = "\$/person")
+    rcons_per_cap_SLRRemainConsumption = Variable(index=[time, region], unit="\$/person") # include?
+    rgdp_per_cap_SLRRemainGDP = Variable(index=[time, region], unit="\$/person")
 
     function run_timestep(p, v, d, t)
 
         for r in d.region
             v.cons_percap_aftercosts[t, r] = p.cons_percap_consumption[t, r] - p.tct_per_cap_totalcostspercap[t, r] - p.act_percap_adaptationcosts[t, r]
-            v.gdp_percap_aftercosts[t,r]=v.cons_percap_aftercosts[t, r]/(1 - p.save_savingsrate/100)
+            v.gdp_percap_aftercosts[t,r] = v.cons_percap_aftercosts[t, r] / (1 - p.save_savingsrate / 100)
 
-            if (p.s_sealevel[t]-p.atl_adjustedtolerablelevelofsealevelrise[t,r]) < 0
+            if (p.s_sealevel[t] - p.atl_adjustedtolerablelevelofsealevelrise[t,r]) < 0
                 v.i_regionalimpactSLR[t,r] = 0
             else
-                v.i_regionalimpactSLR[t,r] = p.s_sealevel[t]-p.atl_adjustedtolerablelevelofsealevelrise[t,r]
+                v.i_regionalimpactSLR[t,r] = p.s_sealevel[t] - p.atl_adjustedtolerablelevelofsealevelrise[t,r]
             end
 
-            v.iref_ImpactatReferenceGDPperCapSLR[t,r]= p.wincf_weightsfactor[r]*((p.W_SatCalibrationSLR + p.iben_SLRInitialBenefit * p.scal_calibrationSLR)*
-                (v.i_regionalimpactSLR[t,r]/p.scal_calibrationSLR)^p.pow_SLRImpactFxnExponent - v.i_regionalimpactSLR[t,r] * p.iben_SLRInitialBenefit)
+            v.iref_ImpactatReferenceGDPperCapSLR[t,r] = p.wincf_weightsfactor[r] * ((p.W_SatCalibrationSLR + p.iben_SLRInitialBenefit * p.scal_calibrationSLR) *
+                (v.i_regionalimpactSLR[t,r] / p.scal_calibrationSLR)^p.pow_SLRImpactFxnExponent - v.i_regionalimpactSLR[t,r] * p.iben_SLRInitialBenefit)
 
-            v.igdp_ImpactatActualGDPperCapSLR[t,r]= v.iref_ImpactatReferenceGDPperCapSLR[t,r]*
-                    (v.gdp_percap_aftercosts[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_SLRIncomeFxnExponent
+            v.igdp_ImpactatActualGDPperCapSLR[t,r] = v.iref_ImpactatReferenceGDPperCapSLR[t,r] *
+                    (v.gdp_percap_aftercosts[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_SLRIncomeFxnExponent
 
             if v.igdp_ImpactatActualGDPperCapSLR[t,r] < p.isatg_impactfxnsaturation
                 v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.igdp_ImpactatActualGDPperCapSLR[t,r]
             else
-                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = p.isatg_impactfxnsaturation+
-                    ((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)*
-                    ((v.igdp_ImpactatActualGDPperCapSLR[t,r]-p.isatg_impactfxnsaturation)/
-                    (((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)+
-                    (v.igdp_ImpactatActualGDPperCapSLR[t,r]- p.isatg_impactfxnsaturation)))
-                end
+                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = p.isatg_impactfxnsaturation +
+                    ((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) *
+                    ((v.igdp_ImpactatActualGDPperCapSLR[t,r] - p.isatg_impactfxnsaturation) /
+                    (((100 - p.save_savingsrate) - p.isatg_impactfxnsaturation) +
+                    (v.igdp_ImpactatActualGDPperCapSLR[t,r] - p.isatg_impactfxnsaturation)))
+            end
             if v.i_regionalimpactSLR[t,r] < p.impmax_maxSLRforadaptpolicySLR[r]
-                v.isat_ImpactinclSaturationandAdaptationSLR[t,r]=v.isat_ImpactinclSaturationandAdaptationSLR[t,r]*(1-p.imp_actualreductionSLR[t,r]/100)
+                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.isat_ImpactinclSaturationandAdaptationSLR[t,r] * (1 - p.imp_actualreductionSLR[t,r] / 100)
             else
-                v.isat_ImpactinclSaturationandAdaptationSLR[t,r]=v.isat_ImpactinclSaturationandAdaptationSLR[t,r]*(1-(p.imp_actualreductionSLR[t,r]/100)* p.impmax_maxSLRforadaptpolicySLR[r] /
+                v.isat_ImpactinclSaturationandAdaptationSLR[t,r] = v.isat_ImpactinclSaturationandAdaptationSLR[t,r] * (1 - (p.imp_actualreductionSLR[t,r] / 100) * p.impmax_maxSLRforadaptpolicySLR[r] /
                     v.i_regionalimpactSLR[t,r])
-                end
+            end
 
-                v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptationSLR[t,r]/100)*v.gdp_percap_aftercosts[t,r]
-                v.rcons_per_cap_SLRRemainConsumption[t,r] = v.cons_percap_aftercosts[t,r] - v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r]
-                v.rgdp_per_cap_SLRRemainGDP[t,r] = v.rcons_per_cap_SLRRemainConsumption[t,r]/(1-p.save_savingsrate/100)
+            v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptationSLR[t,r] / 100) * v.gdp_percap_aftercosts[t,r]
+            v.rcons_per_cap_SLRRemainConsumption[t,r] = v.cons_percap_aftercosts[t,r] - v.isat_per_cap_SLRImpactperCapinclSaturationandAdaptation[t,r]
+            v.rgdp_per_cap_SLRRemainGDP[t,r] = v.rcons_per_cap_SLRRemainConsumption[t,r] / (1 - p.save_savingsrate / 100)
 
         end
 
@@ -93,6 +93,6 @@ end
 function addslrdamages(model::Model)
     SLRDamagescomp  = add_comp!(model, SLRDamages)
 
-    set_param!(model, :SLRDamages, :impmax_maxSLRforadaptpolicySLR, readpagedata(model, "data/impmax_sealevel.csv"))
+set_param!(model, :SLRDamages, :impmax_maxSLRforadaptpolicySLR, readpagedata(model, "data/impmax_sealevel.csv"))
     return SLRDamagescomp
 end
diff --git a/src/components/SeaLevelRise.jl b/src/components/SeaLevelRise.jl
index d77bbea5..c0abffc3 100644
--- a/src/components/SeaLevelRise.jl
+++ b/src/components/SeaLevelRise.jl
@@ -4,36 +4,36 @@
 
   # Parameters
 
-  rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
-  sltemp_SLtemprise=Parameter(unit = "m-degreeC", default=1.7333333333333334)
-  sla_SLbaselinerise=Parameter(unit = "m", default=1.00)
-  sltau_SLresponsetime=Parameter(unit = "years", default=1000.)
-  s0_initialSL=Parameter(unit = "m", default=0.15)
-  y_year=Parameter(index=[time], unit="year")
-  y_year_0=Parameter(unit="year")
+    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
+    sltemp_SLtemprise = Parameter(unit="m-degreeC", default=1.7333333333333334)
+    sla_SLbaselinerise = Parameter(unit="m", default=1.00)
+    sltau_SLresponsetime = Parameter(unit="years", default=1000.)
+    s0_initialSL = Parameter(unit="m", default=0.15)
+    y_year = Parameter(index=[time], unit="year")
+    y_year_0 = Parameter(unit="year")
 
   # Variables
 
-  es_equilibriumSL=Variable(index=[time], unit = "m")
-  s_sealevel=Variable(index=[time], unit = "m")
-  expfs_exponential=Variable(index=[time], unit = "unitless")
-  yp_timestep=Variable(index=[time], unit = "years")
+    es_equilibriumSL = Variable(index=[time], unit="m")
+    s_sealevel = Variable(index=[time], unit="m")
+    expfs_exponential = Variable(index=[time], unit="unitless")
+    yp_timestep = Variable(index=[time], unit="years")
 
 
-  function run_timestep(p, v, d, t)
-    if is_first(t)
-      v.yp_timestep[t]=p.y_year[TimestepIndex(1)] - p.y_year_0
-      v.es_equilibriumSL[t]=p.sltemp_SLtemprise*p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
-      v.expfs_exponential[t]=exp(-v.yp_timestep[t]/p.sltau_SLresponsetime)
-      v.s_sealevel[t]=p.s0_initialSL + (v.es_equilibriumSL[t] - p.s0_initialSL)*(1-v.expfs_exponential[t])
+    function run_timestep(p, v, d, t)
+        if is_first(t)
+            v.yp_timestep[t] = p.y_year[TimestepIndex(1)] - p.y_year_0
+            v.es_equilibriumSL[t] = p.sltemp_SLtemprise * p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
+            v.expfs_exponential[t] = exp(-v.yp_timestep[t] / p.sltau_SLresponsetime)
+            v.s_sealevel[t] = p.s0_initialSL + (v.es_equilibriumSL[t] - p.s0_initialSL) * (1 - v.expfs_exponential[t])
 
-    else
-      v.yp_timestep[t]=p.y_year[t] - p.y_year[t-1]
-      v.es_equilibriumSL[t]=p.sltemp_SLtemprise*p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
-      v.expfs_exponential[t]=exp(-v.yp_timestep[t]/p.sltau_SLresponsetime)
-      v.s_sealevel[t]=v.s_sealevel[t-1] + (v.es_equilibriumSL[t] - v.s_sealevel[t-1])*(1-v.expfs_exponential[t])
+        else
+            v.yp_timestep[t] = p.y_year[t] - p.y_year[t - 1]
+            v.es_equilibriumSL[t] = p.sltemp_SLtemprise * p.rt_g_globaltemperature[t] + p.sla_SLbaselinerise
+            v.expfs_exponential[t] = exp(-v.yp_timestep[t] / p.sltau_SLresponsetime)
+            v.s_sealevel[t] = v.s_sealevel[t - 1] + (v.es_equilibriumSL[t] - v.s_sealevel[t - 1]) * (1 - v.expfs_exponential[t])
 
+        end
     end
-  end
 end
 
diff --git a/src/components/SulphateForcing.jl b/src/components/SulphateForcing.jl
index da24cd11..d212875b 100644
--- a/src/components/SulphateForcing.jl
+++ b/src/components/SulphateForcing.jl
@@ -6,7 +6,7 @@
     se0_sulphateemissionsbase = Parameter(index=[region], unit="TgS/year")
     pse_sulphatevsbase = Parameter(index=[time, region], unit="%")
     se_sulphateemissions = Variable(index=[time, region], unit="TgS/year")
-    area = Parameter(index=[region], unit ="km^2")
+    area = Parameter(index=[region], unit="km^2")
 
     sfx_sulphateflux = Variable(index=[time, region], unit="TgS/km^2/yr")
 
@@ -27,7 +27,7 @@
             # Update for Eq. 18 from Hope (2009) - sulfate radiative forcing effect
             bigSFD0 = p.d_sulphateforcingbase * bigSFX0[rr] / (sum(bigSFX0 .* p.area) / sum(p.area))
             fsd_term = bigSFD0 * v.sfx_sulphateflux[tt,rr] / bigSFX0[rr]
-            fsi_term = p.ind_slopeSEforcing_indirect/log(2) * log((p.nf_naturalsfx[rr] + v.sfx_sulphateflux[tt, rr]) / p.nf_naturalsfx[rr])
+            fsi_term = p.ind_slopeSEforcing_indirect / log(2) * log((p.nf_naturalsfx[rr] + v.sfx_sulphateflux[tt, rr]) / p.nf_naturalsfx[rr])
 
             v.fs_sulphateforcing[tt, rr] = fsd_term + fsi_term
         end
diff --git a/src/components/TotalAbatementCosts.jl b/src/components/TotalAbatementCosts.jl
index 18413f71..74747e36 100644
--- a/src/components/TotalAbatementCosts.jl
+++ b/src/components/TotalAbatementCosts.jl
@@ -7,16 +7,16 @@
     tc_totalcosts_ch4 = Parameter(index=[time, region], unit="\$million")
     tc_totalcosts_n2o = Parameter(index=[time, region], unit="\$million")
     tc_totalcosts_linear = Parameter(index=[time, region], unit="\$million")
-    pop_population = Parameter(index=[time, region], unit= "million person")
+    pop_population = Parameter(index=[time, region], unit="million person")
 
-    tct_totalcosts = Variable(index=[time,region], unit= "\$million")
-    tct_per_cap_totalcostspercap = Variable(index=[time,region], unit= "\$/person")
+    tct_totalcosts = Variable(index=[time,region], unit="\$million")
+    tct_per_cap_totalcostspercap = Variable(index=[time,region], unit="\$/person")
         
     function run_timestep(p, v, d, t)
 
         for r in d.region
             v.tct_totalcosts[t,r] = p.tc_totalcosts_co2[t,r] + p.tc_totalcosts_n2o[t,r] + p.tc_totalcosts_ch4[t,r] + p.tc_totalcosts_linear[t,r]
-            v.tct_per_cap_totalcostspercap[t,r] = v.tct_totalcosts[t,r]/p.pop_population[t,r]
+            v.tct_per_cap_totalcostspercap[t,r] = v.tct_totalcosts[t,r] / p.pop_population[t,r]
         end
     end
 end
diff --git a/src/components/TotalAdaptationCosts.jl b/src/components/TotalAdaptationCosts.jl
index 637f44ed..9a0a5e64 100644
--- a/src/components/TotalAdaptationCosts.jl
+++ b/src/components/TotalAdaptationCosts.jl
@@ -2,8 +2,8 @@
 
 @defcomp TotalAdaptationCosts begin
     region = Index()
-    #Total Adaptation Costs
-    pop_population = Parameter(index=[time, region], unit= "million person")
+    # Total Adaptation Costs
+    pop_population = Parameter(index=[time, region], unit="million person")
     ac_adaptationcosts_economic = Parameter(index=[time, region], unit="\$million")
     ac_adaptationcosts_noneconomic = Parameter(index=[time, region], unit="\$million")
     ac_adaptationcosts_sealevelrise = Parameter(index=[time, region], unit="\$million")
@@ -15,7 +15,7 @@
 
         for r in d.region
             v.act_adaptationcosts_total[t,r] = p.ac_adaptationcosts_economic[t,r] + p.ac_adaptationcosts_sealevelrise[t,r] + p.ac_adaptationcosts_noneconomic[t,r]
-            v.act_percap_adaptationcosts[t,r] = v.act_adaptationcosts_total[t,r]/p.pop_population[t,r]
+            v.act_percap_adaptationcosts[t,r] = v.act_adaptationcosts_total[t,r] / p.pop_population[t,r]
         end
     end
 end
diff --git a/src/components/TotalCosts.jl b/src/components/TotalCosts.jl
index 3fad2570..91d95de1 100644
--- a/src/components/TotalCosts.jl
+++ b/src/components/TotalCosts.jl
@@ -1,30 +1,30 @@
 @defcomp TotalCosts begin
 
     # Per capita costs to aggregate
-    abatement_costs_percap_peryear          = Parameter(index = [time, region], unit = "\$/person") # Total abatement costs per capita per year
-    adaptation_costs_percap_peryear         = Parameter(index = [time, region], unit = "\$/person") # Total adaptation costs per capita per year
-    slr_damages_percap_peryear              = Parameter(index = [time, region], unit = "\$/person") # SLR damages per capita per year
-    market_damages_percap_peryear           = Parameter(index = [time, region], unit = "\$/person") # Market damages per capita per year
-    non_market_damages_percap_peryear       = Parameter(index = [time, region], unit = "\$/person") # Non market damages per capita per year
-    discontinuity_damages_percap_peryear    = Parameter(index = [time, region], unit = "\$/person") # Discontinuity damages per capita per year
+    abatement_costs_percap_peryear          = Parameter(index=[time, region], unit="\$/person") # Total abatement costs per capita per year
+    adaptation_costs_percap_peryear         = Parameter(index=[time, region], unit="\$/person") # Total adaptation costs per capita per year
+    slr_damages_percap_peryear              = Parameter(index=[time, region], unit="\$/person") # SLR damages per capita per year
+    market_damages_percap_peryear           = Parameter(index=[time, region], unit="\$/person") # Market damages per capita per year
+    non_market_damages_percap_peryear       = Parameter(index=[time, region], unit="\$/person") # Non market damages per capita per year
+    discontinuity_damages_percap_peryear    = Parameter(index=[time, region], unit="\$/person") # Discontinuity damages per capita per year
 
     # other parameters 
-    population = Parameter(index = [time, region], unit = "million person")
-    period_length = Parameter(index = [time], unit="year")
+    population = Parameter(index=[time, region], unit="million person")
+    period_length = Parameter(index=[time], unit="year")
 
     # Undiscounted costs to calculate
 
-    total_damages_percap_peryear    = Variable(index = [time, region], unit = "\$/person")  # Undiscounted total damages per capita per year   
-    total_damages_peryear           = Variable(index = [time, region], unit = "\$million")  # Undiscounted total damages per year 
-    total_damages_aggregated        = Variable(index = [time, region], unit = "\$million")  # Undiscounted total damages per period
+    total_damages_percap_peryear    = Variable(index=[time, region], unit="\$/person")  # Undiscounted total damages per capita per year   
+    total_damages_peryear           = Variable(index=[time, region], unit="\$million")  # Undiscounted total damages per year 
+    total_damages_aggregated        = Variable(index=[time, region], unit="\$million")  # Undiscounted total damages per period
     
-    total_abatement_percap_peryear  = Variable(index = [time, region], unit = "\$/person")  # Undiscounted total abatement per capita per year     
-    total_abatement_peryear         = Variable(index = [time, region], unit = "\$million")  # Undiscounted total abatement per year 
-    total_abatement_aggregated      = Variable(index = [time, region], unit = "\$million")  # Undiscounted total abatement per period 
+    total_abatement_percap_peryear  = Variable(index=[time, region], unit="\$/person")  # Undiscounted total abatement per capita per year     
+    total_abatement_peryear         = Variable(index=[time, region], unit="\$million")  # Undiscounted total abatement per year 
+    total_abatement_aggregated      = Variable(index=[time, region], unit="\$million")  # Undiscounted total abatement per period 
     
-    total_costs_percap_peryear      = Variable(index = [time, region], unit = "\$/person")  # Undiscounted total costs per capita per year (damages + abatement)
-    total_costs_peryear             = Variable(index = [time, region], unit = "\$million")  # Undiscounted total costs per year                  
-    total_costs_aggregated          = Variable(index = [time, region], unit = "\$million")  # Undiscounted total costs per period                  
+    total_costs_percap_peryear      = Variable(index=[time, region], unit="\$/person")  # Undiscounted total costs per capita per year (damages + abatement)
+    total_costs_peryear             = Variable(index=[time, region], unit="\$million")  # Undiscounted total costs per year                  
+    total_costs_aggregated          = Variable(index=[time, region], unit="\$million")  # Undiscounted total costs per period                  
 
     function run_timestep(p, v, d, t)
 
diff --git a/src/compute_scc.jl b/src/compute_scc.jl
index 1658424a..8e063841 100644
--- a/src/compute_scc.jl
+++ b/src/compute_scc.jl
@@ -12,7 +12,7 @@ end
 function getperiodlength(year)      # same calculations made for yagg_periodspan in the model
     i = getpageindexfromyear(year)
 
-    if year==page_years[1]
+    if year == page_years[1]
         start_year = page_year_0
     else
         start_year = page_years[i - 1]
@@ -41,10 +41,10 @@ function undiscount_scc(m::Model, year::Int)
 end
 
 @defcomp ExtraEmissions begin
-    e_globalCO2emissions = Parameter(index=[time],unit="Mtonne/year")
+    e_globalCO2emissions = Parameter(index=[time], unit="Mtonne/year")
     pulse_size = Parameter(unit="Mtonne CO2")
     pulse_year = Parameter()
-    e_globalCO2emissions_adjusted = Variable(index=[time],unit="Mtonne/year")
+    e_globalCO2emissions_adjusted = Variable(index=[time], unit="Mtonne/year")
 
     function run_timestep(p, v, d, t)
         if gettime(t) == p.pulse_year
@@ -84,15 +84,15 @@ Optionally providing a `seed` value will set the random seed before running the
 results to be replicated.
 """
 function compute_scc(
-        m::Model = get_model();
-        year::Union{Int, Nothing} = nothing,
-        eta::Union{Float64, Nothing} = nothing,
-        prtp::Union{Float64, Nothing} = nothing,
-        equity_weighting::Bool = true,
-        pulse_size = 100_000.,
+        m::Model=get_model();
+        year::Union{Int,Nothing}=nothing,
+        eta::Union{Float64,Nothing}=nothing,
+        prtp::Union{Float64,Nothing}=nothing,
+        equity_weighting::Bool=true,
+        pulse_size=100_000.,
         n::Union{Int,Nothing}=nothing,
-        trials_output_filename::Union{String, Nothing} = nothing,
-        seed::Union{Int, Nothing} = nothing
+        trials_output_filename::Union{String,Nothing}=nothing,
+        seed::Union{Int,Nothing}=nothing
         )
 
     year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
@@ -124,11 +124,11 @@ function compute_scc(
     
     mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
 
-    if n===nothing
+    if n === nothing
         # Run the "best guess" social cost calculation
         run(mm)
         scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
-    elseif n<1
+    elseif n < 1
         error("Invalid `n` value, only values >=1 allowed.")
     else
         # Run a Monte Carlo simulation
@@ -139,7 +139,7 @@ function compute_scc(
        
         seed !== nothing ? Random.seed!(seed) : nothing
         Mimi.set_payload!(simdef, (Vector{Float64}(undef, n), year))  # pass the year and a vector for storing SCC values to the `run` function
-        si = run(simdef, mm, n, trials_output_filename = trials_output_filename, post_trial_func = _scc_post_trial_func)
+        si = run(simdef, mm, n, trials_output_filename=trials_output_filename, post_trial_func=_scc_post_trial_func)
         scc = Mimi.payload(si)[1]   # collect the values computed during the post-trial function        
     end
 
@@ -174,7 +174,7 @@ If no model is provided, the default model from MimiPAGE2009.get_model() is used
 Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model. 
 If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.    
 """
-function compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing, pulse_size = 100000.)
+function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=100000.)
     year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
     !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing 
 
@@ -197,7 +197,7 @@ function compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = noth
     mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
     scc = mm[:EquityWeighting, :td_totaldiscountedimpacts] / undiscount_scc(mm.base, year)
 
-    return (scc = scc, mm = mm)
+return (scc = scc, mm = mm)
 end
 
 """
@@ -207,7 +207,7 @@ Returns a Mimi MarginalModel where the provided m is the base model, and the mar
 If no Model m is provided, the default model from MimiPAGE2009.get_model() is used as the base model.
 Note that the returned MarginalModel has already been run.
 """
-function get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing, pulse_size = 100000.)
+function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing, pulse_size=100000.)
     year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
     !(year in page_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $page_years.") : nothing
     
diff --git a/src/contrib/taxeffect.jl b/src/contrib/taxeffect.jl
index cbf844dc..3300e4a2 100644
--- a/src/contrib/taxeffect.jl
+++ b/src/contrib/taxeffect.jl
@@ -8,18 +8,18 @@
     taxrate = Parameter(index=[time], unit="\$/tonne")
 
     # From external files
-    e0_baselineemissions = Parameter(index=[region], unit= "Mtonne/year")
+    e0_baselineemissions = Parameter(index=[region], unit="Mtonne/year")
 
     # From the AbatementCostParameters
-    zc_zerocostemissions = Parameter(index=[time, region], unit= "%")
-    q0_absolutecutbacksatnegativecost = Parameter(index=[time, region], unit= "Mtonne")
-    blo = Parameter(index=[time, region], unit = "per Mtonne")
-    alo = Parameter(index=[time, region], unit = "\$/tonne")
-    bhi = Parameter(index=[time, region], unit = "per Mtonne")
-    ahi = Parameter(index=[time, region], unit = "\$/tonne")
+    zc_zerocostemissions = Parameter(index=[time, region], unit="%")
+    q0_absolutecutbacksatnegativecost = Parameter(index=[time, region], unit="Mtonne")
+    blo = Parameter(index=[time, region], unit="per Mtonne")
+    alo = Parameter(index=[time, region], unit="\$/tonne")
+    bhi = Parameter(index=[time, region], unit="per Mtonne")
+    ahi = Parameter(index=[time, region], unit="\$/tonne")
 
     # Outputs to AbatementCosts
-    er_emissionsgrowth = Variable(index=[time,region],unit="%")
+    er_emissionsgrowth = Variable(index=[time,region], unit="%")
 
     function run_timestep(p, v, d, t)
 
@@ -29,16 +29,16 @@
 
             # If cbe_absoluteemissionreductions < q0_absolutecutbacksatnegativecost
             if p.alo[t, rr] + 1 > 0
-                er_lowside = p.zc_zerocostemissions[t, rr] - (log(p.taxrate[t] / p.alo[t, rr] + 1) / p.blo[t, rr] + p.q0_absolutecutbacksatnegativecost[t, rr]) / (p.e0_baselineemissions[rr]/100)
-                cbe_lowside = (p.zc_zerocostemissions[t, rr] - er_lowside) * p.e0_baselineemissions[rr]/100
+                er_lowside = p.zc_zerocostemissions[t, rr] - (log(p.taxrate[t] / p.alo[t, rr] + 1) / p.blo[t, rr] + p.q0_absolutecutbacksatnegativecost[t, rr]) / (p.e0_baselineemissions[rr] / 100)
+                cbe_lowside = (p.zc_zerocostemissions[t, rr] - er_lowside) * p.e0_baselineemissions[rr] / 100
             else
                 cbe_lowside = Inf
             end
 
             # Else
             if p.ahi[t, rr] + 1 > 0
-                er_highside = p.zc_zerocostemissions[t, rr] - (log(p.taxrate[t] / p.ahi[t, rr] + 1) / p.bhi[t, rr] + p.q0_absolutecutbacksatnegativecost[t, rr]) / (p.e0_baselineemissions[rr]/100)
-                cbe_highside = (p.zc_zerocostemissions[t, rr] - er_highside) * p.e0_baselineemissions[rr]/100
+                er_highside = p.zc_zerocostemissions[t, rr] - (log(p.taxrate[t] / p.ahi[t, rr] + 1) / p.bhi[t, rr] + p.q0_absolutecutbacksatnegativecost[t, rr]) / (p.e0_baselineemissions[rr] / 100)
+                cbe_highside = (p.zc_zerocostemissions[t, rr] - er_highside) * p.e0_baselineemissions[rr] / 100
             else
                 cbe_highside = -Inf
             end
@@ -65,14 +65,14 @@ function addtaxdrivengrowth(model::Model, class::Symbol)
     elseif class == :N2O
         setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model, "data/e0_baselineN2Oemissions.csv"))
     elseif class == :Lin
-        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model,"data/e0_baselineLGemissions.csv"))
+        setdistinctparameter(model, componentname, :e0_baselineemissions, readpagedata(model, "data/e0_baselineLGemissions.csv"))
     else
         error("Unknown class of abatement costs.")
     end
 end
 
 @defcomp UniformTaxDrivenGrowth begin
-    region=Index()
+    region = Index()
 
     uniformtax = Parameter(index=[time], unit="\$/tonne")
 
@@ -93,7 +93,7 @@ end
 """Construct a model with a uniform (global and all gases, but time-varying) tax."""
 function getuniformtaxmodel()
     m = Model()
-    set_dimension!(m, :time, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
+set_dimension!(m, :time, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
     set_dimension!(m, :region, ["EU", "USA", "OECD","USSR","China","SEAsia","Africa","LatAmerica"])
 
     buildpage(m)
@@ -106,7 +106,7 @@ function getuniformtaxmodel()
         taxgrowth = Symbol("TaxDrivenGrowth$class")
         abateparams = Symbol("AbatementCostParameters$class")
         connect_param!(m, taxgrowth => :zc_zerocostemissions, abateparams => :zc_zerocostemissions)
-        connect_param!(m, taxgrowth => :q0_absolutecutbacksatnegativecost, abateparams =>:q0_absolutecutbacksatnegativecost)
+        connect_param!(m, taxgrowth => :q0_absolutecutbacksatnegativecost, abateparams => :q0_absolutecutbacksatnegativecost)
         connect_param!(m, taxgrowth => :blo, abateparams => :blo)
         connect_param!(m, taxgrowth => :alo, abateparams => :alo)
         connect_param!(m, taxgrowth => :bhi, abateparams => :bhi)
diff --git a/src/mcs.jl b/src/mcs.jl
index 37f64da0..6c1d789e 100644
--- a/src/mcs.jl
+++ b/src/mcs.jl
@@ -3,7 +3,7 @@ using CSVFiles
 using DataFrames
 
 function getsim()
-        mcs = @defsim begin
+    mcs = @defsim begin
                 
                 ############################################################################
                 # Define random variables (RVs) - for UNSHARED parameters
@@ -12,76 +12,76 @@ function getsim()
                 # each component should have the same value for its save_savingsrate,
                 # so we use an RV here because in the model this is not an explicitly
                 # shared parameter, then assign to components
-                rv(RV_save_savingsrate) = TriangularDist(10, 20, 15)
-                GDP.save_savingsrate = RV_save_savingsrate
-                MarketDamages.save_savingsrate = RV_save_savingsrate
-                NonMarketDamages.save_savingsrate = RV_save_savingsrate
-                SLRDamages.save_savingsrate = RV_save_savingsrate
+        rv(RV_save_savingsrate) = TriangularDist(10, 20, 15)
+        GDP.save_savingsrate = RV_save_savingsrate
+        MarketDamages.save_savingsrate = RV_save_savingsrate
+        NonMarketDamages.save_savingsrate = RV_save_savingsrate
+        SLRDamages.save_savingsrate = RV_save_savingsrate
 
                 # each component should have the same value for its tcal_CalibrationTemp
                 # so we use an RV here because in the model this is not an explicitly
                 # shared parameter, then assign to components
-                rv(RV_tcal_CalibrationTemp) = TriangularDist(2.5, 3.5, 3.)
-                MarketDamages.tcal_CalibrationTemp = RV_tcal_CalibrationTemp
-                NonMarketDamages.tcal_CalibrationTemp = RV_tcal_CalibrationTemp
+        rv(RV_tcal_CalibrationTemp) = TriangularDist(2.5, 3.5, 3.)
+        MarketDamages.tcal_CalibrationTemp = RV_tcal_CalibrationTemp
+        NonMarketDamages.tcal_CalibrationTemp = RV_tcal_CalibrationTemp
 
                 # CO2cycle
-                co2cycle.air_CO2fractioninatm = TriangularDist(57, 67, 62)
-                co2cycle.res_CO2atmlifetime = TriangularDist(50, 100, 70)
-                co2cycle.ccf_CO2feedback = TriangularDist(4, 15, 10)
-                co2cycle.ccfmax_maxCO2feedback = TriangularDist(30, 80, 50)
-                co2cycle.stay_fractionCO2emissionsinatm = TriangularDist(0.25,0.35,0.3)
+        co2cycle.air_CO2fractioninatm = TriangularDist(57, 67, 62)
+        co2cycle.res_CO2atmlifetime = TriangularDist(50, 100, 70)
+        co2cycle.ccf_CO2feedback = TriangularDist(4, 15, 10)
+        co2cycle.ccfmax_maxCO2feedback = TriangularDist(30, 80, 50)
+        co2cycle.stay_fractionCO2emissionsinatm = TriangularDist(0.25, 0.35, 0.3)
                 
                 # SulphateForcing
-                SulphateForcing.d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
-                SulphateForcing.ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)
+        SulphateForcing.d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
+        SulphateForcing.ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)
                 
                 # ClimateTemperature
-                ClimateTemperature.rlo_ratiolandocean = TriangularDist(1.2, 1.6, 1.4)
-                ClimateTemperature.pole_polardifference = TriangularDist(1, 2, 1.5)
-                ClimateTemperature.frt_warminghalflife = TriangularDist(10, 65, 30)
-                ClimateTemperature.tcr_transientresponse = TriangularDist(1, 2.8, 1.3)
+        ClimateTemperature.rlo_ratiolandocean = TriangularDist(1.2, 1.6, 1.4)
+        ClimateTemperature.pole_polardifference = TriangularDist(1, 2, 1.5)
+        ClimateTemperature.frt_warminghalflife = TriangularDist(10, 65, 30)
+        ClimateTemperature.tcr_transientresponse = TriangularDist(1, 2.8, 1.3)
                 
                 # SeaLevelRise
-                SeaLevelRise.s0_initialSL = TriangularDist(0.1, 0.2, 0.15)
-                SeaLevelRise.sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)
-                SeaLevelRise.sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)
-                SeaLevelRise.sltau_SLresponsetime = TriangularDist(500, 1500, 1000)
+        SeaLevelRise.s0_initialSL = TriangularDist(0.1, 0.2, 0.15)
+        SeaLevelRise.sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)
+        SeaLevelRise.sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)
+        SeaLevelRise.sltau_SLresponsetime = TriangularDist(500, 1500, 1000)
                 
                 # GDP
-                GDP.isat0_initialimpactfxnsaturation = TriangularDist(20, 50, 30) 
+        GDP.isat0_initialimpactfxnsaturation = TriangularDist(20, 50, 30) 
                 
                 # MarketDamages
-                MarketDamages.iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
-                MarketDamages.W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
-                MarketDamages.pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
-                MarketDamages.ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)
+        MarketDamages.iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
+        MarketDamages.W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
+        MarketDamages.pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
+        MarketDamages.ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)
 
                 # NonMarketDamages
-                NonMarketDamages.iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
-                NonMarketDamages.w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
-                NonMarketDamages.pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
-                NonMarketDamages.ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)
+        NonMarketDamages.iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
+        NonMarketDamages.w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
+        NonMarketDamages.pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
+        NonMarketDamages.ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)
                 
                 # SLRDamages
-                SLRDamages.scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
-                #SLRDamages.iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
-                SLRDamages.W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
-                SLRDamages.pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
-                SLRDamages.ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)
+        SLRDamages.scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
+                # SLRDamages.iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
+        SLRDamages.W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
+        SLRDamages.pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
+        SLRDamages.ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)
                 
                 # Discontinuity
-                Discontinuity.rand_discontinuity = Uniform(0, 1)
-                Discontinuity.tdis_tolerabilitydisc = TriangularDist(2, 4, 3)
-                Discontinuity.pdis_probability = TriangularDist(10, 30, 20)
-                Discontinuity.wdis_gdplostdisc = TriangularDist(5, 25, 15)
-                Discontinuity.ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
-                Discontinuity.distau_discontinuityexponent = TriangularDist(20, 200, 50)
+        Discontinuity.rand_discontinuity = Uniform(0, 1)
+        Discontinuity.tdis_tolerabilitydisc = TriangularDist(2, 4, 3)
+        Discontinuity.pdis_probability = TriangularDist(10, 30, 20)
+        Discontinuity.wdis_gdplostdisc = TriangularDist(5, 25, 15)
+        Discontinuity.ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
+        Discontinuity.distau_discontinuityexponent = TriangularDist(20, 200, 50)
                 
                 # EquityWeighting
-                EquityWeighting.civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
-                EquityWeighting.ptp_timepreference = TriangularDist(0.1,2,1)
-                EquityWeighting.emuc_utilityconvexity = TriangularDist(0.5,2,1)
+        EquityWeighting.civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
+        EquityWeighting.ptp_timepreference = TriangularDist(0.1, 2, 1)
+        EquityWeighting.emuc_utilityconvexity = TriangularDist(0.5, 2, 1)
                 
                 ############################################################################
                 # Define random variables (RVs) - for SHARED parameters
@@ -89,109 +89,109 @@ function getsim()
 
                 # shared parameter linked to components: MarketDamages, NonMarketDamages, 
                 # SLRDamages, Discountinuity
-                wincf_weightsfactor["USA"] = TriangularDist(.6, 1, .8) 
-                wincf_weightsfactor["OECD"] = TriangularDist(.4, 1.2, .8)
-                wincf_weightsfactor["USSR"] = TriangularDist(.2, .6, .4)
-                wincf_weightsfactor["China"] = TriangularDist(.4, 1.2, .8)
-                wincf_weightsfactor["SEAsia"] = TriangularDist(.4, 1.2, .8)
-                wincf_weightsfactor["Africa"] = TriangularDist(.4, .8, .6)
-                wincf_weightsfactor["LatAmerica"] = TriangularDist(.4, .8, .6)
+        wincf_weightsfactor["USA"] = TriangularDist(.6, 1, .8) 
+        wincf_weightsfactor["OECD"] = TriangularDist(.4, 1.2, .8)
+        wincf_weightsfactor["USSR"] = TriangularDist(.2, .6, .4)
+        wincf_weightsfactor["China"] = TriangularDist(.4, 1.2, .8)
+        wincf_weightsfactor["SEAsia"] = TriangularDist(.4, 1.2, .8)
+        wincf_weightsfactor["Africa"] = TriangularDist(.4, .8, .6)
+        wincf_weightsfactor["LatAmerica"] = TriangularDist(.4, .8, .6)
 
                 # shared parameter linked to components: AdaptationCosts, AbatementCosts
-                automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  
+        automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  
 
-                #the following variables need to be set, but set the same in all 4 abatement cost components
-                #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables
+                # the following variables need to be set, but set the same in all 4 abatement cost components
+                # note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables
                 
-                emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8,1.2,1.0)
-                emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8,1.2,1.0)
-                emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65,1.35,1.0)
-                emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5,1.5,1.0)
-                emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5,1.5,1.0)
-                emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5,1.5,1.0)
-                emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5,1.5,1.0)
+        emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8, 1.2, 1.0)
+        emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8, 1.2, 1.0)
+        emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65, 1.35, 1.0)
+        emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5, 1.5, 1.0)
+        emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5, 1.5, 1.0)
+        emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5, 1.5, 1.0)
+        emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5, 1.5, 1.0)
                 
-                q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75,1.5,1.0)
-                q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75,1.25,1.0)
-                q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4,1.0,0.7)      
-                q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4,1.0,0.7)
-                q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4,1.0,0.7)
-                q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4,1.0,0.7)
-                q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)
+        q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75, 1.5, 1.0)
+        q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75, 1.25, 1.0)
+        q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4, 1.0, 0.7)      
+        q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4, 1.0, 0.7)
+        q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4, 1.0, 0.7)
+        q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4, 1.0, 0.7)
+        q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4, 1.0, 0.7)
                 
-                cmaxf_maxcostfactor["USA"] = TriangularDist(0.8,1.2,1.0)
-                cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0,1.5,1.2)
-                cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4,1.0,0.7)
-                cmaxf_maxcostfactor["China"] = TriangularDist(0.8,1.2,1.0)
-                cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1,1.5,1.2)
-                cmaxf_maxcostfactor["Africa"] = TriangularDist(1,1.5,1.2)
-                cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)
+        cmaxf_maxcostfactor["USA"] = TriangularDist(0.8, 1.2, 1.0)
+        cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0, 1.5, 1.2)
+        cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4, 1.0, 0.7)
+        cmaxf_maxcostfactor["China"] = TriangularDist(0.8, 1.2, 1.0)
+        cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1, 1.5, 1.2)
+        cmaxf_maxcostfactor["Africa"] = TriangularDist(1, 1.5, 1.2)
+        cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4, 1.0, 0.7)
                 
-                cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
-                cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
-                cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
-                cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
-                cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
-                cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
-                cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)
+        cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
+        cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
+        cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
+        cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
+        cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
+        cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
+        cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)
 
                 # Other
-                q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3,1.2,0.7)
-                qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1,1.5,1.3)
-                c0mult_mostnegativecostinfinalyear = TriangularDist(0.5,1.2,0.8)
-                curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25,0.8,0.45)
-                curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1,0.7,0.4)
-                cross_experiencecrossoverratio = TriangularDist(0.1,0.3,0.2)
-                learn_learningrate = TriangularDist(0.05,0.35,0.2)
+        q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3, 1.2, 0.7)
+        qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1, 1.5, 1.3)
+        c0mult_mostnegativecostinfinalyear = TriangularDist(0.5, 1.2, 0.8)
+        curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25, 0.8, 0.45)
+        curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1, 0.7, 0.4)
+        cross_experiencecrossoverratio = TriangularDist(0.1, 0.3, 0.2)
+        learn_learningrate = TriangularDist(0.05, 0.35, 0.2)
 
                 # NOTE: the below can probably be resolved into unique, unshared parameters with the same name
                 # in the new Mimi paradigm of shared and unshared parameters, but for now this will 
                 # continue to work!
 
                 # AbatementCosts
-                AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,75,0)
-                AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-25,100,0)
-                AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)
-                AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)
+        AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50, 75, 0)
+        AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-25, 100, 0)
+        AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50, 50, 0)
+        AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50, 50, 0)
                 
-                AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,40,20)
-                AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
-                AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
-                AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
+        AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 40, 20)
+        AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)
+        AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)
+        AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0, 20, 10)
                 
-                AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
-                AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000,-1000,-4000)
-                AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000,0,-7000)
-                AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
+        AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400, -100, -200)
+        AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000, -1000, -4000)
+        AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000, 0, -7000)
+        AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400, -100, -200)
                 
-                AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
-                AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
-                AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
-                AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
+        AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60, 80, 70)
+        AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35, 70, 50)
+        AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35, 70, 50)
+        AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60, 80, 70)
                 
-                AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,700,400)
-                AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000,10000,6000)
-                AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000,60000,20000)
-                AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,600,300)
+        AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100, 700, 400)
+        AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000, 10000, 6000)
+        AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000, 60000, 20000)
+        AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100, 600, 300)
                 
-                AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000,200000,150000)
-                AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
-                AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30,80,50)
-                AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
+        AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000, 200000, 150000)
+        AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500, 2500, 2000)
+        AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30, 80, 50)
+        AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500, 2500, 2000)
 
                 # AdaptationCosts
-                AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
-                AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
-                AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
-                AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
-                AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
-                AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)
+        AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
+        AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
+        AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
+        AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
+        AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
+        AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)
 
                 ############################################################################
                 # Indicate which parameters to save for each model run
                 ############################################################################
                 
-                save(EquityWeighting.td_totaldiscountedimpacts,
+        save(EquityWeighting.td_totaldiscountedimpacts,
                         EquityWeighting.tpc_totalaggregatedcosts, 
                         EquityWeighting.tac_totaladaptationcosts,
                         EquityWeighting.te_totaleffect,
@@ -204,52 +204,52 @@ function getsim()
                         NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                         Discontinuity.rgdp_per_cap_NonMarketRemainGDP)
 
-                end #def
+    end # def
         return mcs 
 end
 
-#Reformat the RV results into the format used for testing
-function reformat_RV_outputs(samplesize::Int; output_path::String = joinpath(@__DIR__, "../output"))         
+# Reformat the RV results into the format used for testing
+function reformat_RV_outputs(samplesize::Int; output_path::String=joinpath(@__DIR__, "../output"))         
 
-    #create vectors to hold results of Monte Carlo runs
-    td=zeros(samplesize);
-    tpc=zeros(samplesize);
-    tac=zeros(samplesize);
-    te=zeros(samplesize);
-    ft=zeros(samplesize);
-    rt_g=zeros(samplesize);
-    s=zeros(samplesize);
-    c_co2concentration=zeros(samplesize);
-    rgdppercap_slr=zeros(samplesize);
-    rgdppercap_market=zeros(samplesize);
-    rgdppercap_nonmarket=zeros(samplesize);
-    rgdppercap_disc=zeros(samplesize);
+    # create vectors to hold results of Monte Carlo runs
+    td = zeros(samplesize);
+    tpc = zeros(samplesize);
+    tac = zeros(samplesize);
+    te = zeros(samplesize);
+    ft = zeros(samplesize);
+    rt_g = zeros(samplesize);
+    s = zeros(samplesize);
+    c_co2concentration = zeros(samplesize);
+    rgdppercap_slr = zeros(samplesize);
+    rgdppercap_market = zeros(samplesize);
+    rgdppercap_nonmarket = zeros(samplesize);
+    rgdppercap_disc = zeros(samplesize);
 
-    #load raw data
-    #no filter
-    td      = load_RV("EquityWeighting_td_totaldiscountedimpacts", "td_totaldiscountedimpacts"; output_path = output_path)
-    tpc     = load_RV("EquityWeighting_tpc_totalaggregatedcosts", "tpc_totalaggregatedcosts"; output_path = output_path)
-    tac     = load_RV("EquityWeighting_tac_totaladaptationcosts", "tac_totaladaptationcosts"; output_path = output_path)
-    te      = load_RV("EquityWeighting_te_totaleffect", "te_totaleffect"; output_path = output_path)
+    # load raw data
+    # no filter
+    td      = load_RV("EquityWeighting_td_totaldiscountedimpacts", "td_totaldiscountedimpacts"; output_path=output_path)
+    tpc     = load_RV("EquityWeighting_tpc_totalaggregatedcosts", "tpc_totalaggregatedcosts"; output_path=output_path)
+    tac     = load_RV("EquityWeighting_tac_totaladaptationcosts", "tac_totaladaptationcosts"; output_path=output_path)
+    te      = load_RV("EquityWeighting_te_totaleffect", "te_totaleffect"; output_path=output_path)
 
-    #time index
-    c_co2concentration = load_RV("co2cycle_c_CO2concentration", "c_CO2concentration"; output_path = output_path)
-    ft      = load_RV("TotalForcing_ft_totalforcing", "ft_totalforcing"; output_path = output_path)
-    rt_g    = load_RV("ClimateTemperature_rt_g_globaltemperature", "rt_g_globaltemperature"; output_path = output_path)
-    s       = load_RV("SeaLevelRise_s_sealevel", "s_sealevel"; output_path = output_path)
+    # time index
+    c_co2concentration = load_RV("co2cycle_c_CO2concentration", "c_CO2concentration"; output_path=output_path)
+    ft      = load_RV("TotalForcing_ft_totalforcing", "ft_totalforcing"; output_path=output_path)
+    rt_g    = load_RV("ClimateTemperature_rt_g_globaltemperature", "rt_g_globaltemperature"; output_path=output_path)
+    s       = load_RV("SeaLevelRise_s_sealevel", "s_sealevel"; output_path=output_path)
 
-    #region index
-    rgdppercap_slr          = load_RV("SLRDamages_rgdp_per_cap_SLRRemainGDP", "rgdp_per_cap_SLRRemainGDP"; output_path = output_path)
-    rgdppercap_market       = load_RV("MarketDamages_rgdp_per_cap_MarketRemainGDP", "rgdp_per_cap_MarketRemainGDP"; output_path = output_path)
-    rgdppercap_nonmarket    =load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path = output_path)
-    rgdppercap_disc         = load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path = output_path)
+    # region index
+    rgdppercap_slr          = load_RV("SLRDamages_rgdp_per_cap_SLRRemainGDP", "rgdp_per_cap_SLRRemainGDP"; output_path=output_path)
+    rgdppercap_market       = load_RV("MarketDamages_rgdp_per_cap_MarketRemainGDP", "rgdp_per_cap_MarketRemainGDP"; output_path=output_path)
+    rgdppercap_nonmarket    = load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path=output_path)
+    rgdppercap_disc         = load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path=output_path)
 
-    #resave data
-    df=DataFrame(td=td,tpc=tpc,tac=tac,te=te,c_co2concentration=c_co2concentration,ft=ft,rt_g=rt_g,sealevel=s,rgdppercap_slr=rgdppercap_slr,rgdppercap_market=rgdppercap_market,rgdppercap_nonmarket=rgdppercap_nonmarket,rgdppercap_di=rgdppercap_disc)
-    save(joinpath(output_path, "mimipagemontecarlooutput.csv"),df)
+    # resave data
+    df = DataFrame(td=td, tpc=tpc, tac=tac, te=te, c_co2concentration=c_co2concentration, ft=ft, rt_g=rt_g, sealevel=s, rgdppercap_slr=rgdppercap_slr, rgdppercap_market=rgdppercap_market, rgdppercap_nonmarket=rgdppercap_nonmarket, rgdppercap_di=rgdppercap_disc)
+    save(joinpath(output_path, "mimipagemontecarlooutput.csv"), df)
 end
 
-function do_monte_carlo_runs(samplesize::Int; output_dir::String = joinpath(@__DIR__, "../output"))
+function do_monte_carlo_runs(samplesize::Int; output_dir::String=joinpath(@__DIR__, "../output"))
 
         # get simulation
         mcs = getsim()
@@ -259,8 +259,8 @@ function do_monte_carlo_runs(samplesize::Int; output_dir::String = joinpath(@__D
         run(m)
 
         # Run
-        res = run(mcs, m, samplesize; trials_output_filename = joinpath(output_dir, "trialdata.csv"), results_output_dir = output_dir)
+        res = run(mcs, m, samplesize; trials_output_filename=joinpath(output_dir, "trialdata.csv"), results_output_dir=output_dir)
 
         # reformat outputs for testing and analysis
-        MimiPAGE2009.reformat_RV_outputs(samplesize, output_path = output_dir)
+        MimiPAGE2009.reformat_RV_outputs(samplesize, output_path=output_dir)
 end
diff --git a/src/utils/load_parameters.jl b/src/utils/load_parameters.jl
index 765fc0d0..8c174025 100644
--- a/src/utils/load_parameters.jl
+++ b/src/utils/load_parameters.jl
@@ -1,7 +1,7 @@
 using DelimitedFiles
 
 function checkregionorder(model::Model, regions, file)
-    regionaliases = Dict{AbstractString, Vector{AbstractString}}("EU" => [],
+    regionaliases = Dict{AbstractString,Vector{AbstractString}}("EU" => [],
                                                                  "USA" => ["US"],
                                                                  "OECD" => ["OT"],
                                                                  "Africa" => ["AF"],
@@ -41,7 +41,7 @@ function readpagedata(model::Model, filepath::AbstractString)
         # Check that regions are in the right order
         checkregionorder(model, data[1][:, 1], basename(filepath))
 
-        return convert(Vector{Float64},vec(data[1][:, 2]))
+        return convert(Vector{Float64}, vec(data[1][:, 2]))
     elseif firstline == "# Index: time"
         data = readdlm(filepath, ',', header=true, comments=true)
 
@@ -50,7 +50,7 @@ function readpagedata(model::Model, filepath::AbstractString)
 
         return convert(Vector{Float64}, vec(data[1][:, 2]))
     elseif firstline == "# Index: time, region"
-        data = readdlm(filepath, ',', header=true, comments = true)
+        data = readdlm(filepath, ',', header=true, comments=true)
 
         # Check that both dimension match
         checktimeorder(model, data[1][:, 1], basename(filepath))
@@ -63,10 +63,10 @@ function readpagedata(model::Model, filepath::AbstractString)
 end
 
 function load_parameters(model::Model, policy::String="policy-a")
-    parameters = Dict{Any, Any}()
+    parameters = Dict{Any,Any}()
 
     parameter_directory = joinpath(dirname(@__FILE__), "..", "..", "data")
-    for file in filter(q->splitext(q)[2]==".csv", readdir(parameter_directory))
+    for file in filter(q -> splitext(q)[2] == ".csv", readdir(parameter_directory))
         parametername = splitext(file)[1]
 
         if policy != "policy-a" && isfile(joinpath(parameter_directory, policy, file))
diff --git a/src/utils/mctools.jl b/src/utils/mctools.jl
index aced08a7..7a820282 100644
--- a/src/utils/mctools.jl
+++ b/src/utils/mctools.jl
@@ -13,15 +13,15 @@ end
 Load raw RV output into reformat_RV_outputs 
 """
 function load_RV(filename::String, RVname::String; 
-                    output_path::String = joinpath(@__DIR__, "../../output/"), 
-                    time_filter::Int = 2200,
-                    region_filter::String = "LatAmerica")
+                    output_path::String=joinpath(@__DIR__, "../../output/"), 
+                    time_filter::Int=2200,
+                    region_filter::String="LatAmerica")
     
     df = DataFrame(load(joinpath(output_path, "$filename.csv")))
     cols = names(df)
     
-    #apply filters if necessary, currently the function supports a time filter 
-    #of a single time value and a region filter of a single region
+    # apply filters if necessary, currently the function supports a time filter 
+    # of a single time value and a region filter of a single region
     if in("time", cols)
 
         if in("region", cols) 
@@ -29,7 +29,7 @@ function load_RV(filename::String, RVname::String;
                 @where i.time == time_filter
                 @where i.region == region_filter
                 @select i
-                end) |> DataFrame
+            end) |> DataFrame
 
         else
             filtered_df = df |> @query(i, begin
@@ -40,7 +40,7 @@ function load_RV(filename::String, RVname::String;
 
         return filtered_df[!, Symbol(RVname)]
           
-    #no filters applied
+    # no filters applied
     else
         return df[!, Symbol(RVname)]
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 767e2895..58f8fa97 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,45 +11,45 @@ function test_page_model()
     set_dimension!(m, :region, ["EU", "USA", "OECD", "USSR", "China", "SEAsia", "Africa", "LatAmerica"])
     
     return m
- end
+end
 
 @testset "MimiPAGE2009.jl" begin
 
-include("test_climatemodel.jl")
-include("test_AbatementCosts.jl")
-include("test_AdaptationCosts.jl")
-include("test_CH4cycle.jl")
-include("test_CH4emissions.jl")
-include("test_CH4forcing.jl")
-include("test_ClimateTemperature.jl")
-include("test_CO2cycle.jl")
-include("test_CO2emissions.jl")
-include("test_CO2forcing.jl")
-include("test_Discontinuity.jl")
-include("test_EquityWeighting.jl")
-include("test_GDP.jl")
-include("test_LGcycle.jl")
-include("test_LGemissions.jl")
-include("test_LGforcing.jl")
-include("test_loadparameters.jl")
-include("test_mainmodel.jl") 
-include("test_mainmodel_policyb.jl") 
-include("test_MarketDamages.jl")
-include("test_N2Ocycle.jl")
-include("test_N2Oemissions.jl")
-include("test_N2Oforcing.jl")
-include("test_NonMarketDamages.jl")
-include("test_Population.jl")
-include("test_SeaLevelRise.jl")
-include("test_SLRDamages.jl")
-include("test_SulphateForcing.jl") 
-include("test_TotalAbatementCosts.jl")
-include("test_TotalAdaptationCosts.jl")
-include("test_TotalCosts.jl")
-include("test_TotalForcing.jl")
-include("test_mcs.jl")
-include("contrib/test_taxeffect.jl") 
+    include("test_climatemodel.jl")
+    include("test_AbatementCosts.jl")
+    include("test_AdaptationCosts.jl")
+    include("test_CH4cycle.jl")
+    include("test_CH4emissions.jl")
+    include("test_CH4forcing.jl")
+    include("test_ClimateTemperature.jl")
+    include("test_CO2cycle.jl")
+    include("test_CO2emissions.jl")
+    include("test_CO2forcing.jl")
+    include("test_Discontinuity.jl")
+    include("test_EquityWeighting.jl")
+    include("test_GDP.jl")
+    include("test_LGcycle.jl")
+    include("test_LGemissions.jl")
+    include("test_LGforcing.jl")
+    include("test_loadparameters.jl")
+    include("test_mainmodel.jl") 
+    include("test_mainmodel_policyb.jl") 
+    include("test_MarketDamages.jl")
+    include("test_N2Ocycle.jl")
+    include("test_N2Oemissions.jl")
+    include("test_N2Oforcing.jl")
+    include("test_NonMarketDamages.jl")
+    include("test_Population.jl")
+    include("test_SeaLevelRise.jl")
+    include("test_SLRDamages.jl")
+    include("test_SulphateForcing.jl") 
+    include("test_TotalAbatementCosts.jl")
+    include("test_TotalAdaptationCosts.jl")
+    include("test_TotalCosts.jl")
+    include("test_TotalForcing.jl")
+    include("test_mcs.jl")
+    include("contrib/test_taxeffect.jl") 
 
-include("test_standard_api.jl")
+    include("test_standard_api.jl")
 
 end
diff --git a/test/test_AbatementCosts.jl b/test/test_AbatementCosts.jl
index 0d018d08..51535e5f 100644
--- a/test/test_AbatementCosts.jl
+++ b/test/test_AbatementCosts.jl
@@ -22,7 +22,7 @@ for gas in [:CO2, :CH4, :N2O, :Lin]
 
 end
 
-set_param!(m, :yagg, readpagedata(m,"test/validationdata/yagg_periodspan.csv"))
+set_param!(m, :yagg, readpagedata(m, "test/validationdata/yagg_periodspan.csv"))
 
 p = load_parameters(m)
 p["y_year_0"] = 2008.
@@ -36,24 +36,24 @@ run(m)
 @test !isnan(m[:AbatementCostsN2O, :tc_totalcost][10, 5])
 @test !isnan(m[:AbatementCostsLin, :tc_totalcost][10, 5])
 
-#compare output to validation data
-tc_compare_co2=readpagedata(m, "test/validationdata/tc_totalcosts_co2.csv")
-tc_compare_ch4=readpagedata(m, "test/validationdata/tc_totalcosts_ch4.csv")
-tc_compare_n2o=readpagedata(m, "test/validationdata/tc_totalcosts_n2o.csv")
-tc_compare_lin=readpagedata(m, "test/validationdata/tc_totalcosts_linear.csv")
+# compare output to validation data
+tc_compare_co2 = readpagedata(m, "test/validationdata/tc_totalcosts_co2.csv")
+tc_compare_ch4 = readpagedata(m, "test/validationdata/tc_totalcosts_ch4.csv")
+tc_compare_n2o = readpagedata(m, "test/validationdata/tc_totalcosts_n2o.csv")
+tc_compare_lin = readpagedata(m, "test/validationdata/tc_totalcosts_linear.csv")
 
-zc_compare_co2=readpagedata(m, "test/validationdata/zc_zerocostemissionsCO2.csv")
-zc_compare_ch4=readpagedata(m, "test/validationdata/zc_zerocostemissionsCH4.csv")
-zc_compare_n2o=readpagedata(m, "test/validationdata/zc_zerocostemissionsN2O.csv")
-zc_compare_lin=readpagedata(m, "test/validationdata/zc_zerocostemissionsLG.csv")
+zc_compare_co2 = readpagedata(m, "test/validationdata/zc_zerocostemissionsCO2.csv")
+zc_compare_ch4 = readpagedata(m, "test/validationdata/zc_zerocostemissionsCH4.csv")
+zc_compare_n2o = readpagedata(m, "test/validationdata/zc_zerocostemissionsN2O.csv")
+zc_compare_lin = readpagedata(m, "test/validationdata/zc_zerocostemissionsLG.csv")
 
 
-@test m[:AbatementCostsCO2, :tc_totalcost] ≈ tc_compare_co2 rtol=1e-2
-@test m[:AbatementCostsCH4, :tc_totalcost] ≈ tc_compare_ch4 rtol=1e-2
-@test m[:AbatementCostsN2O, :tc_totalcost] ≈ tc_compare_n2o rtol=1e-2
-@test m[:AbatementCostsLin, :tc_totalcost] ≈ tc_compare_lin rtol=1e-2
+@test m[:AbatementCostsCO2, :tc_totalcost] ≈ tc_compare_co2 rtol = 1e-2
+@test m[:AbatementCostsCH4, :tc_totalcost] ≈ tc_compare_ch4 rtol = 1e-2
+@test m[:AbatementCostsN2O, :tc_totalcost] ≈ tc_compare_n2o rtol = 1e-2
+@test m[:AbatementCostsLin, :tc_totalcost] ≈ tc_compare_lin rtol = 1e-2
 
-@test m[:AbatementCostParametersCO2, :zc_zerocostemissions] ≈ zc_compare_co2 rtol=1e-2
-@test m[:AbatementCostParametersCH4, :zc_zerocostemissions] ≈ zc_compare_ch4 rtol=1e-3
-@test m[:AbatementCostParametersN2O, :zc_zerocostemissions] ≈ zc_compare_n2o rtol=1e-3
-@test m[:AbatementCostParametersLin, :zc_zerocostemissions] ≈ zc_compare_lin rtol=1e-3
+@test m[:AbatementCostParametersCO2, :zc_zerocostemissions] ≈ zc_compare_co2 rtol = 1e-2
+@test m[:AbatementCostParametersCH4, :zc_zerocostemissions] ≈ zc_compare_ch4 rtol = 1e-3
+@test m[:AbatementCostParametersN2O, :zc_zerocostemissions] ≈ zc_compare_n2o rtol = 1e-3
+@test m[:AbatementCostParametersLin, :zc_zerocostemissions] ≈ zc_compare_lin rtol = 1e-3
diff --git a/test/test_AdaptationCosts.jl b/test/test_AdaptationCosts.jl
index ffe55d6a..528e5bad 100644
--- a/test/test_AdaptationCosts.jl
+++ b/test/test_AdaptationCosts.jl
@@ -31,11 +31,11 @@ imp_economic_compare = readpagedata(m, "test/validationdata/imp_1.csv")
 acp_economic_compare = readpagedata(m, "test/validationdata/acp_adaptivecostplateau_economic.csv")
 aci_economic_compare = readpagedata(m, "test/validationdata/aci_adaptivecostimpact_economic.csv")
 
-@test autofac ≈ autofac_compare rtol=1e-8
-@test atl_economic ≈ atl_economic_compare rtol=1e-8
-@test imp_economic ≈ imp_economic_compare rtol=1e-8
-@test acp_economic ≈ acp_economic_compare rtol=1e-3
-@test aci_economic ≈ aci_economic_compare rtol=1e-3
+@test autofac ≈ autofac_compare rtol = 1e-8
+@test atl_economic ≈ atl_economic_compare rtol = 1e-8
+@test imp_economic ≈ imp_economic_compare rtol = 1e-8
+@test acp_economic ≈ acp_economic_compare rtol = 1e-3
+@test aci_economic ≈ aci_economic_compare rtol = 1e-3
 
 ac_noneconomic = m[:AdaptiveCostsNonEconomic, :ac_adaptivecosts]
 ac_economic = m[:AdaptiveCostsEconomic, :ac_adaptivecosts]
@@ -45,6 +45,6 @@ ac_noneconomic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts
 ac_economic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv")
 ac_sealevel_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_sealevelrise.csv")
 
-@test ac_noneconomic ≈ ac_noneconomic_compare rtol=1e-2
-@test ac_economic ≈ ac_economic_compare rtol=1e-3
-@test ac_sealevel ≈ ac_sealevel_compare rtol=1e-2
+@test ac_noneconomic ≈ ac_noneconomic_compare rtol = 1e-2
+@test ac_economic ≈ ac_economic_compare rtol = 1e-3
+@test ac_sealevel ≈ ac_sealevel_compare rtol = 1e-2
diff --git a/test/test_CH4cycle.jl b/test/test_CH4cycle.jl
index 06266cdf..a2125e09 100644
--- a/test/test_CH4cycle.jl
+++ b/test/test_CH4cycle.jl
@@ -6,18 +6,18 @@ include("../src/components/CH4cycle.jl")
 
 add_comp!(m, ch4cycle, :ch4cycle)
 
-set_param!(m, :ch4cycle, :e_globalCH4emissions, readpagedata(m,"test/validationdata/e_globalCH4emissions.csv"))
-set_param!(m, :ch4cycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))
+set_param!(m, :ch4cycle, :e_globalCH4emissions, readpagedata(m, "test/validationdata/e_globalCH4emissions.csv"))
+set_param!(m, :ch4cycle, :rtl_g_landtemperature, readpagedata(m, "test/validationdata/rtl_g_landtemperature.csv"))
 set_param!(m,:ch4cycle,:y_year_0,2008.)
 
 p = load_parameters(m)
 p["y_year"] = Mimi.dim_keys(m.md, :time)
 set_leftover_params!(m, p)
 
-#running Model
+# running Model
 run(m)
 
-conc=m[:ch4cycle,  :c_CH4concentration]
-conc_compare=readpagedata(m,"test/validationdata/c_ch4concentration.csv")
+conc = m[:ch4cycle,  :c_CH4concentration]
+conc_compare = readpagedata(m, "test/validationdata/c_ch4concentration.csv")
 
-@test conc ≈ conc_compare rtol=1e-4
+@test conc ≈ conc_compare rtol = 1e-4
diff --git a/test/test_CH4emissions.jl b/test/test_CH4emissions.jl
index d2fc4679..0dcc4439 100644
--- a/test/test_CH4emissions.jl
+++ b/test/test_CH4emissions.jl
@@ -5,16 +5,16 @@ include("../src/components/CH4emissions.jl")
 
 add_comp!(m, ch4emissions)
 
-set_param!(m, :ch4emissions, :e0_baselineCH4emissions, readpagedata(m, "data/e0_baselineCH4emissions.csv")) #PAGE 2009 documentation pp38
+set_param!(m, :ch4emissions, :e0_baselineCH4emissions, readpagedata(m, "data/e0_baselineCH4emissions.csv")) # PAGE 2009 documentation pp38
 set_param!(m, :ch4emissions, :er_CH4emissionsgrowth, readpagedata(m, "data/er_CH4emissionsgrowth.csv"))
 
 ##running Model
 run(m)
 
 # Generated data
-emissions= m[:ch4emissions,  :e_regionalCH4emissions]
+emissions = m[:ch4emissions,  :e_regionalCH4emissions]
 
 # Recorded data
-emissions_compare=readpagedata(m, "test/validationdata/e_regionalCH4emissions.csv")
+emissions_compare = readpagedata(m, "test/validationdata/e_regionalCH4emissions.csv")
 
-@test emissions ≈ emissions_compare rtol=1e-3
+@test emissions ≈ emissions_compare rtol = 1e-3
diff --git a/test/test_CH4forcing.jl b/test/test_CH4forcing.jl
index 00162a3a..8960cc9e 100644
--- a/test/test_CH4forcing.jl
+++ b/test/test_CH4forcing.jl
@@ -6,15 +6,15 @@ include("../src/components/CH4forcing.jl")
 
 add_comp!(m, ch4forcing, :ch4forcing)
 
-set_param!(m, :ch4forcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/c_n2oconcentration.csv"))
-set_param!(m, :ch4forcing, :c_CH4concentration, readpagedata(m,"test/validationdata/c_ch4concentration.csv"))
+set_param!(m, :ch4forcing, :c_N2Oconcentration, readpagedata(m, "test/validationdata/c_n2oconcentration.csv"))
+set_param!(m, :ch4forcing, :c_CH4concentration, readpagedata(m, "test/validationdata/c_ch4concentration.csv"))
 
 ##running Model
 run(m)
 
 @test !isnan(m[:ch4forcing, :f_CH4forcing][10])
 
-forcing=m[:ch4forcing,:f_CH4forcing]
-forcing_compare=readpagedata(m,"test/validationdata/f_ch4forcing.csv")
+forcing = m[:ch4forcing,:f_CH4forcing]
+forcing_compare = readpagedata(m, "test/validationdata/f_ch4forcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_CO2cycle.jl b/test/test_CO2cycle.jl
index 7895cc6c..6d74b80d 100644
--- a/test/test_CO2cycle.jl
+++ b/test/test_CO2cycle.jl
@@ -7,17 +7,17 @@ include("../src/components/CO2cycle.jl")
 add_comp!(m, co2cycle)
 
 set_param!(m, :co2cycle, :e_globalCO2emissions, readpagedata(m, "test/validationdata/e_globalCO2emissions.csv"))
-set_param!(m, :co2cycle, :y_year,[2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.])#real values
-set_param!(m, :co2cycle, :y_year_0,2008.)#real value
+set_param!(m, :co2cycle, :y_year,[2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.])# real values
+set_param!(m, :co2cycle, :y_year_0,2008.)# real value
 set_param!(m, :co2cycle, :rt_g_globaltemperature, readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv"))
 
-p=load_parameters(m)
+p = load_parameters(m)
 set_leftover_params!(m,p) 
 
 run(m)
 
-pop=m[:co2cycle,  :c_CO2concentration]
+pop = m[:co2cycle,  :c_CO2concentration]
 
-pop_compare=readpagedata(m, "test/validationdata/c_co2concentration.csv")
+pop_compare = readpagedata(m, "test/validationdata/c_co2concentration.csv")
 
-@test pop ≈ pop_compare rtol=1e-4
+@test pop ≈ pop_compare rtol = 1e-4
diff --git a/test/test_CO2emissions.jl b/test/test_CO2emissions.jl
index 35040ecc..09cede45 100644
--- a/test/test_CO2emissions.jl
+++ b/test/test_CO2emissions.jl
@@ -6,15 +6,15 @@ include("../src/components/CO2emissions.jl")
 
 add_comp!(m, co2emissions)
 
-set_param!(m, :co2emissions, :e0_baselineCO2emissions, readpagedata(m,"data/e0_baselineCO2emissions.csv"))
+set_param!(m, :co2emissions, :e0_baselineCO2emissions, readpagedata(m, "data/e0_baselineCO2emissions.csv"))
 set_param!(m, :co2emissions, :er_CO2emissionsgrowth, readpagedata(m, "data/er_CO2emissionsgrowth.csv"))
 
 ##running Model
 run(m)
 
-emissions= m[:co2emissions,  :e_regionalCO2emissions]
+emissions = m[:co2emissions,  :e_regionalCO2emissions]
 
 # Recorded data
-emissions_compare=readpagedata(m, "test/validationdata/e_regionalCO2emissions.csv")
+emissions_compare = readpagedata(m, "test/validationdata/e_regionalCO2emissions.csv")
 
-@test emissions ≈ emissions_compare rtol=1e-3
+@test emissions ≈ emissions_compare rtol = 1e-3
diff --git a/test/test_CO2forcing.jl b/test/test_CO2forcing.jl
index 3123339e..3fe144ce 100644
--- a/test/test_CO2forcing.jl
+++ b/test/test_CO2forcing.jl
@@ -6,12 +6,12 @@ include("../src/components/CO2forcing.jl")
 
 add_comp!(m, co2forcing)
 
-set_param!(m, :co2forcing, :c_CO2concentration, readpagedata(m,"test/validationdata/c_co2concentration.csv"))
+set_param!(m, :co2forcing, :c_CO2concentration, readpagedata(m, "test/validationdata/c_co2concentration.csv"))
 
 ##running Model
 run(m)
 
-forcing=m[:co2forcing,:f_CO2forcing]
-forcing_compare=readpagedata(m,"test/validationdata/f_co2forcing.csv")
+forcing = m[:co2forcing,:f_CO2forcing]
+forcing_compare = readpagedata(m, "test/validationdata/f_co2forcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_ClimateTemperature.jl b/test/test_ClimateTemperature.jl
index 91e11d13..af91cc3c 100644
--- a/test/test_ClimateTemperature.jl
+++ b/test/test_ClimateTemperature.jl
@@ -21,4 +21,4 @@ run(m)
 rt_g = m[:ClimateTemperature, :rt_g_globaltemperature]
 rt_g_compare = readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv")
 
-@test rt_g ≈ rt_g_compare rtol=1e-5
+@test rt_g ≈ rt_g_compare rtol = 1e-5
diff --git a/test/test_Discontinuity.jl b/test/test_Discontinuity.jl
index 7b9e717c..51c20565 100644
--- a/test/test_Discontinuity.jl
+++ b/test/test_Discontinuity.jl
@@ -25,9 +25,9 @@ run(m)
 @test !isnan(m[:Discontinuity, :isat_per_cap_DiscImpactperCapinclSaturation][10,8])
 @test !isnan(m[:Discontinuity, :rcons_per_cap_DiscRemainConsumption][10])
 
-#validating - comparison spreadsheet has discontinuity occuring in 2200
-#keep running model until m[:Discontinuity,:occurdis_occurrencedummy] shows discontiuity occuring in 2200
-output=m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
-validation=readpagedata(m,"test/validationdata/rcons_per_cap_DiscRemainConsumption.csv")
+# validating - comparison spreadsheet has discontinuity occuring in 2200
+# keep running model until m[:Discontinuity,:occurdis_occurrencedummy] shows discontiuity occuring in 2200
+output = m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
+validation = readpagedata(m, "test/validationdata/rcons_per_cap_DiscRemainConsumption.csv")
 
-@test output ≈ validation rtol=1e-2
+@test output ≈ validation rtol = 1e-2
diff --git a/test/test_EquityWeighting.jl b/test/test_EquityWeighting.jl
index 6ca281b2..5dfb47cd 100644
--- a/test/test_EquityWeighting.jl
+++ b/test/test_EquityWeighting.jl
@@ -55,20 +55,20 @@ addt_compare = readpagedata(m, "test/validationdata/addt_equityweightedimpact_di
 addt_gt_compare = 204132238.85242900
 te_compare = 213208136.69903600
 
-@test df ≈ df_compare rtol=1e-8
-@test wtct_percap ≈ wtct_percap_compare rtol=1e-7
-@test pct_percap ≈ pct_percap_compare rtol=1e-7
-@test dr ≈ dr_compare rtol=1e-5
-@test dfc ≈ dfc_compare rtol=1e-7
-@test pct ≈ pct_compare rtol=1e-3
-@test pcdt ≈ pcdt_compare rtol=1e-3
-@test wacdt ≈ wacdt_compare rtol=1e-4
+@test df ≈ df_compare rtol = 1e-8
+@test wtct_percap ≈ wtct_percap_compare rtol = 1e-7
+@test pct_percap ≈ pct_percap_compare rtol = 1e-7
+@test dr ≈ dr_compare rtol = 1e-5
+@test dfc ≈ dfc_compare rtol = 1e-7
+@test pct ≈ pct_compare rtol = 1e-3
+@test pcdt ≈ pcdt_compare rtol = 1e-3
+@test wacdt ≈ wacdt_compare rtol = 1e-4
 
-@test aact ≈ aact_compare rtol=1e-3
+@test aact ≈ aact_compare rtol = 1e-3
 
-@test wit ≈ wit_compare rtol=1e-3
-@test addt ≈ addt_compare rtol=1e-2
-@test addt_gt ≈ addt_gt_compare rtol=1e-2
+@test wit ≈ wit_compare rtol = 1e-3
+@test addt ≈ addt_compare rtol = 1e-2
+@test addt_gt ≈ addt_gt_compare rtol = 1e-2
 
-@test te ≈ te_compare rtol=1e-2
+@test te ≈ te_compare rtol = 1e-2
 
diff --git a/test/test_GDP.jl b/test/test_GDP.jl
index 3578b977..5972f5bd 100644
--- a/test/test_GDP.jl
+++ b/test/test_GDP.jl
@@ -11,7 +11,7 @@ gdp[:pop_population] = readpagedata(m, "test/validationdata/pop_population.csv")
 gdp[:y_year] = Mimi.dim_keys(m.md, :time)
 gdp[:y_year_0] = 2008.
 
-p=load_parameters(m)
+p = load_parameters(m)
 set_leftover_params!(m,p)
 
 # run model
@@ -23,7 +23,7 @@ gdp = m[:GDP, :gdp]
 # Recorded data
 gdp_compare = readpagedata(m, "test/validationdata/gdp.csv")
 
-@test gdp ≈ gdp_compare rtol=100
+@test gdp ≈ gdp_compare rtol = 100
 
 cons_percap_consumption_0_compare = readpagedata(m, "test/validationdata/cons_percap_consumption_0.csv")
-@test m[:GDP, :cons_percap_consumption_0] ≈ cons_percap_consumption_0_compare rtol=1e-2
+@test m[:GDP, :cons_percap_consumption_0] ≈ cons_percap_consumption_0_compare rtol = 1e-2
diff --git a/test/test_LGcycle.jl b/test/test_LGcycle.jl
index b4f31dfa..3a0c29e6 100644
--- a/test/test_LGcycle.jl
+++ b/test/test_LGcycle.jl
@@ -6,18 +6,18 @@ include("../src/components/LGcycle.jl")
 
 add_comp!(m, LGcycle)
 
-set_param!(m, :LGcycle, :e_globalLGemissions, readpagedata(m,"test/validationdata/e_globalLGemissions.csv"))
-set_param!(m, :LGcycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) #real value
-set_param!(m, :LGcycle, :y_year_0, 2008.) #real value
-set_param!(m, :LGcycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))
+set_param!(m, :LGcycle, :e_globalLGemissions, readpagedata(m, "test/validationdata/e_globalLGemissions.csv"))
+set_param!(m, :LGcycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) # real value
+set_param!(m, :LGcycle, :y_year_0, 2008.) # real value
+set_param!(m, :LGcycle, :rtl_g_landtemperature, readpagedata(m, "test/validationdata/rtl_g_landtemperature.csv"))
 
-p=load_parameters(m)
-set_leftover_params!(m,p) #important for setting left over component values
+p = load_parameters(m)
+set_leftover_params!(m,p) # important for setting left over component values
 
 # run Model
 run(m)
 
-conc=m[:LGcycle,  :c_LGconcentration]
-conc_compare=readpagedata(m,"test/validationdata/c_LGconcentration.csv")
+conc = m[:LGcycle,  :c_LGconcentration]
+conc_compare = readpagedata(m, "test/validationdata/c_LGconcentration.csv")
 
-@test conc ≈ conc_compare rtol=1e-4
+@test conc ≈ conc_compare rtol = 1e-4
diff --git a/test/test_LGemissions.jl b/test/test_LGemissions.jl
index 8372f423..697a535d 100644
--- a/test/test_LGemissions.jl
+++ b/test/test_LGemissions.jl
@@ -6,13 +6,13 @@ include("../src/components/LGemissions.jl")
 
 add_comp!(m, LGemissions)
 
-set_param!(m, :LGemissions, :e0_baselineLGemissions, readpagedata(m,"data/e0_baselineLGemissions.csv"))
+set_param!(m, :LGemissions, :e0_baselineLGemissions, readpagedata(m, "data/e0_baselineLGemissions.csv"))
 set_param!(m, :LGemissions, :er_LGemissionsgrowth, readpagedata(m, "data/er_LGemissionsgrowth.csv"))
 
 # run Model
 run(m)
 
-emissions= m[:LGemissions,  :e_regionalLGemissions]
-emissions_compare=readpagedata(m, "test/validationdata/e_regionalLGemissions.csv")
+emissions = m[:LGemissions,  :e_regionalLGemissions]
+emissions_compare = readpagedata(m, "test/validationdata/e_regionalLGemissions.csv")
 
-@test emissions ≈ emissions_compare rtol=1e-3
+@test emissions ≈ emissions_compare rtol = 1e-3
diff --git a/test/test_LGforcing.jl b/test/test_LGforcing.jl
index 3886f5db..f8334209 100644
--- a/test/test_LGforcing.jl
+++ b/test/test_LGforcing.jl
@@ -6,12 +6,12 @@ include("../src/components/LGforcing.jl")
 
 add_comp!(m, LGforcing)
 
-set_param!(m, :LGforcing, :c_LGconcentration, readpagedata(m,"test/validationdata/c_LGconcentration.csv"))
+set_param!(m, :LGforcing, :c_LGconcentration, readpagedata(m, "test/validationdata/c_LGconcentration.csv"))
 
 # run Model
 run(m)
 
-forcing=m[:LGforcing,:f_LGforcing]
-forcing_compare=readpagedata(m,"test/validationdata/f_LGforcing.csv")
+forcing = m[:LGforcing,:f_LGforcing]
+forcing_compare = readpagedata(m, "test/validationdata/f_LGforcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_MarketDamages.jl b/test/test_MarketDamages.jl
index fdf24b54..1973289d 100644
--- a/test/test_MarketDamages.jl
+++ b/test/test_MarketDamages.jl
@@ -8,10 +8,10 @@ include("../src/components/MarketDamages.jl")
 marketdamages = addmarketdamages(m)
 
 set_param!(m, :MarketDamages, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/rtl_realizedtemperature.csv"))
-set_param!(m, :MarketDamages, :rcons_per_cap_SLRRemainConsumption, readpagedata(m,"test/validationdata/rcons_per_cap_SLRRemainConsumption.csv"))
-set_param!(m, :MarketDamages, :rgdp_per_cap_SLRRemainGDP, readpagedata(m,"test/validationdata/rgdp_per_cap_SLRRemainGDP.csv"))
-set_param!(m, :MarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m,"test/validationdata/atl_adjustedtolerableleveloftemprise_market.csv"))
-set_param!(m, :MarketDamages, :imp_actualreduction, readpagedata(m,"test/validationdata/imp_actualreduction_market.csv"))
+set_param!(m, :MarketDamages, :rcons_per_cap_SLRRemainConsumption, readpagedata(m, "test/validationdata/rcons_per_cap_SLRRemainConsumption.csv"))
+set_param!(m, :MarketDamages, :rgdp_per_cap_SLRRemainGDP, readpagedata(m, "test/validationdata/rgdp_per_cap_SLRRemainGDP.csv"))
+set_param!(m, :MarketDamages, :atl_adjustedtolerableleveloftemprise, readpagedata(m, "test/validationdata/atl_adjustedtolerableleveloftemprise_market.csv"))
+set_param!(m, :MarketDamages, :imp_actualreduction, readpagedata(m, "test/validationdata/imp_actualreduction_market.csv"))
 set_param!(m, :MarketDamages, :isatg_impactfxnsaturation, 28.333333333333336)
 
 p = load_parameters(m)
@@ -23,8 +23,8 @@ run(m)
 
 rcons_per_cap = m[:MarketDamages, :rcons_per_cap_MarketRemainConsumption]
 rcons_per_cap_compare = readpagedata(m, "test/validationdata/rcons_per_cap_MarketRemainConsumption.csv")
-@test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-1
+@test rcons_per_cap ≈ rcons_per_cap_compare rtol = 1e-1
 
 rgdp_per_cap = m[:MarketDamages, :rgdp_per_cap_MarketRemainGDP]
 rgdp_per_cap_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_MarketRemainGDP.csv")
-@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2
+@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol = 1e-2
diff --git a/test/test_N2Ocycle.jl b/test/test_N2Ocycle.jl
index 7d6f005a..d42404b7 100644
--- a/test/test_N2Ocycle.jl
+++ b/test/test_N2Ocycle.jl
@@ -6,15 +6,15 @@ include("../src/components/N2Ocycle.jl")
 
 add_comp!(m, n2ocycle)
 
-set_param!(m, :n2ocycle, :e_globalN2Oemissions, readpagedata(m,"test/validationdata/e_globalN2Oemissions.csv"))
-set_param!(m, :n2ocycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) #real value
-set_param!(m, :n2ocycle, :y_year_0, 2008.) #real value
-set_param!(m, :n2ocycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))
+set_param!(m, :n2ocycle, :e_globalN2Oemissions, readpagedata(m, "test/validationdata/e_globalN2Oemissions.csv"))
+set_param!(m, :n2ocycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) # real value
+set_param!(m, :n2ocycle, :y_year_0, 2008.) # real value
+set_param!(m, :n2ocycle, :rtl_g_landtemperature, readpagedata(m, "test/validationdata/rtl_g_landtemperature.csv"))
 
 ##running Model
 run(m)
 
-conc=m[:n2ocycle,  :c_N2Oconcentration]
-conc_compare=readpagedata(m,"test/validationdata/c_n2oconcentration.csv")
+conc = m[:n2ocycle,  :c_N2Oconcentration]
+conc_compare = readpagedata(m, "test/validationdata/c_n2oconcentration.csv")
 
-@test conc ≈ conc_compare rtol=1e-4
+@test conc ≈ conc_compare rtol = 1e-4
diff --git a/test/test_N2Oemissions.jl b/test/test_N2Oemissions.jl
index 07707e00..5443c385 100644
--- a/test/test_N2Oemissions.jl
+++ b/test/test_N2Oemissions.jl
@@ -6,15 +6,15 @@ include("../src/components/N2Oemissions.jl")
 
 add_comp!(m, n2oemissions)
 
-set_param!(m, :n2oemissions, :e0_baselineN2Oemissions, readpagedata(m,"data/e0_baselineN2Oemissions.csv"))
+set_param!(m, :n2oemissions, :e0_baselineN2Oemissions, readpagedata(m, "data/e0_baselineN2Oemissions.csv"))
 set_param!(m, :n2oemissions, :er_N2Oemissionsgrowth, readpagedata(m, "data/er_N2Oemissionsgrowth.csv"))
 
 ##running Model
 run(m)
 
 # Generated data
-emissions= m[:n2oemissions,  :e_regionalN2Oemissions]
+emissions = m[:n2oemissions,  :e_regionalN2Oemissions]
 # Recorded data
-emissions_compare=readpagedata(m, "test/validationdata/e_regionalN2Oemissions.csv")
+emissions_compare = readpagedata(m, "test/validationdata/e_regionalN2Oemissions.csv")
 
-@test emissions ≈ emissions_compare rtol=1e-3
+@test emissions ≈ emissions_compare rtol = 1e-3
diff --git a/test/test_N2Oforcing.jl b/test/test_N2Oforcing.jl
index bc2d4f65..a7ee0b1a 100644
--- a/test/test_N2Oforcing.jl
+++ b/test/test_N2Oforcing.jl
@@ -6,13 +6,13 @@ include("../src/components/N2Oforcing.jl")
 
 add_comp!(m, n2oforcing)
 
-set_param!(m, :n2oforcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/c_n2oconcentration.csv"))
-set_param!(m, :n2oforcing, :c_CH4concentration, readpagedata(m,"test/validationdata/c_ch4concentration.csv"))
+set_param!(m, :n2oforcing, :c_N2Oconcentration, readpagedata(m, "test/validationdata/c_n2oconcentration.csv"))
+set_param!(m, :n2oforcing, :c_CH4concentration, readpagedata(m, "test/validationdata/c_ch4concentration.csv"))
 
 ##running Model
 run(m)
 
-forcing=m[:n2oforcing,:f_N2Oforcing]
-forcing_compare=readpagedata(m,"test/validationdata/f_n2oforcing.csv")
+forcing = m[:n2oforcing,:f_N2Oforcing]
+forcing_compare = readpagedata(m, "test/validationdata/f_n2oforcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_NonMarketDamages.jl b/test/test_NonMarketDamages.jl
index e8baa4bd..7f3dc003 100644
--- a/test/test_NonMarketDamages.jl
+++ b/test/test_NonMarketDamages.jl
@@ -27,8 +27,8 @@ run(m)
 
 rcons_per_cap = m[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption]
 rcons_per_cap_compare = readpagedata(m, "test/validationdata/rcons_per_cap_NonMarketRemainConsumption.csv")
-@test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-2
+@test rcons_per_cap ≈ rcons_per_cap_compare rtol = 1e-2
 
 rgdp_per_cap = m[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP]
 rgdp_per_cap_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_NonMarketRemainGDP.csv")
-@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2
+@test rgdp_per_cap ≈ rgdp_per_cap_compare rtol = 1e-2
diff --git a/test/test_Population.jl b/test/test_Population.jl
index 697da7ea..f13ed506 100644
--- a/test/test_Population.jl
+++ b/test/test_Population.jl
@@ -20,4 +20,4 @@ pop = m[:Population, :pop_population]
 # Recorded data
 pop_compare = readpagedata(m, "test/validationdata/pop_population.csv")
 
-@test pop ≈ pop_compare rtol=1e-3
+@test pop ≈ pop_compare rtol = 1e-3
diff --git a/test/test_SLRDamages.jl b/test/test_SLRDamages.jl
index 899b078a..d584f3b4 100644
--- a/test/test_SLRDamages.jl
+++ b/test/test_SLRDamages.jl
@@ -23,28 +23,28 @@ run(m)
 
 cons_percap_aftercosts = m[:SLRDamages, :cons_percap_aftercosts]
 cons_percap_aftercosts_compare = readpagedata(m, "test/validationdata/cons_percap_aftercosts.csv")
-@test ones(10, 8) ≈ cons_percap_aftercosts ./ cons_percap_aftercosts_compare rtol=.001
+@test ones(10, 8) ≈ cons_percap_aftercosts ./ cons_percap_aftercosts_compare rtol = .001
 
 i_regionalimpactSLR = m[:SLRDamages, :i_regionalimpactSLR]
 i_regionalimpactSLR_compare = readpagedata(m, "test/validationdata/i_regionalimpact_SLRise.csv")
-@test i_regionalimpactSLR ≈ i_regionalimpactSLR_compare rtol=.001
+@test i_regionalimpactSLR ≈ i_regionalimpactSLR_compare rtol = .001
 
 iref_ImpactatReferenceGDPperCapSLR = m[:SLRDamages, :iref_ImpactatReferenceGDPperCapSLR]
 iref_ImpactatReferenceGDPperCapSLR_compare = readpagedata(m, "test/validationdata/iref_ImpactatReferenceGDPperCap_sea.csv")
-@test iref_ImpactatReferenceGDPperCapSLR ≈ iref_ImpactatReferenceGDPperCapSLR_compare rtol=.0001
+@test iref_ImpactatReferenceGDPperCapSLR ≈ iref_ImpactatReferenceGDPperCapSLR_compare rtol = .0001
 
 rcons_per_cap_SLRRemainConsumption = m[:SLRDamages, :rcons_per_cap_SLRRemainConsumption]
 rcons_per_cap_SLRRemainConsumption_compare = readpagedata(m, "test/validationdata/rcons_per_cap_SLRRemainConsumption.csv")
-@test rcons_per_cap_SLRRemainConsumption ≈ rcons_per_cap_SLRRemainConsumption_compare rtol=.01
+@test rcons_per_cap_SLRRemainConsumption ≈ rcons_per_cap_SLRRemainConsumption_compare rtol = .01
 
 rgdp_per_cap_SLRRemainGDP = m[:SLRDamages, :rgdp_per_cap_SLRRemainGDP]
 rgdp_per_cap_SLRRemainGDP_compare = readpagedata(m, "test/validationdata/rgdp_per_cap_SLRRemainGDP.csv")
-@test rgdp_per_cap_SLRRemainGDP ≈ rgdp_per_cap_SLRRemainGDP_compare rtol=.01
+@test rgdp_per_cap_SLRRemainGDP ≈ rgdp_per_cap_SLRRemainGDP_compare rtol = .01
 
 igdp_ImpactatActualGDPperCapSLR = m[:SLRDamages, :igdp_ImpactatActualGDPperCapSLR]
 igdp_ImpactatActualGDPperCapSLR_compare = readpagedata(m, "test/validationdata/igdp_ImpactatActualGDPperCap_sea.csv")
-@test igdp_ImpactatActualGDPperCapSLR ≈ igdp_ImpactatActualGDPperCapSLR_compare rtol=.01
+@test igdp_ImpactatActualGDPperCapSLR ≈ igdp_ImpactatActualGDPperCapSLR_compare rtol = .01
 
 isat_ImpactinclSaturationandAdaptationSLR = m[:SLRDamages, :isat_ImpactinclSaturationandAdaptationSLR]
 isat_ImpactinclSaturationandAdaptationSLR_compare = readpagedata(m, "test/validationdata/isat_ImpactinclSaturationandAdaptation_SLRise.csv")
-@test isat_ImpactinclSaturationandAdaptationSLR ≈ isat_ImpactinclSaturationandAdaptationSLR_compare rtol=.001
+@test isat_ImpactinclSaturationandAdaptationSLR ≈ isat_ImpactinclSaturationandAdaptationSLR_compare rtol = .001
diff --git a/test/test_SeaLevelRise.jl b/test/test_SeaLevelRise.jl
index 2c40fc46..b7a9a7d8 100644
--- a/test/test_SeaLevelRise.jl
+++ b/test/test_SeaLevelRise.jl
@@ -17,12 +17,12 @@ run(m)
 
 es_equilibriumSL = m[:SeaLevelRise, :es_equilibriumSL]
 es_equilibriumSL_compare = readpagedata(m, "test/validationdata/es_equilibriumSL.csv")
-@test ones(10) ≈ es_equilibriumSL ./ es_equilibriumSL rtol=.01
+@test ones(10) ≈ es_equilibriumSL ./ es_equilibriumSL rtol = .01
 
 s_sealevel = m[:SeaLevelRise, :s_sealevel]
 s_sealevel_compare = readpagedata(m, "test/validationdata/s_sealevel.csv")
-@test ones(10) ≈ s_sealevel ./ s_sealevel_compare rtol=.01
+@test ones(10) ≈ s_sealevel ./ s_sealevel_compare rtol = .01
 
 expfs_exponential = m[:SeaLevelRise, :expfs_exponential]
 expfs_exponential_compare = readpagedata(m, "test/validationdata/expfs_exponential.csv")
-@test ones(10) ≈ expfs_exponential ./ expfs_exponential_compare rtol=.01
+@test ones(10) ≈ expfs_exponential ./ expfs_exponential_compare rtol = .01
diff --git a/test/test_SulphateForcing.jl b/test/test_SulphateForcing.jl
index 9ace389e..25e8ca3c 100644
--- a/test/test_SulphateForcing.jl
+++ b/test/test_SulphateForcing.jl
@@ -13,7 +13,7 @@ set_leftover_params!(m, p)
 
 run(m)
 
-forcing=m[:SulphateForcing,:fs_sulphateforcing]
-forcing_compare=readpagedata(m,"test/validationdata/fs_sulphateforcing.csv")
+forcing = m[:SulphateForcing,:fs_sulphateforcing]
+forcing_compare = readpagedata(m, "test/validationdata/fs_sulphateforcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_TotalAbatementCosts.jl b/test/test_TotalAbatementCosts.jl
index b9c7ae8f..35f7228e 100644
--- a/test/test_TotalAbatementCosts.jl
+++ b/test/test_TotalAbatementCosts.jl
@@ -23,5 +23,5 @@ abate_cost_per_cap = m[:TotalAbatementCosts, :tct_per_cap_totalcostspercap]
 cost_compare = readpagedata(m, "test/validationdata/tct_totalcosts.csv")
 cost_cap_compare = readpagedata(m, "test/validationdata/tct_per_cap_totalcostspercap.csv")
 
-@test abate_cost ≈ cost_compare rtol=1e-4
-@test abate_cost_per_cap ≈ cost_cap_compare rtol=1e-7
+@test abate_cost ≈ cost_compare rtol = 1e-4
+@test abate_cost_per_cap ≈ cost_cap_compare rtol = 1e-7
diff --git a/test/test_TotalAdaptationCosts.jl b/test/test_TotalAdaptationCosts.jl
index d9d03a8a..7874e187 100644
--- a/test/test_TotalAdaptationCosts.jl
+++ b/test/test_TotalAdaptationCosts.jl
@@ -22,5 +22,5 @@ adapt_cost_per_cap = m[:TotalAdaptationCosts, :act_percap_adaptationcosts]
 cost_compare = readpagedata(m, "test/validationdata/act_adaptationcosts_tot.csv")
 cost_cap_compare = readpagedata(m, "test/validationdata/act_percap_adaptationcosts.csv")
 
-@test adapt_cost ≈ cost_compare rtol=1e-3
-@test adapt_cost_per_cap ≈ cost_cap_compare rtol=1e-6
+@test adapt_cost ≈ cost_compare rtol = 1e-3
+@test adapt_cost_per_cap ≈ cost_cap_compare rtol = 1e-6
diff --git a/test/test_TotalCosts.jl b/test/test_TotalCosts.jl
index 640e2d7b..177c95a7 100644
--- a/test/test_TotalCosts.jl
+++ b/test/test_TotalCosts.jl
@@ -10,8 +10,8 @@ atol = 1e-8
 
 @test all(isapprox.(m[:TotalCosts, :total_costs_percap_peryear],
     m[:EquityWeighting, :cons_percap_aftercosts] - m[:EquityWeighting, :rcons_percap_dis] + m[:EquityWeighting, :act_percap_adaptationcosts] + m[:EquityWeighting, :tct_percap_totalcosts_total],
-    atol = atol))
+    atol=atol))
 
 @test all(isapprox.(m[:TotalCosts, :total_damages_percap_peryear],
     m[:EquityWeighting, :cons_percap_aftercosts] - m[:EquityWeighting, :rcons_percap_dis] + m[:EquityWeighting, :act_percap_adaptationcosts], # without abatement costs, just damages
-    atol = atol))
\ No newline at end of file
+    atol=atol))
\ No newline at end of file
diff --git a/test/test_TotalForcing.jl b/test/test_TotalForcing.jl
index 42f27f56..b8306a82 100644
--- a/test/test_TotalForcing.jl
+++ b/test/test_TotalForcing.jl
@@ -6,15 +6,15 @@ include("../src/components/TotalForcing.jl")
 
 totalforcing = add_comp!(m, TotalForcing)
 
-totalforcing[:f_CO2forcing] = readpagedata(m,"test/validationdata/f_co2forcing.csv")
-totalforcing[:f_CH4forcing] = readpagedata(m,"test/validationdata/f_ch4forcing.csv")
-totalforcing[:f_N2Oforcing] = readpagedata(m,"test/validationdata/f_n2oforcing.csv")
-totalforcing[:f_lineargasforcing] = readpagedata(m,"test/validationdata/f_LGforcing.csv")
-totalforcing[:exf_excessforcing] = readpagedata(m,"data/exf_excessforcing.csv")
+totalforcing[:f_CO2forcing] = readpagedata(m, "test/validationdata/f_co2forcing.csv")
+totalforcing[:f_CH4forcing] = readpagedata(m, "test/validationdata/f_ch4forcing.csv")
+totalforcing[:f_N2Oforcing] = readpagedata(m, "test/validationdata/f_n2oforcing.csv")
+totalforcing[:f_lineargasforcing] = readpagedata(m, "test/validationdata/f_LGforcing.csv")
+totalforcing[:exf_excessforcing] = readpagedata(m, "data/exf_excessforcing.csv")
 
 run(m)
 
-forcing=m[:TotalForcing, :ft_totalforcing]
-forcing_compare=readpagedata(m,"test/validationdata/ft_totalforcing.csv")
+forcing = m[:TotalForcing, :ft_totalforcing]
+forcing_compare = readpagedata(m, "test/validationdata/ft_totalforcing.csv")
 
-@test forcing ≈ forcing_compare rtol=1e-3
+@test forcing ≈ forcing_compare rtol = 1e-3
diff --git a/test/test_climatemodel.jl b/test/test_climatemodel.jl
index 91a1521e..9892de24 100644
--- a/test/test_climatemodel.jl
+++ b/test/test_climatemodel.jl
@@ -6,4 +6,4 @@ include("../src/climate_model.jl")
 rt_g = m[:ClimateTemperature, :rt_g_globaltemperature]
 rt_g_compare = readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv")
 
-@test rt_g ≈ rt_g_compare rtol=1e-4
+@test rt_g ≈ rt_g_compare rtol = 1e-4
diff --git a/test/test_mainmodel.jl b/test/test_mainmodel.jl
index 4a4d5795..6562be1f 100644
--- a/test/test_mainmodel.jl
+++ b/test/test_mainmodel.jl
@@ -7,51 +7,51 @@ while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,0.,0.
     run(m)
 end
 
-#climate component
-temp=m[:ClimateTemperature,:rt_g_globaltemperature]
-temp_compare=readpagedata(m,"test/validationdata/rt_g_globaltemperature.csv")
-@test temp ≈ temp_compare rtol=1e-4
-
-slr=m[:SeaLevelRise,:s_sealevel]
-slr_compare=readpagedata(m,"test/validationdata/s_sealevel.csv")
-@test slr ≈ slr_compare rtol=1e-2
-
-#Socio-Economics
-gdp=m[:GDP,:gdp]
-gdp_compare=readpagedata(m,"test/validationdata/gdp.csv")
-@test gdp ≈ gdp_compare rtol=1
-
-pop=m[:Population,:pop_population]
-pop_compare=readpagedata(m,"test/validationdata/pop_population.csv")
-@test pop ≈ pop_compare rtol=0.001
-
-#Abatement Costs
-abatement=m[:TotalAbatementCosts,:tct_totalcosts]
-abatement_compare=readpagedata(m,"test/validationdata/tct_totalcosts.csv")
-@test abatement ≈ abatement_compare rtol=1e-2
-
-#Adaptation Costs
-adaptation=m[:TotalAdaptationCosts,:act_adaptationcosts_total]
-adaptation_compare=readpagedata(m, "test/validationdata/act_adaptationcosts_tot.csv")
-@test adaptation ≈ adaptation_compare rtol=1e-2
-
-#Damages
-damages=m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
-damages_compare=readpagedata(m,"test/validationdata/rcons_per_cap_DiscRemainConsumption.csv")
-@test damages ≈ damages_compare rtol=10
-#SLR damages
-slrdamages=m[:SLRDamages,:rcons_per_cap_SLRRemainConsumption]
-slrdamages_compare=readpagedata(m, "test/validationdata/rcons_per_cap_SLRRemainConsumption.csv")
-@test slrdamages ≈ slrdamages_compare rtol=0.1
-#Market damages
-mdamages=m[:MarketDamages,:rcons_per_cap_MarketRemainConsumption]
-mdamages_compare=readpagedata(m,"test/validationdata/rcons_per_cap_MarketRemainConsumption.csv")
-@test mdamages ≈ mdamages_compare rtol=1
-#NonMarket Damages
-nmdamages=m[:NonMarketDamages,:rcons_per_cap_NonMarketRemainConsumption]
-nmdamages_compare=readpagedata(m,"test/validationdata/rcons_per_cap_NonMarketRemainConsumption.csv")
-@test nmdamages ≈ nmdamages_compare rtol=1
+# climate component
+temp = m[:ClimateTemperature,:rt_g_globaltemperature]
+temp_compare = readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv")
+@test temp ≈ temp_compare rtol = 1e-4
+
+slr = m[:SeaLevelRise,:s_sealevel]
+slr_compare = readpagedata(m, "test/validationdata/s_sealevel.csv")
+@test slr ≈ slr_compare rtol = 1e-2
+
+# Socio-Economics
+gdp = m[:GDP,:gdp]
+gdp_compare = readpagedata(m, "test/validationdata/gdp.csv")
+@test gdp ≈ gdp_compare rtol = 1
+
+pop = m[:Population,:pop_population]
+pop_compare = readpagedata(m, "test/validationdata/pop_population.csv")
+@test pop ≈ pop_compare rtol = 0.001
+
+# Abatement Costs
+abatement = m[:TotalAbatementCosts,:tct_totalcosts]
+abatement_compare = readpagedata(m, "test/validationdata/tct_totalcosts.csv")
+@test abatement ≈ abatement_compare rtol = 1e-2
+
+# Adaptation Costs
+adaptation = m[:TotalAdaptationCosts,:act_adaptationcosts_total]
+adaptation_compare = readpagedata(m, "test/validationdata/act_adaptationcosts_tot.csv")
+@test adaptation ≈ adaptation_compare rtol = 1e-2
+
+# Damages
+damages = m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
+damages_compare = readpagedata(m, "test/validationdata/rcons_per_cap_DiscRemainConsumption.csv")
+@test damages ≈ damages_compare rtol = 10
+# SLR damages
+slrdamages = m[:SLRDamages,:rcons_per_cap_SLRRemainConsumption]
+slrdamages_compare = readpagedata(m, "test/validationdata/rcons_per_cap_SLRRemainConsumption.csv")
+@test slrdamages ≈ slrdamages_compare rtol = 0.1
+# Market damages
+mdamages = m[:MarketDamages,:rcons_per_cap_MarketRemainConsumption]
+mdamages_compare = readpagedata(m, "test/validationdata/rcons_per_cap_MarketRemainConsumption.csv")
+@test mdamages ≈ mdamages_compare rtol = 1
+# NonMarket Damages
+nmdamages = m[:NonMarketDamages,:rcons_per_cap_NonMarketRemainConsumption]
+nmdamages_compare = readpagedata(m, "test/validationdata/rcons_per_cap_NonMarketRemainConsumption.csv")
+@test nmdamages ≈ nmdamages_compare rtol = 1
 
 te = m[:EquityWeighting, :te_totaleffect]
 te_compare = 213208136.69903600
-@test te ≈ te_compare rtol=1e4
+@test te ≈ te_compare rtol = 1e4
diff --git a/test/test_mainmodel_policyb.jl b/test/test_mainmodel_policyb.jl
index 33de8dd3..5e03ea08 100644
--- a/test/test_mainmodel_policyb.jl
+++ b/test/test_mainmodel_policyb.jl
@@ -12,14 +12,14 @@ end
 
 # Climate component
 temp = m[:ClimateTemperature, :rt_g_globaltemperature]
-temp_compare = readpagedata(m,"test/validationdata/policy-b/rt_g_globaltemperature.csv")
-@test temp ≈ temp_compare rtol=1e-4
+temp_compare = readpagedata(m, "test/validationdata/policy-b/rt_g_globaltemperature.csv")
+@test temp ≈ temp_compare rtol = 1e-4
 
 # Abatement Costs
 abatement = m[:TotalAbatementCosts, :tct_totalcosts]
-abatement_compare = readpagedata(m,"test/validationdata/policy-b/tct_totalcosts.csv")
-@test abatement ≈ abatement_compare rtol=1e-2
+abatement_compare = readpagedata(m, "test/validationdata/policy-b/tct_totalcosts.csv")
+@test abatement ≈ abatement_compare rtol = 1e-2
 
 te = m[:EquityWeighting, :te_totaleffect]
 te_compare = 224812034.49558273
-@test te ≈ te_compare rtol=1e4
+@test te ≈ te_compare rtol = 1e4
diff --git a/test/test_mcs.jl b/test/test_mcs.jl
index cf1fe897..c634912f 100644
--- a/test/test_mcs.jl
+++ b/test/test_mcs.jl
@@ -29,7 +29,7 @@ if regenerate
     MimiPAGE2009.do_monte_carlo_runs(100_000)
     df = DataFrame(load(joinpath(@__DIR__, "../output/mimipagemontecarlooutput.csv")))
 
-     for ii in 1:nrow(compare)
+    for ii in 1:nrow(compare)
         name = Symbol(compare[ii, :Variable_Name])
         if kurtosis(df[name]) > 2.9 # exponential distribution
             if name == :tpc # negative across all quantiles
@@ -62,7 +62,7 @@ for ii in 1:nrow(compare)
 
         expected = transform(compare[ii, Symbol("perc_$(trunc(Int, qval * 100))")])
 
-        #println("$name x $qval: $estimated ≈ $expected rtol=$(ceil(confidence * stderr, -trunc(Int, log10(stderr))))")
-        @test estimated ≈ expected rtol=ceil(confidence * stderr; digits = -trunc(Int, log10(stderr)))
+        # println("$name x $qval: $estimated ≈ $expected rtol=$(ceil(confidence * stderr, -trunc(Int, log10(stderr))))")
+        @test estimated ≈ expected rtol = ceil(confidence * stderr; digits=-trunc(Int, log10(stderr)))
     end
 end
diff --git a/test/test_standard_api.jl b/test/test_standard_api.jl
index dffaadd7..bf9ed29f 100644
--- a/test/test_standard_api.jl
+++ b/test/test_standard_api.jl
@@ -4,52 +4,52 @@ using Mimi
 @testset "Standard API" begin 
 
 # Test that the function does not error and returns a valid value
-scc1 = MimiPAGE2009.compute_scc(year=2020)
-@test scc1 isa Float64
+    scc1 = MimiPAGE2009.compute_scc(year=2020)
+    @test scc1 isa Float64
 
 # Test that a higher discount rate makes a lower scc value
-scc2 = MimiPAGE2009.compute_scc(year=2020, eta=0., prtp=0.03)
-@test scc2 < scc1   
+    scc2 = MimiPAGE2009.compute_scc(year=2020, eta=0., prtp=0.03)
+    @test scc2 < scc1   
 
 # Test with a modified model
-m = MimiPAGE2009.get_model()
-set_param!(m, :tcr_transientresponse, 3)
-scc3 = MimiPAGE2009.compute_scc(m, year=2020)
-@test scc3 > scc1
+    m = MimiPAGE2009.get_model()
+    set_param!(m, :tcr_transientresponse, 3)
+    scc3 = MimiPAGE2009.compute_scc(m, year=2020)
+    @test scc3 > scc1
 
 # Test get_marginal_model
-mm = MimiPAGE2009.get_marginal_model(year = 2040)
-mm[:ClimateTemperature, :rt_realizedtemperature]
+    mm = MimiPAGE2009.get_marginal_model(year=2040)
+    mm[:ClimateTemperature, :rt_realizedtemperature]
 
 # Test compute_scc_mm
-result = MimiPAGE2009.compute_scc_mm(year=2050)
-@test result.scc > scc1
-@test result.mm isa Mimi.MarginalModel
+    result = MimiPAGE2009.compute_scc_mm(year=2050)
+    @test result.scc > scc1
+    @test result.mm isa Mimi.MarginalModel
 
 # Test Monte Carlo SCC support
-sccs1 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=350)
-sccs2 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=350)
-@test sccs1 == sccs2
+    sccs1 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=350)
+    sccs2 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=350)
+    @test sccs1 == sccs2
 
-sccs3 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=351)
-@test sccs3 != sccs1
+    sccs3 = MimiPAGE2009.compute_scc(year=2020, n=10, seed=351)
+    @test sccs3 != sccs1
 
-sccs4 = MimiPAGE2009.compute_scc(year=2020, prtp=0.02, eta=1.5, n=10, seed=350)
-@test sccs4 != sccs1    # test that the user specified discounting scheme overrides the default random variable values
+    sccs4 = MimiPAGE2009.compute_scc(year=2020, prtp=0.02, eta=1.5, n=10, seed=350)
+    @test sccs4 != sccs1    # test that the user specified discounting scheme overrides the default random variable values
 
-#Test equity weighting options
+# Test equity weighting options
 
 # Test that regional equity weighting inreases the SCC
-@test MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = true) >
-        MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = false)
+    @test MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=true) >
+        MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=false)
 
 # Test that when eta==0, equity_weighting does not change the SCC value
-@test MimiPAGE2009.compute_scc(year=2020, eta=0., prtp = 0.03, equity_weighting = true) ≈
-        MimiPAGE2009.compute_scc(year=2020, eta=0., prtp = 0.03, equity_weighting = false) atol=1e3
+    @test MimiPAGE2009.compute_scc(year=2020, eta=0., prtp=0.03, equity_weighting=true) ≈
+        MimiPAGE2009.compute_scc(year=2020, eta=0., prtp=0.03, equity_weighting=false) atol = 1e3
 
 # Test Monte Carlo w/ and w/o equity weighting, with the same seed
-scc10a = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = true, n = 10, seed=350)
-scc10b = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp = 0.01, equity_weighting = false, n = 10, seed=350)
-@test all(scc10a .> scc10b)
+    scc10a = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=true, n=10, seed=350)
+    scc10b = MimiPAGE2009.compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=false, n=10, seed=350)
+    @test all(scc10a .> scc10b)
 
 end
\ No newline at end of file