docs: def surface_to_array(surface):
      """
      :param surface: a Surface object
      :returns: a 3D array of shape (3, N, M)
      """
      if hasattr(surface, 'surfaces'):  # multiplanar surfaces
          n = len(surface.surfaces)
          arr = numpy.zeros((3, n, 4), F32)
          for i, surf in enumerate(surface.surfaces):
              arr[:, i] = surf.mesh.array
          return arr
      mesh = surface.mesh
      if len(mesh.lons.shape) == 1:  # 1D mesh
          shp = (3, 1) + mesh.lons.shape
      else:  # 2D mesh
          shp = (3,) + mesh.lons.shape
      return mesh.array.reshape(shp)