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