def surface_state(num_lat=90,
                  num_lon=None,
                  water_depth=10.,
                  T0=12.,
                  T2=-40.):
    """Sets up a state variable dictionary for a surface model
    (e.g. :class:`~climlab.model.ebm.EBM`) with a uniform slab ocean depth.

    The domain is either 1D (latitude) or 2D (latitude, longitude)
    depending on whether the input argument num_lon is supplied.

    Returns a single state variable `Ts`, the temperature of the surface
    mixed layer (slab ocean).

    The temperature is initialized to a smooth equator-to-pole shape given by

    .. math::

        T(\phi) = T_0 + T_2 P_2(\sin\phi)

    where :math:`\phi` is latitude, and :math:`P_2` is the second Legendre
    polynomial :class:`~climlab.utils.legendre.P2`.

    **Function-call arguments** \n

    :param int num_lat:         number of latitude points [default: 90]
    :param int num_lat:         (optional) number of longitude points [default: None]
    :param float water_depth:   depth of the slab ocean in meters [default: 10.]
    :param float T0:            global-mean initial temperature in :math:`^{\circ} \\textrm{C}` [default: 12.]
    :param float T2:            2nd Legendre coefficient for equator-to-pole gradient in
                                initial temperature, in :math:`^{\circ} \\textrm{C}` [default: -40.]

    :returns:                   dictionary with temperature
                                :class:`~climlab.domain.field.Field`
                                for surface mixed layer ``Ts``
    :rtype:                     dict


    :Example:

        ::

            >>> from climlab.domain import initial
            >>> import numpy as np

            >>> T_dict = initial.surface_state(num_lat=36)

            >>> print np.squeeze(T_dict['Ts'])
            [-27.88584094 -26.97777479 -25.18923361 -22.57456133 -19.21320344
             -15.20729309 -10.67854785  -5.76457135  -0.61467228   4.61467228
               9.76457135  14.67854785  19.20729309  23.21320344  26.57456133
              29.18923361  30.97777479  31.88584094  31.88584094  30.97777479
              29.18923361  26.57456133  23.21320344  19.20729309  14.67854785
               9.76457135   4.61467228  -0.61467228  -5.76457135 -10.67854785
             -15.20729309 -19.21320344 -22.57456133 -25.18923361 -26.97777479
             -27.88584094]

    """
    if num_lon is None:
        sfc = domain.zonal_mean_surface(num_lat=num_lat,
                                        water_depth=water_depth)
    else:
        sfc = domain.surface_2D(num_lat=num_lat,
                                num_lon=num_lon,
                                water_depth=water_depth)
    if 'lon' in sfc.axes:
        lon, lat = np.meshgrid(sfc.axes['lon'].points, sfc.axes['lat'].points)
    else:
        lat = sfc.axes['lat'].points
    sinphi = np.sin(np.deg2rad(lat))
    initial = T0 + T2 * legendre.P2(sinphi)
    Ts = Field(initial, domain=sfc)
    #if num_lon is None:
    #    Ts = Field(initial, domain=sfc)
    #else:
    #    Ts = Field([[initial for k in range(num_lon)]], domain=sfc)
    state = AttrDict()
    state['Ts'] = Ts
    return state
Sets up a state variable dictionary for a surface model
    (e.g. :class:`~climlab.model.ebm.EBM`) with a uniform slab ocean depth.

    The domain is either 1D (latitude) or 2D (latitude, longitude)
    depending on whether the input argument num_lon is supplied.

    Returns a single state variable `Ts`, the temperature of the surface
    mixed layer (slab ocean).

    The temperature is initialized to a smooth equator-to-pole shape given by

    .. math::

        T(\phi) = T_0 + T_2 P_2(\sin\phi)

    where :math:`\phi` is latitude, and :math:`P_2` is the second Legendre
    polynomial :class:`~climlab.utils.legendre.P2`.

    **Function-call arguments** \n

    :param int num_lat:         number of latitude points [default: 90]
    :param int num_lat:         (optional) number of longitude points [default: None]
    :param float water_depth:   depth of the slab ocean in meters [default: 10.]
    :param float T0:            global-mean initial temperature in :math:`^{\circ} \\textrm{C}` [default: 12.]
    :param float T2:            2nd Legendre coefficient for equator-to-pole gradient in
                                initial temperature, in :math:`^{\circ} \\textrm{C}` [default: -40.]

    :returns:                   dictionary with temperature
                                :class:`~climlab.domain.field.Field`
                                for surface mixed layer ``Ts``
    :rtype:                     dict


    :Example:

        ::

            >>> from climlab.domain import initial
            >>> import numpy as np

            >>> T_dict = initial.surface_state(num_lat=36)

            >>> print np.squeeze(T_dict['Ts'])
            [-27.88584094 -26.97777479 -25.18923361 -22.57456133 -19.21320344
             -15.20729309 -10.67854785  -5.76457135  -0.61467228   4.61467228
               9.76457135  14.67854785  19.20729309  23.21320344  26.57456133
              29.18923361  30.97777479  31.88584094  31.88584094  30.97777479
              29.18923361  26.57456133  23.21320344  19.20729309  14.67854785
               9.76457135   4.61467228  -0.61467228  -5.76457135 -10.67854785
             -15.20729309 -19.21320344 -22.57456133 -25.18923361 -26.97777479
             -27.88584094]