diff --git a/src/pgen/disk.hpp b/src/pgen/disk.hpp
index 7c34ea3..65247de 100644
--- a/src/pgen/disk.hpp
+++ b/src/pgen/disk.hpp
@@ -39,6 +39,7 @@
 #include "nbody/nbody.hpp"
 #include "utils/artemis_utils.hpp"
 #include "utils/eos/eos.hpp"
+#include "utils/units.hpp"
 
 using ArtemisUtils::EOS;
 using ArtemisUtils::VI;
@@ -60,6 +61,7 @@ struct DiskParams {
   Real alpha, nu0, nu_indx;
   Real mdot;
   Real temp_soft2;
+  Real kbmu;
   bool do_dust;
   bool nbody_temp;
   bool quiet_start;
@@ -101,7 +103,8 @@ Real TempProfile(struct DiskParams pgen, const Real R, const Real z) {
   const Real H = R * pgen.h0 * std::pow(R / pgen.r0, pgen.flare);
   const Real ir1 = 1.0 / std::sqrt(R * R + pgen.temp_soft2);
   const Real omk2 = SQR(pgen.Omega0) * ir1 * ir1 * ir1;
-  const Real T0 = omk2 * H * H / pgen.Gamma;
+  // c_iso^2 = P/rho = kb/mu T = Omk^2 H^2
+  const Real T0 = omk2 * H * H / (pgen.kbmu * pgen.Gamma);
   return T0 * std::pow(rho / rho0, pgen.Gamma - 1.0);
 }
 
@@ -275,6 +278,9 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) {
     disk_params.l0 = pin->GetOrAddReal("problem", "l0", 0.0);
     disk_params.dust_to_gas = pin->GetOrAddReal("problem", "dust_to_gas", 0.01);
     disk_params.temp_soft2 = pin->GetOrAddReal("problem", "temp_soft", 0.0);
+    const auto mu = gas_pkg->Param<Real>("mu");
+    auto &constants = artemis_pkg->Param<ArtemisUtils::Constants>("constants");
+    disk_params.kbmu = constants.GetKBCode() / (mu * constants.GetAMUCode());
 
     disk_params.do_dust = params.Get<bool>("do_dust");
 
@@ -312,7 +318,7 @@ inline void InitDiskParams(MeshBlock *pmb, ParameterInput *pin) {
         disk_params.nu0 = disk_params.alpha * disk_params.gamma_gas *
                           SQR(disk_params.h0 * disk_params.r0 * disk_params.Omega0);
         disk_params.nu_indx = 1.5 + disk_params.q;
-      } else if (vtype == "powerlaw") {
+      } else if ((vtype == "powerlaw") || (vtype == "constant")) {
         disk_params.nu0 = pin->GetReal("gas/viscosity", "nu");
         disk_params.nu_indx = pin->GetOrAddReal("gas/viscosity", "r_exp", 0.0);
       } else {