Exact two layer solution

Project.toml

@@ -1,6 +1,6 @@
 name = "ClimateMARGO"
 uuid = "d3f62095-a717-45bf-aadc-ac9dfc258fa6"
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

README.md

@@ -32,7 +32,7 @@
 
 The MARGO model is described in full in an accompanying manuscript (not yet peer-reviewed, and currently being revised), available to all at [EarthArXiv.org/5bgyc](https://eartharxiv.org/5bgyc/).
 
-Try out the MARGO model by working through the Binder tutorial right in your browser (you don't need to download anything – just click the [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/hdrake/ClimateMARGO.jl/master?filepath=examples%2Ftutorial.ipynb) button and, once the tutorial loads, click on the code cells and press ``Enter`` to run them).
+Try out the MARGO model by working through the Binder tutorial right in your browser (you don't need to download anything – just click the [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/hdrake/ClimateMARGO.jl/master?filepath=examples%2Ftutorial.ipynb) button and, once the tutorial loads, click on the code cells and press ``Enter`` to run them) or the ClimateWidgets dashboard [![Binder](https://mybinder.org/badge_logo.svg)](https://camo.githubusercontent.com/483bae47a175c24dfbfc57390edd8b6982ac5fb3/68747470733a2f2f6d7962696e6465722e6f72672f62616467655f6c6f676f2e737667).
 
 ClimateMARGO.jl is currently in alpha testing and basic model documentation is slowly being added. Anyone interested in helping develop the model post an Issue here or contact Henri F. Drake directly (henrifdrake `at` gmail.com), until explicit guidelines for contributing to the model are posted at a later date.
 

configurations/default_params.json

@@ -1 +1 @@
-{"name":"default","domain":{"dt":5.0,"present_year":2020.0,"initial_year":2020.0,"final_year":2200.0},"economics":{"E0":100.0,"γ":0.02,"β":0.0022222222222222222,"ρ":0.01,"Finf":8.5,"mitigate_cost":0.034,"remove_cost":13.0,"geoeng_cost":0.046,"adapt_cost":4.5,"mitigate_init":0.1,"remove_init":0.0,"geoeng_init":0.0,"adapt_init":null,"baseline_emissions":[7.5,8.4375,9.375,10.3125,11.25,12.1875,13.125,14.0625,15.0,15.9375,16.875,17.8125,18.75,19.6875,20.625,21.5625,22.5,20.25,18.0,15.75,13.5,11.25,9.0,6.75,4.5,2.25,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],"extra_CO₂":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]},"physics":{"c0":460.0,"T0":1.1,"a":4.977297891066924,"B":1.13,"Cd":106.0,"κ":0.73,"r":0.5}}
+{"name":"default","grid":{"dt":5.0,"present_year":2020.0,"initial_year":2020.0,"final_year":2200.0},"economics":{"E0":100.0,"γ":0.02,"β":0.0022222222222222222,"ρ":0.01,"Finf":8.5,"mitigate_cost":0.034,"remove_cost":13.0,"geoeng_cost":0.046,"adapt_cost":4.5,"mitigate_init":0.1,"remove_init":0.0,"geoeng_init":0.0,"adapt_init":null,"baseline_emissions":[7.5,8.4375,9.375,10.3125,11.25,12.1875,13.125,14.0625,15.0,15.9375,16.875,17.8125,18.75,19.6875,20.625,21.5625,22.5,20.25,18.0,15.75,13.5,11.25,9.0,6.75,4.5,2.25,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],"extra_CO₂":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]},"physics":{"c0":460.0,"T0":1.1,"r":0.5,"a":4.977297891066924,"λ":1.13,"Cu":7.3,"Cd":106.0,"κ":0.73,"ϵ":1.0}}

src/ClimateMARGO.jl

@@ -8,8 +8,8 @@ include("Models/Models.jl")
 include("Utils/Utils.jl")
 include("Diagnostics/Diagnostics.jl")
 include("Optimization/Optimization.jl")
-include("IO/IO.jl")
-include("PolicyResponse/PolicyResponse.jl")
-include("Plotting/Plotting.jl")
+#include("IO/IO.jl")
+#include("PolicyResponse/PolicyResponse.jl")
+#include("Plotting/Plotting.jl")
 
 end

src/Diagnostics/Diagnostics.jl

@@ -4,10 +4,10 @@ using ClimateMARGO.Models
 using ClimateMARGO.Utils
 
 export
-    t, future_mask,
+    t, future_mask, allow_control, deferred,
     ramp_emissions, emissions, effective_emissions,
-    c, F, Tslow, Tfast, T,
-    τd, B, F2x, ECS,
+    c, F, T_mode, T,
+    af, τf, as, τs, calc_λ, F2x, ECS,
     discount, f, E,
     damage, cost, benefit,
     net_benefit, net_present_cost, net_present_benefit

src/Diagnostics/carbon.jl

@@ -3,7 +3,7 @@
 #
 # See below link for 2020 initial condition:
 # https://www.eea.europa.eu/data-and-maps/indicators/atmospheric-greenhouse-gas-concentrations-6/assessment-1
-function ramp_emissions(t, q0::Float64, n::Float64, t1::Float64, t2::Float64)
+function ramp_emissions(t, q0::Real, n::Real, t1::Real, t2::Real)
     t0 = t[1]
     Δt0 = t1 - t0
     Δt1 = t2 - t1
@@ -15,15 +15,15 @@ function ramp_emissions(t, q0::Float64, n::Float64, t1::Float64, t2::Float64)
     q[t .> t2] .= 0.
     return q
 end
-function ramp_emissions(dom::Domain)
-    return ramp_emissions(t(dom), 7.5, 3., 2100., 2150.)
+function ramp_emissions(grid::Grid)
+    return ramp_emissions(t(grid), 7.5, 3., 2100., 2150.)
 end
 
 emissions(q, M) = q .* (1. .- M)
 function emissions(m::ClimateModel; M=false)
     return emissions(
         m.economics.baseline_emissions,
-        m.controls.mitigate .* (1. .- .~future_mask(m) * ~M)
+        m.controls.mitigate .* allow_control(m, M),
     )
 end
 
@@ -34,8 +34,8 @@ function effective_emissions(m; M=false, R=false)
     return effective_emissions(
         m.physics.r,
         m.economics.baseline_emissions,
-        m.controls.mitigate .* (1. .- .~future_mask(m) * ~M),
-        m.controls.remove .* (1. .- .~future_mask(m) * ~R)
+        m.controls.mitigate .* allow_control(m, M),
+        m.controls.remove .* allow_control(m, R)
     )
 end
 
@@ -47,6 +47,6 @@ function c(m; M=false, R=false)
     return c(
         m.physics.c0,
         effective_emissions(m, M=M, R=R),
-        m.domain.dt
+        m.grid.dt
     )
 end

src/Diagnostics/cost_benefit.jl

@@ -4,7 +4,7 @@ E(t, E0, γ) = E0 * (1. .+ γ).^(t .- t[1])
 E(m) = E(t(m), m.economics.E0, m.economics.γ)
 
 discount(t, ρ, tp) = .~future_mask(t, tp) .* (1. .+ ρ) .^ (- (t .- tp))
-discount(m::ClimateModel) = discount(t(m), m.economics.ρ, m.domain.present_year)
+discount(m::ClimateModel) = discount(t(m), m.economics.ρ, m.grid.present_year)
 
 damage(β, E, T, A; discount=1.) = ((1. .- A) .* β .* E .* T.^2) .* discount
 
@@ -63,12 +63,12 @@ function net_present_cost(
     )
     return net_present_cost(
         cost(m, discounting=discounting, M=M, R=R, G=G, A=A),
-        m.domain.dt
+        m.grid.dt
     )
 end
 
 net_present_benefit(net_benefit, dt) = sum(net_benefit*dt)
 net_present_benefit(m::ClimateModel; discounting=true, M=false, R=false, G=false, A=false) = net_present_benefit(
     net_benefit(m, discounting=discounting, M=M, R=R, G=G, A=A),
-    m.domain.dt
+    m.grid.dt
 )

src/Diagnostics/energy_balance.jl

@@ -12,102 +12,101 @@ function F(m; M=false, R=false, G=false, F0=0.)
     return F(
         m.physics.a, m.physics.c0, m.economics.Finf,
         c(m, M=M, R=R),
-        m.controls.geoeng .* (1. .- .~future_mask(m) * ~G),
+        m.controls.geoeng .* allow_control(m, G),
         F0=F0
     )
 end
 
 F2x(a::Float64) = a*log(2)
 F2x(m::ClimateModel) = F2x(m.physics.a)
 
-ECS(a, B) = F2x(a)/B
-ECS(params::ClimateModelParameters) = ECS(params.physics.a, m.physics.B)
-ECS(m::ClimateModel) = ECS(m.physics.a, m.physics.B)
+ECS(a, λ) = F2x(a)/λ
+ECS(params::ClimateModelParameters) = ECS(params.physics.a, m.physics.λ)
+ECS(m::ClimateModel) = ECS(m.physics.a, m.physics.λ)
 
-calc_B(a::Float64, ECS::Float64) = F2x(a)/ECS
-calc_B(params::ClimateModelParameters; ECS=ECS(params)) = calc_B(params.physics.a, ECS)
-calc_B(m::ClimateModel; ECS=ECS(m)) = calc_B(m.physics.a, ECS)
+calc_λ(a::Float64, ECS::Float64) = F2x(a)/ECS
+calc_λ(params::ClimateModelParameters; ECS=ECS(params)) = calc_λ(params.physics.a, ECS)
+calc_λ(m::ClimateModel; ECS=ECS(m)) = calc_λ(m.physics.a, ECS)
 
-τd(Cd, B, κ) = (Cd/B) * (B+κ)/κ
-τd(phys::Physics) = τd(phys.Cd, phys.B, phys.κ)
-τd(m::ClimateModel) = τd(m.physics)
+## Two-layer model
+# Diagnostic constants
+b(λ, κ, ϵ, Cu, Cd) = (λ+ϵ*κ)/Cu + κ/Cd
+b(phys::Physics) = b(phys.λ, phys.κ, phys.ϵ, phys.Cu, phys.Cd)
 
-"""
-    T_fast(F, κ, B; A=0.)
-"""
-T_fast(F, κ, B; A=0.) = sqrt.(1. .- A) .* F/(κ + B)
+δ(b_, λ, κ, Cu, Cd) = b_^2 - 4(λ*κ)/(Cu*Cd)
+δ(phys::Physics) = δ(b(phys), phys.λ, phys.κ, phys.Cu, phys.Cd)
 
-"""
-    T_fast(m::ClimateModel; M=false, R=false, G=false, A=false)
-"""
-T_fast(m::ClimateModel; M=false, R=false, G=false, A=false) = T_fast(
-    F(m, M=M, R=R, G=G),
-    m.physics.κ,
-    m.physics.B,
-    A=m.controls.adapt .* (1. .- .~future_mask(m) * ~A)
-)
+τf(b_, δ_, λ, κ, Cu, Cd) = Cu*Cd/(2*λ*κ)*(b_ - √δ_)
+τf(phys::Physics) = τf(b(phys), δ(phys), phys.λ, phys.κ, phys.Cu, phys.Cd)
+
+τs(b_, δ_, λ, κ, Cu, Cd) = Cu*Cd/(2*λ*κ)*(b_ + √δ_)
+τs(phys::Physics) = τs(b(phys), δ(phys), phys.λ, phys.κ, phys.Cu, phys.Cd)
+
+bϕ(λ, κ, ϵ, Cu, Cd) = (λ+ϵ*κ)/Cu - κ/Cd
+bϕ(phys::Physics) = bϕ(phys.λ, phys.κ, phys.ϵ, phys.Cu, phys.Cd)
+
+ϕf(bϕ_, δ_, κ, ϵ, Cu) = Cu/(2*ϵ*κ)*(bϕ_ - √δ_)
+ϕf(phys::Physics) = ϕf(bϕ(phys), δ(phys), phys.κ, phys.ϵ, phys.Cu)
+
+ϕs(bϕ_, δ_, κ, ϵ, Cu) = Cu/(2*ϵ*κ)*(bϕ_ + √δ_)
+ϕs(phys::Physics) = ϕs(bϕ(phys), δ(phys), phys.κ, phys.ϵ, phys.Cu)
+
+af(ϕs_, ϕf_, τf_, λ, Cu) = ϕs_*τf_ /(Cu*(ϕs_ - ϕf_))*λ
+af(phys::Physics) = af(ϕs(phys), ϕf(phys), τf(phys), phys.λ, phys.Cu)
 
+as(ϕs_, ϕf_, τs_, λ, Cu) = -ϕf_*τs_ /(Cu*(ϕs_ - ϕf_))*λ
+as(phys::Physics) = as(ϕs(phys), ϕf(phys), τs(phys), phys.λ, phys.Cu)
+
+# Diagnostic variables
 """
-    T_slow(F, Cd, κ, B, t, dt; A=0.)
+    T_mode(F, λ, am, τm, t, dt)
 """
-function T_slow(F, Cd, κ, B, t, dt; A=0.)
-    τ = τd(Cd, κ, B)
-    return sqrt.(1. .- A) .* (
-        (κ/B) / (κ + B) *
-        exp.( - (t .- (t[1] - dt)) / τ) .*
-        cumsum( (exp.( (t .- (t[1] - dt)) / τ) / τ) .* F * dt)
+function T_mode(F, λ, am, τm, t, dt)
+    return (am/λ) .* (
+        exp.( - (t .- t[1]) / τm) .*
+        cumsum( exp.( (t .- t[1] ) / τm) .* F * dt ) / τm 
     )
 end
 
 """
-    T_slow(m::ClimateModel; M=false, R=false, G=false, A=false)
-"""
-T_slow(m::ClimateModel; M=false, R=false, G=false, A=false) = T_slow(
-    F(m, M=M, R=R, G=G),
-    m.physics.Cd,
-    m.physics.κ,
-    m.physics.B,
-    t(m),
-    m.domain.dt,
-    A=m.controls.adapt .* (1. .- .~future_mask(m) * ~A)
-)
-
-"""
-    T(T0, F, Cd, κ, B, t, dt; A=0.)
+    T(T0, F, λ, af_, τf_, as_, τs_, t, dt; A=0.)
 
 Returns the sum of the initial, fast mode, and slow mode temperature change.
 
 # Arguments
 - `T0::Float64`: warming relative to pre-industrial [°C]
 - `F::Array{Float64}`: change in radiative forcing since the initial time ``t_{0}`` [W/m``{2}``]
-- `Cd::Float64`: deep ocean heat capacity [W yr m``^{2}`` K``^{-1}``]
-- `κ::Float64`: ocean heat uptake rate [W m``^{2}`` K``^{-1}``]
-- `B::Float64`: feedback parameter [W m``^{2}`` K``^{-1}``]
+- `λ::Float64`: feedback parameter [W m``^{2}`` K``^{-1}``]
+- `af_::Float64`: fast mode fraction [1]
+- `τf_::Float64`: fast mode timescale [years]
+- `as_::Float64`: slow mode fraction [1]
+- `τs_::Float64`: slow mode timescale [years]
 - `t::Array{Float64}`: year [years]
 - `dt::Float64`: timestep [years]
-- `A::Float64`: Adaptation control [fraction]
+- `A::Float64`: adaptation control [fraction]
 
 """
-T(T0, F, Cd, κ, B, t, dt; A=0.) = sqrt.(1. .- A) .* (
-    T0 .+
-    T_fast(F, κ, B) .+
-    T_slow(F, Cd, κ, B, t, dt)
+function T(T0, F, λ, af_, τf_, as_, τs_, t, dt; A=0.)
+    T_fast_mode = T_mode(F, λ, af_, τf_, t, dt)
+    T_slow_mode = T_mode(F, λ, as_, τs_, t, dt)
+    return (T0 .+ T_fast_mode .+ T_slow_mode) .* sqrt.(1. .- A)
+end
+
+T(grid::Grid, phys::Physics, F; A=0.) = T(
+    phys.T0, F, phys.λ, af(phys), τf(phys), as(phys), τs(phys),
+    t(grid), grid.dt,
+    A=A
 )
 
 """
     T(m::ClimateModel; M=false, R=false, G=false, A=false)
-
 Returns the sum of the initial, fast mode, and slow mode temperature change,
 as diagnosed from `m` and modified by the climate controls activated by the
 Boolean kwargs.
 """
 T(m::ClimateModel; M=false, R=false, G=false, A=false) = T(
-    m.physics.T0,
+    m.grid,
+    m.physics,
     F(m, M=M, R=R, G=G),
-    m.physics.Cd,
-    m.physics.κ,
-    m.physics.B,
-    t(m),
-    m.domain.dt,
-    A=m.controls.adapt .* (1. .- .~future_mask(m) * ~A)
-)
+    A=m.controls.adapt .* allow_control(m, A)
+)

src/Diagnostics/utils.jl

@@ -1,7 +1,13 @@
 t(t0, tf, dt) = t0:dt:tf
-t(dom::Domain) = t(dom.initial_year,dom.final_year,dom.dt)
-t(params::ClimateModelParameters) = t(params.domain)
-t(m::ClimateModel) = t(m.domain)
+t(grid::Grid) = t(grid.initial_year,grid.final_year,grid.dt)
+t(params::ClimateModelParameters) = t(params.grid)
+t(m::ClimateModel) = t(m.grid)
 
 future_mask(t, present_year) = t .<= present_year
-future_mask(model) = future_mask(t(model), model.domain.present_year)
+future_mask(m::ClimateModel) = future_mask(t(m), model.grid.present_year)
+
+allow_control(t_, present_year, C) = (1. .- .~future_mask(t_, present_year) * ~C)
+allow_control(grid::Grid, C) = allow_control(t(grid), grid.present_year, C)
+allow_control(m::ClimateModel, C) = allow_control(m.grid, C)
+
+deferred(arr; n=1) = [zeros(n); arr[1:end-n]]

src/Models/Models.jl

@@ -1,32 +1,32 @@
 module Models
 
-export Domain, Physics, Controls, Economics, ClimateModelParameters, ClimateModel
+export Grid, Physics, Controls, Economics, ClimateModelParameters, ClimateModel
 
-include("domain.jl")
+include("grid.jl")
 include("physics.jl")
 include("controls.jl")
 include("economics.jl")
 
 """
-    ClimateModelParameters(name, domain::Domain, economics::Economics, physics::Physics)
+    ClimateModelParameters(name, grid::Grid, economics::Economics, physics::Physics)
 
 Create a named instance of the MARGO climate model parameters, which include
 economic input parameters (`economics`), physical climate parameters (`physics`),
-and climate control policies (`controls`) on some spatial-temporal grid (`domain`).
+and climate control policies (`controls`) on some spatial-temporal grid (`grid`).
 
 Use these to construct a [`ClimateModel`](@ref), which also contains the optimized 
 controls.
 """
 mutable struct ClimateModelParameters
     name::String
-    domain::Domain
+    grid::Grid
     economics::Economics
     physics::Physics
 end
 
 mutable struct ClimateModel
     name::String
-    domain::Domain
+    grid::Grid
     economics::Economics
     physics::Physics
     controls::Controls
@@ -42,14 +42,14 @@ the optimized [`Controls`](@ref). These can be computed using
 """
 ClimateModel(params::ClimateModelParameters, controls::Controls) = ClimateModel(
     params.name,
-    params.domain,
+    params.grid,
     params.economics,
     params.physics,
     controls
 )
 function ClimateModel(params::ClimateModelParameters)
-    dom = params.domain
-    t = collect(dom.initial_year:dom.dt:dom.final_year)
+    grid = params.grid
+    t = collect(grid.initial_year:grid.dt:grid.final_year)
     return ClimateModel(
         params,
         Controls(