def METHOD_NAME(var_name, var_file, lsm_name="landseamask", lsm_file=None):
    if lsm_file is None:
        lsm_file = var_file
    var_DataArray = xr.open_dataset(var_file)
    lsm_DataArray = xr.open_dataset(lsm_file)
    var, lsm = var_DataArray[var_name], lsm_DataArray[lsm_name]
    assert var.dims[-1] == "landpoint"
    if var.ndim == 1:
        var_unpacked = np.empty((lsm.data.shape))
    elif var.ndim == 2:
        var_unpacked = np.empty((var.shape[0],) + lsm.data.shape)
    elif var.ndim == 3:
        var_unpacked = np.empty(
            (
                var.shape[0],
                var.shape[1],
            )
            + lsm.data.shape
        )
    else:
        print(
            "Opps, your array has more dimensions than unpack_jsbach.py knows how to handle! Goodbye!"
        )
        sys.exit()
    var_unpacked[:] = np.nan
    mask = lsm.data.astype(bool)
    if var.ndim == 1:
        np.place(var_unpacked[:, :], mask, var[:])
    elif var.ndim == 2:
        for d1 in range(var.shape[0]):
            np.place(var_unpacked[d1, :, :], mask, var[d1, :])
    elif var.ndim == 3:
        for d1 in range(var.shape[0]):
            for d2 in range(var.shape[1]):
                np.place(var_unpacked[d1, d2, :, :], mask, var[d1, d2, :])
    return var_unpacked