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
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)