From d37c6cf4d3b73ac95c84f536cac0631834f7dd98 Mon Sep 17 00:00:00 2001
From: Liam Megill <liam.megill@dlr.de>
Date: Thu, 24 Oct 2024 11:02:46 +0200
Subject: [PATCH] Include PMrel capability

---
 example/example.toml     | 3 ++-
 openairclim/calc_cont.py | 9 ++++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/example/example.toml b/example/example.toml
index 15d70ef..942922e 100644
--- a/example/example.toml
+++ b/example/example.toml
@@ -10,7 +10,7 @@ inv = ["CO2", "H2O", "NOx", "distance"]
 nox = "NO"
 # Output / response species
 # possible values: "CO2", "H2O"
-out = ["CO2", "H2O", "CH4", "cont"]
+out = ["CO2", "H2O", "cont"]
 
 # Emission inventories                                                    
 [inventories]
@@ -80,6 +80,7 @@ CH4.rf.method = "Etminan_2016"
 cont.response_grid = "cont"
 cont.G_comp = 0.04  # conventional: 0.04; hydrogen: 0.12
 cont.eff_fac = 1.0  # 
+cont.PMrel = 1.0  # relative PM emissions compared to kerosene (1.0)
 
 # Temperature options
 [temperature]
diff --git a/openairclim/calc_cont.py b/openairclim/calc_cont.py
index 37806b7..0c66320 100644
--- a/openairclim/calc_cont.py
+++ b/openairclim/calc_cont.py
@@ -233,10 +233,17 @@ def calc_cont_RF(config, cccov_tot_dict, inv_dict):
             between the simulation start and end years.
     """
     
+    # calculate RF factor due to PM reduction, from AirClim 2.1
+    PMrel = config["responses"]["cont"]["PMrel"]
+    if PMrel >= 0.033:
+        PM_factor = 0.92 * np.arctan(1.902 * PMrel ** 0.74)
+    else:
+        PM_factor = 0.92 * np.arctan(1.902 * 0.033 ** 0.74)
+    
     # calculate contrail RF
     cont_RF_at_inv = []  # RF at inventory years
     for year, cccov_tot in cccov_tot_dict.items():
-        cont_RF = np.sum(14.9 * cccov_tot)
+        cont_RF = 14.9 * np.sum(cccov_tot) * PM_factor
         cont_RF_at_inv.append(cont_RF) 
     
     # interpolate RF to all simulation years