Merge pull request #104 from brian-rose/land

Improve and standardize surface flux processes for land modeling

climlab/__init__.py

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

climlab/convection/emanuel_convection.py

@@ -98,7 +98,8 @@ class EmanuelConvection(TimeDependentProcess):
     Diagnostics computed:
 
         - CBMF (cloud base mass flux in kg/m2/s) -- this is actually stored internally and used as input for subsequent timesteps
-        - PRECIP (convective precipitation rate in mm/day)
+        - precipitation (convective precipitation rate in kg/m2/s or mm/s)
+        - relative_humidity (dimensionless)
 
         :Example:
 
@@ -108,6 +109,7 @@ class EmanuelConvection(TimeDependentProcess):
             This example also demonstrates *asynchronous coupling*:
             the radiation uses a longer timestep than the other model components::
 
+                import numpy as np
                 import climlab
                 from climlab import constants as const
                 # Temperatures in a single column
@@ -169,11 +171,17 @@ def __init__(self,
         super(EmanuelConvection, self).__init__(**kwargs)
         self.time_type = 'explicit'
         #  Define inputs and diagnostics
-        #surface_shape = self.Tatm[...,0].shape
-        #  For some strange reason self.Tatm is breaking tests under Python 3.5 in some configurations
         surface_shape = self.state['Tatm'][...,0].shape
-        self.add_diagnostic('CBMF', np.zeros(surface_shape))  # cloud base mass flux
-        self.add_diagnostic('PRECIP', np.zeros(surface_shape)) # Precip rate (mm/day)
+        #  Hack to handle single column and multicolumn
+        if surface_shape is ():
+            init = np.atleast_1d(np.zeros(surface_shape))
+            self.multidim=False
+        else:
+            init = np.zeros(surface_shape)[...,np.newaxis]
+            self.multidim=True
+        self.add_diagnostic('CBMF', init*0.)  # cloud base mass flux
+        self.add_diagnostic('precipitation', init*0.) # Precip rate (kg/m2/s)
+        self.add_diagnostic('relative_humidity', 0*self.Tatm)
         self.add_input('MINORIG', MINORIG)
         self.add_input('ELCRIT', ELCRIT)
         self.add_input('TLCRIT', TLCRIT)
@@ -236,6 +244,12 @@ def _compute(self):
         if 'V' in self.state:
             tendencies['V'] = _convect_to_climlab(FV) * np.ones_like(self.state['V'])
         self.CBMF = CBMFnew
-        self.PRECIP = PRECIP
+        #  Need to convert from mm/day to mm/s or kg/m2/s
+        #  Hack to handle single column and multicolumn
+        if self.multidim:
+            self.precipitation[:,0] = _convect_to_climlab(PRECIP)/const.seconds_per_day
+        else:
+            self.precipitation[:] = _convect_to_climlab(PRECIP)/const.seconds_per_day
         self.IFLAG = IFLAG
+        self.relative_humidity[:] = self.q / qsat(self.Tatm,self.lev)
         return tendencies

climlab/surface/__init__.py

@@ -1,4 +1,11 @@
-'''Modules for surface processes in climlab.'''
+'''Modules for surface processes in climlab.
+
+``climlab.surface.turbulent`` contains ``SensibleHeatFlux`` and ``LatentHeatFlux``
+processes that implement bulk aerodynamic fluxes for surface energy and water exchange.
+
+``climlab.surface.albedo`` contains processes for surface albedo
+and interactive ice and snow lines for energy balance models.
+'''
 from __future__ import absolute_import
 from .turbulent import SensibleHeatFlux, LatentHeatFlux
 from .albedo import ConstantAlbedo, P2Albedo, Iceline, StepFunctionAlbedo

climlab/surface/turbulent.py

@@ -1,3 +1,59 @@
+'''Processes for surface turbulent heat and moisture fluxes
+
+:class:`~climlab.surface.SensibleHeatFlux` and
+:class:`~climlab.surface.LatentHeatFlux` implement standard bulk formulae
+for the turbulent heat fluxes, assuming that the heating or moistening
+occurs in the lowest atmospheric model level.
+
+        :Example:
+
+            Here is an example of setting up a single-column
+            Radiative-Convective model with interactive water vapor
+            and surface latent and sensible heat fluxes.
+
+            This example also demonstrates *asynchronous coupling*:
+            the radiation uses a longer timestep than the other model components::
+
+                import numpy as np
+                import climlab
+                from climlab import constants as const
+                # Temperatures in a single column
+                full_state = climlab.column_state(num_lev=30, water_depth=2.5)
+                temperature_state = {'Tatm':full_state.Tatm,'Ts':full_state.Ts}
+                #  Initialize a nearly dry column (small background stratospheric humidity)
+                q = np.ones_like(full_state.Tatm) * 5.E-6
+                #  Add specific_humidity to the state dictionary
+                full_state['q'] = q
+                #  ASYNCHRONOUS COUPLING -- the radiation uses a much longer timestep
+                #  The top-level model
+                model = climlab.TimeDependentProcess(state=full_state,
+                                              timestep=const.seconds_per_hour)
+                #  Radiation coupled to water vapor
+                rad = climlab.radiation.RRTMG(state=temperature_state,
+                                              specific_humidity=full_state.q,
+                                              albedo=0.3,
+                                              timestep=const.seconds_per_day
+                                              )
+                #  Convection scheme -- water vapor is a state variable
+                conv = climlab.convection.EmanuelConvection(state=full_state,
+                                              timestep=const.seconds_per_hour)
+                #  Surface heat flux processes
+                shf = climlab.surface.SensibleHeatFlux(state=temperature_state, Cd=0.5E-3,
+                                              timestep=const.seconds_per_hour)
+                lhf = climlab.surface.LatentHeatFlux(state=full_state, Cd=0.5E-3,
+                                              timestep=const.seconds_per_hour)
+                #  Couple all the submodels together
+                model.add_subprocess('Radiation', rad)
+                model.add_subprocess('Convection', conv)
+                model.add_subprocess('SHF', shf)
+                model.add_subprocess('LHF', lhf)
+                print(model)
+
+                #  Run the model
+                model.integrate_years(1)
+                #  Check for energy balance
+                print(model.ASR - model.OLR)
+'''
 from __future__ import division
 import numpy as np
 from climlab.utils.thermo import qsat
@@ -6,10 +62,12 @@
 from climlab.domain.field import Field
 
 
-class SurfaceFlux(EnergyBudget):
-    def __init__(self, Cd=3E-3, **kwargs):
-        super(SurfaceFlux, self).__init__(**kwargs)
+class _SurfaceFlux(EnergyBudget):
+    '''Abstract parent class for SensibleHeatFlux and LatentHeatFlux'''
+    def __init__(self, Cd=3E-3, resistance=1., **kwargs):
+        super(_SurfaceFlux, self).__init__(**kwargs)
         self.Cd = Cd
+        self.add_input('resistance', resistance)
         self.heating_rate['Tatm'] = np.zeros_like(self.Tatm)
         #  fixed wind speed (for now)
         self.add_input('U', 5. * np.ones_like(self.Ts))
@@ -27,10 +85,36 @@ def _air_density(self, Ta):
         return self.ps * const.mb_to_Pa / const.Rd / Ta
 
 
-class SensibleHeatFlux(SurfaceFlux):
+class SensibleHeatFlux(_SurfaceFlux):
+    r'''Surface turbulent sensible heat flux implemented through a bulk aerodynamic formula.
+
+    The flux is computed from
+
+    .. math::
+
+        SH = r ~ c_p ~\rho ~ C_D ~ U \left( T_s - T_a \right)
+
+
+    where:
+
+    - :math:`c_p` and :math:`\rho` are the specific heat and density of air
+    - :math:`C_D` is a drag coefficient (stored as ``self.Cd``, default value is 3E-3)
+    - :math:`U` is the near-surface wind speed, stored as ``self.U``, default value is 5 m/s
+    - :math:`r` is an optional resistance parameter (stored as ``self.resistance``, default value = 1)
+
+    The surface temperature :math:`T_s` is taken directly from ``self.state['Ts']``,
+    while the near-surface air temperature :math:`T_a` is taken as the lowest model
+    level in ``self.state['Tatm']``
+
+    Diagnostic quantity ``self.SHF`` gives the sensible heat flux in W/m2.
+
+    Temperature tendencies associated with this flux are computed for
+    ``Ts`` and for the lowest model level in ``Tatm``. All other tendencies
+    (including air temperature tendencies at other levels) are set to zero.
+    '''
     def __init__(self, Cd=3E-3, **kwargs):
         super(SensibleHeatFlux, self).__init__(Cd=Cd, **kwargs)
-        self.add_diagnostic('SHF')
+        self.add_diagnostic('SHF', 0.*self.Ts)
 
     def _compute_flux(self):
         # this ensure same dimensions as Ts
@@ -40,14 +124,51 @@ def _compute_flux(self):
         DeltaT = Ts - Ta
         rho = self._air_density(Ta)
         #  flux from bulk formula
-        self._flux = const.cp * rho * self.Cd * self.U * DeltaT
+        self._flux = self.resistance * const.cp * rho * self.Cd * self.U * DeltaT
         self.SHF = self._flux
 
 
-class LatentHeatFlux(SurfaceFlux):
+
+class LatentHeatFlux(_SurfaceFlux):
+    r'''Surface turbulent latent heat flux implemented through a bulk aerodynamic formula.
+
+    The flux is computed from
+
+    .. math::
+
+        LH = r ~ L ~\rho ~ C_D ~ U \left( q_s - q_a \right)
+
+
+    where:
+
+    - :math:`L` and :math:`\rho` are the latent heat of vaporization and density of air
+    - :math:`C_D` is a drag coefficient (stored as ``self.Cd``, default value is 3E-3)
+    - :math:`U` is the near-surface wind speed, stored as ``self.U``, default value is 5 m/s
+    - :math:`r` is an optional resistance parameter (stored as ``self.resistance``, default value = 1)
+
+    The surface specific humidity :math:`q_s` is computed as the saturation specific
+    humidity at the surface temperature ``self.state['Ts']`` and surface pressure
+    ``self.ps``, while the near-surface specific humidity :math:`q_a` is taken as the lowest model
+    level in the field ``self.q`` (which must be provided either as a state variable or as input).
+
+    Two diagnostics are computed:
+
+    - ``self.LHF`` gives the sensible heat flux in W/m2.
+    - ``self.evaporation`` gives the evaporation rate in kg/m2/s (or mm/s)
+
+    How the tendencies are computed depends on whether specific humidity ``q``
+    is a state variable (i.e. is present in ``self.state``):
+
+    - If ``q`` is in ``self.state`` then the evaporation determines the specific humidity tendency ``self.tendencies['q']``. The water vapor is added to the lowest model level only. Evaporation cools the surface through the surface tendency ``self.tendencies['Ts']``. Air temperature tendencies are zero everywhere.
+    - If ``q`` is not in ``self.state`` then we compute an equivalent air temperature tendency for the lowest model layer instead of a specific humidity tendency (i.e. the latent heat flux is applied in the same way as a sensible heat flux).
+
+    This process does not apply a tendency to the surface water amount.
+    In the absence of other water processes this implies an infinite water source at the surface (slab ocean).
+    '''
     def __init__(self, Cd=3E-3, **kwargs):
         super(LatentHeatFlux, self).__init__(Cd=Cd, **kwargs)
-        self.add_diagnostic('LHF')
+        self.add_diagnostic('LHF', 0.*self.Ts)
+        self.add_diagnostic('evaporation', 0.*self.Ts)  # in kg/m2/s or mm/s
 
     def _compute_flux(self):
         #  specific humidity at lowest model level
@@ -58,8 +179,10 @@ def _compute_flux(self):
         Deltaq = Field(qs - q, domain=self.Ts.domain)
         rho = self._air_density(Ta)
         #  flux from bulk formula
-        self._flux = const.Lhvap * rho * self.Cd * self.U * Deltaq
-        self.LHF = self._flux
+        self._flux = self.resistance * const.Lhvap * rho * self.Cd * self.U * Deltaq
+        self.LHF[:] = self._flux
+        # evporation rate, convert from W/m2 to kg/m2/s (or mm/s)
+        self.evaporation[:] = self.LHF/const.Lhvap
 
     def _compute(self):
         '''Overides the _compute method of EnergyBudget'''
@@ -71,6 +194,6 @@ def _compute(self):
             Pa_per_hPa = 100.
             air_mass_per_area = self.Tatm.domain.lev.delta[...,-1] * Pa_per_hPa / const.g
             specific_humidity_tendency = 0.*self.q
-            specific_humidity_tendency[...,-1] = self.LHF/const.Lhvap / air_mass_per_area
+            specific_humidity_tendency[...,-1,np.newaxis] = self.LHF/const.Lhvap / air_mass_per_area
             tendencies['q'] = specific_humidity_tendency
         return tendencies

conda-recipe/bld.bat

@@ -6,5 +6,5 @@ echo compiler=mingw32 >> "%CFG%"
 echo [build_ext] >> "%CFG%"
 echo compiler=mingw32 >> "%CFG%"
 
-"%PYTHON%" setup.py install --single-version-externally-managed --record=record.txt
+"%PYTHON%" -m pip install . --no-deps -vv
 if errorlevel 1 exit 1

conda-recipe/meta.yaml

@@ -1,4 +1,4 @@
-{% set version = "0.7.2.dev3" %}
+{% set version = "0.7.2.dev4" %}
 
 package:
   name: climlab
@@ -9,7 +9,7 @@ source:
 
 build:
   skip: True  # [win32 or (win and py27)]
-  number: 0
+  number: 2
 
 requirements:
   build:

setup.py

@@ -1,7 +1,7 @@
 import os, sys
 import textwrap
 
-VERSION = '0.7.2.dev3'
+VERSION = '0.7.2.dev4'
 
 # BEFORE importing setuptools, remove MANIFEST. Otherwise it may not be
 # properly updated when the contents of directories change (true for distutils,