<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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]