def METHOD_NAME(lat, lon, latvar, lonvar, n=1):
    """Find the nearest grid point to a given lat/lon pair.
    Parameters
    ----------
    lat : float
        Latitude value at which to find the nearest grid point.
    lon : float
        Longitude value at which to find the nearest grid point.
    latvar : xarray.DataArray
        DataArray corresponding to latitude variable.
    lonVar : xarray.DataArray
        DataArray corresponding to longitude variable.
    n : int, optional
        Number of nearest grid points to return. Default is to return the
        single closest grid point.
    Returns
    -------
    iy, ix, dist_sq
        A tuple of numpy arrays:
        - ``iy``: the y indices of the nearest grid points
        - ``ix``: the x indices of the nearest grid points
        - dist_sq: squared distance
    """
    # Note the use of the squeeze method: it removes single-dimensional entries
    # from the shape of an array. For example, in the GIOPS mesh file the
    # longitude of the U velocity points is defined as an array with shape
    # (1, 1, 1021, 1442). The squeeze method converts this into the equivalent
    # array with shape (1021, 1442).
    latvar = latvar.squeeze()
    lonvar = lonvar.squeeze()
    rad_factor = pi / 180.0
    latvals = latvar[:] * rad_factor
    lonvals = lonvar[:] * rad_factor
    clat, clon = np.cos(latvals), np.cos(lonvals)
    slat, slon = np.sin(latvals), np.sin(lonvals)
    if latvar.ndim == 1:
        # If latitude and longitude are 1D arrays (as is the case with the
        # GIOPS forecast data currently pulled from datamart), then we need to
        # handle this situation in a special manner. The clat array will be of
        # some size m, say, and the clon array will be of size n. By virtue of
        # being defined with different dimensions, the product of these two
        # arrays will be of size (m, n) because xarray will automatically
        # broadcast the arrays so that the multiplication makes sense to do.
        # Thus, the array calculated from
        #
        #   np.ravel(clat * clon)
        #
        # will be of size mn. However, the array
        #
        #   np.ravel(slat)
        #
        # will be of size m and this will cause the KDTree() call to fail. To
        # resolve this issue, we broadcast slat to the appropriate size and
        # shape.
        shape = (slat.size, slon.size)
        slat = np.broadcast_to(slat.values[:, np.newaxis], shape)
    else:
        shape = latvar.shape
    triples = np.array(
        [np.ravel(clat * clon), np.ravel(clat * slon), np.ravel(slat)]
    ).transpose()
    kdt = KDTree(triples)
    dist_sq, iy, ix = _find_index(lat, lon, kdt, shape, n)
    # The results returned from _find_index are two-dimensional arrays (if
    # n > 1) because it can handle the case of finding indices closest to
    # multiple lat/lon locations (i.e., where lat and lon are arrays, not
    # scalars). Currently, this function is intended only for a single lat/lon,
    # so we redefine the results as one-dimensional arrays.
    if n > 1:
        return iy, ix, dist_sq
    else:
        return int(iy), int(ix), dist_sq