def rho_moist(T,p,q):
e = vapor_pressure_from_specific_humidity(p,q)<EOL>w = mixing_ratio_from_vapor_pressure(p,e)<EOL>Tv = virtual_temperature_from_mixing_ratio(T,w)<EOL>return p*<NUM_LIT>/ Rd/Tv<EOL>
Density of moist air.
    T is air temperature (K)
    p is air pressure (hPa)
    q is specific humidity (hPa)

    returns density in kg/m3