Add Brunt-Vaisala freq into Smag model

[AMLattanzi]

This PR brings in the Brunt-Vaisala scaling for the Smagorinsky model. It includes dry and moist computations of the fb scaling parameter.

NOTE: This will likely modify all the Smagorinsky regtests and the vertical gradients are not aware of the boundary conditions at z_hi or z_lo.

See the follwowing sources:
Sec 2.7.2 (14) - (15)
Approximation in (36)

[AMLattanzi]

This has stability issues with the MOST boundary condition. It is very easy for the fb parameter to become 0 and thus the mu_turb at the bottom surface is 0 and the MOST BC crashes.

[AMLattanzi]

To avoid the stability issue with MOST, I have added Qfx MultiFabs like what was done for Hfx with theta so we can use the erf.use_explicit_most pathway. This directly prescribes the fluxes and stores them. I have also taken the N^2 calculation from the WRF manual (below) and made this into an inline GPU function in the Utils header.

Finally, some non-negligible errors were found with the MOST boundary condition. I reformulated it to compute the covariance term directly (e.g. rho * <\phi'w'> = -K * d\phi/dz). This is consistent with what is done for momentum in defining \tau_{ij}. The scalar quantities now use the negative sign in their definition of the diffusive flux. There could be lingering issues though if some negative sign(s) were not caught...