docs: def zonal_mean_column(num_lat=90, num_lev=30, water_depth=10., lat=None,
                        lev=None, **kwargs):
      """Creates two Domains with one water cell, a latitude axis and
      a level/height axis.
  
          * SlabOcean:    one water cell and a latitude axis above
            (similar to :func:`zonal_mean_surface`)
          * Atmosphere: a latitude axis and a level/height axis (two dimensional)
  
  
      **Function-call argument** \n
  
      :param int num_lat:         number of latitude points on the axis
                                  [default: 90]
      :param int num_lev:         number of pressure levels
                                  (evenly spaced from surface to TOA) [default: 30]
      :param float water_depth:   depth of the water cell (slab ocean) [default: 10.]
      :param lat:                 specification for latitude axis (optional)
      :type lat:                  :class:`~climlab.domain.axis.Axis` or latitude array
      :param lev:                 specification for height axis (optional)
      :type lev:                  :class:`~climlab.domain.axis.Axis` or pressure array
      :raises: :exc:`ValueError`  if `lat` is given but neither Axis nor latitude array.
      :raises: :exc:`ValueError`  if `lev` is given but neither Axis nor pressure array.
      :returns:                   a list of 2 Domain objects (slab ocean, atmosphere)
      :rtype:                     :py:class:`list` of :class:`SlabOcean`, :class:`Atmosphere`
  
      :Example:
  
          ::
  
              >>> from climlab import domain
              >>> sfc, atm = domain.zonal_mean_column(num_lat=36,num_lev=10)
  
              >>> print sfc
              climlab Domain object with domain_type=ocean and shape=(36, 1)
  
              >>> print atm
              climlab Domain object with domain_type=atm and shape=(36, 10)
  
  
      """
      if lat is None:
          latax = Axis(axis_type='lat', num_points=num_lat)
      elif isinstance(lat, Axis):
          latax = lat
      else:
          try:
              latax = Axis(axis_type='lat', points=lat)
          except:
              raise ValueError('lat must be Axis object or latitude array')
      if lev is None:
          levax = Axis(axis_type='lev', num_points=num_lev)
      elif isinstance(lev, Axis):
          levax = lev
      else:
          try:
              levax = Axis(axis_type='lev', points=lev)
          except:
              raise ValueError('lev must be Axis object or pressure array')
  
      depthax = Axis(axis_type='depth', bounds=[water_depth, 0.])
      #axes = {'depth': depthax, 'lat': latax, 'lev': levax}
      slab = SlabOcean(axes={'lat':latax, 'depth':depthax}, **kwargs)
      atm = Atmosphere(axes={'lat':latax, 'lev':levax}, **kwargs)
      return slab, atm