def load_dataframe(
    pfile, df, func='mean', units='unknown',
    dims=None, coords=False, keys=None
):
    """
    Arguments
    ---------
    pfile : PseudoNetCDFFile
        must have dimensions ('TSTEP', 'LAY', 'ROW', 'COL')
    df : pandas.DataFrame
        must have columns matching dimensions
    func : str
        Used to aggregate values within a unique coordinate
    units : str
        Default value for units
    dims : tuple or None
        If None, defaults to ('TSTEP', 'LAY', 'ROW', 'COL'), otherwise
        tuple of strings indicating dimensions for new variables
    coords : bool
        If True, copy dimension variables as IN_{dim}
    keys : list or None
        If None, all columns will be exported. If a list, only those
        keys will be output.

    Returns
    -------
    pfile : PseudoNetCDFFile
    """
    incols = df.columns
    # dims = pfile.dimensions
    dims = ('TSTEP', 'LAY', 'ROW', 'COL')
    gkeys = [k for k in incols if k in pfile.dimensions]
    dfgb = df.groupby(gkeys, as_index=False)
    gdf = getattr(dfgb, func)()
    vals = {k: gdf[k].values for k in gkeys}
    slices = {}
    for k, v in vals.items():
        if k in pfile.variables:
            slices[k] = pfile.val2idx(k, v)
        else:
            slices[k] = v.astype('i')

    slicer = tuple([
        slices.get(k, slice(None))
        for k in dims
    ])

    if keys is None:
        keys = incols

    for outkey in keys:
        outval = gdf[outkey]
        if outkey in gkeys:
            if not coords:
                continue

            outkey = f'IN_{outkey}'
        var = pfile.createVariable(
            outkey, outval.dtype.char, dims,
            units=units, long_name=outkey.ljust(16), var_desc=outkey.ljust(80)
        )
        var[slicer] = outval

    return pfile