Add damages for crops from methane emissions

data/MarketDamageAQ/Historical_MethaneEmissions.csv

@@ -0,0 +1,3 @@
+1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019
+,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+244.7851084,247.877912,251.2872804,257.3561574,258.9272937,261.8797995,264.9630002,267.5363652,269.0503088,273.2023324,271.8763774,267.5901506,268.9638341,269.9928233,274.3478234,276.9162625,280.9257592,283.0158911,288.7016316,294.7454026,294.8209681,293.1377146,293.0369675,292.352545,294.3197125,298.134385,301.2144647,300.8369132,297.3285842,298.536472,302.4784253,302.6381616,302.4973569,311.5903345,318.6138763,324.8513638,331.3708279,334.6652641,340.3105459,339.3839151,346.3654889,354.6388557,358.3110531,359.8600998,362.4142898,363.7954805,365.3850233,370.4264365,376.7292251,379.6269702

src/components/MarketDamageAQ_AsthmaERVisits.jl

@@ -0,0 +1,76 @@
+methane_on_AsthmaERVisits = myloadcsv("data/MarketDamageAQ/Methane_AsthmaERVisits.csv")
+historical_emissions = myloadcsv("data/MarketDamageAQ/Historical_MethaneEmissions.csv")
+
+function find_model_years(uu, model_years)
+    for i in 1:(length(model_years) - 1)
+        if model_years[i] <= uu && uu <= model_years[i+1]
+            return i, i+1
+        end
+    end
+    error("Year not found between model years.")
+end
+
+function addMarketDamageAQ_AsthmaERVisits(model::Model)
+    marketdamageaq_asthmaervisitscomp = add_comp!(model, MarketDamageAQ_AsthmaERVisits)
+    formatteddata = zeros(dim_count(model, :country), 50)
+    for t in 1:50
+        formatteddata[:, t] = readcountrydata_i_const(model, methane_on_AsthmaERVisits, :ISO3, Symbol(t))
+    end
+
+    # The year can be passed through model parameters
+    marketdamageaq_asthmaervisitscomp[:y_year_0] = 2015.
+    marketdamageaq_asthmaervisitscomp[:y_year] = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
+
+    marketdamageaq_asthmaervisitscomp[:asthma_er_visits_per_mton_ch4] = formatteddata
+    return marketdamageaq_asthmaervisitscomp
+end
+
+@defcomp MarketDamageAQ_AsthmaERVisits begin
+    country = Index()
+
+    # Incoming parameter: Global CH4 emissions from ch4emissions component
+    global_ch4_emissions = Parameter(index=[time], unit="Mtonne/year")
+
+    # Impact parameters: value per Mt methane emitted in year 1, year 2, ..., year 50
+    asthma_er_visits_per_mton_ch4 = Parameter(index=[country, 50], unit="\$/Mtonne")
+
+    # Component variable to store total AsthmaERVisits yield value
+    total_asthma_er_visits = Variable(index=[time, country], unit="\$")
+
+    # Define the y_year-0 parameter
+    y_year_0 = Parameter(unit="year")
+
+    # Define the y_year parameter (for interpolation)
+    y_year = Parameter(index=[10], unit="year")  # 10 time points
+
+
+    function run_timestep(p, v, d, t)
+        for c in d.country
+            v.total_asthma_er_visits[t, c] = 0
+
+            for tt in 1:50
+                uu = gettime(t) - tt + 1  # Current Year
+
+                # Using Historical Data
+                if uu >= 1970 && uu <= 2019
+                    v.total_asthma_er_visits[t, c] += historical_emissions[2, string(uu)] * p.asthma_er_visits_per_mton_ch4[c, tt]
+                elseif uu >= 2020 && uu <= 2300
+                    if uu in p.y_year
+                        idx = findfirst(x -> x == uu, p.y_year)
+                        println(idx)
+                        e = p.global_ch4_emissions[TimestepIndex(idx)]
+                        v.total_asthma_er_visits[t, c] += e * p.asthma_er_visits_per_mton_ch4[c, tt]
+                    else
+                        i1, i2 = find_model_years(uu, p.y_year)
+                        y1, y2 = p.y_year[i1], p.y_year[i2]
+                        e1 = p.global_ch4_emissions[TimestepIndex(i1)]
+                        e2 = p.global_ch4_emissions[TimestepIndex(i2)]
+                        interpolated_emission = e1 + (e2 - e1) * (uu - y1) / (y2 - y1)
+                        v.total_asthma_er_visits[t, c] += interpolated_emission * p.asthma_er_visits_per_mton_ch4[c, tt]
+                    end
+                end
+            end
+        end
+    end
+end
+

src/components/MarketDamageAQ_CropLoss.jl

@@ -0,0 +1,76 @@
+methane_on_crops = myloadcsv("data/MarketDamageAQ/Methane_CropLoss.csv")
+historical_emissions = myloadcsv("data/MarketDamageAQ/Historical_MethaneEmissions.csv")
+
+function find_model_years(uu, model_years)
+    for i in 1:(length(model_years) - 1)
+        if model_years[i] <= uu && uu <= model_years[i+1]
+            return i, i+1
+        end
+    end
+    error("Year not found between model years.")
+end
+
+function addMarketDamageAQ_CropLoss(model::Model)
+    marketdamageaq_croplosscomp = add_comp!(model, MarketDamageAQ_CropLoss)
+    formatteddata = zeros(dim_count(model, :country), 50)
+    for t in 1:50
+        formatteddata[:, t] = readcountrydata_i_const(model, methane_on_crops, :ISO3, Symbol(t))
+    end
+
+    # The year can be passed through model parameters
+    marketdamageaq_croplosscomp[:y_year_0] = 2015.
+    marketdamageaq_croplosscomp[:y_year] = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
+
+    marketdamageaq_croplosscomp[:crop_yield_value_per_mton_ch4] = formatteddata
+    return marketdamageaq_croplosscomp
+end
+
+@defcomp MarketDamageAQ_CropLoss begin
+    country = Index()
+
+    # Incoming parameter: Global CH4 emissions from ch4emissions component
+    global_ch4_emissions = Parameter(index=[time], unit="Mtonne/year")
+
+    # Impact parameters: value per Mt methane emitted in year 1, year 2, ..., year 50
+    crop_yield_value_per_mton_ch4 = Parameter(index=[country, 50], unit="\$/Mtonne")
+
+    # Component variable to store total crop yield value
+    total_crop_yield_value = Variable(index=[time, country], unit="\$")
+
+    # Define the y_year-0 parameter
+    y_year_0 = Parameter(unit="year")
+
+    # Define the y_year parameter (for interpolation)
+    y_year = Parameter(index=[10], unit="year")  # 10 time points
+
+
+    function run_timestep(p, v, d, t)
+        for c in d.country
+            v.total_crop_yield_value[t, c] = 0
+
+            for tt in 1:50
+                uu = gettime(t) - tt + 1  # Current Year
+
+                # Using Historical Data
+                if uu >= 1970 && uu <= 2019
+                    v.total_crop_yield_value[t, c] += historical_emissions[2, string(uu)] * p.crop_yield_value_per_mton_ch4[c, tt]
+                elseif uu >= 2020 && uu <= 2300
+                    if uu in p.y_year
+                        idx = findfirst(x -> x == uu, p.y_year)
+                        println(idx)
+                        e = p.global_ch4_emissions[TimestepIndex(idx)]
+                        v.total_crop_yield_value[t, c] += e * p.crop_yield_value_per_mton_ch4[c, tt]
+                    else
+                        i1, i2 = find_model_years(uu, p.y_year)
+                        y1, y2 = p.y_year[i1], p.y_year[i2]
+                        e1 = p.global_ch4_emissions[TimestepIndex(i1)]
+                        e2 = p.global_ch4_emissions[TimestepIndex(i2)]
+                        interpolated_emission = e1 + (e2 - e1) * (uu - y1) / (y2 - y1)
+                        v.total_crop_yield_value[t, c] += interpolated_emission * p.crop_yield_value_per_mton_ch4[c, tt]
+                    end
+                end
+            end
+        end
+    end
+end
+

src/components/MarketDamageAQ_LostWorkHours.jl

@@ -0,0 +1,76 @@
+methane_on_lostworkhours = myloadcsv("data/MarketDamageAQ/Methane_LostWorkHours.csv")
+historical_emissions = myloadcsv("data/MarketDamageAQ/Historical_MethaneEmissions.csv")
+
+function find_model_years(uu, model_years)
+    for i in 1:(length(model_years) - 1)
+        if model_years[i] <= uu && uu <= model_years[i+1]
+            return i, i+1
+        end
+    end
+    error("Year not found between model years.")
+end
+
+function addMarketDamageAQ_LostWorkHours(model::Model)
+    marketdamageaq_lostworkhours = add_comp!(model, MarketDamageAQ_LostWorkHours)
+    formatteddata = zeros(dim_count(model, :country), 50)
+    for t in 1:50
+        formatteddata[:, t] = readcountrydata_i_const(model, methane_on_lostworkhours, :ISO3, Symbol(t))
+    end
+
+    # The year can be passed through model parameters
+    marketdamageaq_lostworkhours[:y_year_0] = 2015.
+    marketdamageaq_lostworkhours[:y_year] = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
+
+    marketdamageaq_lostworkhours[:lost_work_per_mton_ch4] = formatteddata
+    return marketdamageaq_lostworkhours
+end
+
+@defcomp MarketDamageAQ_LostWorkHours begin
+    country = Index()
+
+    # Incoming parameter: Global CH4 emissions from ch4emissions component
+    global_ch4_emissions = Parameter(index=[time], unit="Mtonne/year")
+
+    # Impact parameters: value per Mt methane emitted in year 1, year 2, ..., year 50
+    lost_work_per_mton_ch4 = Parameter(index=[country, 50], unit="\$/Mtonne")
+
+    # Component variable to store total lost work hours yield value
+    total_methane_on_lostworkhours = Variable(index=[time, country], unit="\$")
+
+    # Define the y_year-0 parameter
+    y_year_0 = Parameter(unit="year")
+
+    # Define the y_year parameter (for interpolation)
+    y_year = Parameter(index=[10], unit="year")  # 10 time points
+
+
+    function run_timestep(p, v, d, t)
+        for c in d.country
+            v.total_methane_on_lostworkhours[t, c] = 0
+
+            for tt in 1:50
+                uu = gettime(t) - tt + 1  # Current Year
+
+                # Using Historical Data
+                if uu >= 1970 && uu <= 2019
+                    v.total_methane_on_lostworkhours[t, c] += historical_emissions[2, string(uu)] * p.lost_work_per_mton_ch4[c, tt]
+                elseif uu >= 2020 && uu <= 2300
+                    if uu in p.y_year
+                        idx = findfirst(x -> x == uu, p.y_year)
+                        println(idx)
+                        e = p.global_ch4_emissions[TimestepIndex(idx)]
+                        v.total_methane_on_lostworkhours[t, c] += e * p.lost_work_per_mton_ch4[c, tt]
+                    else
+                        i1, i2 = find_model_years(uu, p.y_year)
+                        y1, y2 = p.y_year[i1], p.y_year[i2]
+                        e1 = p.global_ch4_emissions[TimestepIndex(i1)]
+                        e2 = p.global_ch4_emissions[TimestepIndex(i2)]
+                        interpolated_emission = e1 + (e2 - e1) * (uu - y1) / (y2 - y1)
+                        v.total_methane_on_lostworkhours[t, c] += interpolated_emission * p.lost_work_per_mton_ch4[c, tt]
+                    end
+                end
+            end
+        end
+    end
+end
+

src/components/MarketDamageAQ_RespiratoryAdmissions.jl

@@ -0,0 +1,76 @@
+methane_on_respiratoryadmissions = myloadcsv("data/MarketDamageAQ/Methane_RespiratoryAdmissions.csv")
+historical_emissions = myloadcsv("data/MarketDamageAQ/Historical_MethaneEmissions.csv")
+
+function find_model_years(uu, model_years)
+    for i in 1:(length(model_years) - 1)
+        if model_years[i] <= uu && uu <= model_years[i+1]
+            return i, i+1
+        end
+    end
+    error("Year not found between model years.")
+end
+
+function addMarketDamageAQ_RespiratoryAdmissions(model::Model)
+    marketdamageaq_respiratoryadmissionscomp = add_comp!(model, MarketDamageAQ_RespiratoryAdmissions)
+    formatteddata = zeros(dim_count(model, :country), 50)
+    for t in 1:50
+        formatteddata[:, t] = readcountrydata_i_const(model, methane_on_respiratoryadmissions, :ISO3, Symbol(t))
+    end
+
+    # The year can be passed through model parameters
+    marketdamageaq_respiratoryadmissionscomp[:y_year_0] = 2015.
+    marketdamageaq_respiratoryadmissionscomp[:y_year] = [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
+
+    marketdamageaq_respiratoryadmissionscomp[:respiratory_admissions_per_mton_ch4] = formatteddata
+    return marketdamageaq_respiratoryadmissionscomp
+end
+
+@defcomp MarketDamageAQ_RespiratoryAdmissions begin
+    country = Index()
+
+    # Incoming parameter: Global CH4 emissions from ch4emissions component
+    global_ch4_emissions = Parameter(index=[time], unit="Mtonne/year")
+
+    # Impact parameters: value per Mt methane emitted in year 1, year 2, ..., year 50
+    respiratory_admissions_per_mton_ch4 = Parameter(index=[country, 50], unit="\$/Mtonne")
+
+    # Component variable to store total respiratory admissions value
+    methane_on_respiratoryadmissions = Variable(index=[time, country], unit="\$")
+
+    # Define the y_year-0 parameter
+    y_year_0 = Parameter(unit="year")
+
+    # Define the y_year parameter (for interpolation)
+    y_year = Parameter(index=[10], unit="year")  # 10 time points
+
+
+    function run_timestep(p, v, d, t)
+        for c in d.country
+            v.methane_on_respiratoryadmissions[t, c] = 0
+
+            for tt in 1:50
+                uu = gettime(t) - tt + 1  # Current Year
+
+                # Using Historical Data
+                if uu >= 1970 && uu <= 2019
+                    v.methane_on_respiratoryadmissions[t, c] += historical_emissions[2, string(uu)] * p.respiratory_admissions_per_mton_ch4[c, tt]
+                elseif uu >= 2020 && uu <= 2300
+                    if uu in p.y_year
+                        idx = findfirst(x -> x == uu, p.y_year)
+                        println(idx)
+                        e = p.global_ch4_emissions[TimestepIndex(idx)]
+                        v.methane_on_respiratoryadmissions[t, c] += e * p.respiratory_admissions_per_mton_ch4[c, tt]
+                    else
+                        i1, i2 = find_model_years(uu, p.y_year)
+                        y1, y2 = p.y_year[i1], p.y_year[i2]
+                        e1 = p.global_ch4_emissions[TimestepIndex(i1)]
+                        e2 = p.global_ch4_emissions[TimestepIndex(i2)]
+                        interpolated_emission = e1 + (e2 - e1) * (uu - y1) / (y2 - y1)
+                        v.methane_on_respiratoryadmissions[t, c] += interpolated_emission * p.respiratory_admissions_per_mton_ch4[c, tt]
+                    end
+                end
+            end
+        end
+    end
+end
+

src/main_model.jl

@@ -50,6 +50,11 @@ include("components/EquityWeighting.jl")
 include("components/PermafrostSiBCASA.jl")
 include("components/PermafrostJULES.jl")
 include("components/PermafrostTotal.jl")
+include("components/MarketDamageAQ_CropLoss.jl")
+include("components/MarketDamageAQ_LostWorkHours.jl")
+include("components/MarketDamageAQ_AsthmaERVisits.jl")
+include("components/MarketDamageAQ_RespiratoryAdmissions.jl")
+
 
 include("models/main_model_def.jl")
 

src/models/main_model_def.jl

@@ -84,14 +84,24 @@ function buildpage(m::Model, scenario::String, use_permafrost::Bool=true, use_se
     nonmarketdamages = addnonmarketdamages(m)
     discontinuity = adddiscontinuity(m)
 
-    # Total costs component
-    add_comp!(m, TotalCosts)
+    # Add MarketDamageAQ_CropLoss (Country version)
+    marketdamageaq_croplosscomp = addMarketDamageAQ_CropLoss(m)
+
+    # Add MarketDamageAQ_LostWorkHours (Country version)
+    marketdamageaq_lostworkhours = addMarketDamageAQ_LostWorkHours(m)
+
+    # Add MarketDamageAQ_AsthmaERVisits (Country version)
+    marketdamageaq_asthmaervisitscomp = addMarketDamageAQ_AsthmaERVisits(m)
 
-    # Equity weighting and Total Costs
+    # Add MarketDamageAQ_RespiratoryAdmissions (Country version)
+    marketdamageaq_respiratoryadmissionscomp = addMarketDamageAQ_RespiratoryAdmissions(m)
+
+
+    add_comp!(m, TotalCosts)
     countrylevelnpv = addcountrylevelnpv(m)
     equityweighting = addequityweighting(m)
 
-    # connect parameters together
+    #continue
     connect_param!(m, :GlobalTemperature => :fant_anthroforcing, :TotalForcing => :fant_anthroforcing)
     regtemp[:rt_g_globaltemperature] = glotemp[:rt_g_globaltemperature]
 
@@ -221,6 +231,11 @@ function buildpage(m::Model, scenario::String, use_permafrost::Bool=true, use_se
     connect_param!(m, :MarketDamagesBurke => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)
     connect_param!(m, :MarketDamagesBurke => :gdp, :GDP => :gdp)
     connect_param!(m, :MarketDamagesBurke => :pop_population, :Population => :pop_population)
+
+    connect_param!(m, :MarketDamageAQ_CropLoss => :global_ch4_emissions, :ch4emissions => :e_globalCH4emissions)
+    connect_param!(m, :MarketDamageAQ_LostWorkHours => :global_ch4_emissions, :ch4emissions => :e_globalCH4emissions)
+    connect_param!(m, :MarketDamageAQ_AsthmaERVisits => :global_ch4_emissions, :ch4emissions => :e_globalCH4emissions)
+    connect_param!(m, :MarketDamageAQ_RespiratoryAdmissions => :global_ch4_emissions, :ch4emissions => :e_globalCH4emissions)
 
     connect_param!(m, :NonMarketDamages => :rtl_realizedtemperature_change, :RegionTemperature => :rtl_realizedtemperature_change)
     connect_param!(m, :NonMarketDamages => :rtl_g_landtemperature, :RegionTemperature => :rtl_g_landtemperature)