docs: 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