pure real(wp) function density(s) class(State),intent(in) :: s density = s%pressure * s%mean_molecular_weight / (Rstar * s%temperature) end function density