From 35c1381a3f45d2a25b8fb42a55bcc50a2adeb7cc Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 07:48:22 -0500 Subject: [PATCH 1/9] Initial carbon setter --- beam_carbon/beam.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index cdcf8a7..e0553b7 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -41,6 +41,8 @@ def __init__(self, emissions=None, time_step=1, intervals=10): self._A = None self._B = None + self._initial_carbon = np.array([808.9, 725., 35641.]) + self._linear_temperature = False @property @@ -48,7 +50,11 @@ def initial_carbon(self): """Values for initial carbon in atmosphere, upper and lower oceans in GtC. Default values are from 2005. """ - return np.array([808.9, 725., 35641.]) + return self._initial_carbon + + @initial_carbon.setter + def initial_carbon(self, value): + self._initial_carbon = value @property def transfer_matrix(self): From 1a01b7290c3dbc1afdddbb9a02ebf7cfdb8c4403 Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 07:49:02 -0500 Subject: [PATCH 2/9] Carbon mass setter --- beam_carbon/beam.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index e0553b7..f8ab644 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -42,6 +42,7 @@ def __init__(self, emissions=None, time_step=1, intervals=10): self._B = None self._initial_carbon = np.array([808.9, 725., 35641.]) + self._carbon_mass = None self._linear_temperature = False @@ -56,6 +57,16 @@ def initial_carbon(self): def initial_carbon(self, value): self._initial_carbon = value + @property + def carbon_mass(self): + if self._carbon_mass is None: + self._carbon_mass = self.initial_carbon.copy() + return self._carbon_mass + + @carbon_mass.setter + def carbon_mass(self, value): + self._carbon_mass = value + @property def transfer_matrix(self): """3 by 3 matrix of transfer coefficients for carbon cycle. From 6036bdc4810d68e943728dc85fcd87fff4d8af94 Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 07:50:32 -0500 Subject: [PATCH 3/9] kh, k1, k2 calibrate at t0 if unspecified --- beam_carbon/beam.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index f8ab644..fe04523 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -138,6 +138,8 @@ def delta(self): def k_h(self): """CO2 solubility. """ + if self._k_h is None: + self._k_h = self.get_kh(0) return self._k_h @k_h.setter @@ -148,6 +150,8 @@ def k_h(self, value): def k_1(self): """First dissacoiation constant. """ + if self._k_1 is None: + self._k_1 = self.get_k1(0) return self._k_1 @k_1.setter @@ -158,6 +162,8 @@ def k_1(self, value): def k_2(self): """Second dissacoiation constant. """ + if self._k_2 is None: + self._k_2 = self.get_k2(0) return self._k_2 @k_2.setter From 9d03896316db9db98ef7654a574be759ff30b103 Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 07:51:50 -0500 Subject: [PATCH 4/9] Alkalinity setter --- beam_carbon/beam.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index fe04523..ee83cd8 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -40,6 +40,7 @@ def __init__(self, emissions=None, time_step=1, intervals=10): self._k_h = 1.23e3 self._A = None self._B = None + self._Alk = 767. self._initial_carbon = np.array([808.9, 725., 35641.]) self._carbon_mass = None @@ -186,7 +187,11 @@ def OM(self): def Alk(self): """Alkalinity in GtC. """ - return 767. + return self._Alk + + @Alk.setter + def Alk(self, value): + self._Alk = value @property def A(self): From 8f5be3435cbf9d9f1c6f48ae03aabab175e14124 Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 08:01:20 -0500 Subject: [PATCH 5/9] Calibrate ks at t0 from Tupper(0) --- beam_carbon/beam.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index ee83cd8..f82baab 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -140,7 +140,7 @@ def k_h(self): """CO2 solubility. """ if self._k_h is None: - self._k_h = self.get_kh(0) + self._k_h = self.get_kh(self.temperature.initial_temp[1]) return self._k_h @k_h.setter @@ -152,7 +152,7 @@ def k_1(self): """First dissacoiation constant. """ if self._k_1 is None: - self._k_1 = self.get_k1(0) + self._k_1 = self.get_k1(self.temperature.initial_temp[1]) return self._k_1 @k_1.setter @@ -164,7 +164,7 @@ def k_2(self): """Second dissacoiation constant. """ if self._k_2 is None: - self._k_2 = self.get_k2(0) + self._k_2 = self.get_k2(self.temperature.initial_temp[1]) return self._k_2 @k_2.setter @@ -469,9 +469,16 @@ def write_beam(output, csv=None): if __name__ == '__main__': b = BEAMCarbon() b.time_step = 10. - b.intervals = 20 + b.intervals = 200 N = 100 - b.emissions = np.concatenate((10. * np.exp(-np.arange(N) / 40), np.zeros(100-N))) + b.emissions = np.array([ + # 7.10, 7.97, + 9.58, 12.25, 14.72, 16.07, 17.43, 19.16, 20.89, 23.22, 26.15, 29.09 + ]) + # df = pd.DataFrame.from_csv('webDICE-CSV.csv', header=-1, index_col=0) + # b.emissions = np.array(df.ix['emissions_total', :]) + # b.emissions = np.concatenate((10. * np.exp(-np.arange(N) / 40), np.zeros(100-N))) b.temperature_dependent = False b.linear_temperature = False - r = b.run() \ No newline at end of file + r = b.run() + print(r) \ No newline at end of file From e75b0226c81554f5502e4479acc88f7f3dbb205d Mon Sep 17 00:00:00 2001 From: njmattes Date: Tue, 10 May 2016 08:06:25 -0500 Subject: [PATCH 6/9] Edits for clarity --- beam_carbon/beam.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index f82baab..f516e7e 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -353,6 +353,11 @@ def get_k2(self, temp_ocean): def run(self): N = self.n * self.intervals + self.carbon_mass = self.initial_carbon.copy() + total_carbon = 0 + emissions = np.zeros(3) + temp_atmosphere, temp_ocean = self.temperature.initial_temp + output = np.tile(np.concatenate(( self.initial_carbon, self.temperature.initial_temp, @@ -367,10 +372,6 @@ def run(self): 'cumulative', 'A', 'B'], columns=np.arange(self.n + 1) * self.time_step, ) - carbon_mass = self.initial_carbon.copy() - total_carbon = 0 - emissions = np.zeros(3) - temp_atmosphere, temp_ocean = self.temperature.initial_temp for i in xrange(N): @@ -379,10 +380,17 @@ def run(self): if i % self.intervals == 0 and self.temperature_dependent: # First interval in time step self.temp_calibrate(temp_ocean) - h = self.get_H(carbon_mass[1], re_solve=False) + h = self.get_H(self.carbon_mass[1], re_solve=False) self.B = self.get_B(h) - if i % self.intervals == 0: # First interval in time step + emissions[0] = self.emissions[_i] * self.time_step / self.intervals + total_carbon += emissions[0] + + self.carbon_mass += ( + (self.transfer_matrix * self.carbon_mass + / self.intervals).sum(axis=1) + emissions) + + if (i + 1) % self.intervals == 0: emissions[0] = self.emissions[_i] * self.time_step total_carbon += emissions[0] @@ -390,25 +398,21 @@ def run(self): ta = temp_atmosphere temp_atmosphere = self.temperature.temp_atmosphere( index=_i, temp_atmosphere=ta, - temp_ocean=temp_ocean, mass_atmosphere=carbon_mass[0], + temp_ocean=temp_ocean, mass_atmosphere=self.carbon_mass[0], carbon=total_carbon, initial_carbon=self.initial_carbon, phi11=self.transfer_matrix[0][0], phi21=self.transfer_matrix[1][0]) temp_ocean = self.temperature.temp_ocean( ta, temp_ocean) - carbon_mass += (((self.transfer_matrix * carbon_mass) / - self.intervals).sum(axis=1) + - emissions / self.intervals) - - if (i + 1) % self.intervals == 0: output.iloc[:, _i + 1] = ( - np.concatenate((carbon_mass.copy(), + np.concatenate((self.carbon_mass.copy(), np.array([temp_atmosphere, temp_ocean]), np.array([self.transfer_matrix[0][1], self.transfer_matrix[1][1], total_carbon, self.A, self.B])))) + return output From 2295efd2883279163db38d90253c8bce600d69c7 Mon Sep 17 00:00:00 2001 From: Nathan Matteson Date: Tue, 10 May 2016 11:00:05 -0500 Subject: [PATCH 7/9] Cleanup --- beam_carbon/beam.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index f516e7e..1af8e98 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -8,7 +8,7 @@ from beam_carbon.temperature import DICETemperature, LinearTemperature -__version__ = '0.1' +__version__ = '0.3' class BEAMCarbon(object): @@ -18,11 +18,13 @@ def __init__(self, emissions=None, time_step=1, intervals=10): """BEAMCarbon init Args: - :param emissions: Array of annual emissions in GtC, beginning in 2005 + :param emissions: Array of annual emissions in GtC, beginning + in 2005 :type emissions: list :param time_step: Time between emissions values in years :type time_step: float - :param intervals: Nuber of times to calculate BEAM carbon per timestep + :param intervals: Number of times to calculate BEAM carbon + per timestep :type intervals: int """ self._temperature_dependent = True @@ -149,7 +151,7 @@ def k_h(self, value): @property def k_1(self): - """First dissacoiation constant. + """First dissociation constant. """ if self._k_1 is None: self._k_1 = self.get_k1(self.temperature.initial_temp[1]) @@ -161,7 +163,7 @@ def k_1(self, value): @property def k_2(self): - """Second dissacoiation constant. + """Second dissociation constant. """ if self._k_2 is None: self._k_2 = self.get_k2(self.temperature.initial_temp[1]) @@ -228,7 +230,7 @@ def salinity(self): @property def temperature_dependent(self): - """Switch for calculating temperature-dependent paramters k_1, + """Switch for calculating temperature-dependent parameters k_1, k_2, and k_h. """ return self._temperature_dependent @@ -240,9 +242,11 @@ def linear_temperature(self): @linear_temperature.setter def linear_temperature(self, value): if value: - self.temperature = LinearTemperature(self.time_step, self.intervals, self.n) + self.temperature = LinearTemperature(self.time_step, + self.intervals, self.n) else: - self.temperature = DICETemperature(self.time_step, self.intervals, self.n) + self.temperature = DICETemperature(self.time_step, + self.intervals, self.n) self._linear_temperature = value @temperature_dependent.setter @@ -364,7 +368,7 @@ def run(self): np.array([ self.transfer_matrix[0][1], self.transfer_matrix[1][1]]), np.zeros(3), - )).reshape((10, 1)).copy(), self.n + 1) + )).reshape((10, 1)).copy(), (self.n + 1, )) output = pd.DataFrame( output, index=['mass_atmosphere', 'mass_upper', 'mass_lower', @@ -377,7 +381,7 @@ def run(self): _i = int(floor(i / self.intervals)) # time_step - if i % self.intervals == 0 and self.temperature_dependent: # First interval in time step + if i % self.intervals == 0 and self.temperature_dependent: self.temp_calibrate(temp_ocean) h = self.get_H(self.carbon_mass[1], re_solve=False) @@ -387,8 +391,9 @@ def run(self): total_carbon += emissions[0] self.carbon_mass += ( - (self.transfer_matrix * self.carbon_mass - / self.intervals).sum(axis=1) + emissions) + (self.transfer_matrix * + np.divide(self.carbon_mass, self.intervals)).sum(axis=1) + + emissions) if (i + 1) % self.intervals == 0: From 894b1d4c28f5a6cde6a69a8c97671ab8b29324d8 Mon Sep 17 00:00:00 2001 From: Nathan Matteson Date: Tue, 10 May 2016 11:06:53 -0500 Subject: [PATCH 8/9] Use numpy roots --- beam_carbon/beam.py | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index 1af8e98..54ff880 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -283,7 +283,7 @@ def get_A(self): """ return self.k_h * self.AM / (self.OM / (self.delta + 1)) - def get_H(self, mass_upper, re_solve=False): + def get_H(self, mass_upper): """Solve for H+, the concentration of hydrogen ions (the (pH) of seawater). @@ -292,21 +292,10 @@ def get_H(self, mass_upper, re_solve=False): :return: H :rtype: float """ - if re_solve: - h = symbols('h') - a = mass_upper / self.Alk - f = Eq( - (h**2 + self.k_1 * h + self.k_1 * self.k_2) / self.k_1, - a * (h + 2 * self.k_2) - ) - return max(solve(f, h)) - return ( - -self.k_1 * (self.Alk - mass_upper) / (2 * self.Alk) + - sqrt(self.k_1 * (self.Alk ** 2 * self.k_1 - - 4 * self.Alk ** 2 * self.k_2 - - 2 * self.Alk * self.k_1 * mass_upper + - 8 * self.Alk * self.k_2 * mass_upper + - self.k_1 * mass_upper ** 2)) / (2 * self.Alk)) + p0 = 1 + p1 = (self.k_1 - mass_upper * self.k_1 / self.Alk) + p2 = (1 - 2 * mass_upper / self.Alk) * self.k_1 * self.k_2 + return max(np.roots([p0, p1, p2])) def get_kh(self, temp_ocean): """Calculate temperature dependent k_h From d655ab10fa19dc324c72c417139db257b46c3581 Mon Sep 17 00:00:00 2001 From: Nathan Matteson Date: Tue, 10 May 2016 11:09:37 -0500 Subject: [PATCH 9/9] Separate pk1 and pk2 methods. Cleanup. --- beam_carbon/beam.py | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/beam_carbon/beam.py b/beam_carbon/beam.py index 54ff880..926f862 100644 --- a/beam_carbon/beam.py +++ b/beam_carbon/beam.py @@ -1,10 +1,9 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- from __future__ import division -from math import floor, sqrt +from math import floor import numpy as np import pandas as pd -from sympy import symbols, solve, Eq from beam_carbon.temperature import DICETemperature, LinearTemperature @@ -314,6 +313,18 @@ def get_kh(self, temp_ocean): self.A = kh * self.AM / (self.OM / (self.delta + 1.)) return kh + def get_pk1(self, t): + return ( + -13.721 + 0.031334 * t + 3235.76 / t + 1.3e-5 * self.salinity * t - + 0.1031 * self.salinity ** 0.5) + + def get_pk2(self, t): + return ( + 5371.96 + 1.671221 * t + 0.22913 * self.salinity + + 18.3802 * np.log10(self.salinity)) - (128375.28 / t + + 2194.30 * np.log10(t) + 8.0944e-4 * self.salinity * t + + 5617.11 * np.log10(self.salinity) / t) + 2.136 * self.salinity / t + def get_k1(self, temp_ocean): """Calculate temperature dependent k_1 @@ -322,11 +333,7 @@ def get_k1(self, temp_ocean): :return: k_1 :rtype: float """ - t = 283.15 + temp_ocean - pk1 = ( - -13.721 + 0.031334 * t + 3235.76 / t + 1.3e-5 * self.salinity * t - - 0.1031 * self.salinity ** 0.5) - return 10 ** -pk1 + return 10 ** -self.get_pk1(283.15 + temp_ocean) def get_k2(self, temp_ocean): """Calculate temperature dependent k_2 @@ -336,13 +343,7 @@ def get_k2(self, temp_ocean): :return: k_2 :rtype: float """ - t = 283.15 + temp_ocean - pk2 = ( - 5371.96 + 1.671221 * t + 0.22913 * self.salinity + - 18.3802 * np.log10(self.salinity)) - (128375.28 / t + - 2194.30 * np.log10(t) + 8.0944e-4 * self.salinity * t + - 5617.11 * np.log10(self.salinity) / t) + 2.136 * self.salinity / t - return 10 ** -pk2 + return 10 ** -self.get_pk2(283.15 + temp_ocean) def run(self): N = self.n * self.intervals @@ -373,7 +374,7 @@ def run(self): if i % self.intervals == 0 and self.temperature_dependent: self.temp_calibrate(temp_ocean) - h = self.get_H(self.carbon_mass[1], re_solve=False) + h = self.get_H(self.carbon_mass[1]) self.B = self.get_B(h) emissions[0] = self.emissions[_i] * self.time_step / self.intervals