def modis_bs_alb_mcd43b1(fiso, fvol, fgeo, sza):
    """
    Calculate black-sky albedo from MODIS BRDF parameters.

    Parameters:
    fiso, fvol, fgeo, sza : numpy arrays of the same dimension (sza in radians)

    Returns:
    bs_alb : numpy array
    """
    sza2 = np.power(sza, 2)
    sza3 = sza2 * sza
    bs_alb = (fiso + 
              fvol * (-0.007574 - 0.070987 * sza2 + 0.307588 * sza3) + 
              fgeo * (-1.284909 - 0.166314 * sza2 + 0.041840 * sza3))
    return bs_alb