Merge pull request #7 from brian-rose/newapi

version 0.3.0

.gitignore

@@ -1,5 +1,10 @@
 # Compiled python modules.
 *.pyc
+# Compiled fortran extensions
+*.so
+*.so.dSYM
+# f2py signature files
+*.pyf
 
 # Setuptools distribution folder.
 /dist/

README.rst

@@ -24,35 +24,49 @@ Installation
 About climlab
 --------------
 ``climlab`` is a flexible engine for process-oriented climate modeling.
-It is based on a very general concept of a model as a collection of individual, 
+It is based on a very general concept of a model as a collection of individual,
 interacting processes. ``climlab`` defines a base class called ``Process``, which
-can contain an arbitrarily complex tree of sub-processes (each also some 
-sub-class of ``Process``). Every climate process (radiative, dynamical, 
+can contain an arbitrarily complex tree of sub-processes (each also some
+sub-class of ``Process``). Every climate process (radiative, dynamical,
 physical, turbulent, convective, chemical, etc.) can be simulated as a stand-alone
-process model given appropriate input, or as a sub-process of a more complex model. 
+process model given appropriate input, or as a sub-process of a more complex model.
 New classes of model can easily be defined and run interactively by putting together an
 appropriate collection of sub-processes.
 
-Most of the actual computation uses vectorized ``numpy`` array functions. 
-It should run out-of-the-box on a standard scientific Python distribution.
-Future versions of ``climlab`` will provide hooks to compiled Fortran code for 
-more numerically intensive processes.
+Most of the actual computation for simpler model components use vectorized
+``numpy`` array functions. It should run out-of-the-box on a standard scientific
+Python distribution, such as ``Anaconda`` or ``Enthought Canopy``.
 
-Currently, ``climlab`` has out-of-the-box support and documented examples for 
+New in version 0.3, ``climlab`` now includes Python wrappers for more
+numerically intensive processes implemented in Fortran code (specifically the
+CAM3 radiation module). These require a Fortran compiler on your system,
+but otherwise have no other library dependencies.  ``climlab`` uses a compile-on-demand
+strategy. The compiler is invoked automatically as necessary when a new process
+in created by the user.
 
-- 1D grey-radiation and radiative-convective single column models
+Currently, ``climlab`` has out-of-the-box support and documented examples for
+
+- 1D radiative and radiative-convective single column models, with various radiation schemes:
+    - Grey Gas
+    - Simplified band-averaged models (4 bands each in longwave and shortwave)
+    - One GCM-level radiation module (CAM3)
 - 1D diffusive energy balance models
 - Seasonal and steady-state models
-- orbital / insolation calculations.
+- Arbitrary combinations of the above, for example:
+    - 2D latitude-pressure models with radiation, horizontal diffusion, and fixed relative humidity
+- orbital / insolation calculations
+- boundary layer sensible and latent heat fluxes
+
 
 Example usage
 ------------------
-The directory ``climlab/courseware/`` contains a collection of IPython notebooks (*.ipynb)
-used for teaching some basics of climate science, 
+The directory ``climlab/courseware/`` contains a collection of IPython / Jupyter
+notebooks (*.ipynb) used for teaching some basics of climate science,
 and documenting use of the ``climlab`` package.
 These are self-describing, and should all run out-of-the-box once the package is installed, e.g:
 
-``ipython notebook Insolation.ipynb``
+``jupyter notebook Insolation.ipynb``
+
 
 History
 ----------------------
@@ -61,11 +75,25 @@ in support of an undergraduate course at the University at Albany.
 See the original course webpage at
 http://www.atmos.albany.edu/facstaff/brose/classes/ENV480_Spring2014/
 
-The package and its API was completely redesigned around a truly object-oriented 
+The package and its API was completely redesigned around a truly object-oriented
 modeling framework in January 2015.
 
-It will be used extensively for a graduate-level climate modeling course in Spring 2015:
+It was used extensively for a graduate-level climate modeling course in Spring 2015:
 http://www.atmos.albany.edu/facstaff/brose/classes/ATM623_Spring2015/
+Many more examples are found in the online lecture notes for that course:
+http://nbviewer.jupyter.org/github/brian-rose/ClimateModeling_courseware/blob/master/index.ipynb
+
+Version 0.3 was released in February 2016. It includes many internal changes and
+some backwards-incompatible changes (hopefully simplifications) to the public API.
+It also includes the CAM3 radiation module.
+
+
+Contact and Bug Reports
+----------------------
+Users are strongly encouraged to submit bug reports and feature requests on
+github at
+https://github.com/brian-rose/climlab
+
 
 License
 ---------------

climlab/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '0.2.13'
+__version__ = '0.3.0'
 
 #  This list defines all the modules that will be loaded if a user invokes
 #   from climLab import *
@@ -10,7 +10,7 @@
 #           "column", "convadj"]
 
 #from climlab import radiation
-# this should ensure that we can still import constants.py as climlab.constants           
+# this should ensure that we can still import constants.py as climlab.constants
 from climlab.utils import constants
 from climlab.utils import thermo, legendre
 # some more useful shorcuts
@@ -20,6 +20,7 @@
 from climlab.domain import domain
 from climlab.domain.field import Field, global_mean
 from climlab.domain.axis import Axis
+from climlab.domain.initial import column_state, surface_state
 from climlab.process.process import Process, process_like, get_axes
 from climlab.process.time_dependent_process import TimeDependentProcess
 from climlab.process.implicit import ImplicitProcess

climlab/convection/convadj.py

@@ -7,7 +7,7 @@
 class ConvectiveAdjustment(TimeDependentProcess):
     '''Convective adjustment process
     Instantly returns column to neutral lapse rate
-    
+
     Adjustment includes the surface IF 'Ts' is included in the state
     dictionary. Otherwise only the atmopsheric temperature is adjusted.'''
     def __init__(self, adj_lapse_rate=None, **kwargs):
@@ -20,9 +20,9 @@ def __init__(self, adj_lapse_rate=None, **kwargs):
         patm = self.lev
         c_atm = self.Tatm.domain.heat_capacity
         if 'Ts' in self.state:
-            c_sfc = self.Ts.domain.heat_capacity       
-            self.pnew = np.flipud(np.append(np.flipud(patm), const.ps))
-            self.cnew = np.flipud(np.append(np.flipud(c_atm), c_sfc))
+            c_sfc = self.Ts.domain.heat_capacity
+            self.pnew = np.append(patm, const.ps)
+            self.cnew = np.append(c_atm, c_sfc)
         else:
             self.pnew = patm
             self.cnew = c_atm
@@ -36,27 +36,42 @@ def adj_lapse_rate(self, lapserate):
         else:
             self._adj_lapse_rate = lapserate
         self.param['adj_lapse_rate'] = self._adj_lapse_rate
-    
-    def compute(self):
+
+    def _compute(self):
         #lapse_rate = self.param['adj_lapse_rate']
         if self.adj_lapse_rate is None:
-            self.adjusted_state = self.state
+            self.adjustment = self.state * 0.
         else:
             #  For now, let's assume that the vertical axis is the last axis
             unstable_Tatm = self.Tatm
             if 'Ts' in self.state:
                 unstable_Ts = np.atleast_1d(self.Ts)
-                Tcol = np.concatenate((unstable_Ts, unstable_Tatm),axis=-1)
+                #Tcol = np.concatenate((unstable_Ts, unstable_Tatm),axis=-1)
+                Tcol = np.concatenate((unstable_Tatm, unstable_Ts),axis=-1)
             else:
                 Tcol = unstable_Tatm
-            Tadj = convective_adjustment_direct(self.pnew, Tcol, self.cnew, lapserate=self.adj_lapse_rate)            
-            if 'Ts' in self.state:            
-                Ts = Field(Tadj[...,0], domain=self.Ts.domain)
-                Tatm = Field(Tadj[...,1:], domain=self.Tatm.domain)
+            #  convective adjustment routine expect reversered vertical axis
+            pflip = self.pnew[..., ::-1]
+            Tflip = Tcol[..., ::-1]
+            cflip = self.cnew[..., ::-1]
+            Tadj_flip = convective_adjustment_direct(pflip, Tflip, cflip, lapserate=self.adj_lapse_rate)
+            Tadj = Tadj_flip[..., ::-1]
+            if 'Ts' in self.state:
+                Ts = Field(Tadj[...,-1], domain=self.Ts.domain)
+                Tatm = Field(Tadj[...,:-1], domain=self.Tatm.domain)
                 self.adjustment['Ts'] = Ts - self.Ts
             else:
                 Tatm = Field(Tadj, domain=self.Tatm.domain)
             self.adjustment['Tatm'] = Tatm - self.Tatm
+        # # express the adjustment (already accounting for the finite time step)
+        # #  as a tendency per unit time, so that it can be applied along with explicit
+        # tendencies = {}
+        # for name, adj in self.adjustment.iteritems():
+        #     tendencies[name] = adj / self.param['timestep']
+        # return tendencies
+        #  go back to just returning the adjustment, independent of timestep
+        #  because the parent process might have set a different timestep!
+        return self.adjustment
 
 
 def convective_adjustment_direct(p, T, c, lapserate=6.5):
@@ -162,7 +177,7 @@ def Akamaev_adjustment(theta, q, beta, n_k, theta_k, s_k, t_k):
         theta_k[k-1] = thistheta
         # back to step 2
 
-               
+
 #
 #    for l in range(2, L+1):  # l = 2,3,4,...,L
 #        n = 1
@@ -225,7 +240,7 @@ def Akamaev_adjustment(theta, q, beta, n_k, theta_k, s_k, t_k):
         n = n_k[k-1]
         thistheta = theta_k[k-1]
         # back to step 8
-    
+
 #    for l in range(L, 0, -1):
 #        if n > 1:
 #        while n>1:

climlab/domain/__init__.py

@@ -1 +1 @@
-__all__ = ['axis', 'domain']
+__all__ = ['axis', 'domain', 'field', 'initial']

climlab/domain/axis.py

@@ -89,9 +89,11 @@ def __init__(self, axis_type='abstract', num_points=10, points=None, bounds=None
         self.num_points = num_points
         self.units = defaultUnits[axis_type]
         # pressure axis should decrease from surface to TOA
-        if axis_type is 'lev':
-            points = np.flipud(points)
-            bounds = np.flipud(bounds)
+        #  NO! Now define the lowest (near-to-surface) element as lev[-1]
+        #   and the nearest to space as lev[0]
+        #if axis_type is 'lev':
+        #    points = np.flipud(points)
+        #    bounds = np.flipud(bounds)
         self.points = points
         self.bounds = bounds
         self.delta = np.abs(np.diff(self.bounds))

climlab/domain/initial.py

@@ -0,0 +1,57 @@
+'''Convenience routines for setting up initial conditions.'''
+import numpy as np
+from climlab.domain import domain
+from climlab.domain.field import Field
+from climlab.utils.attr_dict import AttrDict
+from climlab.utils import legendre
+
+
+def column_state(num_lev=30,
+                 num_lat=1,
+                 lev=None,
+                 lat=None,
+                 water_depth=1.0):
+    '''Set up a state variable dictionary consisting of temperatures
+    for atmospheric column (`Tatm`) and surface mixed layer (`Ts`).
+    '''
+    if lat is not None:
+        num_lat = np.array(lat).size
+    if lev is not None:
+        num_lev = np.array(lev).size
+
+    if num_lat is 1:
+        sfc, atm = domain.single_column(water_depth=water_depth,
+                                        num_lev=num_lev,
+                                        lev=lev)
+    else:
+        sfc, atm = domain.zonal_mean_column(water_depth=water_depth,
+                                            num_lev=num_lev,
+                                            lev=lev,
+                                            num_lat=num_lat,
+                                            lat=lat)
+    num_lev = atm.lev.num_points
+    Ts = Field(288.*np.ones(sfc.shape), domain=sfc)
+    Tinitial = np.tile(np.linspace(200., 288.-10., num_lev), sfc.shape)
+    Tatm = Field(Tinitial, domain=atm)
+    state = AttrDict()
+    state['Ts'] = Ts
+    state['Tatm'] = Tatm
+    return state
+
+
+def surface_state(num_lat=90,
+                  water_depth=10.,
+                  T0=12.,
+                  T2=-40.):
+    '''Set up a state variable dictionary for a zonal-mean surface model
+    (e.g. basic EBM). There is a single state variable `Ts`, the temperature
+    of the surface mixed layer.
+    '''
+    sfc = domain.zonal_mean_surface(num_lat=num_lat,
+                                    water_depth=water_depth)
+    sinphi = np.sin(np.deg2rad(sfc.axes['lat'].points))
+    initial = T0 + T2 * legendre.P2(sinphi)
+    Ts = Field(initial, domain=sfc)
+    state = AttrDict()
+    state['Ts'] = Ts
+    return state