diff --git a/.gitcustom/git-done b/.gitcustom/git-done
index c834451e..76a02546 100755
--- a/.gitcustom/git-done
+++ b/.gitcustom/git-done
@@ -5,4 +5,4 @@ currentBranch=$(git symbolic-ref --short -q HEAD) # Getting the current branch
 git checkout main
 git fetch
 git rebase origin/main
-git branch -D currentBranch
+git branch -D $currentBranch
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 3c2728e6..77d90d35 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -144,9 +144,6 @@ jobs:
     runs-on: [self-hosted, linux, 4core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -175,9 +172,6 @@ jobs:
     runs-on: [self-hosted, Linux, 8core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -190,6 +184,8 @@ jobs:
         run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j8
 
       - name: run simulation
+        env:
+          PATH: $PATH:/usr/lib64/mpich/bin
         run: cd ${{github.workspace}}/build/bin && mpiexec -n 8 ./robert.release -i robert2d.inp
 
       - name: combile outputs
@@ -206,9 +202,6 @@ jobs:
     runs-on: [self-hosted, Linux, 8core]
     needs: build-non-hydrostatic
     steps:
-      - name: source
-        run: source ~/.bashrc
-
       - uses: actions/checkout@v3
 
       - name: create build directory
@@ -221,6 +214,8 @@ jobs:
         run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j8
 
       - name: run simulation
+        env:
+          PATH: $PATH:/usr/lib64/mpich/bin
         run: cd ${{github.workspace}}/build/bin && mpiexec -n 8 ./robert.release -i robert3d.inp
 
       - name: combile outputs
diff --git a/examples/2023-jupiter-mwr-eq/juno.cpp b/examples/2023-jupiter-mwr-eq/juno.cpp
index ca14f05b..51c54bea 100644
--- a/examples/2023-jupiter-mwr-eq/juno.cpp
+++ b/examples/2023-jupiter-mwr-eq/juno.cpp
@@ -147,6 +147,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
   Real cp = gamma/(gamma - 1.)*Rd;
 
   // index
+  auto pindex = IndexMap::GetInstance();
   int iH2O = pindex->GetVaporId("H2O");
   int iNH3 = pindex->GetVaporId("NH3");
 
@@ -218,24 +219,6 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
       }
     }
 
-  // set tracers, electron and Na
-  int ielec = pindex->GetTracerId("e-");
-  int iNa = pindex->GetTracerId("Na");
-  auto ptracer = pimpl->ptracer;
-
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j)
-      for (int i = is; i <= ie; ++i) {
-        Real temp = pthermo->GetTemp(this,k,j,i);
-        Real pH2S = xH2S*phydro->w(IPR,k,j,i);
-        Real pNa = xNa*phydro->w(IPR,k,j,i);
-        Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
-        pNa = std::min(svp, pNa);
-
-        ptracer->u(iNa,k,j,i) = pNa/(Constants::kBoltz*temp);
-        ptracer->u(ielec,k,j,i) = saha_ionization_electron_density(temp, ptracer->u(iNa,k,j,i), 5.14);
-      }
-
   // if requested
   // modify atmospheric stability
   Real adlnTdlnP = pin->GetOrAddReal("problem", "adlnTdlnP", 0.);
@@ -275,6 +258,24 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
         }
   }
 
+  // set tracers, electron and Na
+  int ielec = pindex->GetTracerId("e-");
+  int iNa = pindex->GetTracerId("Na");
+  auto ptracer = pimpl->ptracer;
+
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        Real temp = pthermo->GetTemp(this,k,j,i);
+        Real pH2S = xH2S*phydro->w(IPR,k,j,i);
+        Real pNa = xNa*phydro->w(IPR,k,j,i);
+        Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
+        pNa = std::min(svp, pNa);
+
+        ptracer->u(iNa,k,j,i) = pNa/(Constants::kBoltz*temp);
+        ptracer->u(ielec,k,j,i) = saha_ionization_electron_density(temp, ptracer->u(iNa,k,j,i), 5.14);
+      }
+
   // primitive to conserved conversion
   peos->PrimitiveToConserved(phydro->w, pfield->bcc,
     phydro->u, pcoord, is, ie, js, je, ks, ke);
diff --git a/examples/2023-jupiter-mwr-eq/juno.inp b/examples/2023-jupiter-mwr-eq/juno.inp
index b2aa7815..fc61e7aa 100644
--- a/examples/2023-jupiter-mwr-eq/juno.inp
+++ b/examples/2023-jupiter-mwr-eq/juno.inp
@@ -116,7 +116,6 @@ xNa         = 3.97E-6
 xKCl        = 0.
 metallicity = -2
 
-#rdlnTdlnP   = 1.015
-rdlnTdlnP   = 0.985
-rdlnTdlnP.pmin = 1.E5
-rdlnTdlnP.pmax = 20.E5
+adlnTdlnP   = 0.015
+adlnTdlnP.pmin = 1.E5
+adlnTdlnP.pmax = 20.E5
diff --git a/src/harp/absorber.cpp b/src/harp/absorber.cpp
index 57e14cbe..64a365ed 100644
--- a/src/harp/absorber.cpp
+++ b/src/harp/absorber.cpp
@@ -1,6 +1,5 @@
 // athena
 #include <athena/athena.hpp>
-#include <athena/mesh/mesh.hpp>
 
 // canoe
 #include <index_map.hpp>
@@ -17,15 +16,16 @@ Absorber::Absorber(std::string name) : name_(name) {
   app->Log("Create Absorber " + name_);
 }
 
-Absorber::Absorber(MeshBlock* pmb, std::string name,
-                   std::vector<std::string> species,
+Absorber::Absorber(std::string name, std::vector<std::string> species,
                    std::map<std::string, Real> params)
     : name_(name), params_(params) {
   Application::Logger app("harp");
   app->Log("Create Absorber " + name_);
 
+  auto pindex = IndexMap::GetInstance();
+
   for (auto s : species) {
-    imols_.push_back(pmb->pindex->GetSpeciesId(s));
+    imols_.push_back(pindex->GetSpeciesId(s));
   }
 
   app->Log("Dependent species ids", imols_);
diff --git a/src/harp/absorber.hpp b/src/harp/absorber.hpp
index d538e5ee..c225ac3d 100644
--- a/src/harp/absorber.hpp
+++ b/src/harp/absorber.hpp
@@ -24,7 +24,6 @@
 #include <athena/athena.hpp>
 
 class Variable;
-class MeshBlock;
 
 using ParameterMap = std::map<std::string, Real>;
 
@@ -32,7 +31,7 @@ class Absorber {
  public:
   explicit Absorber(std::string name);
 
-  Absorber(MeshBlock* pmb, std::string name, std::vector<std::string> species,
+  Absorber(std::string name, std::vector<std::string> species,
            ParameterMap params);
 
   std::string GetName() const { return name_; }
diff --git a/src/harp/add_absorber_giants.cpp b/src/harp/add_absorber_giants.cpp
index a150fa19..103b26bf 100644
--- a/src/harp/add_absorber_giants.cpp
+++ b/src/harp/add_absorber_giants.cpp
@@ -45,32 +45,32 @@ void RadiationBand::addAbsorberGiants(ParameterInput *pin, YAML::Node &node) {
   if (type_ == "radio") {
     if (aname == "NH3") {
       auto ab = std::make_unique<gp::MwrAbsorberNH3>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "H2O") {
       auto ab = std::make_unique<gp::MwrAbsorberH2O>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "H2S") {
       auto ab = std::make_unique<gp::MwrAbsorberH2S>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "PH3") {
       auto ab = std::make_unique<gp::MwrAbsorberPH3>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "CIA") {
       auto ab = std::make_unique<gp::MwrAbsorberCIA>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else if (aname == "Electron") {
       auto ab = std::make_unique<gp::MwrAbsorberElectron>(
-          pmy_block_, species, ToParameterMap(node["parameters"]));
+          species, ToParameterMap(node["parameters"]));
 
       absorbers_.push_back(std::move(ab));
     } else {
diff --git a/src/impl.cpp b/src/impl.cpp
index 7abb8851..49e21b26 100644
--- a/src/impl.cpp
+++ b/src/impl.cpp
@@ -24,10 +24,14 @@
 
 // canoe
 #include "impl.hpp"
+#include "index_map.hpp"
 
 MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
   du.NewAthenaArray(NHYDRO, pmb->ncells3, pmb->ncells2, pmb->ncells1);
 
+  // index map
+  IndexMap::InitFromAthenaInput(pin);
+
   // thermodynamics
   Thermodynamics::InitFromAthenaInput(pin);
 
@@ -64,96 +68,89 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
 #endif  // HYDROSTATIC
 }
 
-MeshBlock::Impl::~Impl() { Thermodynamics::Destroy(); }
+MeshBlock::Impl::~Impl() {
+  Thermodynamics::Destroy();
+  IndexMap::Destroy();
+}
 
-void MeshBlock::Impl::GatherFromPrimitive(Variable *var, int k, int j, int i) const {
-  auto& pmb = pmy_block_;
+void MeshBlock::Impl::GatherFromPrimitive(Variable *var, int k, int j,
+                                          int i) const {
+  auto &pmb = pmy_block_;
   auto mytype = var->GetType();
   var->SetType(Variable::Type::MassFrac);
 
-  for (int n = 0; n < NHYDRO; ++n)
-    var->w[n] = pmb->phydro->w(n, k, j, i);
+  for (int n = 0; n < NHYDRO; ++n) var->w[n] = pmb->phydro->w(n, k, j, i);
 
-  for (int n = 0; n < NCLOUD; ++n)
-    var->c[n] = pcloud->w(n, k, j, i);
+  for (int n = 0; n < NCLOUD; ++n) var->c[n] = pcloud->w(n, k, j, i);
 
-  for (int n = 0; n < NTRACER; ++n)
-    var->x[n] = ptracer->w(n, k, j, i);
+  for (int n = 0; n < NTRACER; ++n) var->x[n] = ptracer->w(n, k, j, i);
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    var->q[n] = pchem->w(n, k, j, i);
+  for (int n = 0; n < NCHEMISTRY; ++n) var->q[n] = pchem->w(n, k, j, i);
 
   var->ConvertTo(mytype);
 }
 
-void MeshBlock::Impl::GatherFromConserved(Variable *var, int k, int j, int i) const {
-  auto& pmb = pmy_block_;
+void MeshBlock::Impl::GatherFromConserved(Variable *var, int k, int j,
+                                          int i) const {
+  auto &pmb = pmy_block_;
   auto mytype = var->GetType();
   var->SetType(Variable::Type::MassConc);
 
-  for (int n = 0; n < NHYDRO; ++n)
-    var->w[n] = pmb->phydro->u(n, k, j, i);
+  for (int n = 0; n < NHYDRO; ++n) var->w[n] = pmb->phydro->u(n, k, j, i);
 
-  for (int n = 0; n < NCLOUD; ++n)
-    var->c[n] = pcloud->u(n, k, j, i);
+  for (int n = 0; n < NCLOUD; ++n) var->c[n] = pcloud->u(n, k, j, i);
 
-  for (int n = 0; n < NTRACER; ++n)
-    var->x[n] = ptracer->u(n, k, j, i);
+  for (int n = 0; n < NTRACER; ++n) var->x[n] = ptracer->u(n, k, j, i);
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    var->q[n] = pchem->u(n, k, j, i);
+  for (int n = 0; n < NCHEMISTRY; ++n) var->q[n] = pchem->u(n, k, j, i);
 
   var->ConvertTo(mytype);
 }
 
-void MeshBlock::Impl::DistributeToPrimitive(Variable const& var_in, int k, int j, int i) {
+void MeshBlock::Impl::DistributeToPrimitive(Variable const &var_in, int k,
+                                            int j, int i) {
   Variable *var;
-  auto& pmb = pmy_block_;
+  auto &pmb = pmy_block_;
 
   if (var_in.GetType() != Variable::Type::MassFrac) {
     var = new Variable(var_in);
     var->ConvertToMassFraction();
   } else {
-    var = const_cast<Variable*>(&var_in);
+    var = const_cast<Variable *>(&var_in);
   }
 
   for (int n = 0; n < NHYDRO; ++n) pmb->phydro->w(n, k, j, i) = var->w[n];
 
-  for (int n = 0; n < NCLOUD; ++n)
-    pcloud->w(n, k, j, i) = var->c[n];
+  for (int n = 0; n < NCLOUD; ++n) pcloud->w(n, k, j, i) = var->c[n];
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    pchem->w(n, k, j, i) = var->q[n];
+  for (int n = 0; n < NCHEMISTRY; ++n) pchem->w(n, k, j, i) = var->q[n];
 
-  for (int n = 0; n < NTRACER; ++n)
-    ptracer->w(n, k, j, i) = var->x[n];
+  for (int n = 0; n < NTRACER; ++n) ptracer->w(n, k, j, i) = var->x[n];
 
   if (var_in.GetType() != Variable::Type::MassFrac) {
     delete var;
   }
 }
 
-void MeshBlock::Impl::DistributeToConserved(Variable const& var_in, int k, int j, int i) {
+void MeshBlock::Impl::DistributeToConserved(Variable const &var_in, int k,
+                                            int j, int i) {
   Variable *var;
-  auto& pmb = pmy_block_;
+  auto &pmb = pmy_block_;
 
   if (var_in.GetType() != Variable::Type::MassConc) {
     var = new Variable(var_in);
     var->ConvertToMassConcentration();
   } else {
-    var = const_cast<Variable*>(&var_in);
+    var = const_cast<Variable *>(&var_in);
   }
 
   for (int n = 0; n < NHYDRO; ++n) pmb->phydro->u(n, k, j, i) = var->w[n];
 
-  for (int n = 0; n < NCLOUD; ++n)
-    pcloud->u(n, k, j, i) = var->c[n];
+  for (int n = 0; n < NCLOUD; ++n) pcloud->u(n, k, j, i) = var->c[n];
 
-  for (int n = 0; n < NCHEMISTRY; ++n)
-    pchem->u(n, k, j, i) = var->q[n];
+  for (int n = 0; n < NCHEMISTRY; ++n) pchem->u(n, k, j, i) = var->q[n];
 
-  for (int n = 0; n < NTRACER; ++n)
-    ptracer->u(n, k, j, i) = var->x[n];
+  for (int n = 0; n < NTRACER; ++n) ptracer->u(n, k, j, i) = var->x[n];
 
   if (var_in.GetType() != Variable::Type::MassConc) {
     delete var;
diff --git a/src/index_map.cpp b/src/index_map.cpp
index 250f47b2..c99327c5 100644
--- a/src/index_map.cpp
+++ b/src/index_map.cpp
@@ -1,3 +1,7 @@
+// C/C++
+#include <mutex>
+#include <string>
+
 // canoe
 #include <configure.hpp>
 
@@ -6,6 +10,7 @@
 #include <athena/parameter_input.hpp>
 
 // application
+#include <application/application.hpp>
 #include <application/exceptions.hpp>
 
 // utils
@@ -14,8 +19,37 @@
 // canoe
 #include "index_map.hpp"
 
-MeshBlock::IndexMap::IndexMap(MeshBlock *pmb, ParameterInput *pin)
-    : pmy_block_(pmb) {
+static std::mutex imap_mutex;
+
+IndexMap::~IndexMap() {
+  Application::Logger app("canoe");
+  app->Log("Destroy IndexMap");
+}
+
+IndexMap const* IndexMap::GetInstance() {
+  // RAII
+  std::unique_lock<std::mutex> lock(imap_mutex);
+
+  if (myindex_map_ == nullptr) {
+    myindex_map_ = new IndexMap();
+  }
+
+  return myindex_map_;
+}
+
+IndexMap const* IndexMap::InitFromAthenaInput(ParameterInput* pin) {
+  Application::Logger app("canoe");
+  app->Log("Initialize IndexMap");
+
+  myindex_map_ = new IndexMap();
+
+  auto& vapor_index_map = myindex_map_->vapor_index_map_;
+  auto& cloud_index_map = myindex_map_->cloud_index_map_;
+  auto& chemistry_index_map = myindex_map_->chemistry_index_map_;
+  auto& tracer_index_map = myindex_map_->tracer_index_map_;
+  auto& static_index_map = myindex_map_->static_index_map_;
+  auto& particle_index_map = myindex_map_->particle_index_map_;
+
   // vapor id
   std::string str = pin->GetOrAddString("species", "vapor", "");
   std::vector<std::string> names = Vectorize<std::string>(str.c_str(), " ,");
@@ -25,7 +59,7 @@ MeshBlock::IndexMap::IndexMap(MeshBlock *pmb, ParameterInput *pin)
   }
 
   for (size_t i = 0; i < names.size(); ++i) {
-    vapor_index_map_[names[i]] = 1 + i;
+    vapor_index_map[names[i]] = 1 + i;
   }
 
   // cloud id
@@ -37,32 +71,32 @@ MeshBlock::IndexMap::IndexMap(MeshBlock *pmb, ParameterInput *pin)
   }
 
   for (size_t i = 0; i < names.size(); ++i) {
-    cloud_index_map_[names[i]] = i;
+    cloud_index_map[names[i]] = i;
   }
 
-  // tracer id
-  str = pin->GetOrAddString("species", "tracer", "");
+  // chemistry id
+  str = pin->GetOrAddString("species", "chemistry", "");
   names = Vectorize<std::string>(str.c_str(), " ,");
 
-  if (names.size() > NTRACER) {
-    throw ValueError("IndexMap", "Number of tracers", NTRACER, names.size());
+  if (names.size() > NCHEMISTRY) {
+    throw ValueError("IndexMap", "Number of chemistry", NCHEMISTRY,
+                     names.size());
   }
 
   for (size_t i = 0; i < names.size(); ++i) {
-    tracer_index_map_[names[i]] = i;
+    chemistry_index_map[names[i]] = i;
   }
 
-  // chemistry id
-  str = pin->GetOrAddString("species", "chemistry", "");
+  // tracer id
+  str = pin->GetOrAddString("species", "tracer", "");
   names = Vectorize<std::string>(str.c_str(), " ,");
 
-  if (names.size() > NCHEMISTRY) {
-    throw ValueError("IndexMap", "Number of chemistry", NCHEMISTRY,
-                     names.size());
+  if (names.size() > NTRACER) {
+    throw ValueError("IndexMap", "Number of tracers", NTRACER, names.size());
   }
 
   for (size_t i = 0; i < names.size(); ++i) {
-    chemistry_index_map_[names[i]] = i;
+    tracer_index_map[names[i]] = i;
   }
 
   // static variable id
@@ -75,18 +109,20 @@ MeshBlock::IndexMap::IndexMap(MeshBlock *pmb, ParameterInput *pin)
   }
 
   for (size_t i = 0; i < names.size(); ++i) {
-    static_index_map_[names[i]] = i;
+    static_index_map[names[i]] = i;
   }
 
   // particle id
   str = pin->GetOrAddString("species", "particle", "");
   names = Vectorize<std::string>(str.c_str(), " ,");
   for (size_t i = 0; i < names.size(); ++i) {
-    particle_index_map_[names[i]] = i;
+    particle_index_map[names[i]] = i;
   }
+
+  return myindex_map_;
 }
 
-size_t MeshBlock::IndexMap::GetSpeciesId(std::string category_name) const {
+size_t IndexMap::GetSpeciesId(std::string category_name) const {
   std::string delimiter = ".";
 
   // Find the position of the delimiter
@@ -113,3 +149,14 @@ size_t MeshBlock::IndexMap::GetSpeciesId(std::string category_name) const {
     throw NotFoundError("GetSpeciesId", "Category " + category);
   }
 }
+
+void IndexMap::Destroy() {
+  std::unique_lock<std::mutex> lock(imap_mutex);
+
+  if (IndexMap::myindex_map_ != nullptr) {
+    delete IndexMap::myindex_map_;
+    IndexMap::myindex_map_ = nullptr;
+  }
+}
+
+IndexMap* IndexMap::myindex_map_ = nullptr;
diff --git a/src/index_map.hpp b/src/index_map.hpp
index d25b2ebc..34742a93 100644
--- a/src/index_map.hpp
+++ b/src/index_map.hpp
@@ -5,15 +5,21 @@
 #include <map>
 #include <string>
 
-// athena
-#include <athena/mesh/mesh.hpp>
-
 class ParameterInput;
 
-class MeshBlock::IndexMap {
+class IndexMap {
+ protected:
+  //! Protected ctor access thru static member function Instance
+  IndexMap() {}
+
  public:
-  IndexMap(MeshBlock *pmb, ParameterInput *pin);
-  ~IndexMap() {}
+  ~IndexMap();
+
+  static IndexMap const *GetInstance();
+
+  static IndexMap const *InitFromAthenaInput(ParameterInput *pin);
+
+  static void Destroy();
 
   size_t GetVaporId(std::string name) const {
     return vapor_index_map_.at(name);
@@ -23,25 +29,26 @@ class MeshBlock::IndexMap {
     return cloud_index_map_.at(name);
   }
 
-  size_t GetTracerId(std::string name) const {
-    return tracer_index_map_.at(name);
-  }
-
   size_t GetChemistryId(std::string name) const {
     return chemistry_index_map_.at(name);
   }
 
+  size_t GetTracerId(std::string name) const {
+    return tracer_index_map_.at(name);
+  }
+
   size_t GetSpeciesId(std::string category_name) const;
 
  private:
-  MeshBlock *pmy_block_;
-
   std::map<std::string, size_t> vapor_index_map_;
-  std::map<std::string, size_t> tracer_index_map_;
   std::map<std::string, size_t> cloud_index_map_;
   std::map<std::string, size_t> chemistry_index_map_;
-  std::map<std::string, size_t> particle_index_map_;
+  std::map<std::string, size_t> tracer_index_map_;
   std::map<std::string, size_t> static_index_map_;
+  std::map<std::string, size_t> particle_index_map_;
+
+  //! Pointer to the single Application instance
+  static IndexMap *myindex_map_;
 };
 
 #endif  // SRC_INDEX_MAP_HPP_
diff --git a/src/mesh_setup.cpp b/src/mesh_setup.cpp
index f1461733..35938ed1 100644
--- a/src/mesh_setup.cpp
+++ b/src/mesh_setup.cpp
@@ -97,7 +97,6 @@ void mesh_setup(ParameterInput*& pinput, Mesh*& pmesh) {
   // set up additional components
   for (int b = 0; b < pmesh->nblocal; ++b) {
     MeshBlock* pmb = pmesh->my_blocks(b);
-    pmb->pindex = std::make_shared<MeshBlock::IndexMap>(pmb, pinput);
     pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pinput);
   }
 
diff --git a/src/opacity/Giants/microwave/mwr_absorber_cia.cpp b/src/opacity/Giants/microwave/mwr_absorber_cia.cpp
index 434e62f3..46bac955 100644
--- a/src/opacity/Giants/microwave/mwr_absorber_cia.cpp
+++ b/src/opacity/Giants/microwave/mwr_absorber_cia.cpp
@@ -1,6 +1,3 @@
-// athena
-#include <athena/mesh/mesh.hpp>
-
 // canoe
 #include <configure.hpp>
 #include <variable.hpp>
@@ -14,9 +11,9 @@
 
 namespace GiantPlanets {
 
-MwrAbsorberCIA::MwrAbsorberCIA(MeshBlock* pmb, std::vector<std::string> species,
+MwrAbsorberCIA::MwrAbsorberCIA(std::vector<std::string> species,
                                ParameterMap params)
-    : Absorber(pmb, "radio-CIA", species, params) {
+    : Absorber("radio-CIA", species, params) {
   if (!params_.count("xHe")) {
     throw NotFoundError("MwrAbsorberCIA", "parameter 'xHe'");
   }
diff --git a/src/opacity/Giants/microwave/mwr_absorber_h2o.cpp b/src/opacity/Giants/microwave/mwr_absorber_h2o.cpp
index f5d6e822..aa0ffaba 100644
--- a/src/opacity/Giants/microwave/mwr_absorber_h2o.cpp
+++ b/src/opacity/Giants/microwave/mwr_absorber_h2o.cpp
@@ -1,6 +1,3 @@
-// athena
-#include <athena/mesh/mesh.hpp>
-
 // canoe
 #include <configure.hpp>
 #include <variable.hpp>
@@ -14,9 +11,9 @@
 
 namespace GiantPlanets {
 
-MwrAbsorberH2O::MwrAbsorberH2O(MeshBlock* pmb, std::vector<std::string> species,
+MwrAbsorberH2O::MwrAbsorberH2O(std::vector<std::string> species,
                                ParameterMap params)
-    : Absorber(pmb, "radio-H2O", species, params) {
+    : Absorber("radio-H2O", species, params) {
   if (!params_.count("xHe")) {
     throw NotFoundError("MwrAbsorberH2O", "parameter 'xHe'");
   }
diff --git a/src/opacity/Giants/microwave/mwr_absorber_nh3.cpp b/src/opacity/Giants/microwave/mwr_absorber_nh3.cpp
index c009039c..a2a951c7 100644
--- a/src/opacity/Giants/microwave/mwr_absorber_nh3.cpp
+++ b/src/opacity/Giants/microwave/mwr_absorber_nh3.cpp
@@ -2,9 +2,6 @@
 #include <configure.hpp>
 #include <variable.hpp>
 
-// athena
-#include <athena/mesh/mesh.hpp>
-
 // climath
 #include <climath/interpolation.h>
 
@@ -17,9 +14,9 @@
 
 namespace GiantPlanets {
 
-MwrAbsorberNH3::MwrAbsorberNH3(MeshBlock* pmb, std::vector<std::string> species,
+MwrAbsorberNH3::MwrAbsorberNH3(std::vector<std::string> species,
                                ParameterMap params)
-    : Absorber(pmb, "radio-NH3", species, params) {
+    : Absorber("radio-NH3", species, params) {
   if (!params_.count("xHe")) {
     throw NotFoundError("MwrAbsorberNH3", "parameter 'xHe'");
   }
diff --git a/src/opacity/Giants/microwave/mwr_absorbers.hpp b/src/opacity/Giants/microwave/mwr_absorbers.hpp
index 714b1c56..9743b02d 100644
--- a/src/opacity/Giants/microwave/mwr_absorbers.hpp
+++ b/src/opacity/Giants/microwave/mwr_absorbers.hpp
@@ -26,8 +26,7 @@ namespace GiantPlanets {
 
 class MwrAbsorberCIA : public Absorber {
  public:
-  MwrAbsorberCIA(MeshBlock* pmb, std::vector<std::string> species,
-                 ParameterMap params);
+  MwrAbsorberCIA(std::vector<std::string> species, ParameterMap params);
 
   Real GetAttenuation(Real wave1, Real wave2, Variable const& var) const;
 
@@ -37,8 +36,7 @@ class MwrAbsorberCIA : public Absorber {
 
 class MwrAbsorberNH3 : public Absorber {
  public:
-  MwrAbsorberNH3(MeshBlock* pmb, std::vector<std::string> species,
-                 ParameterMap params);
+  MwrAbsorberNH3(std::vector<std::string> species, ParameterMap params);
 
   MwrAbsorberNH3& SetModelHanley() {
     model_name_ = "Hanley09";
@@ -66,9 +64,8 @@ class MwrAbsorberNH3 : public Absorber {
 
 class MwrAbsorberPH3 : public Absorber {
  public:
-  MwrAbsorberPH3(MeshBlock* pmb, std::vector<std::string> species,
-                 ParameterMap params)
-      : Absorber(pmb, "radio-PH3", species, params) {}
+  MwrAbsorberPH3(std::vector<std::string> species, ParameterMap params)
+      : Absorber("radio-PH3", species, params) {}
 
   MwrAbsorberPH3& SetModelRadtran() {
     model_name_ = "Radtran";
@@ -85,8 +82,7 @@ class MwrAbsorberPH3 : public Absorber {
 class MwrAbsorberH2O : public Absorber {
  public:
   // TODO(cli) check Karpowics model
-  MwrAbsorberH2O(MeshBlock* pmb, std::vector<std::string> species,
-                 ParameterMap params);
+  MwrAbsorberH2O(std::vector<std::string> species, ParameterMap params);
 
   MwrAbsorberH2O& SetModeldeBoer() {
     model_name_ = "deBoer";
@@ -110,18 +106,16 @@ class MwrAbsorberH2O : public Absorber {
 
 class MwrAbsorberH2S : public Absorber {
  public:
-  MwrAbsorberH2S(MeshBlock* pmb, std::vector<std::string> species,
-                 ParameterMap params)
-      : Absorber(pmb, "radio-H2S", species, params) {}
+  MwrAbsorberH2S(std::vector<std::string> species, ParameterMap params)
+      : Absorber("radio-H2S", species, params) {}
 
   Real GetAttenuation(Real wave1, Real wave2, Variable const& var) const;
 };
 
 class MwrAbsorberElectron : public Absorber {
  public:
-  MwrAbsorberElectron(MeshBlock* pmb, std::vector<std::string> species,
-                      ParameterMap params)
-      : Absorber(pmb, "radio-Electron", species, params) {}
+  MwrAbsorberElectron(std::vector<std::string> species, ParameterMap params)
+      : Absorber("radio-Electron", species, params) {}
 
   MwrAbsorberElectron& SetModelAppletonHartree() {
     model_name_ = "AppletonHartree";
diff --git a/src/snap/thermodynamics/enroll_vapor_functions.cpp b/src/snap/thermodynamics/enroll_vapor_functions.cpp
index d8d6800d..20c6643b 100644
--- a/src/snap/thermodynamics/enroll_vapor_functions.cpp
+++ b/src/snap/thermodynamics/enroll_vapor_functions.cpp
@@ -8,6 +8,7 @@
 
 // canoe
 #include <configure.hpp>
+#include <index_map.hpp>
 #include <variable.hpp>
 
 // snap
@@ -41,11 +42,11 @@ void Thermodynamics::enrollVaporFunctionsGiants() {
   enrollVaporFunctionsJupiterJuno();
 }
 
-// TODO(cli): This is a temporary solution.
-// Make IndexMap a singleton class can solve this problem.
 void Thermodynamics::enrollVaporFunctionsJupiterJuno() {
+  auto pindex = IndexMap::GetInstance();
+
   // water svp:
-  int iH2O = 1;
+  int iH2O = pindex->GetVaporId("H2O");
   for (int n = 0; n < cloud_index_set_[iH2O].size(); ++n) {
     int j = cloud_index_set_[iH2O][n];
     if (n == 0) {
@@ -58,7 +59,7 @@ void Thermodynamics::enrollVaporFunctionsJupiterJuno() {
   }
 
   // ammonia svp:
-  int iNH3 = 2;
+  int iNH3 = pindex->GetVaporId("NH3");
   for (int n = 0; n < cloud_index_set_[iNH3].size(); ++n) {
     int j = cloud_index_set_[iNH3][n];
     if (n == 0) {
diff --git a/src/snap/thermodynamics/thermodynamics.hpp b/src/snap/thermodynamics/thermodynamics.hpp
index 3d22f5c9..d56ee37a 100644
--- a/src/snap/thermodynamics/thermodynamics.hpp
+++ b/src/snap/thermodynamics/thermodynamics.hpp
@@ -70,7 +70,7 @@ class Thermodynamics {
   };
 
   // member functions
-  virtual ~Thermodynamics();
+  ~Thermodynamics();
 
   //! Return a pointer to the one and only instance of Thermodynamics
   static Thermodynamics const *GetInstance();
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3dbef685..37df1d4e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -42,13 +42,13 @@ enable_testing()
 add_test_case(test_weno.cpp)
 add_test_case(test_absorber.cpp)
 add_test_case(test_yaml_read.cpp)
-add_test_case(test_variable.cpp)
 
 if (${NVAPOR} EQUAL 0)
   add_test_case(test_mesh.cpp)
 endif()
 
 if (${NVAPOR} EQUAL 2)
+  add_test_case(test_variable.cpp)
   add_test_case(test_radiation.cpp)
   add_test_case(test_microwave_opacity.cpp)
   add_test_case(test_thermodynamics.cpp)
diff --git a/tests/test_mesh.cpp b/tests/test_mesh.cpp
index c11217eb..a324e9d7 100644
--- a/tests/test_mesh.cpp
+++ b/tests/test_mesh.cpp
@@ -9,7 +9,6 @@
 
 // canoe
 #include <impl.hpp>
-#include <index_map.hpp>
 
 class TestMesh : public testing::Test {
  protected:
@@ -33,7 +32,6 @@ class TestMesh : public testing::Test {
     // set up components
     for (int b = 0; b < pmesh->nblocal; ++b) {
       MeshBlock *pmb = pmesh->my_blocks(b);
-      pmb->pindex = std::make_shared<MeshBlock::IndexMap>(pmb, pinput);
       pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pinput);
     }
 
diff --git a/tests/test_microwave_opacity.cpp b/tests/test_microwave_opacity.cpp
index ef955f97..215e7692 100644
--- a/tests/test_microwave_opacity.cpp
+++ b/tests/test_microwave_opacity.cpp
@@ -9,7 +9,6 @@
 
 // snap
 #include <impl.hpp>
-#include <index_map.hpp>
 #include <variable.hpp>
 
 // harp
@@ -42,7 +41,6 @@ class TestMicrowaveOpacity : public testing::Test {
     // set up components
     for (int b = 0; b < pmesh->nblocal; ++b) {
       MeshBlock *pmb = pmesh->my_blocks(b);
-      pmb->pindex = std::make_shared<MeshBlock::IndexMap>(pmb, pinput);
       pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pinput);
     }
 
diff --git a/tests/test_radiation.cpp b/tests/test_radiation.cpp
index abe05877..e60128f0 100644
--- a/tests/test_radiation.cpp
+++ b/tests/test_radiation.cpp
@@ -38,7 +38,6 @@ class TestRadiation : public testing::Test {
     // set up components
     for (int b = 0; b < pmesh->nblocal; ++b) {
       MeshBlock *pmb = pmesh->my_blocks(b);
-      pmb->pindex = std::make_shared<MeshBlock::IndexMap>(pmb, pinput);
       pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pinput);
     }
 
@@ -55,9 +54,9 @@ class TestRadiation : public testing::Test {
 };
 
 TEST_F(TestRadiation, Species) {
-  MeshBlock *pmb = pmesh->my_blocks(0);
-  EXPECT_EQ(pmb->pindex->GetSpeciesId("vapor.H2O"), 1);
-  EXPECT_EQ(pmb->pindex->GetSpeciesId("vapor.NH3"), 2);
+  auto pindex = IndexMap::GetInstance();
+  EXPECT_EQ(pindex->GetSpeciesId("vapor.H2O"), 1);
+  EXPECT_EQ(pindex->GetSpeciesId("vapor.NH3"), 2);
 }
 
 TEST_F(TestRadiation, Radiation) {
diff --git a/tests/test_thermodynamics.cpp b/tests/test_thermodynamics.cpp
index 88bb713d..17e3b02a 100644
--- a/tests/test_thermodynamics.cpp
+++ b/tests/test_thermodynamics.cpp
@@ -27,6 +27,7 @@ class TestThermodynamics : public testing::Test {
     pinput->LoadFromFile(infile);
     infile.Close();
 
+    IndexMap::InitFromAthenaInput(pinput);
     Thermodynamics::InitFromAthenaInput(pinput);
   }
 
@@ -35,6 +36,7 @@ class TestThermodynamics : public testing::Test {
     // ok to through exceptions from here if need be
 
     Thermodynamics::Destroy();
+    IndexMap::Destroy();
     delete pinput;
   }
 };
diff --git a/tests/test_variable.cpp b/tests/test_variable.cpp
index c8d13d9c..6d2cc378 100644
--- a/tests/test_variable.cpp
+++ b/tests/test_variable.cpp
@@ -30,6 +30,7 @@ class TestVariable : public testing::Test {
     pinput->LoadFromFile(infile);
     infile.Close();
 
+    IndexMap::InitFromAthenaInput(pinput);
     Thermodynamics::InitFromAthenaInput(pinput);
 
     a.SetType(Variable::Type::MoleFrac);
@@ -49,6 +50,7 @@ class TestVariable : public testing::Test {
     // ok to through exceptions from here if need be
     //
     Thermodynamics::Destroy();
+    IndexMap::Destroy();
     delete pinput;
   }
 };