Merge pull request #1 from m-kreuzer/changes_moritz

Changes moritz

climlab/dynamics/__init__.py

@@ -0,0 +1 @@
+from budyko_transport import BudykoTransport

climlab/dynamics/budyko_transport.py

@@ -0,0 +1,29 @@
+from climlab.process.energy_budget import EnergyBudget
+from climlab.domain.field import global_mean
+
+
+class BudykoTransport(EnergyBudget):
+    ''' module which calculates the ebm heat transport as the difference 
+    between the local temperature and the global mean temperature.
+    
+    parameters: b   - budyko transport parameter
+    
+    implemented by Moritz Kreuzer
+    '''
+    def __init__(self, b=3.81, **kwargs):
+        super(BudykoTransport, self).__init__(**kwargs)
+        self.b = b        
+
+    @property
+    def b(self):
+        return self._b
+    @b.setter
+    def b(self, value):
+        self._b = value
+        self.param['b'] = value
+
+    def _compute_heating_rates(self):
+        '''Compute energy flux convergences to get heating rates in W / m**2.'''
+        for varname, value in self.state.iteritems():
+            self.heating_rate[varname] = - self.b * (value - global_mean(value))
+

climlab/model/ebm.py

@@ -32,13 +32,15 @@ def __init__(self,
                  a2=0.078,
                  ai=0.62,
                  timestep=const.seconds_per_year/90.,
+                 T_init_0 = 12.,
+                 T_init_P2 = -40.,
                  **kwargs):
         super(EBM, self).__init__(timestep=timestep, **kwargs)
         if not self.domains and not self.state:  # no state vars or domains yet
             sfc = domain.zonal_mean_surface(num_lat=num_lat,
                                             water_depth=water_depth)
             lat = sfc.axes['lat'].points
-            initial = 12. - 40. * legendre.P2(np.sin(np.deg2rad(lat)))
+            initial = T_init_0 + T_init_P2 * legendre.P2(np.sin(np.deg2rad(lat)))
             self.set_state('Ts', Field(initial, domain=sfc))
         self.param['S0'] = S0
         self.param['A'] = A

climlab/process/time_dependent_process.py

@@ -155,3 +155,18 @@ def integrate_days(self, days=1.0, verbose=True):
         '''Timestep the model forward a specified number of days.'''
         years = days / const.days_per_year
         self.integrate_years(years=years, verbose=verbose)
+
+    def integrate_converge(self, crit=1e-4, verbose=True):
+        '''integrate until solution is converging
+        param:  crit    - exit criteria for difference of iterated solutions
+        '''
+        for varname, value in self.state.iteritems():
+            value_old = copy.deepcopy(value)
+            self.integrate_years(1,verbose=False)
+            while np.max(np.abs(value_old-value)) > crit : 
+                value_old = copy.deepcopy(value)                
+                self.integrate_years(1,verbose=False)              
+        if verbose == True:   
+            print("Total elapsed time is %s years." 
+                  % str(self.time['days_elapsed']/const.days_per_year))
+

climlab/radiation/AplusBT.py

@@ -19,6 +19,8 @@
 print olr.state
 '''
 from climlab.process.energy_budget import EnergyBudget
+from climlab.utils import constants as const
+import numpy as np
 
 
 class AplusBT(EnergyBudget):
@@ -55,3 +57,38 @@ def _compute_heating_rates(self):
         self.emission()
         for varname, value in self.state.iteritems():
             self.heating_rate[varname] = -self.OLR
+
+
+class AplusBT_CO2(EnergyBudget):
+    '''longwave radiation module considering CO2 concentration
+    see Caldeira & Kasting [1992] for further reading.
+    
+    parameter:  CO2  - CO2 concentration in atmosphere (ppm)
+    
+    implemented by Moritz Kreuzer'''
+    def __init__(self, CO2=300, **kwargs):
+        super(AplusBT_CO2, self).__init__(**kwargs)
+        self.CO2 = CO2
+
+    @property
+    def CO2(self):
+        return self._CO2
+    @CO2.setter
+    def CO2(self, value):
+        self._CO2 = value
+        self.param['CO2'] = value
+
+    def emission(self):
+        l = np.log(self.CO2/300.)
+        A = -326.400 + 9.16100*l - 3.16400*l**2 + 0.546800*l**3
+        B =    1.953 - 0.04866*l + 0.01309*l**2 - 0.002577*l**3
+        for varname, value in self.state.iteritems():
+            flux = A + B * (value + const.tempCtoK)
+            self.OLR = flux
+            self.diagnostics['OLR'] = self.OLR
+
+    def _compute_heating_rates(self):
+        '''Compute energy flux convergences to get heating rates in W / m**2.'''
+        self.emission()
+        for varname, value in self.state.iteritems():
+            self.heating_rate[varname] = -self.OLR

climlab/radiation/Boltzmann.py

@@ -0,0 +1,45 @@
+from climlab import constants as const
+from climlab.process.energy_budget import EnergyBudget
+
+
+class Boltzmann(EnergyBudget):
+    '''a radiation subprocess which computes longwave radiation
+    with the Stefan-Boltzmann law for black/grey body radiation
+    
+    parameter:  eps - emissivity
+                tau - transmissivity
+    
+    implemented by Moritz Kreuzer    
+    '''
+    def __init__(self, eps= 0.65, tau=0.95, **kwargs):
+        super(Boltzmann, self).__init__(**kwargs)
+        self.eps = eps
+        self.tau = tau
+
+    @property
+    def eps(self):
+        return self._eps
+    @eps.setter
+    def eps(self, value):
+        self._eps = value
+        self.param['eps'] = value
+
+    @property
+    def tau(self):
+        return self._tau
+    @tau.setter
+    def tau(self, value):
+        self._tau = value
+        self.param['tau'] = value
+
+    def emission(self):
+        for varname, value in self.state.iteritems():
+            flux = self.eps * self.tau * const.sigma * (value + const.tempCtoK)**4.
+            self.OLR = flux
+            self.diagnostics['OLR'] = self.OLR
+
+    def _compute_heating_rates(self):
+        '''Compute energy flux convergences to get heating rates in W / m**2.'''
+        self.emission()
+        for varname, value in self.state.iteritems():
+            self.heating_rate[varname] = -self.OLR

climlab/radiation/__init__.py

@@ -1,4 +1,6 @@
 from AplusBT import AplusBT
+from AplusBT import AplusBT_CO2
+from Boltzmann import Boltzmann
 from insolation import FixedInsolation, P2Insolation, AnnualMeanInsolation, DailyInsolation
 from radiation import Radiation
 #from three_band import ThreeBandSW