<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def interpolate_to_isosurface(level_var, interp_var, level, **kwargs):
 r"""Linear interpolation of a variable to a given vertical level from given values. This function assumes that highest vertical level (lowest pressure) is zeroth index. A classic use of this function would be to compute the potential temperature on the dynamic tropopause (2 PVU surface). Parameters level_var: array_like (P, M, N) Level values in 3D grid on common vertical coordinate (e.g., PV values on isobaric levels). Assumes height dimension is highest to lowest in atmosphere. interp_var: array_like (P, M, N) Variable on 3D grid with same vertical coordinate as level_var to interpolate to given level (e.g., potential temperature on isobaric levels) level: int or float Desired interpolated level (e.g., 2 PVU surface) Other Parameters bottom_up_search : bool, optional Controls whether to search for levels bottom-up, or top-down. Defaults to True, which is bottom-up search. Returns ------- interp_level: (M, N) ndarray The interpolated variable (e.g., potential temperature) on the desired level (e.g., 2 PVU surface) Notes ----- This function implements a linear interpolation to estimate values on a given surface. The prototypical example is interpolation of potential temperature to the dynamic tropopause (e.g., 2 PVU surface) """

    # Change when Python 2.7 no longer supported
    # Pull out keyword arguments
    bottom_up_search = kwargs.pop('bottom_up_search', True)

    # Find index values above and below desired interpolated surface value
    above, below, good = metpy.calc.find_bounding_indices(level_var, [level], axis=0,
                                                          from_below=bottom_up_search)

    # Linear interpolation of variable to interpolated surface value
    interp_level = (((level - level_var[above]) / (level_var[below] - level_var[above]))
                    * (interp_var[below] - interp_var[above])) + interp_var[above]

    # Handle missing values and instances where no values for surface exist above and below
    interp_level[~good] = np.nan
    minvar = (np.min(level_var, axis=0) >= level)
    maxvar = (np.max(level_var, axis=0) <= level)
    interp_level[0][minvar] = interp_var[-1][minvar]
    interp_level[0][maxvar] = interp_var[0][maxvar]
    return interp_level.squeeze()