diff --git a/CMakeLists.txt b/CMakeLists.txt
index 14dd980..16fe2f9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -92,6 +92,7 @@ include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/fmt.cmake)
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/gtest.cmake)
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/yamlpp.cmake)
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/eigen.cmake)
+include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/libarchive.cmake)
 
 if (CMAKE_CUDA_COMPILER)
   find_package(CUDAToolkit REQUIRED)
diff --git a/cmake/libarchive.cmake b/cmake/libarchive.cmake
new file mode 100644
index 0000000..8fbb746
--- /dev/null
+++ b/cmake/libarchive.cmake
@@ -0,0 +1,19 @@
+include(FetchContent)
+
+# ---- libarchive options (tune as needed) ----
+set(ENABLE_WERROR OFF CACHE BOOL "" FORCE)
+set(ENABLE_TEST OFF CACHE BOOL "" FORCE)
+
+set(ENABLE_OPENSSL OFF CACHE BOOL "" FORCE)
+set(ENABLE_LZMA OFF CACHE BOOL "" FORCE)
+set(ENABLE_ZSTD OFF CACHE BOOL "" FORCE)
+set(ENABLE_BZip2 OFF CACHE BOOL "" FORCE)
+
+set(PACKAGE_NAME libarchive)
+set(REPO_URL "https://github.com/libarchive/libarchive")
+set(REPO_TAG "v3.8.5")
+
+add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "" ON)
+set(LIBARCHIVE_INCLUDE_DIR
+   "${${PACKAGE_NAME}_SOURCE_DIR}/libarchive"
+   CACHE PATH "libarchive include directory")
diff --git a/examples/run_hydro.cpp b/examples/run_hydro.cpp
index 00005d0..2f2ecd4 100644
--- a/examples/run_hydro.cpp
+++ b/examples/run_hydro.cpp
@@ -135,7 +135,7 @@ int main(int argc, char **argv) {
   // initialize
   std::map<std::string, torch::Tensor> vars;
   vars["hydro_w"] = w;
-  block->initialize(vars);
+  double current_time = block->initialize(vars, cli->restart_filename);
 
   // user output variables
   // (1) total precipitable mass fraction [kg/kg]
@@ -152,8 +152,9 @@ int main(int argc, char **argv) {
   kinet->to(device);
 
   // time loop
-  double current_time = 0.;
-  block->make_outputs(vars, current_time);
+  if (cli->restart_filename == nullptr) {
+    block->make_outputs(vars, current_time);
+  }
 
   while (!block->pintg->stop(block->cycle++, current_time)) {
     auto dt = block->max_time_step(vars);
diff --git a/examples/straka.cpp b/examples/straka.cpp
index 7de6ae9..4b7bbed 100644
--- a/examples/straka.cpp
+++ b/examples/straka.cpp
@@ -74,7 +74,8 @@ int main(int argc, char** argv) {
   // initialize
   std::map<std::string, torch::Tensor> vars;
   vars["hydro_w"] = w;
-  block->initialize(vars);
+  char const* restart = nullptr;
+  double current_time = block->initialize(vars, restart);
 
   block->user_output_callback = [Rd, cp, p0](Variables const& vars) {
     auto w = vars.at("hydro_w");
@@ -86,8 +87,9 @@ int main(int argc, char** argv) {
     return out;
   };
 
-  double current_time = 0.;
-  block->make_outputs(vars, current_time);
+  if (restart == nullptr) {
+    block->make_outputs(vars, current_time);
+  }
 
   while (!block->pintg->stop(block->cycle++, current_time)) {
     auto dt = block->max_time_step(vars);
diff --git a/pyproject.toml b/pyproject.toml
index 5f2643f..9d0e9a0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,7 +5,7 @@ requires = [
   "build",
   "torch>=2.7.0,<=2.7.1",
   "ninja",
-  "kintera>=1.2.6",
+  "kintera>=1.2.9",
   ]
 build-backend = "setuptools.build_meta"
 
@@ -21,7 +21,7 @@ readme = "README.md"
 requires-python = ">=3.9"
 dependencies = [
   "numpy",
-  "kintera>=1.2.6",
+  "kintera>=1.2.9",
 ]
 classifiers = [
   "Development Status :: 3 - Alpha",
diff --git a/python/csrc/pymesh.cpp b/python/csrc/pymesh.cpp
index 69f65cf..0d5c8b0 100644
--- a/python/csrc/pymesh.cpp
+++ b/python/csrc/pymesh.cpp
@@ -130,12 +130,19 @@ void bind_mesh(py::module &m) {
           py::arg("offset"), py::arg("exterior") = true,
           py::arg("extend_x1") = 0, py::arg("extend_x2") = 0,
           py::arg("extend_x3") = 0)
-      //.def("initialize", &snap::MeshBlockImpl::initialize)
-      .def("initialize",
-           [](snap::MeshBlockImpl &self, snap::Variables &vars) {
-             auto time = self.initialize(vars);
-             return std::make_pair(vars, time);
-           })
+      .def(
+          "initialize",
+          [](snap::MeshBlockImpl &self, snap::Variables &vars,
+             std::string restart_file) {
+            double time;
+            if (restart_file.empty()) {
+              time = self.initialize(vars);
+            } else {
+              time = self.initialize(vars, restart_file.c_str());
+            }
+            return std::make_pair(vars, time);
+          },
+          py::arg("vars"), py::arg("restart_file") = "")
       .def("print_cycle_info", &snap::MeshBlockImpl::print_cycle_info)
       .def("finalize", &snap::MeshBlockImpl::finalize)
       .def("device", &snap::MeshBlockImpl::device)
diff --git a/python/csrc/snapy.cpp b/python/csrc/snapy.cpp
index 4d2baf8..761e4a3 100644
--- a/python/csrc/snapy.cpp
+++ b/python/csrc/snapy.cpp
@@ -8,6 +8,7 @@
 #include <snap/snap.h>
 
 #include <snap/input/parameter_input.hpp>
+#include <snap/input/read_restart_file.hpp>
 
 // python
 #include "pyoptions.hpp"
@@ -39,6 +40,14 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
   m.attr("kIPR") = (int)snap::IPR;
   m.attr("kICY") = (int)snap::ICY;
 
+  m.def(
+      "load_restart",
+      [](snap::Variables &vars, std::string const &path) {
+        snap::load_restart(vars, path);
+        return vars;
+      },
+      py::arg("vars"), py::arg("path"));
+
   bind_layout(m);
   bind_bc(m);
   bind_coord(m);
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index dd52e2b..384b51a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -19,6 +19,7 @@ include_directories(
   ${CMAKE_BINARY_DIR}
   ${NETCDF_INCLUDES}
   ${PNETCDF_INCLUDE_DIR}
+  SYSTEM ${LIBARCHIVE_INCLUDE_DIR}
   SYSTEM ${TORCH_API_INCLUDE_DIR}
   SYSTEM ${TORCH_INCLUDE_DIR}
   SYSTEM ${KINTERA_INCLUDE_DIR}
@@ -73,6 +74,7 @@ target_link_libraries(${namel}_${buildl}
     ${TORCH_LIBRARY}
     ${TORCH_CPU_LIBRARY}
     ${C10_LIBRARY}
+    archive
     fmt::fmt
     yaml-cpp::yaml-cpp
     )
diff --git a/src/eos/ideal_moist.hpp b/src/eos/ideal_moist.hpp
index 112fecc..b4afed6 100644
--- a/src/eos/ideal_moist.hpp
+++ b/src/eos/ideal_moist.hpp
@@ -40,7 +40,7 @@ class IdealMoistImpl final : public torch::nn::Cloneable<IdealMoistImpl>,
   torch::Tensor compute(std::string ab,
                         std::vector<torch::Tensor> const& args) override;
 
-  //! \brief Inverse of the mean molecular weight
+  //! \brief Virtual factor
   /*!
    *! Eq.16 in Li2019
    *! $ \frac{R}{R_d} = \frac{\mu_d}{\mu}$
diff --git a/src/input/read_restart_file.cpp b/src/input/read_restart_file.cpp
new file mode 100644
index 0000000..4953c0f
--- /dev/null
+++ b/src/input/read_restart_file.cpp
@@ -0,0 +1,274 @@
+// C/C++
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <optional>
+#include <random>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <vector>
+
+// archive
+#include <archive.h>
+#include <archive_entry.h>
+
+// kintera
+#include <kintera/utils/serialize.hpp>
+
+// snap
+#include <snap/layout/layout.hpp>
+
+namespace fs = std::filesystem;
+
+namespace snap {
+
+// -------------------------
+// Small helpers
+// -------------------------
+
+struct RestartFields {
+  std::string basename;
+  std::string blockid;
+  std::string filenumber;
+};
+
+static RestartFields parse_part_filename(const std::string& name) {
+  constexpr std::string_view suffix = ".part";
+
+  if (name.size() <= suffix.size() ||
+      name.compare(name.size() - suffix.size(), suffix.size(), suffix) != 0) {
+    throw std::invalid_argument("filename does not end with .part");
+  }
+
+  // Strip ".part"
+  const std::string_view core(name.data(), name.size() - suffix.size());
+
+  // Find last two dots
+  const size_t dot2 = core.rfind('.');
+  if (dot2 == std::string::npos) {
+    throw std::invalid_argument("filename missing filenumber field");
+  }
+
+  const size_t dot1 = core.rfind('.', dot2 - 1);
+  if (dot1 == std::string::npos) {
+    throw std::invalid_argument("filename missing block_id field");
+  }
+
+  RestartFields out;
+  out.basename = std::string(core.substr(0, dot1));
+  out.blockid = std::string(core.substr(dot1 + 1, dot2 - dot1 - 1));
+  out.filenumber = std::string(core.substr(dot2 + 1));
+
+  if (out.basename.empty() || out.blockid.empty() || out.filenumber.empty()) {
+    throw std::invalid_argument("one or more filename fields are empty");
+  }
+
+  return out;
+}
+
+static std::string dtype_to_string(const at::ScalarType t) {
+  // at::toString exists in many builds; this is safe enough.
+  return std::string(at::toString(t));
+}
+
+static std::string device_to_string(const at::Device& d) {
+  std::ostringstream oss;
+  oss << d;
+  return oss.str();
+}
+
+static std::string shape_to_string(const at::Tensor& t) {
+  std::ostringstream oss;
+  oss << "(";
+  for (int64_t i = 0; i < t.dim(); ++i) {
+    oss << t.size(i);
+    if (i + 1 < t.dim()) oss << ", ";
+  }
+  oss << ")";
+  return oss.str();
+}
+
+static bool ends_with(std::string const& s, std::string const& suffix) {
+  return s.size() >= suffix.size() &&
+         s.compare(s.size() - suffix.size(), suffix.size(), suffix) == 0;
+}
+
+static bool is_tar_archive(std::string const& path) {
+  if (!fs::is_regular_file(path)) return false;
+
+  struct archive* ar = archive_read_new();
+  if (!ar) return false;
+
+  archive_read_support_filter_all(ar);
+  archive_read_support_format_all(ar);
+
+  // Try opening as an archive; if it succeeds, treat as tar-like.
+  int r = archive_read_open_filename(ar, path.c_str(), 10240);
+  if (r != ARCHIVE_OK) {
+    archive_read_free(ar);
+    return false;
+  }
+
+  // Some files might be recognized as other archive formats too; in practice
+  // this matches Python's "is_tarfile" intent well.
+  archive_read_close(ar);
+  archive_read_free(ar);
+  return true;
+}
+
+// Create a unique temp file path; not bulletproof, but good enough
+static fs::path make_temp_path(std::string_view suffix) {
+  fs::path dir = fs::temp_directory_path();
+
+  std::random_device rd;
+  std::mt19937_64 gen(rd());
+  std::uniform_int_distribution<uint64_t> dis;
+
+  for (int tries = 0; tries < 20; ++tries) {
+    uint64_t r = dis(gen);
+    std::ostringstream name;
+    name << "tmp_" << std::hex << r << suffix;
+    fs::path p = dir / name.str();
+    if (!fs::exists(p)) return p;
+  }
+
+  // Fallback (very unlikely to collide)
+  return dir / ("tmp_fallback" + std::string(suffix));
+}
+
+static void load_pt_from_tar(Variables& vars, struct archive* ar,
+                             struct archive_entry* entry) {
+  const char* name_c = archive_entry_pathname(entry);
+  std::string member_name =
+      name_c ? std::string{name_c} : std::string{"<unknown>"};
+
+  // Extract this entry into a temporary file (TorchScript loader prefers
+  // real/seekable file)
+  fs::path tmp_path = make_temp_path(".part");
+  std::ofstream out(tmp_path, std::ios::binary);
+  if (!out) {
+    std::cerr << "\n=== " << member_name << " ===\n";
+    std::cerr << "  ERROR: could not create temp file: " << tmp_path.string()
+              << "\n";
+    // Must still consume/skip entry data:
+    archive_read_data_skip(ar);
+    return;
+  }
+
+  std::vector<char> buf(1 << 20);
+  while (true) {
+    la_ssize_t n = archive_read_data(ar, buf.data(), buf.size());
+    if (n == 0) break;  // end of this entry
+    if (n < 0) {
+      std::cerr << "\n=== " << member_name << " ===\n";
+      std::cerr << "  ERROR: could not extract file from tar: "
+                << archive_error_string(ar) << "\n";
+      out.close();
+      std::error_code ec;
+      fs::remove(tmp_path, ec);
+      return;
+    }
+    out.write(buf.data(), static_cast<std::streamsize>(n));
+    if (!out) {
+      std::cerr << "\n=== " << member_name << " ===\n";
+      std::cerr << "  ERROR: failed writing temp file\n";
+      out.close();
+      std::error_code ec;
+      fs::remove(tmp_path, ec);
+      return;
+    }
+  }
+
+  out.flush();
+  out.close();
+
+  // load the extracted .part
+  kintera::load_tensors(vars, tmp_path.string());
+
+  // remove empty tensors (if any)
+  for (auto it = vars.begin(); it != vars.end();) {
+    if (!it->second.defined() || it->second.numel() == 0) {
+      it = vars.erase(it);
+    } else {
+      ++it;
+    }
+  }
+
+  // Cleanup
+  std::error_code ec;
+  fs::remove(tmp_path, ec);
+}
+
+void load_restart(Variables& vars, std::string const& path) {
+  // Dispatch based on whether `path` is a .part file or a tar archive.
+  if (is_tar_archive(path)) {
+    struct archive* ar = archive_read_new();
+    if (!ar) {
+      std::cerr << path << ": failed to allocate archive reader\n";
+      return;
+    }
+
+    archive_read_support_filter_all(ar);
+    archive_read_support_format_all(ar);
+
+    int r = archive_read_open_filename(ar, path.c_str(), 10240);
+    if (r != ARCHIVE_OK) {
+      std::cerr << path
+                << ": failed to open archive: " << archive_error_string(ar)
+                << "\n";
+      archive_read_free(ar);
+      return;
+    }
+
+    bool found_part = false;
+
+    struct archive_entry* entry = nullptr;
+    while ((r = archive_read_next_header(ar, &entry)) == ARCHIVE_OK) {
+      const char* name_c = archive_entry_pathname(entry);
+      std::string name = name_c ? std::string{name_c} : std::string{};
+
+      if (ends_with(name, ".part")) {
+        found_part = true;
+        auto out = parse_part_filename(name);
+        // find the block rank number after "block"
+        int rank = std::stoi(out.blockid.substr(5, out.blockid.size() - 5));
+        int my_rank = get_rank();
+        if (rank != my_rank) {
+          // Not for this rank; skip
+          archive_read_data_skip(ar);
+        } else {
+          load_pt_from_tar(vars, ar, entry);
+          return;
+        }
+
+        // Note: consume the entry data (via archive_read_data*)
+        // or skip it, otherwise the next header read will misbehave.
+      } else {
+        // Skip non-.part entries quickly
+        archive_read_data_skip(ar);
+      }
+    }
+
+    if (!found_part) {
+      std::cerr << path << ": no .part files found in tar archive\n";
+    }
+
+    if (r != ARCHIVE_EOF && r != ARCHIVE_OK) {
+      std::cerr << path
+                << ": error while reading archive: " << archive_error_string(ar)
+                << "\n";
+    }
+
+    archive_read_close(ar);
+    archive_read_free(ar);
+  } else {
+    // Treat as a single .part TorchScript file
+    std::cout << "single .part file detected\n";
+    kintera::load_tensors(vars, path);
+  }
+}
+
+}  // namespace snap
diff --git a/src/input/read_restart_file.hpp b/src/input/read_restart_file.hpp
index 89a71b2..2341c1b 100644
--- a/src/input/read_restart_file.hpp
+++ b/src/input/read_restart_file.hpp
@@ -5,22 +5,14 @@
 #include <cstring>
 #include <string>
 
-// snap
-#include <snap/mesh/meshblock.hpp>
-
-// kintera
-#include <kintera/utils/serialize.hpp>
+// torch
+#include <torch/script.h>
+#include <torch/torch.h>
 
 namespace snap {
 
-void set_hydro_interior(MeshBlockImpl *block, torch::Tensor &hydro_w,
-                        Variables &in_vars) {
-  auto interior = block->part({0, 0, 0});
-  hydro_w.index(interior)[IDN] = in_vars["rho"];
-  hydro_w.index(interior)[IVX] = in_vars["vel1"];
-  hydro_w.index(interior)[IVY] = in_vars["vel2"];
-  hydro_w.index(interior)[IVZ] = in_vars["vel3"];
-  hydro_w.index(interior)[IPR] = in_vars["press"];
-}
+using Variables = std::map<std::string, torch::Tensor>;
+
+void load_restart(Variables& vars, std::string const& path);
 
 }  // namespace snap
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index b41cc16..bde508e 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -5,10 +5,8 @@
 #include <limits>
 #include <mutex>
 
-// kintera
-#include <kintera/utils/serialize.hpp>
-
 // snap
+#include <snap/input/read_restart_file.hpp>
 #include <snap/output/output_formats.hpp>
 #include <snap/utils/log.hpp>
 #include <snap/utils/signal_handler.hpp>
@@ -268,7 +266,7 @@ std::vector<torch::indexing::TensorIndex> MeshBlockImpl::part(
   }
 }
 
-double MeshBlockImpl::initialize(Variables& vars) {
+double MeshBlockImpl::initialize(Variables& vars, char const* restart_file) {
   /*c10d::BarrierOptions op;
   op.device_ids = {options->layout()->local_rank()};
   _playout->pg->barrier(op)->wait();*/
@@ -277,8 +275,8 @@ double MeshBlockImpl::initialize(Variables& vars) {
   //// ------------ (1) Set up a signal handler ------------ ////
   SignalHandler::GetInstance();
 
-  if (pintg->options->restart() != "") {
-    return _init_from_restart(vars);
+  if (restart_file != nullptr) {
+    return _init_from_restart(vars, std::string(restart_file));
   }
 
   BoundaryFuncOptions bops;
@@ -828,46 +826,56 @@ torch::Device MeshBlockImpl::device() const {
   }
 }
 
-double MeshBlockImpl::_init_from_restart(Variables& vars) {
-  // create filename: <file_basename>.<block_id>.<fid>.restart
-  std::string fid = pintg->options->restart();
-
-  std::string fname;
-  char bid[12];
-  snprintf(bid, sizeof(bid), "block%d", options->layout()->rank());
-
-  fname.append(options->basename());
-  fname.append(".");
-  fname.append(bid);
-  fname.append(".");
-  fname.append(fid);
-  fname.append(".restart");
+double MeshBlockImpl::_init_from_restart(Variables& vars, std::string fname) {
+  // timing data
+  vars["last_time"] = torch::Tensor();
+  vars["last_cycle"] = torch::Tensor();
+  vars["file_number"] = torch::Tensor();
+  vars["next_time"] = torch::Tensor();
 
-  // load variables from disk
-  kintera::load_tensors(vars, fname);
+  // hydro data
+  vars["hydro_u"] = torch::Tensor();
+  vars["hydro_w"] = torch::Tensor();
 
-  // load auxiliary timing data from disk
-  Variables timing_vars;
+  // solid data
+  vars["solid"] = torch::Tensor();
+  vars["fill_solid_hydro_w"] = torch::Tensor();
+  vars["fill_solid_hydro_u"] = torch::Tensor();
 
-  timing_vars["last_time"] = torch::Tensor();
-  timing_vars["last_cycle"] = torch::Tensor();
-  timing_vars["file_number"] = torch::Tensor();
-  timing_vars["next_time"] = torch::Tensor();
+  // scalar data
+  if (pscalar->nvar() > 0) {
+    vars["scalar_s"] = torch::Tensor();
+    vars["scalar_r"] = torch::Tensor();
+  }
 
-  kintera::load_tensors(timing_vars, fname);
+  load_restart(vars, fname);
 
-  cycle = timing_vars.at("last_cycle").item<int64_t>();
+  cycle = vars.at("last_cycle").item<int64_t>() - 1;
   for (int n = 0; n < output_types.size(); ++n) {
-    output_types[n]->file_number =
-        timing_vars.at("file_number")[n].item<int64_t>();
-    output_types[n]->next_time = timing_vars.at("next_time")[n].item<double>();
+    output_types[n]->file_number = vars.at("file_number")[n].item<int64_t>();
+    output_types[n]->next_time = vars.at("next_time")[n].item<double>();
   }
 
   // start timing
   _time_start = clock();
   _cycle_start = cycle;
 
-  return timing_vars.at("last_time").item<double>();
+  auto current_time = vars.at("last_time").item<double>();
+
+  // move to device
+  for (auto& [name, tensor] : vars) {
+    if (tensor.defined()) {
+      tensor = tensor.to(device());
+    }
+  }
+
+  // remove timing data
+  vars.erase("last_time");
+  vars.erase("last_cycle");
+  vars.erase("file_number");
+  vars.erase("next_time");
+
+  return current_time;
 }
 
 }  // namespace snap
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index eb045d5..b6ab09a 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -126,7 +126,7 @@ class MeshBlockImpl : public torch::nn::Cloneable<MeshBlockImpl> {
    * \param vars: variables to initialize
    * \return: initial simulation time
    */
-  double initialize(Variables& vars);
+  double initialize(Variables& vars, char const* restart_file = nullptr);
 
   //! compute the maximum allowable time step
   /*!
@@ -177,9 +177,10 @@ class MeshBlockImpl : public torch::nn::Cloneable<MeshBlockImpl> {
   //! initialize from restart file
   /*!
    * \param vars: variables to initialize
-   * \return: simulation time from the restart file
+   * \param fname: restart filename
+   * \return: last simulation time from the restart file
    */
-  double _init_from_restart(Variables& vars);
+  double _init_from_restart(Variables& vars, std::string fname);
 
  private:
   //! clock and cycle at time start
diff --git a/src/output/load_diag_output_data.cpp b/src/output/load_diag_output_data.cpp
index 96dfac1..65e61d5 100644
--- a/src/output/load_diag_output_data.cpp
+++ b/src/output/load_diag_output_data.cpp
@@ -30,8 +30,10 @@ void OutputType::loadDiagOutputData(MeshBlockImpl* pmb, Variables const& vars) {
 
     int ny = thermo_y->options->species().size() - 1;
     auto temp = peos->compute("W->T", {w});
+    auto dens = w[IDN];
     auto pres = w[IPR];
-    auto xfrac = thermo_y->compute("Y->X", {w.narrow(0, ICY, ny)});
+    auto yfrac = w.narrow(0, ICY, ny);
+    auto xfrac = thermo_y->compute("Y->X", {yfrac});
 
     // mole concentration [mol/m^3]
     auto conc = thermo_x->compute("TPX->V", {temp, pres, xfrac});
@@ -57,6 +59,11 @@ void OutputType::loadDiagOutputData(MeshBlockImpl* pmb, Variables const& vars) {
     // potential temperature [K]
     auto theta = (entropy_vol / cp_vol).exp();
 
+    // virtual potential temperature [K]
+    auto Rgas = kintera::constants::Rgas / mu;
+    auto feps = pres / (dens * Rgas * temp);
+    auto theta_v = theta * feps;
+
     // relative humidity
     auto rh = kintera::relative_humidity(temp, conc, thermo_x->stoich,
                                          thermo_x->options->nucleation());
@@ -77,6 +84,14 @@ void OutputType::loadDiagOutputData(MeshBlockImpl* pmb, Variables const& vars) {
     AppendOutputDataNode(pod);
     num_vars_ += 1;
 
+    // virtual potential temperature
+    pod = new OutputData;
+    pod->type = "SCALARS";
+    pod->name = "theta_v";
+    pod->data.CopyFromTensor(theta);
+    AppendOutputDataNode(pod);
+    num_vars_ += 1;
+
     // entropy
     pod = new OutputData;
     pod->type = "SCALARS";