Skip to content

Commit

Permalink
Format files using DocumentFormat
Browse files Browse the repository at this point in the history
  • Loading branch information
github-actions[bot] authored May 5, 2021
1 parent bba29d5 commit cbe8336
Show file tree
Hide file tree
Showing 28 changed files with 701 additions and 701 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)
20 changes: 10 additions & 10 deletions examples/other_ghg.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
using MimiIWG

#------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Deterministic single runs of SC-CH4 and SC-N2O
#------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

# DICE
scch4 = MimiIWG.compute_scc(DICE, USG1, gas=:CH4, year=2020, discount=0.025)
Expand All @@ -16,19 +16,19 @@ scn2o = MimiIWG.compute_scc(FUND, USG1, gas=:N2O, year=2020, discount=0.025)
scch4 = MimiIWG.compute_scc(PAGE, USG1, gas=:CH4, year=2020, discount=0.025)
scn2o = MimiIWG.compute_scc(PAGE, USG1, gas=:N2O, year=2020, discount=0.025)

#------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Full Monte Carlo simulations for SC-CH4 and SC-N2O for each model
#------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

# 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")
2 changes: 1 addition & 1 deletion src/MimiIWG.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ using MimiDICE2010
using MimiFUND # pinned to version 3.8 in package registration Compat.toml
using MimiPAGE2009
using StatsBase
using XLSX: readxlsx
using XLSX:readxlsx

export DICE, FUND, PAGE, # export the enumerated model_choice options
USG1, USG2, USG3, USG4, USG5 # export the enumerated scenario_choice options
Expand Down
40 changes: 20 additions & 20 deletions src/components/IWG_DICE_ScenarioChoice.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,36 +6,36 @@
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
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)
if is_first(t)
# Get the specified scenario
scenario_num = p.scenario_num
if ! (scenario_num in d.scenarios)
error("Invalid :scenario_num in :IWGScenarioChoice component: $scenario_num. :scenario_num must be in $(d.scenarios).")
else
scenario_num = p.scenario_num
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]
end
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
end
32 changes: 16 additions & 16 deletions src/components/IWG_DICE_climatedynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,28 +20,28 @@
function run_timestep(p, v, d, t)

# TA_eq added by EPA "to avoid over- and under-shoots to eq T when CS is very low
if is_first(t)
v.TA_eq[t] = 0.083
else
v.TA_eq[t] = p.FORC[t] * p.t2xco2 / p.fco22x
end
if is_first(t)
v.TA_eq[t] = 0.083
else
v.TA_eq[t] = p.FORC[t] * p.t2xco2 / p.fco22x
end

# Define function for TATM
if is_first(t)
v.TATM[t] = p.tatm0
elseif p.t2xco2 < 0.5
if is_first(t)
v.TATM[t] = p.tatm0
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])))
end
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])
end

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])
end

end

end
34 changes: 17 additions & 17 deletions src/components/IWG_DICE_co2cycle.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,28 +25,28 @@
function run_timestep(p, v, d, t)

# Define function for MU
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
end
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
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
end
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
end

# Define function for MAT
if is_first(t)
v.MAT[t] = p.mat0
if is_first(t)
v.MAT[t] = p.mat0
# 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]
end

end
if !is_last(t)
v.MAT[t + 1] = v.MAT[t] * p.b11 + v.MU[t] * p.b21 + p.E[t]
end

end

end
24 changes: 12 additions & 12 deletions src/components/IWG_DICE_neteconomy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,30 +24,30 @@
function run_timestep(p, v, d, t)

# Define function for YNET
v.YNET[t] = p.YGROSS[t] / (1 + p.DAMFRAC[t])
v.YNET[t] = p.YGROSS[t] / (1 + p.DAMFRAC[t])

# Define function for ABATECOST
v.ABATECOST[t] = p.YGROSS[t] * p.cost1[t] * (p.MIU[t]^p.expcost2) * (p.partfract[t]^(1 - p.expcost2))
v.ABATECOST[t] = p.YGROSS[t] * p.cost1[t] * (p.MIU[t]^p.expcost2) * (p.partfract[t]^(1 - p.expcost2))

# Define function for Y
v.Y[t] = v.YNET[t] - v.ABATECOST[t]
v.Y[t] = v.YNET[t] - v.ABATECOST[t]

# Define function for I
v.I[t] = p.S[t] * v.Y[t]
v.I[t] = p.S[t] * v.Y[t]

# Define function for C
v.C[t] = v.Y[t] - v.I[t]
v.C[t] = v.Y[t] - v.I[t]

# Define function for CPC
v.CPC[t] = 1000 * v.C[t] / p.l[t]
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]
else
v.CPRICE[t] = p.pbacktime[t] * 1000 * p.MIU[t] ^ (p.expcost2 - 1)
end

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)
end

end

end
14 changes: 7 additions & 7 deletions src/components/IWG_DICE_radiativeforcing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,14 @@

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
end
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
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
78 changes: 39 additions & 39 deletions src/components/IWG_DICE_simple_gas_cycle.jl
Original file line number Diff line number Diff line change
@@ -1,64 +1,64 @@
@defcomp IWG_DICE_simple_gas_cycle begin

# Output: forcing
F_CH4 = Variable(index = [decades]) # Path of CH4 forcing [W/m^2] - decadal time step
F_N2O = Variable(index = [decades]) # Path of N2O forcing [W/m^2] - decadal time step
F_CH4 = Variable(index=[decades]) # Path of CH4 forcing [W/m^2] - decadal time step
F_N2O = Variable(index=[decades]) # Path of N2O forcing [W/m^2] - decadal time step

# Input: emissions
E_CH4A = Parameter(index = [time]) # Annual CH4 emissions
E_N2OA = Parameter(index = [time]) # Annual N2O emissions
E_CH4A = Parameter(index=[time]) # Annual CH4 emissions
E_N2OA = Parameter(index=[time]) # Annual N2O emissions

# Atmospheric concentration dynamics coefficients
delta_ch4 = Parameter(default = 1/12) # Annual rate of decay of CH4 (IPCC FAR)
gamma_ch4 = Parameter(default = 0.3597) # Mass to volume conversion factor [ppb CH4/Mt CH4]
alpha_ch4 = Parameter(default = 0.036) # Radiative forcing parameter (TAR page 358)
delta_n2o = Parameter(default = 1/114) # Rate of decay of N2O (IPCC FAR)
gamma_n2o = Parameter(default = 0.2079) # Mass to volume conversion factor [ppb N2O/Mt N]
alpha_n2o = Parameter(default = 0.12) # Radiative forcing parameter
ipcc_adj = Parameter(default = 1.4) # IPCC AR4 adjustment for tropospheric ozone effect (25%) and stratospheric water vapor effect (15%)

M_pre = Parameter(default = 700) # Pre-industrial CH4 concentrations [ppb](TAR page 358)
N_pre = Parameter(default = 270) # Pre-industrial N2O concentrations [ppb](TAR page 358)
delta_ch4 = Parameter(default=1 / 12) # Annual rate of decay of CH4 (IPCC FAR)
gamma_ch4 = Parameter(default=0.3597) # Mass to volume conversion factor [ppb CH4/Mt CH4]
alpha_ch4 = Parameter(default=0.036) # Radiative forcing parameter (TAR page 358)
delta_n2o = Parameter(default=1 / 114) # Rate of decay of N2O (IPCC FAR)
gamma_n2o = Parameter(default=0.2079) # Mass to volume conversion factor [ppb N2O/Mt N]
alpha_n2o = Parameter(default=0.12) # Radiative forcing parameter
ipcc_adj = Parameter(default=1.4) # IPCC AR4 adjustment for tropospheric ozone effect (25%) and stratospheric water vapor effect (15%)

M_pre = Parameter(default=700) # Pre-industrial CH4 concentrations [ppb](TAR page 358)
N_pre = Parameter(default=270) # Pre-industrial N2O concentrations [ppb](TAR page 358)

M_AA_2005 = Parameter(default = 1774) # 2005 concentration of CH4
N_AA_2005 = Parameter(default = 319) # 2005 concentration of N2O
M_AA_2005 = Parameter(default=1774) # 2005 concentration of CH4
N_AA_2005 = Parameter(default=319) # 2005 concentration of N2O

# Other intermediate variables to calculate
M_AA = Variable(index = [time]) # Atmospheric CH4 concentration (annual)
F_CH4A = Variable(index = [time]) # Contribution of atmospheric CH4 to radiative forcing (annual)
N_AA = Variable(index = [time]) # Atmospheric N2O concentration (annual)
F_N2OA = Variable(index = [time]) # Contribution of atmospheric N2O to radiative forcing (annual)
M_AA = Variable(index=[time]) # Atmospheric CH4 concentration (annual)
F_CH4A = Variable(index=[time]) # Contribution of atmospheric CH4 to radiative forcing (annual)
N_AA = Variable(index=[time]) # Atmospheric N2O concentration (annual)
F_N2OA = Variable(index=[time]) # Contribution of atmospheric N2O to radiative forcing (annual)
pre_f = Variable() # pre-industrial interaction effect

function run_timestep(p, v, d, t)

function f(M_A, N_A)
function f(M_A, N_A)
# calculate the interaction effect on radiative forcing
return 0.47 * log(1 + 2.01 * 10 ^ -5 * (M_A * N_A) ^ 0.75 + 5.31 * 10^-15 * M_A * (M_A * N_A) ^ 1.52)
end
return 0.47 * log(1 + 2.01 * 10^-5 * (M_A * N_A)^0.75 + 5.31 * 10^-15 * M_A * (M_A * N_A)^1.52)
end

# Calculate the annual atmospheric concentrations
if is_first(t)
v.M_AA[t] = p.M_AA_2005
v.N_AA[t] = p.N_AA_2005

v.pre_f = f(p.M_pre, p.N_pre) # only need to calculate this once; used in each timestep below
else
v.M_AA[t] = (1 - p.delta_ch4) * v.M_AA[t - 1] + p.delta_ch4 * p.M_pre + p.gamma_ch4 * p.E_CH4A[t - 1]
v.N_AA[t] = (1 - p.delta_n2o) * v.N_AA[t - 1] + p.delta_n2o * p.N_pre + p.gamma_n2o * (28/44) * p.E_N2OA[t - 1]
end
if is_first(t)
v.M_AA[t] = p.M_AA_2005
v.N_AA[t] = p.N_AA_2005

v.pre_f = f(p.M_pre, p.N_pre) # only need to calculate this once; used in each timestep below
else
v.M_AA[t] = (1 - p.delta_ch4) * v.M_AA[t - 1] + p.delta_ch4 * p.M_pre + p.gamma_ch4 * p.E_CH4A[t - 1]
v.N_AA[t] = (1 - p.delta_n2o) * v.N_AA[t - 1] + p.delta_n2o * p.N_pre + p.gamma_n2o * (28 / 44) * p.E_N2OA[t - 1]
end

# Calculate the annual forcing
v.F_CH4A[t] = p.ipcc_adj * (p.alpha_ch4 * (sqrt(v.M_AA[t]) - sqrt(p.M_pre)) - (f(v.M_AA[t], p.N_pre) - v.pre_f))
v.F_N2OA[t] = p.alpha_n2o * (sqrt(v.N_AA[t]) - sqrt(p.N_pre)) - (f(p.M_pre, v.N_AA[t]) - v.pre_f)
v.F_CH4A[t] = p.ipcc_adj * (p.alpha_ch4 * (sqrt(v.M_AA[t]) - sqrt(p.M_pre)) - (f(v.M_AA[t], p.N_pre) - v.pre_f))
v.F_N2OA[t] = p.alpha_n2o * (sqrt(v.N_AA[t]) - sqrt(p.N_pre)) - (f(p.M_pre, v.N_AA[t]) - v.pre_f)

# calculate the decadal forcing at the end
if is_last(t)
v.F_CH4[:] = mean(reshape(v.F_CH4A[:], 10, length(d.decades)), dims=1)
v.F_N2O[:] = mean(reshape(v.F_N2OA[:], 10, length(d.decades)), dims=1)
end

if is_last(t)
v.F_CH4[:] = mean(reshape(v.F_CH4A[:], 10, length(d.decades)), dims=1)
v.F_N2O[:] = mean(reshape(v.F_N2OA[:], 10, length(d.decades)), dims=1)
end

end
end


Expand Down
Loading

0 comments on commit cbe8336

Please sign in to comment.