Skip to content

Commit

Permalink
Format files using DocumentFormat
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidanthoff
davidanthoff authored Oct 14, 2022
1 parent 25caa2e commit a10f6fc
Show file tree
Hide file tree
Showing 6 changed files with 89 additions and 89 deletions.
24 changes: 12 additions & 12 deletions src/MimiMAGICC.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,26 +9,26 @@ include("components/rf_ch4.jl")
include("components/rf_ch4h2o.jl")
include("components/rf_o3.jl")

function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2300)
function get_model(; rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2300)

# ---------------------------------------------
# Load and clean up necessary data.
# ---------------------------------------------

# Load RCP emissions and concentration scenario values (RCP options = "RCP26" or "RCP85").
rcp_emissions = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario*"_EMISSIONS.csv"), skiplines_begin=36))
rcp_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario*"_CONCENTRATIONS.csv"), skiplines_begin=37))
rcp_emissions = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario * "_EMISSIONS.csv"), skiplines_begin=36))
rcp_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario * "_CONCENTRATIONS.csv"), skiplines_begin=37))

# Load temperature scenario (mean response from SNEASY across 100,000 posterior parameter samples).
rcp_temperature = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario*"_temperature_sneasy_1765_2300.csv")))
rcp_temperature = DataFrame(load(joinpath(@__DIR__, "..", "data", rcp_scenario * "_temperature_sneasy_1765_2300.csv")))

# Load fossil CO₂ emissions scenario to calculate historic tropospheric O₃ forcing (this scenario is specific to how MAGICC calculates this forcing).
foss_hist_for_O₃ = DataFrame(load(joinpath(@__DIR__, "..", "data", "FOSS_HIST_magicc.csv")))
foss_hist_for_O₃ = DataFrame(load(joinpath(@__DIR__, "..", "data", "FOSS_HIST_magicc.csv")))

# Find indices for start and end years to crop scenarios to correct time frame.
start_index, end_index = findall((in)([start_year, end_year]), rcp_emissions.YEARS)
rcp_emissions = rcp_emissions[start_index:end_index, :]
rcp_temperature = rcp_temperature[start_index:end_index, :]
rcp_emissions = rcp_emissions[start_index:end_index, :]
rcp_temperature = rcp_temperature[start_index:end_index, :]
foss_hist_for_O₃ = foss_hist_for_O₃[start_index:end_index, :]

# Set initial CH₄ and N₂O concentrations to RCP 1765 values.
Expand Down Expand Up @@ -64,8 +64,8 @@ function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year
set_param!(m, :ch4_cycle, :AVOC, -0.000315)
set_param!(m, :ch4_cycle, :TAUINIT, 7.73)
set_param!(m, :ch4_cycle, :SCH4, -0.32)
set_param!(m, :ch4_cycle, :TAUSOIL, 160.)
set_param!(m, :ch4_cycle, :TAUSTRAT, 120.)
set_param!(m, :ch4_cycle, :TAUSOIL, 160.0)
set_param!(m, :ch4_cycle, :TAUSTRAT, 120.0)
set_param!(m, :ch4_cycle, :GAM, -1.0)
set_param!(m, :ch4_cycle, :CH4_natural, 266.5)
set_param!(m, :ch4_cycle, :fffrac, 0.18)
Expand Down Expand Up @@ -97,11 +97,11 @@ function get_model(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year
# ---------------------------------------------
# Create connections between Mimi components.
# ---------------------------------------------
connect_param!(m, :rf_ch4, :CH₄, :ch4_cycle, :CH₄)
connect_param!(m, :rf_ch4, :CH₄, :ch4_cycle, :CH₄)
connect_param!(m, :rf_ch4h2o, :CH₄, :ch4_cycle, :CH₄)
connect_param!(m, :rf_o3, :CH₄, :ch4_cycle, :CH₄)
connect_param!(m, :rf_o3, :CH₄, :ch4_cycle, :CH₄)

return(m)
return (m)
end

end #module
82 changes: 41 additions & 41 deletions src/components/ch4_cycle.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,34 +4,34 @@

@defcomp ch4_cycle begin

BBCH4 = Parameter() # Tg to ppb unit conversion between emissions and concentrations Note:(1 Teragram = 1 Megatonne).
ANOX = Parameter() # Hydroxyl lifetime coefficient for nitrogen oxides.
ACO = Parameter() # Hydroxyl lifetime coefficient for carbon monoxide.
AVOC = Parameter() # Hydroxyl lifetime coefficient for non-methane volatile organic compounds.
SCH4 = Parameter() # Hydroxyl lifetime coefficient for methane.
TAUINIT = Parameter() # Initial methane tropospheric lifetime due to hydroxyl abundance (years).
TAUSOIL = Parameter() # Methane soil sink lifetime (years).
TAUSTRAT = Parameter() # Methane stratospheric lifetime (years).
GAM = Parameter() # Hydroxyl lifetime parameter.
CH4_natural = Parameter() # Natural methane emissions (Mt yr⁻¹).
fffrac = Parameter() # Fossil fuel fraction of total methane emissions.
CH₄_0 = Parameter() # Atmospheric methane pre-industrial concentration (ppb).
temperature = Parameter(index=[time]) # Global surface temperature anomaly (°C).
CH4_emissions = Parameter(index=[time]) # Global anthropogenic methane emissions (Mt yr⁻¹).
NOx_emissions = Parameter(index=[time]) # Global nitrogen oxides emissions (Mt yr⁻¹).
CO_emissions = Parameter(index=[time]) # Global carbon monoxide emissions (Mt yr⁻¹).
NMVOC_emissions = Parameter(index=[time]) # Global non-methane volatile organic compound emissions (Mt yr⁻¹).

TAU_OH = Variable(index=[time]) # Methane tropospheric lifetime due to hydroxyl abundance (years).
CH₄ = Variable(index=[time]) # Atmospheric methane concetration for current period (ppb).
emeth = Variable(index=[time]) # Methane that has been oxidized to carbon dioxide.
BBCH4 = Parameter() # Tg to ppb unit conversion between emissions and concentrations Note:(1 Teragram = 1 Megatonne).
ANOX = Parameter() # Hydroxyl lifetime coefficient for nitrogen oxides.
ACO = Parameter() # Hydroxyl lifetime coefficient for carbon monoxide.
AVOC = Parameter() # Hydroxyl lifetime coefficient for non-methane volatile organic compounds.
SCH4 = Parameter() # Hydroxyl lifetime coefficient for methane.
TAUINIT = Parameter() # Initial methane tropospheric lifetime due to hydroxyl abundance (years).
TAUSOIL = Parameter() # Methane soil sink lifetime (years).
TAUSTRAT = Parameter() # Methane stratospheric lifetime (years).
GAM = Parameter() # Hydroxyl lifetime parameter.
CH4_natural = Parameter() # Natural methane emissions (Mt yr⁻¹).
fffrac = Parameter() # Fossil fuel fraction of total methane emissions.
CH₄_0 = Parameter() # Atmospheric methane pre-industrial concentration (ppb).
temperature = Parameter(index=[time]) # Global surface temperature anomaly (°C).
CH4_emissions = Parameter(index=[time]) # Global anthropogenic methane emissions (Mt yr⁻¹).
NOx_emissions = Parameter(index=[time]) # Global nitrogen oxides emissions (Mt yr⁻¹).
CO_emissions = Parameter(index=[time]) # Global carbon monoxide emissions (Mt yr⁻¹).
NMVOC_emissions = Parameter(index=[time]) # Global non-methane volatile organic compound emissions (Mt yr⁻¹).

TAU_OH = Variable(index=[time]) # Methane tropospheric lifetime due to hydroxyl abundance (years).
CH₄ = Variable(index=[time]) # Atmospheric methane concetration for current period (ppb).
emeth = Variable(index=[time]) # Methane that has been oxidized to carbon dioxide.


function run_timestep(p, v, d, t)

# Set initial conditions. NOTE: The CH₄ cycle equations require temperature values from the previous two timesteps (t-1 and t-2), so the cycle switches
# on starting in timestep 3. Further, timestep 1 CH₄ concentrations do not affect the results below (but timestep 2 CH₄ concentrations do). We
# therefore treat period 2 CH₄ concentrations as an initial (and uncertain) condition.
# on starting in timestep 3. Further, timestep 1 CH₄ concentrations do not affect the results below (but timestep 2 CH₄ concentrations do). We
# therefore treat period 2 CH₄ concentrations as an initial (and uncertain) condition.
if t <= TimestepIndex(2)
v.CH₄[TimestepIndex(2)] = p.CH₄_0
v.emeth[t] = 0.0
Expand All @@ -48,14 +48,14 @@

# Calculate change in emissions for NOx, CO, and NMVOC (relative to period 3 when CH₄ cycle model engages due to t-1 and t-2 terms above).
DENOX = p.NOx_emissions[t] - p.NOx_emissions[TimestepIndex(3)]
DECO = p.CO_emissions[t] - p.CO_emissions[TimestepIndex(3)]
DECO = p.CO_emissions[t] - p.CO_emissions[TimestepIndex(3)]
DEVOC = p.NMVOC_emissions[t] - p.NMVOC_emissions[TimestepIndex(3)]

# Add natural CH₄ emissions to anthropogenic CH₄ emissions.
Total_CH4emissions = p.CH4_emissions[t] + p.CH4_natural

# Calcualte TAUOTHER term (based on stratospheric and soil sink lifetimes).
TAUOTHER = 1.0 / (1.0/p.TAUSOIL + 1.0/p.TAUSTRAT)
TAUOTHER = 1.0 / (1.0 / p.TAUSOIL + 1.0 / p.TAUSTRAT)

##########################################################
# Iterate to Calculate OH lifetime and CH₄ concentration.
Expand All @@ -69,7 +69,7 @@
B00 = p.CH₄_0 * p.BBCH4

# Account for other gases on OH abundance and tropospheric methane lifetime.
AAA = exp(p.GAM * (p.ANOX*DENOX + p.ACO*DECO + p.AVOC*DEVOC))
AAA = exp(p.GAM * (p.ANOX * DENOX + p.ACO * DECO + p.AVOC * DEVOC))

# Multiply two parameters together.
X = p.GAM * p.SCH4
Expand All @@ -83,39 +83,39 @@
# FIRST ITERATION
#----------------------------------------
BBAR = B
TAUBAR = U * (BBAR/B00) ^ X
TAUBAR = p.TAUINIT / (p.TAUINIT/TAUBAR + 0.0316 * D_Temp)
DB1 = Total_CH4emissions - (BBAR/TAUBAR) - (BBAR/TAUOTHER)
TAUBAR = U * (BBAR / B00)^X
TAUBAR = p.TAUINIT / (p.TAUINIT / TAUBAR + 0.0316 * D_Temp)
DB1 = Total_CH4emissions - (BBAR / TAUBAR) - (BBAR / TAUOTHER)
B1 = B + DB1

#----------------------------------------
# SECOND ITERATION
#----------------------------------------
BBAR = (B + B1) / 2.0
TAUBAR = U * (BBAR/B00) ^ X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB1/B)
TAUBAR = p.TAUINIT / (p.TAUINIT/TAUBAR + 0.0316 * D_Temp)
DB2 = Total_CH4emissions - (BBAR/TAUBAR) - (BBAR/TAUOTHER)
TAUBAR = U * (BBAR / B00)^X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB1 / B)
TAUBAR = p.TAUINIT / (p.TAUINIT / TAUBAR + 0.0316 * D_Temp)
DB2 = Total_CH4emissions - (BBAR / TAUBAR) - (BBAR / TAUOTHER)
B2 = B + DB2

#----------------------------------------
# THIRD ITERATION
#----------------------------------------
BBAR = (B + B2) / 2.0
TAUBAR = U * (BBAR/B00) ^ X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB2/B)
TAUBAR = p.TAUINIT / (p.TAUINIT/TAUBAR + 0.0316 * D_Temp)
DB3 = Total_CH4emissions - (BBAR/TAUBAR) - (BBAR/TAUOTHER)
TAUBAR = U * (BBAR / B00)^X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB2 / B)
TAUBAR = p.TAUINIT / (p.TAUINIT / TAUBAR + 0.0316 * D_Temp)
DB3 = Total_CH4emissions - (BBAR / TAUBAR) - (BBAR / TAUOTHER)
B3 = B + DB3

#----------------------------------------
# FOURTH ITERATION
#----------------------------------------
BBAR = (B + B3) / 2.0
TAUBAR = U * (BBAR/B00) ^ X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB3/B)
TAUBAR = p.TAUINIT / (p.TAUINIT/TAUBAR + 0.0316 * D_Temp)
DB4 = Total_CH4emissions - (BBAR/TAUBAR) - (BBAR/TAUOTHER)
TAUBAR = U * (BBAR / B00)^X
TAUBAR = TAUBAR * (1.0 - 0.5 * X * DB3 / B)
TAUBAR = p.TAUINIT / (p.TAUINIT / TAUBAR + 0.0316 * D_Temp)
DB4 = Total_CH4emissions - (BBAR / TAUBAR) - (BBAR / TAUOTHER)
B4 = B + DB4

# Calculate CH₄ tropospheric lifetime.
Expand Down
12 changes: 6 additions & 6 deletions src/components/rf_ch4.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,24 +4,24 @@

# Create a function 'interact' to account for the overlap in methane and nitrous oxide in their radiative absoprtion bands.
function interact(M, N)
d = 1.0 + 0.636 * ((M * N)/(10^6))^0.75 + 0.007 * (M/(10^3)) * ((M * N)/(10^6))^1.52
d = 1.0 + 0.636 * ((M * N) / (10^6))^0.75 + 0.007 * (M / (10^3)) * ((M * N) / (10^6))^1.52
return 0.47 * log(d)
end

@defcomp rf_ch4 begin
CH₄_0 = Parameter() # Pre-industrial atmospheric methane concentration (ppb).
N₂O_0 = Parameter() # Pre-industrial atmospheric nitrous oxide concentration (ppb).
CH₄_0 = Parameter() # Pre-industrial atmospheric methane concentration (ppb).
N₂O_0 = Parameter() # Pre-industrial atmospheric nitrous oxide concentration (ppb).
scale_CH₄ = Parameter() # Scaling factor for uncertainty in methane radiative forcing.
CH₄ = Parameter(index=[time]) # Atmospheric methane concentration (ppb).
CH₄ = Parameter(index=[time]) # Atmospheric methane concentration (ppb).

QMeth = Variable(index=[time]) # Direct forcing from atmospheric methane concentrations (Wm⁻²).
QMeth = Variable(index=[time]) # Direct forcing from atmospheric methane concentrations (Wm⁻²).


function run_timestep(p, v, d, t)

if is_first(t)
# Set initial forcings to 0.
v.QMeth[t] = 0.0
v.QMeth[t] = 0.0
else
# Direct radiative forcing from atmospheric CH₄ concentrations.
v.QMeth[t] = (0.036 * (sqrt(p.CH₄[t]) - sqrt(p.CH₄_0)) - (interact(p.CH₄[t], p.N₂O_0) - interact(p.CH₄_0, p.N₂O_0))) * p.scale_CH₄
Expand Down
8 changes: 4 additions & 4 deletions src/components/rf_ch4h2o.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,11 @@

@defcomp rf_ch4h2o begin

CH₄_0 = Parameter() # Pre-industrial atmospheric methane concentration (ppb).
STRATH2O = Parameter() # Share of direct methane radiative forcing that represents stratospheric water vapor forcing from methane oxidation.
CH₄ = Parameter(index=[time]) # Atmospheric methane concentration (ppb).
CH₄_0 = Parameter() # Pre-industrial atmospheric methane concentration (ppb).
STRATH2O = Parameter() # Share of direct methane radiative forcing that represents stratospheric water vapor forcing from methane oxidation.
CH₄ = Parameter(index=[time]) # Atmospheric methane concentration (ppb).

QCH4H2O = Variable(index=[time]) # Radiative forcing for stratoshperic water vapor from methane oxidation (Wm⁻²). rf_O₃ = Variable(index=[time]) # Total radiative forcing from tropospheric ozone (Wm⁻²).
QCH4H2O = Variable(index=[time]) # Radiative forcing for stratoshperic water vapor from methane oxidation (Wm⁻²). rf_O₃ = Variable(index=[time]) # Total radiative forcing from tropospheric ozone (Wm⁻²).


function run_timestep(p, v, d, t)
Expand Down
Loading

0 comments on commit a10f6fc

Please sign in to comment.