Can now pass albedo parameters when initializing EBMs. Seasonal EBM c…

…an have step function albedo if ai is given as input argument.

climlab/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '0.2.6'
+__version__ = '0.2.7'
 
 #  This list defines all the modules that will be loaded if a user invokes
 #   from climLab import *

climlab/model/ebm.py

@@ -24,8 +24,11 @@ def __init__(self,
                  A=210.,
                  B=2.,
                  D=0.555,  # in W / m^2 / degC, same as B
-                 Tf=-10.0,
                  water_depth=10.0,
+                 Tf=-10.,
+                 a0=0.3,
+                 a2=0.078,
+                 ai=0.62,
                  timestep=const.seconds_per_year/90.,
                  **kwargs):
         super(EBM, self).__init__(timestep=timestep, **kwargs)
@@ -40,6 +43,9 @@ def __init__(self,
         self.param['D'] = D
         self.param['Tf'] = Tf
         self.param['water_depth'] = water_depth
+        self.param['a0'] = a0
+        self.param['a2'] = a2
+        self.param['ai'] = ai
         # create sub-models
         self.add_subprocess('LW', AplusBT(state=self.state, **self.param))
         self.add_subprocess('insolation',
@@ -62,7 +68,8 @@ def _compute_heating_rates(self):
         ASR = (1-albedo) * insolation
         self.heating_rate['Ts'] = ASR
         self.diagnostics['ASR'] = ASR
-        self.diagnostics['net_radiation'] = ASR - self.subprocess['LW'].diagnostics['OLR']
+        self.diagnostics['net_radiation'] = (ASR - 
+                                    self.subprocess['LW'].diagnostics['OLR'])
 
     def global_mean_temperature(self):
         '''Convenience method to compute the global mean surface temperature.'''
@@ -73,19 +80,33 @@ class EBM_annual(EBM):
     def __init__(self, **kwargs):
         super(EBM_annual, self).__init__(**kwargs)
         sfc = self.domains['Ts']
-        self.add_subprocess('insolation', AnnualMeanInsolation(domains=sfc, **self.param))
+        self.add_subprocess('insolation',
+                            AnnualMeanInsolation(domains=sfc, **self.param))
 
 
 class EBM_seasonal(EBM):
-    def __init__(self, **kwargs):
-        super(EBM_seasonal, self).__init__(**kwargs)
+    def __init__(self, a0=0.33, a2=0.25, ai=None, **kwargs):
+        '''This EBM uses realistic daily insolation.
+        If ai is not given, the model will not have an albedo feedback.'''
+        super(EBM_seasonal, self).__init__(a0=a0, a2=a2, ai=ai, **kwargs)
         sfc = self.domains['Ts']
-        self.add_subprocess('insolation', DailyInsolation(domains=sfc, **self.param))
-        # By default this EBM does not have an albedo feedback
-        self.add_subprocess('albedo', albedo.P2Albedo(domains=sfc, **self.param))
+        self.add_subprocess('insolation',
+                            DailyInsolation(domains=sfc, **self.param))
+        self.param['a0'] = a0
+        self.param['a2'] = a2
+        if ai is None:
+            # No albedo feedback
+            # Remove unused parameters here for clarity
+            _ = self.param.pop('ai')
+            _ = self.param.pop('Tf')
+            self.add_subprocess('albedo',
+                            albedo.P2Albedo(domains=sfc, **self.param))
+        else:
+            self.param['ai'] = ai
+            self.add_subprocess('albedo',
+                    albedo.StepFunctionAlbedo(state=self.state, **self.param))
+
 
-
-
 # an EBM that computes degree-days has an additional state variable.
 #  Need to implement that
 #  could make a good working example to document creating a new model class

climlab/surface/albedo.py

@@ -8,28 +8,58 @@ class ConstantAlbedo(DiagnosticProcess):
     def __init__(self, albedo=0.33, **kwargs):
         '''Uniform prescribed albedo.'''
         super(ConstantAlbedo, self).__init__(**kwargs)
-        self.param['albedo'] = albedo
-        #lat = self.domains['default'].axes['lat'].points
-        #albedo = np.ones_like(lat)
-        # make sure that the diagnostic has the correct field dimensions.
+        self.albedo = albedo
+
+    @property
+    def albedo(self):
+        return self._albedo
+    @albedo.setter
+    def albedo(self, value):
+        self._albedo = value
+        self.param['albedo'] = value
+        self._compute_fixed()
+    def _compute_fixed(self):
+        '''Recompute any fixed quantities after a change in parameters'''
         dom = self.domains['default']
-        self.diagnostics['albedo'] = Field(albedo, domain=dom)
+        self.diagnostics['albedo'] = Field(self.albedo, domain=dom)
 
 
 class P2Albedo(DiagnosticProcess):
     def __init__(self, a0=0.33, a2=0.25, **kwargs):
         '''Second order Legendre polynomial formula for surface albedo.'''
         super(P2Albedo, self).__init__(**kwargs)
-        self.param['a0'] = a0
-        self.param['a2'] = a2
-        lat = self.domains['default'].axes['lat'].points
-        phi = np.deg2rad(lat)
-        albedo = a0 + a2 * P2(np.sin(phi))
+        self.a0 = a0
+        self.a2 = a2
+
+    @property
+    def a0(self):
+        return self._a0
+    @a0.setter
+    def a0(self, value):
+        self._a0 = value
+        self.param['a0'] = value
+        self._compute_fixed()
+    @property
+    def a2(self):
+        return self._a2
+    @a2.setter
+    def a2(self, value):
+        self._a2 = value
+        self.param['a2'] = value
+        self._compute_fixed()
+    def _compute_fixed(self):
+        '''Recompute any fixed quantities after a change in parameters'''
+        phi = np.deg2rad(self.lat)
+        try:
+            albedo = self.a0 + self.a2 * P2(np.sin(phi))
+        except:
+            albedo = np.zeros_like(phi)
         # make sure that the diagnostic has the correct field dimensions.
         dom = self.domains['default']
         self.diagnostics['albedo'] = Field(albedo, domain=dom)
 
 
+
 class Iceline(DiagnosticProcess):
     def __init__(self, Tf=-10., **kwargs):
         super(DiagnosticProcess, self).__init__(**kwargs)

courseware/Seasonal cycle and heat capacity.ipynb

@@ -2143,7 +2143,7 @@
        "IAiCIAiCIHQK/x/ANJOX2G1x1AAAAABJRU5ErkJggg==\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x110d0fb10>"
+       "<matplotlib.figure.Figure at 0x1031bf590>"
       ]
      },
      "metadata": {},
@@ -2491,7 +2491,7 @@
        "/w/01pa89OTHqAAAAABJRU5ErkJggg==\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x111c1d150>"
+       "<matplotlib.figure.Figure at 0x111269ed0>"
       ]
      },
      "metadata": {},
@@ -2896,7 +2896,7 @@
        "BkIIISDFQAghBKQYCCGEgBQDIYQQAP4f9Tu1o0zeP3cAAAAASUVORK5CYII=\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x111c74850>"
+       "<matplotlib.figure.Figure at 0x110260a10>"
       ]
      },
      "metadata": {},
@@ -3647,7 +3647,7 @@
        "AAAwMOMQAAAAwMC+ACcmkcGlouVRAAAAAElFTkSuQmCC\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x113434ad0>"
+       "<matplotlib.figure.Figure at 0x113357cd0>"
       ]
      },
      "metadata": {},
@@ -4405,7 +4405,7 @@
        "MUmSJEmSJE0dB8UkSZIkSZI0dRwUkyRJkiRJ0tT5T2TNK8weWnjAAAAAAElFTkSuQmCC\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x113f67b10>"
+       "<matplotlib.figure.Figure at 0x113fc4190>"
       ]
      },
      "metadata": {},
@@ -5033,7 +5033,7 @@
        "mZmZmZmZdeSJIzMzMzMzMzMz68gTR2ZmZmZmZmZm1tHPACiZDBjutGQGAAAAAElFTkSuQmCC\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x113f0d0d0>"
+       "<matplotlib.figure.Figure at 0x111f1cdd0>"
       ]
      },
      "metadata": {},
@@ -5443,7 +5443,7 @@
        "RU5ErkJggg==\n"
       ],
       "text/plain": [
-       "<matplotlib.figure.Figure at 0x112e83310>"
+       "<matplotlib.figure.Figure at 0x11332d110>"
       ]
      },
      "metadata": {},