Bulk Richardson Number Calculation

[GGDLW]

I am wanting to make sure I calculate everything correctly. It is to my understanding that the BRN formula is as so:
BRN= CAPE/(1/2 (U ̅-U_0 )^2 ) with the denominator being the shear term, which is half of the the density-weighted mean wind speed of the 0-6000m minus the 0-500m mean wind speed squared. Would this calculation be correct?

prof = mpcalc.parcel_profile(p, T[0], Td[0])
cape, cin = mpcalc.cape_cin(p, T, Td, prof)
ubshr, vbshr = mpcalc.bulk_shear(p[1:], wind_u[1:], wind_v[1:], height=z[1:], depth=6 * units.km)
bshear_6km = mpcalc.wind_speed(ubshr, vbshr)
ubshr500, vbshr500 = mpcalc.bulk_shear(p[1:],wind_u[1:], wind_v[1:], height=z[1:], depth=0.5 * units.km)
bshear_500m = mpcalc.wind_speed(ubshr500, vbshr500)
brnshear = bshear_6km - bshear_500m
brn = cape.magnitude / (0.5 * brnshear.magnitude**2)

[sgdecker]

I believe you should subtract the wind components first (e.g., ubshr - ubshr500), and then do the sum of squares, as in SHARPpy.