Merge pull request #72 from brian-rose/moistadiabat

Convective adjustment to the moist adiabat

climlab/__init__.py

@@ -7,7 +7,7 @@
 has been documented for a comprehensive understanding and traceability.
 '''
 
-__version__ = '0.6.5'
+__version__ = '0.6.6.dev0'
 
 # this should ensure that we can still import constants.py as climlab.constants
 from .utils import constants

climlab/convection/akmaev_adjustment.py

@@ -23,10 +23,9 @@ def convective_adjustment_direct(p, T, c, lapserate=6.5):
     # largely follows notation and algorithm in Akmaev (1991) MWR
     alpha = const.Rd / const.g * lapserate / 1.E3 # same dimensions as lapserate
     L = p.size
-    #  Here const.ps = 1000 hPa is just used as a reference pressure
-    #  to compute potential temperature, so is valid even for different
-    #  surface pressures
-    Pi = (p[:]/const.ps)**alpha  # will need to modify to allow variable lapse rates
+    ###  now handles variable lapse rate
+    pextended = np.insert(p,0,const.ps)  # prepend const.ps = 1000 hPa as ref pressure to compute potential temperature
+    Pi = np.cumprod((p / pextended[:-1])**alpha)  # Akmaev's equation 14 recurrence formula
     beta = 1./Pi
     theta = T * beta
     q = Pi * c

climlab/convection/convadj.py

@@ -2,20 +2,44 @@
 from builtins import range
 import numpy as np
 from climlab import constants as const
+from climlab.utils.thermo import rho_moist, pseudoadiabat
 from climlab.process.time_dependent_process import TimeDependentProcess
 from climlab.domain.field import Field
 from .akmaev_adjustment import convective_adjustment_direct
 
 
 class ConvectiveAdjustment(TimeDependentProcess):
-    '''Convective adjustment process
-    
-    Instantly returns column to neutral lapse rate
+    '''Hard Convective Adjustment to a prescribed lapse rate.
 
-    Adjustment includes the surface IF 'Ts' is included in the state
-    dictionary. Otherwise only the atmopsheric temperature is adjusted.
+    This process computes the instantaneous adjustment to conservatively
+    remove any instabilities in each column.
 
-    Implements the conservative adjustment algorithm from Akmaev (1991) Monthly Weather Review.
+    Instability is defined as a temperature decrease with height that exceeds
+    the prescribed critical lapse rate. This critical rate is set by input argument
+    ``adj_lapse_rate``, which can be either a numerical or string value.
+
+    Numerical values for ``adj_lapse_rate`` are given in units of K / km. Both
+    array and scalar values are valid. For scalar values, the assumption is that
+    the critical lapse rate is the same at every level.
+
+    If an array is given, it is assumed to represent the in-situ critical lapse
+    rate (in K/km) at every grid point.
+
+    Alternatively, string arguments can be given as follows:
+
+        - ``'DALR'`` or ``'dry adiabat'``: critical lapse rate is set to g/cp = 9.8 K / km
+        - ``'MALR'`` or ``'moist adiabat'`` or ``'pseudoadiabat'``: critical lapse rate follows the in-situ moist pseudoadiabat at every level
+
+    Adjustment includes the surface if ``'Ts'`` is included in the ``state``
+    dictionary. This implicitly accounts for turbulent surface fluxes.
+    Otherwise only the atmospheric temperature is adjusted.
+
+    If ``adj_lapse_rate`` is an array, its size must match the number of vertical
+    levels of the adjustment. This is number of pressure levels if the surface is
+    not adjusted, or number of pressure levels + 1 if the surface is adjusted.
+
+    This process implements the conservative adjustment algorithm described in
+    Akmaev (1991) Monthly Weather Review.
     '''
     def __init__(self, adj_lapse_rate=None, **kwargs):
         super(ConvectiveAdjustment, self).__init__(**kwargs)
@@ -24,49 +48,65 @@ def __init__(self, adj_lapse_rate=None, **kwargs):
         self.param['adj_lapse_rate'] = adj_lapse_rate
         self.time_type = 'adjustment'
         self.adjustment = {}
+    @property
+    def pcol(self):
         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.append(patm, const.ps)
             #  surface pressure should correspond to model domain!
             ps = self.lev_bounds[-1]
-            self.pnew = np.append(patm, ps)
-            self.cnew = np.append(c_atm, c_sfc)
+            return np.append(patm, ps)
         else:
-            self.pnew = patm
-            self.cnew = c_atm
+            return patm
+    @property
+    def ccol(self):
+        c_atm = self.Tatm.domain.heat_capacity
+        if 'Ts' in self.state:
+            c_sfc = self.Ts.domain.heat_capacity
+            return np.append(c_atm, c_sfc)
+        else:
+            return c_atm
+    @property
+    def Tcol(self):
+        #  For now, let's assume that the vertical axis is the last axis
+        Tatm = self.Tatm
+        if 'Ts' in self.state:
+            Ts = np.atleast_1d(self.Ts)
+            return np.concatenate((Tatm, Ts),axis=-1)
+        else:
+            return Tatm
     @property
     def adj_lapse_rate(self):
-        return self._adj_lapse_rate
+        lapserate = self._adj_lapse_rate
+        if type(lapserate) is str:
+            if lapserate in ['DALR', 'dry adiabat']:
+                return const.g / const.cp * 1.E3
+            elif lapserate in ['MALR', 'moist adiabat', 'pseudoadiabat']:
+                # critical lapse rate at each level is set by pseudoadiabat
+                dTdp = pseudoadiabat(self.Tcol,self.pcol) / 100.  # K / Pa
+                #  Could include water vapor effect on density here ...
+                #  Replace Tcol with virtual temperature
+                rho = self.pcol*100./const.Rd/self.Tcol  # in kg/m**3
+                return dTdp * const.g * rho * 1000.  # K / km
+            else:
+                raise ValueError('adj_lapse_rate must be either numeric or any of \'DALR\', \'dry adiabat\', \'MALR\', \'moist adiabat\', \'pseudoadiabat\'.')
+        else:
+            return lapserate
     @adj_lapse_rate.setter
     def adj_lapse_rate(self, lapserate):
-        if lapserate is 'DALR':
-            self._adj_lapse_rate = const.g / const.cp * 1.E3
-        else:
-            self._adj_lapse_rate = lapserate
-        self.param['adj_lapse_rate'] = self._adj_lapse_rate
+        self._adj_lapse_rate = lapserate
+        self.param['adj_lapse_rate'] = lapserate
 
     def _compute(self):
-        #lapse_rate = self.param['adj_lapse_rate']
         if self.adj_lapse_rate is None:
-            #self.adjustment = self.state * 0.
             self.adjustment['Ts'] = self.Ts * 0.
             self.adjustment['Tatm'] = self.Tatm * 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_Tatm, unstable_Ts),axis=-1)
-            else:
-                Tcol = unstable_Tatm
             #  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)
+            pflip = self.pcol[..., ::-1]
+            Tflip = self.Tcol[..., ::-1]
+            cflip = self.ccol[..., ::-1]
+            lapseflip = np.atleast_1d(self.adj_lapse_rate)[..., ::-1]
+            Tadj_flip = convective_adjustment_direct(pflip, Tflip, cflip, lapserate=lapseflip)
             Tadj = Tadj_flip[..., ::-1]
             if 'Ts' in self.state:
                 Ts = Field(Tadj[...,-1], domain=self.Ts.domain)

climlab/radiation/rrtm/__init__.py

@@ -16,22 +16,18 @@
             alb = 0.25
             #  State variables (Air and surface temperature)
             state = climlab.column_state(num_lev=30)
-            #  Parent model process
-            rcm = climlab.TimeDependentProcess(state=state)
             #  Fixed relative humidity
-            h2o = climlab.radiation.ManabeWaterVapor(state=state)
+            h2o = climlab.radiation.ManabeWaterVapor(name='WaterVapor', state=state)
             #  Couple water vapor to radiation
-            rad = climlab.radiation.RRTMG(state=state, specific_humidity=h2o.q, albedo=alb)
+            rad = climlab.radiation.RRTMG(name='Radiation', state=state, specific_humidity=h2o.q, albedo=alb)
             #  Convective adjustment
-            conv = climlab.convection.ConvectiveAdjustment(state=state, adj_lapse_rate=6.5)
+            conv = climlab.convection.ConvectiveAdjustment(name='Convection', state=state, adj_lapse_rate=6.5)
             #  Couple everything together
-            rcm.add_subprocess('Radiation', rad)
-            rcm.add_subprocess('WaterVapor', h2o)
-            rcm.add_subprocess('Convection', conv)
+            rcm = climlab.couple([rad,h2o,conv], name='Radiative-Convective Model')
             #  Run the model
             rcm.integrate_years(1)
             #  Check for energy balance
-            print rcm.ASR - rcm.OLR
+            print(rcm.ASR - rcm.OLR)
 
 '''
 from __future__ import division, absolute_import

climlab/tests/test_cam3rad.py

@@ -20,7 +20,7 @@ def rcm():
     # CAM3 radiation with default parameters and interactive water vapor
     rad = climlab.radiation.CAM3(state=state, albedo=alb, specific_humidity=h2o.q, name='Radiation')
     # Couple the models
-    rcm = h2o + convadj + rad
+    rcm = climlab.couple([h2o,convadj,rad], name='RCM')
     return rcm
 
 @pytest.mark.fast

climlab/tests/test_rcm.py

@@ -0,0 +1,40 @@
+from __future__ import division
+import numpy as np
+import climlab
+import pytest
+
+@pytest.fixture()
+def rcm():
+    # initial state (temperatures)
+    state = climlab.column_state(num_lev=40, num_lat=1, water_depth=5.)
+    ## Create individual physical process models:
+    #  fixed relative humidity
+    h2o = climlab.radiation.ManabeWaterVapor(state=state, name='H2O')
+    #  Hard convective adjustment
+    convadj = climlab.convection.ConvectiveAdjustment(state=state, name='ConvectiveAdjustment',
+                                                      adj_lapse_rate=6.5)
+    # RRTMG radiation with default parameters and interactive water vapor
+    rad = climlab.radiation.RRTMG(state=state, albedo=0.2, specific_humidity=h2o.q, name='Radiation')
+    # Couple the models
+    rcm = climlab.couple([h2o,convadj,rad], name='RCM')
+    return rcm
+
+@pytest.mark.fast
+def test_convective_adjustment(rcm):
+    rcm.step_forward()
+    #  test non-scalar critical lapse rate
+    num_lev = rcm.lev.size
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = np.linspace(5., 8., num_lev+1)
+    rcm.step_forward()
+    #  Test two flags for dry adiabatic adjustment
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = 'DALR'
+    rcm.step_forward()
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = 'dry adiabat'
+    rcm.step_forward()
+    #  test pseudoadiabatic critical lapse rate
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = 'pseudoadiabat'
+    rcm.step_forward()
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = 'MALR'
+    rcm.step_forward()
+    rcm.subprocess['ConvectiveAdjustment'].adj_lapse_rate = 'moist adiabat'
+    rcm.step_forward()

climlab/tests/test_thermo.py

@@ -4,6 +4,8 @@
 from climlab.utils import thermo
 import pytest
 
+
+@pytest.mark.fast
 def test_thermo():
     '''Basic single value tests for the thermodynamic routines.'''
     assert np.isclose(thermo.potential_temperature(250., 500.), 304.783)
@@ -16,11 +18,12 @@ def test_thermo():
 
     assert np.isclose(thermo.qsat(300., 1000.), 0.02227839)
 
-    assert np.isclose(thermo.estimated_inversion_strength(300., 290.), 2.58025)
-    assert np.isclose(thermo.EIS(300., 290.), 2.58025)
+    assert np.isclose(thermo.estimated_inversion_strength(300., 290.), 5.3605345)
+    assert np.isclose(thermo.EIS(300., 290.), 5.3605345)
 
     assert np.isclose(thermo.blackbody_emission(300.), 459.3)
 
+@pytest.mark.fast
 def test_thermo_domain():
     '''Can we call qsat, etc on a multi-dim climlab state temperature object?'''
     state = climlab.column_state(num_lev = 30, num_lat=3)

climlab/utils/constants.py

@@ -11,7 +11,7 @@
 """
 
 from __future__ import division
-import numpy as np
+from numpy import pi
 
 a = 6.373E6      # Radius of Earth (m)
 Lhvap = 2.5E6    # Latent heat of vaporization (J / kg)
@@ -23,14 +23,13 @@
 Rv = 461.5       # gas constant for water vapor (J / kg / K)
 cpv = 1875.   # specific heat at constant pressure for water vapor (J / kg / K)
 eps = Rd / Rv
-Omega = 2 * np.math.pi / 24. / 3600.  # Earth's rotation rate, (s**(-1))
+Omega = 2 * pi / 24. / 3600.  # Earth's rotation rate, (s**(-1))
 g = 9.8          # gravitational acceleration (m / s**2)
 kBoltzmann = 1.3806488E-23  # the Boltzmann constant (J / K)
 c_light = 2.99792458E8   # speed of light (m/s)
 hPlanck = 6.62606957E-34  # Planck's constant (J s)
-# sigma = 5.67E-8  # Stefan-Boltzmann constant (W / m**2 / K**4)
-#  sigma derived from fundamental constants
-sigma = (2*np.pi**5 * kBoltzmann**4) / (15 * c_light**2 * hPlanck**3)
+#  Stefan-Boltzmann constant (W / m**2 / K**4) derived from fundamental constants
+sigma = (2*pi**5 * kBoltzmann**4) / (15 * c_light**2 * hPlanck**3)
 
 S0 = 1365.2       # solar constant (W / m**2)
 # value is consistent with Trenberth and Fasullo, Surveys of Geophysics 2012
@@ -66,7 +65,7 @@
 hours_per_month = hours_per_year / months_per_year
 days_per_month = days_per_year / months_per_year
 
-area_earth = 4 * np.math.pi * a**2
+area_earth = 4 * pi * a**2
 
 # present-day orbital parameters, in the same format generated by orbital.py
 orb_present = {'ecc': 0.017236, 'long_peri': 281.37, 'obliquity': 23.446}