Skip to content

Commit

Permalink
Format files using DocumentFormat
Browse files Browse the repository at this point in the history
  • Loading branch information
lrennels authored Oct 17, 2022
1 parent c58314f commit 56db778
Show file tree
Hide file tree
Showing 31 changed files with 943 additions and 947 deletions.
6 changes: 3 additions & 3 deletions examples/main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,10 @@ N = 10000
seed = 350

Random.seed!(seed)
MimiIWG.run_scc_mcs(DICE, gas=:CO2, trials = N)
MimiIWG.run_scc_mcs(DICE, gas=:CO2, trials=N)

Random.seed!(seed)
MimiIWG.run_scc_mcs(FUND, gas=:CO2, trials = N)
MimiIWG.run_scc_mcs(FUND, gas=:CO2, trials=N)

Random.seed!(seed)
MimiIWG.run_scc_mcs(PAGE, gas=:CO2, trials = N)
MimiIWG.run_scc_mcs(PAGE, gas=:CO2, trials=N)
12 changes: 6 additions & 6 deletions examples/other_ghg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,13 @@ scn2o = MimiIWG.compute_scc(PAGE, USG1, gas=:N2O, year=2020, discount=0.025)

# DICE
# (will run for all 3 IWG discount rates if none are specified)
MimiIWG.run_scc_mcs(DICE, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(DICE, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(DICE, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir="output")
MimiIWG.run_scc_mcs(DICE, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir="output")

# FUND
MimiIWG.run_scc_mcs(FUND, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(FUND, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(FUND, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir="output")
MimiIWG.run_scc_mcs(FUND, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir="output")

# PAGE
MimiIWG.run_scc_mcs(PAGE, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(PAGE, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir = "output")
MimiIWG.run_scc_mcs(PAGE, gas=:CH4, trials=10_000, perturbation_years=[2020], output_dir="output")
MimiIWG.run_scc_mcs(PAGE, gas=:N2O, trials=10_000, perturbation_years=[2020], output_dir="output")
6 changes: 3 additions & 3 deletions src/MimiIWG.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ using Interpolations
using Mimi
using MimiDICE2010
using MimiFUND # pinned to version 3.8 in package registration Compat.toml
using MimiPAGE2009
using MimiPAGE2009
using StatsBase
using XLSX: readxlsx, openxlsx, readdata
using CSVFiles
Expand All @@ -17,7 +17,7 @@ using Query
using Statistics

export DICE, FUND, PAGE, # export the enumerated model_choice options
USG1, USG2, USG3, USG4, USG5 # export the enumerated scenario_choice options
USG1, USG2, USG3, USG4, USG5 # export the enumerated scenario_choice options

# General constants and functions
include("core/constants.jl")
Expand Down Expand Up @@ -54,4 +54,4 @@ include("montecarlo/PAGE_mcs.jl")
include("montecarlo/run_scc_mcs.jl")


end
end
32 changes: 16 additions & 16 deletions src/components/IWG_DICE_ScenarioChoice.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,35 +6,35 @@
scenarios = Index()

# Variables (each one has its value set for the chosen scenario in the first timestep)
l = Variable(index = [time]) # Population
E = Variable(index = [time]) # Total CO2 emissions
forcoth = Variable(index = [time]) # other forcing
al = Variable(index = [time]) # total factor productivity
k0 = Variable() # initial capital stock
l = Variable(index=[time]) # Population
E = Variable(index=[time]) # Total CO2 emissions
forcoth = Variable(index=[time]) # other forcing
al = Variable(index=[time]) # total factor productivity
k0 = Variable() # initial capital stock

# The number for which scenario to use
scenario_num = Parameter{Integer}()

# Parameters (each one holds all five scenarios)
l_all = Parameter(index = [time, scenarios])
E_all = Parameter(index = [time, scenarios])
forcoth_all = Parameter(index = [time, scenarios])
al_all = Parameter(index = [time, scenarios])
k0_all = Parameter(index = [scenarios])
l_all = Parameter(index=[time, scenarios])
E_all = Parameter(index=[time, scenarios])
forcoth_all = Parameter(index=[time, scenarios])
al_all = Parameter(index=[time, scenarios])
k0_all = Parameter(index=[scenarios])

function run_timestep(p, v, d, t)
if is_first(t)
# Get the specified scenario
scenario_num = p.scenario_num
if ! (scenario_num in d.scenarios)
if !(scenario_num in d.scenarios)
error("Invalid :scenario_num in :IWGScenarioChoice component: $scenario_num. :scenario_num must be in $(d.scenarios).")
else
# Copy over all of the values for that scenario
v.l[:] = p.l_all[:, scenario_num]
v.E[:] = p.E_all[:, scenario_num]
v.forcoth[:] = p.forcoth_all[:, scenario_num]
v.al[:] = p.al_all[:, scenario_num]
v.k0 = p.k0_all[scenario_num]
v.l[:] = p.l_all[:, scenario_num]
v.E[:] = p.E_all[:, scenario_num]
v.forcoth[:] = p.forcoth_all[:, scenario_num]
v.al[:] = p.al_all[:, scenario_num]
v.k0 = p.k0_all[scenario_num]
end
end
end
Expand Down
22 changes: 11 additions & 11 deletions src/components/IWG_DICE_climatedynamics.jl
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
# Since the IWG sampled from the Roe & Baker distribution for equilibrium climate sensitivity (t2xco2), they added in the catch case for values less than 0.5 in the code below.
@defcomp IWG_DICE_climatedynamics begin

TATM = Variable(index=[time]) # Increase in temperature of atmosphere (degrees C from 1900)
TOCEAN = Variable(index=[time]) # Increase in temperature of lower oceans (degrees C from 1900)
TA_eq = Variable(index=[time])

FORC = Parameter(index=[time]) # Increase in radiative forcing (watts per m2 from 1900)
fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)
t2xco2 = Parameter(default=3) # Equilibrium temp impact (oC per doubling CO2)
tatm0 = Parameter() # Initial atmospheric temp change (C from 1900) in 2005
TATM = Variable(index=[time]) # Increase in temperature of atmosphere (degrees C from 1900)
TOCEAN = Variable(index=[time]) # Increase in temperature of lower oceans (degrees C from 1900)
TA_eq = Variable(index=[time])

FORC = Parameter(index=[time]) # Increase in radiative forcing (watts per m2 from 1900)
fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)
t2xco2 = Parameter(default=3) # Equilibrium temp impact (oC per doubling CO2)
tatm0 = Parameter() # Initial atmospheric temp change (C from 1900) in 2005
tocean0 = Parameter() # Initial lower stratum temp change (C from 1900)

# Transient TSC Correction ("Speed of Adjustment Parameter")
Expand All @@ -32,16 +32,16 @@
elseif p.t2xco2 < 0.5
v.TATM[t] = v.TA_eq[t]
else
v.TATM[t] = v.TATM[t - 1] + p.c1 * ((p.FORC[t] - (p.fco22x/p.t2xco2) * v.TATM[t - 1]) - (p.c3 * (v.TATM[t - 1] - v.TOCEAN[t - 1])))
v.TATM[t] = v.TATM[t-1] + p.c1 * ((p.FORC[t] - (p.fco22x / p.t2xco2) * v.TATM[t-1]) - (p.c3 * (v.TATM[t-1] - v.TOCEAN[t-1])))
end

# Define function for TOCEAN
if is_first(t)
v.TOCEAN[t] = p.tocean0
else
v.TOCEAN[t] = v.TOCEAN[t - 1] + p.c4 * (v.TATM[t - 1] - v.TOCEAN[t - 1])
v.TOCEAN[t] = v.TOCEAN[t-1] + p.c4 * (v.TATM[t-1] - v.TOCEAN[t-1])
end

end

end
end
36 changes: 18 additions & 18 deletions src/components/IWG_DICE_co2cycle.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,23 +3,23 @@
# - E is in units of GtCO2 per decade instead of per year (so it doesn't get mulitiplied by 10)
@defcomp IWG_DICE_co2cycle begin

MAT = Variable(index=[time]) # Carbon concentration increase in atmosphere (GtC from 1750)
ML = Variable(index=[time]) # Carbon concentration increase in lower oceans (GtC from 1750)
MU = Variable(index=[time]) # Carbon concentration increase in shallow oceans (GtC from 1750)
MAT = Variable(index=[time]) # Carbon concentration increase in atmosphere (GtC from 1750)
ML = Variable(index=[time]) # Carbon concentration increase in lower oceans (GtC from 1750)
MU = Variable(index=[time]) # Carbon concentration increase in shallow oceans (GtC from 1750)

E = Parameter(index=[time]) # Total CO2 emissions (GtCO2 per decade)
mat0 = Parameter() # Initial Concentration in atmosphere 2010 (GtC)
ml0 = Parameter() # Initial Concentration in lower strata (GtC)
mu0 = Parameter() # Initial Concentration in upper strata (GtC)
E = Parameter(index=[time]) # Total CO2 emissions (GtCO2 per decade)
mat0 = Parameter() # Initial Concentration in atmosphere 2010 (GtC)
ml0 = Parameter() # Initial Concentration in lower strata (GtC)
mu0 = Parameter() # Initial Concentration in upper strata (GtC)

# Parameters for long-run consistency of carbon cycle
b11 = Parameter() # Carbon cycle transition matrix atmosphere to atmosphere
b12 = Parameter() # Carbon cycle transition matrix atmosphere to shallow ocean
b21 = Parameter() # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
b22 = Parameter() # Carbon cycle transition matrix shallow ocean to shallow oceans
b23 = Parameter() # Carbon cycle transition matrix shallow to deep ocean
b32 = Parameter() # Carbon cycle transition matrix deep ocean to shallow ocean
b33 = Parameter() # Carbon cycle transition matrix deep ocean to deep oceans
b11 = Parameter() # Carbon cycle transition matrix atmosphere to atmosphere
b12 = Parameter() # Carbon cycle transition matrix atmosphere to shallow ocean
b21 = Parameter() # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
b22 = Parameter() # Carbon cycle transition matrix shallow ocean to shallow oceans
b23 = Parameter() # Carbon cycle transition matrix shallow to deep ocean
b32 = Parameter() # Carbon cycle transition matrix deep ocean to shallow ocean
b33 = Parameter() # Carbon cycle transition matrix deep ocean to deep oceans


function run_timestep(p, v, d, t)
Expand All @@ -28,14 +28,14 @@
if is_first(t)
v.MU[t] = p.mu0
else
v.MU[t] = v.MAT[t - 1] * p.b12 + v.MU[t - 1] * p.b22 + v.ML[t - 1] * p.b32
v.MU[t] = v.MAT[t-1] * p.b12 + v.MU[t-1] * p.b22 + v.ML[t-1] * p.b32
end

# Define function for ML
if is_first(t)
v.ML[t] = p.ml0
else
v.ML[t] = v.ML[t - 1] * p.b33 + v.MU[t - 1] * p.b23
v.ML[t] = v.ML[t-1] * p.b33 + v.MU[t-1] * p.b23
end

# Define function for MAT
Expand All @@ -44,9 +44,9 @@
# and also calculate MAT[2] below
end
if !is_last(t)
v.MAT[t + 1] = v.MAT[t] * p.b11 + v.MU[t] * p.b21 + p.E[t]
v.MAT[t+1] = v.MAT[t] * p.b11 + v.MU[t] * p.b21 + p.E[t]
end

end

end
end
40 changes: 20 additions & 20 deletions src/components/IWG_DICE_neteconomy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,23 @@
# In the original component, there are different values for the first and last timestep calculations of C.
@defcomp IWG_DICE_neteconomy begin

ABATECOST = Variable(index=[time]) # Cost of emissions reductions (trillions 2005 USD per year)
C = Variable(index=[time]) # Consumption (trillions 2005 US dollars per year)
CPC = Variable(index=[time]) # Per capita consumption (thousands 2005 USD per year)
CPRICE = Variable(index=[time]) # Carbon price (2005$ per ton of CO2)
I = Variable(index=[time]) # Investment (trillions 2005 USD per year)
Y = Variable(index=[time]) # Gross world product net of abatement and damages (trillions 2005 USD per year)
YNET = Variable(index=[time]) # Output net of damages equation (trillions 2005 USD per year)

cost1 = Parameter(index=[time]) # Abatement cost function coefficient
DAMFRAC = Parameter(index=[time]) # Damages (fraction of gross output)
l = Parameter(index=[time]) # Level of population and labor
MIU = Parameter(index=[time]) # Emission control rate GHGs
partfract = Parameter(index=[time]) # Fraction of emissions in control regime
pbacktime = Parameter(index=[time]) # Backstop price
S = Parameter(index=[time]) # Gross savings rate as fraction of gross world product
YGROSS = Parameter(index=[time]) # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
expcost2 = Parameter() # Exponent of control cost function
ABATECOST = Variable(index=[time]) # Cost of emissions reductions (trillions 2005 USD per year)
C = Variable(index=[time]) # Consumption (trillions 2005 US dollars per year)
CPC = Variable(index=[time]) # Per capita consumption (thousands 2005 USD per year)
CPRICE = Variable(index=[time]) # Carbon price (2005$ per ton of CO2)
I = Variable(index=[time]) # Investment (trillions 2005 USD per year)
Y = Variable(index=[time]) # Gross world product net of abatement and damages (trillions 2005 USD per year)
YNET = Variable(index=[time]) # Output net of damages equation (trillions 2005 USD per year)

cost1 = Parameter(index=[time]) # Abatement cost function coefficient
DAMFRAC = Parameter(index=[time]) # Damages (fraction of gross output)
l = Parameter(index=[time]) # Level of population and labor
MIU = Parameter(index=[time]) # Emission control rate GHGs
partfract = Parameter(index=[time]) # Fraction of emissions in control regime
pbacktime = Parameter(index=[time]) # Backstop price
S = Parameter(index=[time]) # Gross savings rate as fraction of gross world product
YGROSS = Parameter(index=[time]) # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
expcost2 = Parameter() # Exponent of control cost function


function run_timestep(p, v, d, t)
Expand All @@ -42,10 +42,10 @@
v.CPC[t] = 1000 * v.C[t] / p.l[t]

# Define function for CPRICE
if t.t==26
v.CPRICE[t] = v.CPRICE[t - 1]
if t.t == 26
v.CPRICE[t] = v.CPRICE[t-1]
else
v.CPRICE[t] = p.pbacktime[t] * 1000 * p.MIU[t] ^ (p.expcost2 - 1)
v.CPRICE[t] = p.pbacktime[t] * 1000 * p.MIU[t]^(p.expcost2 - 1)
end

end
Expand Down
18 changes: 9 additions & 9 deletions src/components/IWG_DICE_radiativeforcing.jl
Original file line number Diff line number Diff line change
@@ -1,23 +1,23 @@
# Difference from the original DICE radiativeforcing component:
# - In the final timestep, the IWG calculates MAT_avg differently
@defcomp IWG_DICE_radiativeforcing begin
FORC = Variable(index=[time]) # Increase in radiative forcing (watts per m2 from 1900)
FORC = Variable(index=[time]) # Increase in radiative forcing (watts per m2 from 1900)

forcoth = Parameter(index=[time]) # Exogenous forcing for other greenhouse gases
MAT = Parameter(index=[time]) # Carbon concentration increase in atmosphere (GtC from 1750)
fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)
forcoth = Parameter(index=[time]) # Exogenous forcing for other greenhouse gases
MAT = Parameter(index=[time]) # Carbon concentration increase in atmosphere (GtC from 1750)
fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)


function run_timestep(p, v, d, t)

if !is_last(t)
MAT_avg = (p.MAT[t] + p.MAT[t + 1]) / 2
else
MAT_avg = 0.9796 * (p.MAT[t - 1] + p.MAT[t]) / 2
MAT_avg = (p.MAT[t] + p.MAT[t+1]) / 2
else
MAT_avg = 0.9796 * (p.MAT[t-1] + p.MAT[t]) / 2
end

v.FORC[t] = p.fco22x * (log((MAT_avg + 0.000001) / 596.4 ) / log(2)) + p.forcoth[t]
v.FORC[t] = p.fco22x * (log((MAT_avg + 0.000001) / 596.4) / log(2)) + p.forcoth[t]

end

end
end
Loading

0 comments on commit 56db778

Please sign in to comment.