Small function and saved constants for convective velocity

The code replicates convection velocity in several modules as well as having 0.333_r8 hardcoded as a power used in six different modules in the code. Having a function that can be used allows changes to it to be maintained across the code base, and do a better job of documenting what's happening, and ensuring the code is consistent throughout.

```
biogeochem/DustEmisLeung2023.F90: u_sd_slt(p) = wnd_frc_slt * (12_r8 - 0.5_r8 * stblty(p))**0.333_r8
biogeochem/DustEmisLeung2023.F90: u_sd_slt(p) = wnd_frc_slt * (0.001_r8)**0.333_r8 ! should have used 0 theoretically; used 0.001 here to avoid undefined values
biogeophys/BareGroundFluxesMod.F90: wc = beta(c)*(-grav*ustar(p)*thvstar*zii(c)/thv(c))**0.333_r8
biogeophys/CanopyFluxesMod.F90: wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
biogeophys/FrictionVelocityMod.F90: + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8))
biogeophys/FrictionVelocityMod.F90: + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8)) )
biogeophys/FrictionVelocityMod.F90: + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8)) )
biogeophys/FrictionVelocityMod.F90: + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
biogeophys/FrictionVelocityMod.F90: + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
biogeophys/FrictionVelocityMod.F90: + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
biogeophys/FrictionVelocityMod.F90: + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
biogeophys/FrictionVelocityMod.F90: + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
biogeophys/LakeFluxesMod.F90: wc = beta1*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
biogeophys/UrbanFluxesMod.F90: wc = beta(l)*(-grav*ustar(l)*thvstar*zii(l)/thv_g(l))**0.333_r8
```