def calc_overturning_stf(ds,grid,doFlip=True):
    """
    Only for simple domains, compute meridional overturning streamfunction

    Parameters
    ----------
    ds          :   xarray Dataset from MITgcm output, via 
                    e.g. xmitgcm.open_mdsdataset
                    must contain 'V' or 'VVELMASS' fields
    grid        :   xgcm grid object defined via xgcm.Grid(ds)

    doFlip      :   if true, compute by accumulating from bottom to top

    Output
    ------
    ov_stf      :   xarray DataArray containing 2D field with 
                    overturning streamfunction in Sv above
    """

    # Grab the right velocity field from dataset
    if 'V' in ds.keys():
        vstr = 'V'
    elif 'VVELMASS' in ds.keys():
        vstr = 'VVELMASS'
    else:
        raise TypeError('Could not find recognizable velocity field in input dataset')

    # Compute volumetric transport
    v_trsp = ds[vstr] * ds['dxG'] * ds['drF']
    if vstr != 'VVELMASS':
        print(f' *** Multiplying {vstr} by hFacS***')
        v_trsp = v_trsp * ds['hFacS']

    v_trsp = v_trsp.sum(dim=['XC'])

    # flip dim, accumulate in vertical, flip back
    if doFlip:
        v_trsp = v_trsp.isel(Z=slice(None,None,-1))

    ov_stf = grid.cumsum(v_trsp,'Z',boundary='fill')

    if doFlip:
        ov_stf = -ov_stf.isel(Zl=slice(None,None,-1))

    # Convert m/s to Sv
    ov_stf = ov_stf * 10**-6

    return ov_stf