NewFUND with 2-Box Climate, Updated Energy, CO2 Fertilization

In [1]:
using Mimi
using MimiFUND  #downloaded 2020-04-05
modelStartYear = 1950
finalYear = 2200
finalYearPlot = 2200
inflator = 1.5359  #$2018 by GDP deflator
m = MimiFUND.get_model();
In [2]:
@defcomp climatedynamics2 begin
    # 2-box climate model
    # Total radiative forcing
    radforc = Parameter(index=[time])
    # Average global temperature
    temp = Variable(index=[time])
    # Average global deep ocean temperature
    tempocean = Variable(index=[time])
    # mixed layer depth 70 m    
    # Deep layer depth 1200 m
    deeplayer = Parameter()
    # Heat capacity 4180000.0 J/C/m3        
    # CO2 forcing 3.7 W/m2    
    climatesensitivity = Parameter()    
    fpb = Variable()
    fpg = Variable()
    heatcapdeep = Variable()
    dTocean = Variable(index=[time])
    dTemp = Variable(index=[time])
    #
        function init(p, v, d)
        # ocean fraction = 70%
        # heat capacity mixed layer =0.7 * 70 * 4180000 = 204820000.0
        # heat capacity deep layer
        v.heatcapdeep = 0.7 * p.deeplayer * 4180000.0
        #timestep year in seconds
        #year = 60 * 60 * 24 * 365 = 31536000
        # feedback parameter B
            if p.climatesensitivity < 0.4
            v.fpb = 9.25
            else
            v.fpb = 3.7 / p.climatesensitivity
            end
        # deep layer feedback parameter g
        v.fpg = v.fpb / (0.78 - 0.063 * p.climatesensitivity) - v.fpb
        end  
        function run_timestep(p, v, d, t)        
            if is_first(t)
                if p.climatesensitivity > 2.0 
                v.temp[t] = 0.2 
                elseif p.climatesensitivity >= 0.4
                v.temp[t] = 0.2 * p.climatesensitivity /2.0
                else
                v.dTemp[t] = 0.0
                error("climate sensitivity must be zero or >= 0.4")
                end
            v.tempocean[t] = v.temp[t] /8.0
            else     
            v.dTocean[t] = v.fpg * (v.temp[t-1] - v.tempocean[t-1]) / v.heatcapdeep * 31536000.
            v.tempocean[t] = v.tempocean[t-1] + v.dTocean[t]
                if p.climatesensitivity < 0.4
                v.dTemp[t] = 0.0
                else
                v.dTemp[t] = (p.radforc[t] - v.fpb * v.temp[t-1] - v.fpg * (v.temp[t-1] - v.tempocean[t-1]))/ 204820000.0 * 31536000
                end
            v.temp[t] = v.temp[t - 1] + v.dTemp[t]
            end 
        end
end
In [3]:
replace!(m, :climatedynamics => climatedynamics2)
set_param!(m, :deeplayer, 1200)
In [4]:
@defcomp impactcooling2 begin
    # cooling impacts based on absolute regional temperatures at population centroid and emperical changes
    # from FUND-energy-impacts-per-GDP-vs-TempV2.xlsx  V2 has increases space cooling for commercial buildings
    regions = Index()
    cooling = Variable(index=[time,regions])    
    population = Parameter(index=[time,regions])    
    income = Parameter(index=[time,regions])
    heel = Parameter(default = 0.8)    
    cumaeei = Parameter(index=[time,regions])
    rtcfpop = Parameter(index=[regions])
    regtemp1950 = Parameter(index=[regions])
    temp = Parameter(index=[time,regions])    
    cool1950 = Parameter(index=[regions])
    regtemp = Variable(index=[time,regions])
    alpha1950 = Variable(index=[regions])
    alpha = Variable(index=[time,regions])
    ypc =  Variable(index=[time,regions])
    # 
        function init(p, v, d)
            for r in d.regions
                v.alpha1950[r] = 0.44107 .* p.regtemp1950[r] ^2 + 2.71998 .* p.regtemp1950[r] .+ -22.49838
            end            
        end
    function run_timestep(p, v, d, t)  
        if is_first(t)
            for r in d.regions
            v.regtemp[t, r] = p.regtemp1950[r] 
            v.cooling[t, r] = p.cool1950[r]
            end
        else 
            for r in d.regions
            v.ypc[t, r] = p.income[t, r] / p.population[t, r]  
            v.regtemp[t, r] = v.regtemp[t-1, r] + (p.temp[t,r] - p.temp[t-1,r]) .* p.rtcfpop[r]            
            v.alpha[t, r] =  0.44107 .* v.regtemp[t, r] ^2 + 2.71998 .* v.regtemp[t, r] .+ -22.49838                       
            v.cooling[t, r] = v.cooling[TimestepIndex(1), r] .+ (v.alpha1950[r] * p.population[TimestepIndex(1), r] .- v.alpha[t, r] .* p.population[t, r]) .* (v.ypc[t, r] ./ 38.769)^p.heel .* (p.cumaeei[t, r] ./ 0.92108) ./ 1000 
            end
        end                    
    end
end
In [5]:
replace!(m, :impactcooling => impactcooling2)
set_param!(m, :rtcfpop, [1.117,1.290,1.091,1.052,1.034,1.173,1.214,1.194,0.859,0.907,0.715,0.856,1.142,1.192,0.872,0.752])
set_param!(m, :regtemp1950, [10.981,5.551,4.527,12.065,13.427,4.598,3.425,14.612,25.916,23.676,25.520,26.736,14.664,14.958,25.209,26.699])
set_param!(m, :cool1950, [-0.37,-0.02,-0.81,-0.02,0.00,-0.02,-0.13,-0.02,-0.01,-0.05,-0.02,-0.03,-0.08,-0.03,-0.07,0.00])
disconnect_param!(m, :impactcooling, :henl)
disconnect_param!(m, :impactcooling, :hebm)
disconnect_param!(m, :impactcooling, :gdp90)
disconnect_param!(m, :impactcooling, :pop90)
Out[5]:
true
In [6]:
@defcomp impactheating2 begin
    # heating impacts based on absolute regional temperatures at population centroid and emperical changes
    #File \My Files\Climate\FUND\PaperKBG\FUND-energy-impacts-per-GDP-vs-Temp.xlsx
    regions = Index()
    heating = Variable(index=[time,regions])    
    population = Parameter(index=[time,regions])    
    income = Parameter(index=[time,regions])
    heel = Parameter(default = 0.8)    
    cumaeei = Parameter(index=[time,regions])
    rtcfpop = Parameter(index=[regions])
    regtemp1950 = Parameter(index=[regions])
    temp = Parameter(index=[time,regions])    
    heat1950 = Parameter(index=[regions])
    regtemp = Variable(index=[time,regions])
    alpha1950 = Variable(index=[regions])
    alpha = Variable(index=[time,regions])
    ypc =  Variable(index=[time,regions])   
    # 
    function init(p, v, d)
        for r in d.regions
            if p.regtemp1950[r] > 26
            v.alpha1950[r] = 0
            else
            v.alpha1950[r] = 0.231297 .* p.regtemp1950[r] ^2 + -25.69714 .* p.regtemp1950[r] .+ 511.770 
            end
        end
    end
    function run_timestep(p, v, d, t)  
        if is_first(t)
            for r in d.regions
            v.regtemp[t, r] = p.regtemp1950[r] 
            v.heating[t, r] = p.heat1950[r]
            end
        else           
            for r in d.regions
            v.ypc[t, r] = p.income[t, r] / p.population[t, r]  #$k/capaita
            v.regtemp[t, r] = v.regtemp[t-1, r] + (p.temp[t,r] - p.temp[t-1,r]) .* p.rtcfpop[r]            
                if v.regtemp[t, r] > 26
                   v.alpha[t, r] = 0
                else
                   v.alpha[t, r] =  0.231297 .* v.regtemp[t, r] ^2 + -25.69714 .* v.regtemp[t, r] .+ 511.770                         
                end
             v.heating[t, r] = v.heating[TimestepIndex(1), r] .+ (v.alpha1950[r] .* p.population[TimestepIndex(1), r] .- v.alpha[t, r] .* p.population[t, r]) .* (v.ypc[t, r] ./ 38.769)^p.heel .* (p.cumaeei[t, r] ./ 0.92108) ./ 1000 
            end
        end                    
    end
end
In [7]:
replace!(m, :impactheating => impactheating2)
set_param!(m, :heat1950, [2.10,0.10,1.48,0.44,0.05,0.12,0.89,0.44,0.04,0.22,0.06,0.07,0.31,0.0,0.0,0.01])
set_param!(m, :rtcfpop, [1.117,1.290,1.091,1.052,1.034,1.173,1.214,1.194,0.859,0.907,0.715,0.856,1.142,1.192,0.872,0.752])
set_param!(m, :regtemp1950, [10.981,5.551,4.527,12.065,13.427,4.598,3.425,14.612,25.916,23.676,25.520,26.736,14.664,14.958,25.209,26.699])
disconnect_param!(m, :impactheating, :henl)
disconnect_param!(m, :impactheating, :hebm)
disconnect_param!(m, :impactheating, :gdp90)
disconnect_param!(m, :impactheating, :pop90)
Out[7]:
true
In [8]:
#Updated Energy Only
DiscRate = 0.03
AGcbm = [0.089, 0.0402, 0.1541, 0.2319, 0.1048, 0.0952, 0.0671, 0.0943, 0.1641, 0.0596, 0.058, 0.0845, 0.1921, 0.0727, 0.0505, 0.2377]
set_param!(m, :agcbm, AGcbm)
ECS = zeros(9)
ECS[1] = 0.6
Scc = zeros(9)
set_param!(m, :climatesensitivity, ECS[1])
scc = MimiFUND.compute_scco2(m, year = 2020, eta = 0., prtp = DiscRate, equity_weights = false)*inflator
Scc[1] = scc
    for i = 2:9
    ECS[i] = 0.4 + i * 0.2 
    update_param!(m, :climatesensitivity, ECS[i])
    scc = MimiFUND.compute_scco2(m, year = 2020, eta = 0., prtp = DiscRate, equity_weights = false)*inflator
    Scc[i] = scc
    end
println("ECS C   Scc     Probability")
println("Deg.C   USD/tCO2")
Probability = zeros(9)
Probability = [0.09595,0.21974,0.26946,0.20133,0.11202,0.05394,0.02410,0.01157,0.01189]
for i = 1:9
    println(round(ECS[i];digits=1),"     ",round(Scc[i];digits=2),"    ",Probability[i])
end
ECS C   Scc     Probability
Deg.C   USD/tCO2
0.6     -9.25    0.09595
0.8     -8.94    0.21974
1.0     -8.52    0.26946
1.2     -7.98    0.20133
1.4     -7.34    0.11202
1.6     -6.62    0.05394
1.8     -5.82    0.0241
2.0     -4.95    0.01157
2.2     -4.03    0.01189
In [9]:
TotScc = 0.0
for x = 1:9
    TotScc += Probability[x] * Scc[x]
end
println("TotScc =  ",round(TotScc;digits=4))
TotScc =  -8.1787
In [10]:
update_param!(m, :climatesensitivity, 1.1)
scc = MimiFUND.compute_scco2(m, year = 2020, eta = 0., prtp = DiscRate, equity_weights = false)*inflator
Out[10]:
-8.259669233912017
In [11]:
# K Dayaratna, R. McKitrick & P. Michaels, 2020
AGcbm = [0.089, 0.0402, 0.1541, 0.2319, 0.1048, 0.0952, 0.0671, 0.0943, 0.1641, 0.0596, 0.058, 0.0845, 0.1921, 0.0727, 0.0505, 0.2377]
AGcbm_new = AGcbm .* 1.3
set_param!(m, :agcbm, AGcbm_new)
run(m)
AGcbm = m[:impactagriculture, :agcbm]
for j = 1:16
print(round(AGcbm[j];digits=6),"  ")
end
0.1157  0.05226  0.20033  0.30147  0.13624  0.12376  0.08723  0.12259  0.21333  0.07748  0.0754  0.10985  0.24973  0.09451  0.06565  0.30901  
In [12]:
#Updated Energy and CO2 Fertilization
    for i = 1:9
    ECS[i] = 0.4 + i * 0.2 
    update_param!(m, :climatesensitivity, ECS[i])
    scc = MimiFUND.compute_scco2(m, year = 2020, eta = 0., prtp = DiscRate, equity_weights = false)*inflator
    Scc[i] = scc
    end
println("ECS C   Scc     Probability")
println("Deg.C   USD/tCO2")
Probability = zeros(9)
Probability = [0.09595,0.21974,0.26946,0.20133,0.11202,0.05394,0.02410,0.01157,0.01189]
for i = 1:9
    println(round(ECS[i];digits=1),"     ",round(Scc[i];digits=2),"    ",Probability[i])
end
ECS C   Scc     Probability
Deg.C   USD/tCO2
0.6     -12.34    0.09595
0.8     -12.05    0.21974
1.0     -11.82    0.26946
1.2     -11.12    0.20133
1.4     -10.49    0.11202
1.6     -9.77    0.05394
1.8     -8.97    0.0241
2.0     -8.1    0.01157
2.2     -7.17    0.01189
In [13]:
TotScc = 0.0
for x = 1:9
    TotScc += Probability[x] * Scc[x]
end
println("TotScc =  ",round(TotScc;digits=4))
TotScc =  -11.3541
In [14]:
# Use Ramsey Discounting equivalent to USA 3%. 0.03 = US growth rate/capita 2020 X eta + prtp = 0.0171 X eta + 0.01
# eta = (0.03 - 0.01)/0.0171 = 1.1696  prtp = 0.010
# For discount rate 5%, 0.05 = 0.0171 X eta + prtp.  eta = (0.05 - 0.015)/0.0171 = 2.0468  prtp = 0.015
if DiscRate == 0.03
    Eta = 1.1696
    Prtp = 0.010
elseif DiscRate == 0.05
    Eta = 2.0468
    Prtp = 0.015
    else error("select DiscRate 0.03 or 0.05")
end
    for i = 1:9
    ECS[i] = 0.4 + i * 0.2 
    update_param!(m, :climatesensitivity, ECS[i])
    scc = MimiFUND.compute_scco2(m, year = 2020, eta = Eta, prtp = Prtp, equity_weights = false)*inflator
    Scc[i] = scc
    end
println("ECS C   Scc     Probability")
println("Deg.C   USD/tCO2")
Probability = zeros(9)
Probability = [0.09595,0.21974,0.26946,0.20133,0.11202,0.05394,0.02410,0.01157,0.01189]
for i = 1:9
    println(round(ECS[i];digits=1),"     ",round(Scc[i];digits=2),"    ",Probability[i])
end
ECS C   Scc     Probability
Deg.C   USD/tCO2
0.6     -10.72    0.09595
0.8     -10.43    0.21974
1.0     -10.24    0.26946
1.2     -9.55    0.20133
1.4     -8.97    0.11202
1.6     -8.31    0.05394
1.8     -7.58    0.0241
2.0     -6.8    0.01157
2.2     -5.96    0.01189
In [15]:
TotScc = 0.0
for x = 1:9
    TotScc += Probability[x] * Scc[x]
end
println("TotScc =  ",round(TotScc;digits=4))
TotScc =  -9.7874
In [16]:
update_param!(m, :climatesensitivity, 1.1)
run(m)
In [17]:
Temp = zeros(finalYear - modelStartYear + 1)
for i = 1:251
Temp[i] = m[:climatedynamics, :temp][i]
    print(round(Temp[i];digits=7),"  ") 
    end
0.11  0.1099354  0.1118745  0.1144678  0.1170608  0.1195494  0.1218261  0.1239217  0.1263324  0.1293047  0.1335134  0.1370941  0.1402199  0.1438243  0.1476131  0.1516419  0.1558533  0.1598504  0.1638905  0.1686417  0.1746508  0.1808338  0.1855197  0.189324  0.1930552  0.1966815  0.2002299  0.2039224  0.208352  0.2134355  0.217745  0.2214433  0.2255851  0.229871  0.2345535  0.2397608  0.2451295  0.2503445  0.2557473  0.2616956  0.2683007  0.2751355  0.2816562  0.2884387  0.2957495  0.3034116  0.3110302  0.3166913  0.3205001  0.3240797  0.3273817  0.3316147  0.337384  0.3436434  0.3501588  0.3576063  0.3667744  0.3771295  0.387907  0.3988859  0.4094548  0.4190098  0.4279379  0.4368333  0.4458472  0.4550328  0.4643697  0.4738519  0.4834735  0.493234  0.5032793  0.5137931  0.5246891  0.5357782  0.547013  0.5583759  0.5698647  0.581461  0.5931494  0.604927  0.6170216  0.6297976  0.6430788  0.6565685  0.6701843  0.6838772  0.6976913  0.711604  0.7256143  0.7397189  0.7538187  0.7678674  0.7819083  0.7960135  0.8102104  0.8244908  0.8388689  0.8533003  0.8677839  0.8823224  0.8966925  0.9107223  0.9245586  0.9384306  0.9523999  0.9664835  0.9806576  0.9948932  1.0091857  1.0235453  1.0379033  1.0522282  1.0665357  1.08088  1.0952731  1.1097002  1.1242056  1.1388027  1.1534962  1.1682856  1.1831317  1.1979944  1.2128712  1.2277926  1.2427675  1.25779  1.2728893  1.2880145  1.3031564  1.3183091  1.3335564  1.3491829  1.3650885  1.3811108  1.3972016  1.4133065  1.4294804  1.4456853  1.4618969  1.4781218  1.4941743  1.5099645  1.52556  1.5411044  1.5566383  1.5721257  1.5876494  1.6031628  1.6186494  1.6341058  1.6493601  1.664326  1.679091  1.6937808  1.7084242  1.7230231  1.7375734  1.7520604  1.7664804  1.7808183  1.7950769  1.8092418  1.8233132  1.8372795  1.851139  1.8648859  1.878517  1.8920227  1.9054018  1.918646  1.9317497  1.9447162  1.9575365  1.9702104  1.9827292  1.9950912  2.0072917  2.0193303  2.0312037  2.0429084  2.0544443  2.0658061  2.0769909  2.0879997  2.0988278  2.1094766  2.1199476  2.1302364  2.1403409  2.1502568  2.1599884  2.1695392  2.1789029  2.1880821  2.1970784  2.2058876  2.2145173  2.2229548  2.2312082  2.2392819  2.2471777  2.2548894  2.2624192  2.2697718  2.2769508  2.2839622  2.2908034  2.2974799  2.3039965  2.3103452  2.31654  2.3225797  2.3284658  2.3342028  2.3397934  2.3452369  2.3505397  2.3557008  2.3607257  2.3656159  2.3703737  2.3749964  2.379499  2.383878  2.3881373  2.3922752  2.3962981  2.4002078  2.4040065  2.4076988  2.4112867  2.4147755  2.4181629  2.4214563  2.4246552  2.427757  2.4307712  2.4336966  2.4365392  2.4393005  2.4419827  2.4445907  2.4471178  2.449574  2.451959  2.4542789  2.4565316  2.4587215  2.4608463  2.4629139  2.4649261  
In [18]:
REGdrycost = zeros(finalYear - modelStartYear + 1,16)
REGprotcost = zeros(finalYear - modelStartYear + 1,16)
REGentercost = zeros(finalYear - modelStartYear + 1,16)
REGleavecost = zeros(finalYear - modelStartYear + 1,16)
REGwetcost = zeros(finalYear - modelStartYear + 1,16)
REGsealevel = zeros(finalYear - modelStartYear + 1,16)
REGwater = zeros(finalYear - modelStartYear + 1,16)
REGforest = zeros(finalYear - modelStartYear + 1,16)
REGspecies = zeros(finalYear - modelStartYear + 1,16)
REGheating = zeros(finalYear - modelStartYear + 1,16)
REGcooling = zeros(finalYear - modelStartYear + 1,16)
REGenergy = zeros(finalYear - modelStartYear + 1,16)
REGagcost = zeros(finalYear - modelStartYear + 1,16)
REGhurrdam = zeros(finalYear - modelStartYear + 1,16)
REGetstorms = zeros(finalYear - modelStartYear + 1,16)
REGstorms = zeros(finalYear - modelStartYear + 1,16)
REGdeadcost = zeros(finalYear - modelStartYear + 1,16)
REGmorbcost = zeros(finalYear - modelStartYear + 1,16)
REGhealth = zeros(finalYear - modelStartYear + 1,16)
REGecosystems = zeros(finalYear - modelStartYear + 1,16)
REGeloss = zeros(finalYear - modelStartYear + 1,16)
REGsloss = zeros(finalYear - modelStartYear + 1,16)
REGloss = zeros(finalYear - modelStartYear + 1,16)
REGincome = zeros(finalYear - modelStartYear + 1,16)
for i = 2:(finalYear - modelStartYear + 1) 
    for j = 1:16
      REGdrycost[i,j]  = -m[:impactaggregation, :drycost][i,j] * inflator 
      REGprotcost[i,j] = -m[:impactaggregation, :protcost][i,j] * inflator
      REGentercost[i,j] = -m[:impactaggregation, :entercost][i,j] * inflator 
      REGleavecost[i,j] = -m[:impactaggregation, :leavecost][i,j] * inflator
      REGwetcost[i,j] = -m[:impactaggregation, :wetcost][i,j] * inflator
      REGloss[i,j] = -m[:impactaggregation, :loss][i,j]/1000000000 *inflator
      REGwater[i,j] = m[:impactaggregation, :water][i,j] * inflator 
      REGforest[i,j] = m[:impactaggregation, :forests][i,j] * inflator
      REGspecies[i,j] = -m[:impactaggregation, :species][i,j] * inflator  
      REGheating[i,j] = m[:impactaggregation, :heating][i,j] * inflator
      REGcooling[i,j] = m[:impactaggregation, :cooling][i,j] * inflator 
      REGagcost[i,j] = m[:impactaggregation, :agcost][i,j] * inflator
      REGhurrdam[i,j] = -m[:impactaggregation, :hurrdam][i,j] * inflator
      REGetstorms[i,j] = -m[:impactaggregation, :extratropicalstormsdam][i,j] * inflator  
      REGdeadcost[i,j] = -m[:impactaggregation, :deadcost][i,j] * inflator 
      REGmorbcost[i,j] = -m[:impactaggregation, :morbcost][i,j] * inflator
      REGeloss[i,j] = -m[:impactaggregation, :eloss][i,j] * inflator        
      REGsloss[i,j] = -m[:impactaggregation, :sloss][i,j] * inflator
      REGincome[i,j] = m[:socioeconomic, :income][i,j] * inflator
    end
end
REGsealevel = (REGdrycost .+ REGprotcost .+ REGentercost .+ REGleavecost .+ REGwetcost)
REGhealth = (REGdeadcost .+ REGmorbcost)
REGstorms = (REGhurrdam .+ REGetstorms)
REGenergy = (REGheating .+ REGcooling)
REGecosystems = (REGforest .+ REGspecies);
In [19]:
GLdrycost = zeros(finalYear - modelStartYear + 1)
GLprotcost = zeros(finalYear - modelStartYear + 1)
GLentercost = zeros(finalYear - modelStartYear + 1)
GLleavecost = zeros(finalYear - modelStartYear + 1)
GLwetcost = zeros(finalYear - modelStartYear + 1)
GLsealevel = zeros(finalYear - modelStartYear + 1)
GLcomp_SeaLevel = zeros(finalYear - modelStartYear + 1,5)
GLwater = zeros(finalYear - modelStartYear + 1)
GLforest = zeros(finalYear - modelStartYear + 1)
GLspecies = zeros(finalYear - modelStartYear + 1)
GLheating = zeros(finalYear - modelStartYear + 1)
GLcooling = zeros(finalYear - modelStartYear + 1)
GLenergy = zeros(finalYear - modelStartYear + 1)
GLagcost = zeros(finalYear - modelStartYear + 1)
GLhurrdam = zeros(finalYear - modelStartYear + 1)
GLetstorms = zeros(finalYear - modelStartYear + 1)
GLstorms = zeros(finalYear - modelStartYear + 1)
GLdeadcost = zeros(finalYear - modelStartYear + 1)
GLmorbcost = zeros(finalYear - modelStartYear + 1)
GLhealth = zeros(finalYear - modelStartYear + 1)
GLecosystems = zeros(finalYear - modelStartYear + 1)
GLeloss = zeros(finalYear - modelStartYear + 1) 
GLsloss = zeros(finalYear - modelStartYear + 1)
GLincome = zeros(finalYear - modelStartYear + 1) 
for i = 2:(finalYear - modelStartYear + 1) 
    for j = 1:16
 GLdrycost[i] += REGdrycost[i,j]
 GLprotcost[i] += REGprotcost[i,j]
 GLentercost[i] += REGentercost[i,j]   
 GLleavecost[i] += REGleavecost[i,j] 
 GLwetcost[i] += REGwetcost[i,j]
 GLsealevel[i] += REGsealevel[1,j]       
 GLwater[i] += REGwater[i,j]       
 GLforest[i] += REGforest[i,j]   
 GLspecies[i] += REGspecies[i,j]
 GLheating[i] += REGheating[i,j]   
 GLcooling[i] += REGcooling[i,j]  
 GLagcost[i] += REGagcost[i,j]
 GLhurrdam[i] += REGhurrdam[i,j]
 GLetstorms[i] += REGetstorms[i,j]   
 GLdeadcost[i] += REGdeadcost[i,j]
 GLmorbcost[i] += REGmorbcost[i,j]
 GLeloss[i]   += REGeloss[i,j]     
 GLsloss[i]   += REGsloss[i,j]       
    end
GLincome[i] = m[:socioeconomic, :globalincome][i]  *inflator 
end
GLsealevel = (GLdrycost .+ GLprotcost .+ GLentercost .+ GLleavecost .+ GLwetcost)
GLhealth = (GLdeadcost .+ GLmorbcost)
GLstorms = (GLhurrdam .+ GLetstorms)
GLenergy = (GLheating .+ GLcooling)
GLecosystems = (GLforest .+ GLspecies)
GLcomp = zeros(finalYear - modelStartYear + 1,7)
for i = 1:(finalYear - modelStartYear + 1)
    GLcomp[i,1] = GLstorms[i]
    GLcomp[i,2] = GLagcost[i]
    GLcomp[i,3] = GLwater[i]
    GLcomp[i,4] = GLsealevel[i]
    GLcomp[i,5] = GLhealth[i]
    GLcomp[i,6] = GLenergy[i]
    GLcomp[i,7] = GLecosystems[i]
end
In [20]:
for j = 1:7
println(round(GLcomp[251,j];digits=3))
end
-37.487
15437.068
-569.191
-123.972
-98.65
-1325.565
-232.872
In [21]:
using Plots
Plots.plot(modelStartYear:finalYear, GLcomp/1000, grid=true, legend=true, lw=3, 
    label=["Storms" "Agriculture" "Water" "Sea Level" "Health" "Energy"  "Ecosystems"],
    xlims=(1950,finalYearPlot), ylims=(-2,16), xticks=(1950:25:finalYearPlot), yticks=(-2:2:16),
    title="Global Cost Impacts by Component; ECS = 1.1 °C", xlabel="year", ylabel="trillions of 2018 USD")
Out[21]:
In [22]:
GLcompFileName = "NewGLcompECS1.1.csv"
# Open files for writing. They will appear in the same directory as this .ipynb file
GLloss = zeros(finalYear - modelStartYear + 1)
GLlossch = zeros(finalYear - modelStartYear + 1)
GLlossch = GLeloss .+ GLsloss
file1 = open(GLcompFileName, "w") 
header1 = "Global Impacts by Component 2018 USD billions with ECS = 1.1 C"
header2 = "Year,  Storms,  Agriculture,  Water,  Sea Level,  Health , Energy, Ecosystems, Total, Total ch, GDP"
println(file1, header1)
println(file1, header2)
for i = 2:251 
    for j = 1:7
    GLloss[i] += GLcomp[i,j]
    end
end
for i = 2:251
    year = (1949 + i)
    print(file1, year); print(file1, ", ")
    for j = 1:7        
        print(file1, GLcomp[i, j]); print(file1, ", ")             
    end
    print(file1, GLloss[i]); print(file1, ", "); print(file1, GLlossch[i]); print(file1, ", ")
    println(file1, GLincome[i])
end
close(file1)
In [23]:
ecs = m[:climatedynamics, :climatesensitivity]
print(ecs)
1.1
In [24]:
USCooing = m[:impactcooling, :cooling][251,1] * inflator
Out[24]:
-57.38926681347387
In [25]:
USHeating = m[:impactheating, :heating][2,1] * inflator
Out[25]:
2.8248911223740176
In [26]:
GLstormsGDP = zeros(finalYear - modelStartYear + 1)
GLagcostGDP = zeros(finalYear - modelStartYear + 1)
GLwaterGDP = zeros(finalYear - modelStartYear + 1)
GLsealevelGDP = zeros(finalYear - modelStartYear + 1)
GLhealthGDP = zeros(finalYear - modelStartYear + 1)
GLenergyGDP = zeros(finalYear - modelStartYear + 1)
GLecosystemsGDP = zeros(finalYear - modelStartYear + 1)
GLcompGDP = zeros(finalYear - modelStartYear + 1,7)
GLstormsGDP = GLstorms ./ GLincome .*100
GLagcostGDP = GLagcost ./ GLincome .*100
GLwaterGDP = GLwater ./ GLincome .*100
GLsealevelGDP = GLsealevel ./ GLincome .*100
GLhealthGDP = GLhealth ./ GLincome .*100
GLenergyGDP = GLenergy ./ GLincome .*100
GLecosystemsGDP = GLecosystems ./ GLincome .*100
for i = 2:(finalYear - modelStartYear + 1)
    GLcompGDP[i,1] = GLstormsGDP[i]
    GLcompGDP[i,2] = GLagcostGDP[i]
    GLcompGDP[i,3] = GLwaterGDP[i]
    GLcompGDP[i,4] = GLsealevelGDP[i]
    GLcompGDP[i,5] = GLhealthGDP[i]
    GLcompGDP[i,6] = GLenergyGDP[i]
    GLcompGDP[i,7] = GLecosystemsGDP[i]
end
In [27]:
Plots.plot(modelStartYear:finalYear, GLcompGDP, grid=true, legend=true, lw=3, 
    label=["Storms" "Agriculture" "Water" "Sea Level" "Health" "Energy"  "Ecosystems"],
    xlims=(1950,finalYearPlot), ylims=(-0.6,2), xticks=(1950:25:finalYearPlot), yticks=(-0.6:0.2:2),
    title="Global Cost Impacts Per GDP by Component; ECS = 1.1 °C", xlabel="year", ylabel="Percent of GDP")
Out[27]:
In [28]:
CDNdrycost = zeros(finalYear - modelStartYear + 1)
CDNprotcost = zeros(finalYear - modelStartYear + 1)
CDNentercost = zeros(finalYear - modelStartYear + 1)
CDNleavecost = zeros(finalYear - modelStartYear + 1)
CDNwetcost = zeros(finalYear - modelStartYear + 1)
CDNSeaLevel = zeros(finalYear - modelStartYear + 1)
CDNcomp_SeaLevel = zeros(finalYear - modelStartYear + 1,5)
CDNwater = zeros(finalYear - modelStartYear + 1)
CDNforest = zeros(finalYear - modelStartYear + 1)
CDNspecies = zeros(finalYear - modelStartYear + 1)
CDNheating = zeros(finalYear - modelStartYear + 1)
CDNcooling = zeros(finalYear - modelStartYear + 1)
CDNenergy = zeros(finalYear - modelStartYear + 1)
CDNagcost = zeros(finalYear - modelStartYear + 1)
CDNhurrdam = zeros(finalYear - modelStartYear + 1)
CDNetstorms = zeros(finalYear - modelStartYear + 1)
CDNstorms = zeros(finalYear - modelStartYear + 1)
CDNdeadcost = zeros(finalYear - modelStartYear + 1)
CDNmorbcost = zeros(finalYear - modelStartYear + 1)
CDNhealth = zeros(finalYear - modelStartYear + 1)
CDNecosystems = zeros(finalYear - modelStartYear + 1)
CDNeloss = zeros(finalYear - modelStartYear + 1)
CDNsloss = zeros(finalYear - modelStartYear + 1)
CDNincome = zeros(finalYear - modelStartYear + 1)
for i = 1:(finalYear - modelStartYear + 1) 
    CDNdrycost[i] = REGdrycost[i,2]
    CDNprotcost[i] = REGprotcost[i,2]
    CDNentercost[i] = REGentercost[i,2]
    CDNleavecost[i] = REGleavecost[i,2]
    CDNwetcost[i] = REGwetcost[i,2]
    CDNeloss[i] = REGeloss[i,2]
    CDNsloss[i] = REGsloss[i,2]
    CDNcomp_SeaLevel[i,1] = CDNdrycost[i]
    CDNcomp_SeaLevel[i,2] = CDNprotcost[i]
    CDNcomp_SeaLevel[i,3] = CDNentercost[i]
    CDNcomp_SeaLevel[i,4] = CDNleavecost[i]
    CDNcomp_SeaLevel[i,5] = CDNwetcost[i]
    CDNdeadcost[i] = REGdeadcost[i,2]
    CDNmorbcost[i] = REGmorbcost[i,2] 
    CDNhurrdam[i] = REGhurrdam[i,2] 
    CDNetstorms[i] = REGetstorms[i,2]
    CDNheating[i] = REGheating[i,2]
    CDNcooling[i] = REGcooling[i,2]
    CDNwater[i] = REGwater[i,2]  
    CDNforest[i] = REGforest[i,2]
    CDNspecies[i] = REGspecies[i,2]
    CDNagcost[i] = REGagcost[i,2] 
    CDNincome[i] = REGincome[i,2]
end
CDNsealevel = (CDNdrycost .+ CDNprotcost .+ CDNentercost .+ CDNleavecost .+ CDNwetcost)
CDNhealth = (CDNdeadcost .+ CDNmorbcost)
CDNstorms = (CDNhurrdam .+ CDNetstorms)
CDNenergy = (CDNheating .+ CDNcooling)
CDNecosystems = (CDNforest .+ CDNspecies)
CDNcomp = zeros(finalYear - modelStartYear + 1,7)
for i = 2:(finalYear - modelStartYear + 1)
    CDNcomp[i,1] = CDNstorms[i]
    CDNcomp[i,2] = CDNagcost[i]
    CDNcomp[i,3] = CDNwater[i]
    CDNcomp[i,4] = CDNsealevel[i]
    CDNcomp[i,5] = CDNhealth[i]
    CDNcomp[i,6] = CDNenergy[i]
    CDNcomp[i,7] = CDNecosystems[i]
end
In [29]:
print(151+1949,"  ",round(CDNcooling[151];digits=3))
2100  -1.45
In [30]:
for i = 1:10
    print(i+1949,"  ",round(CDNcooling[i];digits=3),"  ",round(CDNheating[i];digits=3),"  ",round(CDNenergy[i];digits=3))
    println()
end
1950  0.0  0.0  0.0
1951  -0.032  0.099  0.068
1952  -0.033  0.041  0.008
1953  -0.034  -0.021  -0.055
1954  -0.035  -0.087  -0.122
1955  -0.037  -0.158  -0.195
1956  -0.038  -0.234  -0.272
1957  -0.039  -0.315  -0.354
1958  -0.041  -0.401  -0.442
1959  -0.043  -0.493  -0.536
In [31]:
Plots.plot(modelStartYear:finalYear, CDNcomp, grid=true, legend=true, lw=3, 
    label=["Storms" "Agriculture" "Water" "Sea Level" "Health" "Energy"  "Ecosystems"],
    xlims=(1950,finalYearPlot), ylims=(-20,40), xticks=(1950:25:finalYearPlot), yticks=(-20:10:40),
    title="Canada Cost Impacts by Component; ECS = 1.1 °C", xlabel="year", ylabel="billions of 2018 USD")
Out[31]:
In [32]:
CDNcompFileName = "NewCDNcompECS1.1.csv"
# Open files for writing. They will appear in the same directory as this .ipynb file
CDNloss = zeros(finalYear - modelStartYear + 1)
CDNlossch = zeros(finalYear - modelStartYear + 1)
CDNlossch = CDNeloss .+ CDNsloss
file2 = open(CDNcompFileName, "w") 
header1 = "Canadian Impacts by Component 2018 USD billions with ECS = 1.1 C"
header2 = "Year,  Storms,  Agriculture,  Water,  Sea Level,  Health , Energy, Ecosystems, Total, Total ch, GDP"
println(file2, header1)
println(file2, header2)
for i = 2:251 
    CDNloss[i] = REGloss[i,2]
    year = (1949 + i)
    print(file2, year); print(file2, ", ")
    for j = 1:7        
        print(file2, CDNcomp[i, j]); print(file2, ", ")             
    end
    print(file2, CDNloss[i],", ",CDNlossch[i],", ")
    println(file2, CDNincome[i])
end
close(file2)
In [33]:
CDNstormsGDP = zeros(finalYear - modelStartYear + 1)
CDNagcostGDP = zeros(finalYear - modelStartYear + 1)
CDNwaterGDP = zeros(finalYear - modelStartYear + 1)
CDNsealevelGDP = zeros(finalYear - modelStartYear + 1)
CDNhealthGDP = zeros(finalYear - modelStartYear + 1)
CDNenergyGDP = zeros(finalYear - modelStartYear + 1)
CDNecosystemsGDP = zeros(finalYear - modelStartYear + 1)
CDNcompGDP = zeros(finalYear - modelStartYear + 1,7)
CDNstormsGDP = CDNstorms ./ CDNincome .*100
CDNagcostGDP = CDNagcost ./ CDNincome .*100
CDNwaterGDP = CDNwater ./ CDNincome .*100
CDNsealevelGDP = CDNsealevel ./ CDNincome .*100
CDNhealthGDP = CDNhealth ./ CDNincome .*100
CDNenergyGDP = CDNenergy ./ CDNincome .*100
CDNecosystemsGDP = CDNecosystems ./ CDNincome .*100
for i = 2:(finalYear - modelStartYear + 1)
    CDNcompGDP[i,1] = CDNstormsGDP[i]
    CDNcompGDP[i,2] = CDNagcostGDP[i]
    CDNcompGDP[i,3] = CDNwaterGDP[i]
    CDNcompGDP[i,4] = CDNsealevelGDP[i]
    CDNcompGDP[i,5] = CDNhealthGDP[i]
    CDNcompGDP[i,6] = CDNenergyGDP[i]
    CDNcompGDP[i,7] = CDNecosystemsGDP[i]
end
In [34]:
Plots.plot(modelStartYear:finalYear, CDNcompGDP, grid=true, legend=true, lw=3, 
    label=["Storms" "Agriculture" "Water" "Sea Level" "Health" "Energy"  "Ecosystems"],
    xlims=(1950,finalYearPlot), ylims=(-0.7,0.6), xticks=(1950:25:finalYearPlot), yticks=(-0.7:0.1:0.6),
    title="Canadian Cost Impacts Per GDP by Component; ECS = 1.1 °C", xlabel="year", ylabel="Percent of GDP")
Out[34]: