[BUG] Failed to calculate BRN when I use Mixed layer parcel

[alment-Li]

When I try to calculate BRN using the Mixed layer parcel, the returned values are all masked. After checking the source code, I found that the logic defining wind shear in the bulk_rich function seems problematic. Why are different wind shear stratifications used for air parcels with flags 1–3 and those with flags ≥ 4? It does not meet the definition of BRN specified in the article "Close Proximity Soundings within Supercell Environments Obtained from the Rapid Update Cycle".

def bulk_rich(prof, pcl):
'''
Calculates the Bulk Richardson Number for a given parcel.

    Parameters
    ----------
    prof : profile object
        Profile object
    pcl : parcel object
        Parcel object

    Returns
    -------
    Bulk Richardson Number : number

    '''

# Make sure parcel is initialized
if not utils.QC(pcl.lplvals):
    pbot = ma.masked
elif pcl.lplvals.flag > 0 and pcl.lplvals.flag <= 4:
    ptop = interp.pres(prof, interp.to_msl(prof, 6000.))
    pbot = prof.pres[prof.sfc]
else:
    h0 = interp.hght(prof, pcl.pres)
    try:
        pbot = interp.pres(prof, h0-500.)
    except:
        pbot = ma.masked
    if utils.QC(pbot): pbot = prof.pres[prof.sfc]
    h1 = interp.hght(prof, pbot)
    ptop = interp.pres(prof, h1+6000.)

if not utils.QC(pbot) or not utils.QC(ptop):
    pcl.brnshear = ma.masked
    pcl.brn = ma.masked
    pcl.brnu = ma.masked
    pcl.brnv = ma.masked
    return pcl

# Calculate the lowest 500m mean wind
p = interp.pres(prof, interp.hght(prof, pbot)+500.)
#print(p, pbot)
mnlu, mnlv = winds.mean_wind(prof, pbot, p)

# Calculate the 6000m mean wind
mnuu, mnuv = winds.mean_wind(prof, pbot, ptop)

# Make sure CAPE and Shear are available
if not utils.QC(pcl.bplus) or not utils.QC(mnlu) or not utils.QC(mnuu):
    pcl.brnshear = ma.masked
    pcl.brnu = ma.masked
    pcl.brnv = ma.masked
    pcl.brn = ma.masked
    return pcl

# Calculate shear between levels
dx = mnuu - mnlu
dy = mnuv - mnlv
pcl.brnu = dx
pcl.brnv = dy
pcl.brnshear = utils.KTS2MS(utils.mag(dx, dy))
pcl.brnshear = pcl.brnshear**2 / 2.
pcl.brn = pcl.bplus / pcl.brnshear
return pcl