def calc_baro_stf(ds,grid):
    """
    Compute barotropic streamfunction as

        u_bt = - d\psi / dy
        \psi = \int_{-90}^{90} u_bt dy

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

    Output
    ------
    baro_stf    :   xarray DataArray containing 2D field with 
                    barotropic streamfunction in Sv, \psi above
    """


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

    # Define barotropic velocity as vertically integrated velocity
    if ustr == 'UVELMASS':
        u_bt = (ds[ustr] * ds['drF']).sum(dim='Z')
    else:
        u_bt = (ds[ustr] * ds['hFacW'] * ds['drF']).sum(dim='Z')

    u_bt = u_bt * ds['dyG']

    # Integrate in Y
    baro_stf = grid.cumsum(-u_bt,'Y',boundary='fill',fill_value=0.)

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

    return baro_stf