docs: def find_icelines(self):
          """Finds iceline according to the surface temperature.
  
          This method is called by the private function
          :func:`~climlab.surface.albedo.Iceline._compute`
          and updates following attributes according to the freezing temperature
          ``self.param['Tf']`` and the surface temperature ``self.param['Ts']``:
  
          **Object attributes** \n
  
          :ivar Field noice:      a Field of booleans which are ``True`` where
                                  :math:`T_s \\ge T_f`
          :ivar Field ice:        a Field of booleans which are ``True`` where
                                  :math:`T_s < T_f`
          :ivar array icelat:     an array with two elements indicating the
                                  ice-edge latitudes
          :ivar float ice_area:   fractional area covered by ice (0 - 1)
          :ivar dict diagnostics: keys ``'icelat'`` and ``'ice_area'`` are updated
  
          """
          Tf = self.param['Tf']
          Ts = self.state['Ts']
          lat_bounds = self.domains['Ts'].axes['lat'].bounds
          self.noice = np.where(Ts >= Tf, True, False)
          self.ice = np.where(Ts < Tf, True, False)
          #  Ice cover in fractional area
          self.ice_area = global_mean(self.ice * np.ones_like(self.Ts))
          #  Express ice cover in terms of ice edge latitudes
          if self.ice.all():
              # 100% ice cover
              self.icelat = np.array([-0., 0.])
          elif self.noice.all():
              # zero ice cover
              self.icelat = np.array([-90., 90.])
          else:  # there is some ice edge
              # Taking np.diff of a boolean array gives True at the boundaries between True and False
              boundary_indices = np.where(np.diff(self.ice.squeeze()))[0]+1
              # check for asymmetry case: [-90,x] or [x,90]
              #  -> boundary_indices hold only one value for icelat
              if boundary_indices.size == 1:
                  if self.ice[0] == True:   # case: [x,90]
                      # extend indice array by missing value for northpole
                      boundary_indices = np.append(boundary_indices, self.ice.size)
                  elif  self.ice[-1] == True:   # case: [-90,x]
                      # extend indice array by missing value for northpole
                      boundary_indices = np.insert(boundary_indices,0 ,0)
              # check for asymmetry case: [-90,x] or [x,90]
              #  -> boundary_indices hold only one value for icelat
              if boundary_indices.size == 1:
                  if self.ice[0] == True:   # case: [x,90]
                      # extend indice array by missing value for northpole
                      boundary_indices = np.append(boundary_indices, self.ice.size)
                  elif  self.ice[-1] == True:   # case: [-90,x]
                      # extend indice array by missing value for northpole
                      boundary_indices = np.insert(boundary_indices,0 ,0)
              self.icelat = lat_bounds[boundary_indices]