From e2ae7e7c861d4d0e816ea5c3b175e6d737f97d4c Mon Sep 17 00:00:00 2001
From: Cheng Li <69489965+chengcli@users.noreply.github.com>
Date: Tue, 31 Oct 2023 00:47:25 -0400
Subject: [PATCH] Add Jupiter CRM (#103)

add jupiter_crm to test case
---
 CMakeLists.txt                                |   8 +
 ...{sedimentation.cmake => jupiter_crm.cmake} |   6 +-
 cmake/macros/macro_setup_problem.cmake        |  23 ++
 cmake/macros/macro_setup_test.cmake           |  31 +++
 examples/2019-Li-snap/CMakeLists.txt          |  60 +----
 .../{bryan_pgen.cpp => bryan.cpp}             |   3 +
 ...yan_mods.cpp => bryan_vapor_functions.hpp} |  10 +-
 examples/2019-Li-snap/jupiter2d.inp           |  99 +++++++++
 examples/2019-Li-snap/jupiter_crm.cpp         | 209 ++++++++++++++++++
 examples/2019-Li-snap/sedimentation.cpp       |   3 +
 examples/2019-Li-snap/water_ammonia.yaml      |   5 +-
 examples/2023-jupiter-mwr-eq/CMakeLists.txt   |  39 +---
 .../{juno_mwr_pgen.cpp => juno_mwr.cpp}       |   3 +
 .../{juno_mwr_mods.cpp => juno_mwr_specs.hpp} |  89 ++++----
 examples/2024-Li-plume/plume_mods.cpp         |  10 +-
 src/CMakeLists.txt                            |   2 +
 src/microphysics/kessler94.cpp                |   7 +-
 src/microphysics/microphysical_schemes.hpp    |  62 ++++--
 src/microphysics/microphysics.cpp             |  19 +-
 src/microphysics/microphysics.hpp             |  43 +++-
 src/snap/eos/ideal_moist_hydro.cpp            |   2 +-
 src/snap/riemann/lmars.cpp                    |  14 +-
 .../thermodynamics/enroll_vapor_functions.cpp |  24 +-
 src/snap/thermodynamics/thermodynamics.cpp    |   2 +-
 src/snap/thermodynamics/thermodynamics.hpp    |  47 ++--
 src/special/CMakeLists.txt                    |   1 +
 .../special/giants_add_absorber_v1.hpp        |   0
 .../giants_enroll_vapor_functions_v1.hpp      |  86 +++----
 .../special}/jupiter_opacity.yaml             |   0
 {examples => src/special}/mars_absorber.cpp   |   0
 src/tasklist/implicit_hydro_tasks.cpp         |   2 +-
 tests/CMakeLists.txt                          | 106 ++-------
 tests/test_air_parcel.cpp                     |   4 +
 tests/test_ammonium_hydrosulfide.cpp          |   4 +
 tests/test_impl.cpp                           |   4 +
 tests/test_microphysics.cpp                   |   4 +
 tests/test_microwave_opacity.cpp              |   4 +
 tests/test_moist_adiabat.cpp                  |  10 +-
 tests/test_radiation.cpp                      |   4 +
 tests/test_thermodynamics.cpp                 |   4 +
 40 files changed, 685 insertions(+), 368 deletions(-)
 rename cmake/{sedimentation.cmake => jupiter_crm.cmake} (81%)
 create mode 100644 cmake/macros/macro_setup_problem.cmake
 create mode 100644 cmake/macros/macro_setup_test.cmake
 rename examples/2019-Li-snap/{bryan_pgen.cpp => bryan.cpp} (98%)
 rename examples/2019-Li-snap/{bryan_mods.cpp => bryan_vapor_functions.hpp} (69%)
 create mode 100644 examples/2019-Li-snap/jupiter2d.inp
 create mode 100644 examples/2019-Li-snap/jupiter_crm.cpp
 rename examples/2023-jupiter-mwr-eq/{juno_mwr_pgen.cpp => juno_mwr.cpp} (99%)
 rename examples/2023-jupiter-mwr-eq/{juno_mwr_mods.cpp => juno_mwr_specs.hpp} (77%)
 create mode 100644 src/special/CMakeLists.txt
 rename examples/giants_add_absorber.cpp => src/special/giants_add_absorber_v1.hpp (100%)
 rename examples/giants_enroll_vapor_functions.cpp => src/special/giants_enroll_vapor_functions_v1.hpp (52%)
 rename {examples => src/special}/jupiter_opacity.yaml (100%)
 rename {examples => src/special}/mars_absorber.cpp (100%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ceb6b502..8a47281b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -77,6 +77,14 @@ endif()
 
 ## 6. set up unit tests and examples
 message(STATUS "6. Set up unit tests")
+
+# load all macros
+FILE(GLOB _macro_files "${CMAKE_SOURCE_DIR}/cmake/macros/*.cmake")
+FOREACH(_file ${_macro_files})
+  MESSAGE(STATUS "Include ${_file}")
+  INCLUDE(${_file})
+ENDFOREACH()
+
 add_subdirectory(examples)
 add_subdirectory(tests)
 
diff --git a/cmake/sedimentation.cmake b/cmake/jupiter_crm.cmake
similarity index 81%
rename from cmake/sedimentation.cmake
rename to cmake/jupiter_crm.cmake
index 5413a859..10848f4c 100644
--- a/cmake/sedimentation.cmake
+++ b/cmake/jupiter_crm.cmake
@@ -1,4 +1,4 @@
-# configure file for test sedimentation
+# configure file for test jupiter crm
 
 macro(SET_IF_EMPTY _variable)
   if("${${_variable}}" STREQUAL "")
@@ -10,8 +10,8 @@ endmacro()
 set_if_empty(NUMBER_GHOST_CELLS 3)
 
 # canoe configure
-set(NVAPOR 2)
-set(NCLOUD 4)
+set(NVAPOR 1)
+set(NCLOUD 2)
 set(NPHASE_LEGACY 3)
 set(NETCDF ON)
 set(PNETCDF ON)
diff --git a/cmake/macros/macro_setup_problem.cmake b/cmake/macros/macro_setup_problem.cmake
new file mode 100644
index 00000000..571a0b75
--- /dev/null
+++ b/cmake/macros/macro_setup_problem.cmake
@@ -0,0 +1,23 @@
+# A small macro used for setting up the build of a problem.
+#
+# Usage: setup_problem(name)
+
+string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
+string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+
+macro(setup_problem namel)
+  add_executable(${namel}.${buildl} ${namel}.cpp
+                                    ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
+
+  set_target_properties(
+    ${namel}.${buildl}
+    PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+               COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}})
+
+  target_link_libraries(${namel}.${buildl} ${CANOE_LIBRARY_${buildu}})
+
+  target_include_directories(
+    ${namel}.${buildl}
+    PRIVATE ${CMAKE_BINARY_DIR} ${CANOE_INCLUDE_DIR} ${EIGEN3_INCLUDE_DIR}
+            ${MPI_CXX_INCLUDE_PATH} ${NETCDF_INCLUDES} ${PNETCDF_INCLUDE_DIR})
+endmacro()
diff --git a/cmake/macros/macro_setup_test.cmake b/cmake/macros/macro_setup_test.cmake
new file mode 100644
index 00000000..768d4c23
--- /dev/null
+++ b/cmake/macros/macro_setup_test.cmake
@@ -0,0 +1,31 @@
+# A small macro used for setting up the build of a test.
+#
+# Usage: setup_test(name)
+
+macro(setup_test namel)
+  string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
+  string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+
+  add_executable(${namel}.${buildl} ${namel}.cpp globals.cpp)
+
+  set_target_properties(${namel}.${buildl}
+                        PROPERTIES COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}})
+
+  target_include_directories(
+    ${namel}.${buildl}
+    PRIVATE $<$<BOOL:${PVFMM}>:${PVFMM_SOURCE_DIR}/SCTL/include>
+            ${CMAKE_BINARY_DIR}
+            ${CANOE_INCLUDE_DIR}
+            ${EIGEN3_INCLUDE_DIR}
+            ${MPI_CXX_INCLUDE_PATH}
+            ${NETCDF_INCLUDES}
+            ${PNETCDF_INCLUDE_DIR}
+            ${OpenMP_CXX_INCLUDE_DIR}
+            ${FFTW_INCLUDE_DIRS})
+
+  target_link_libraries(
+    ${namel}.${buildl} gtest_main $<$<BOOL:${PVFMM}>:pvfmmStatic>
+    ${CANOE_LIBRARY_${buildu}})
+
+  add_test(NAME ${namel}.${buildl} COMMAND ${namel}.${buildl})
+endmacro()
diff --git a/examples/2019-Li-snap/CMakeLists.txt b/examples/2019-Li-snap/CMakeLists.txt
index 101cb26c..de44b373 100644
--- a/examples/2019-Li-snap/CMakeLists.txt
+++ b/examples/2019-Li-snap/CMakeLists.txt
@@ -1,56 +1,20 @@
-string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
-string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+# ============================================
+# Examples published in Li et al. (2019), ApJS
+# ============================================
 
-macro(setup_problem namel)
-  set_target_properties(
-      ${namel}.${buildl}
-      PROPERTIES
-      RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
-      COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
-      )
+# 1. Compile straka, robert and bryan
+setup_problem(straka)
+setup_problem(robert)
+setup_problem(bryan)
 
-  target_link_libraries(${namel}.${buildl}
-      ${CANOE_LIBRARY_${buildu}}
-      )
-
-  target_include_directories(${namel}.${buildl}
-      PRIVATE
-      ${CMAKE_BINARY_DIR}
-      ${CANOE_INCLUDE_DIR}
-      ${EIGEN3_INCLUDE_DIR}
-      ${MPI_CXX_INCLUDE_PATH}
-      ${NETCDF_INCLUDES}
-      ${PNETCDF_INCLUDE_DIR}
-      )
-endmacro()
-
-# 1. Compile straka and robert
-foreach(namel straka robert)
-  add_executable(${namel}.${buildl} ${namel}.cpp
-    ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
-  setup_problem(${namel})
-endforeach()
-
-# 2. Compile bryan
-if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 1 AND ${NPHASE_LEGACY} EQUAL 2)
-  set(namel bryan)
-
-  add_executable(${namel}.${buildl} ${namel}_pgen.cpp ${namel}_mods.cpp
-    ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
-
-  setup_problem(${namel})
-endif()
-
-# 3. Compile sedimentation
+# 2. Compile sedimentation
 if (${NVAPOR} EQUAL 2 AND ${NCLOUD} EQUAL 4 AND ${NPHASE_LEGACY} EQUAL 3)
-  set(namel sedimentation)
-
-  add_executable(${namel}.${buildl} ${namel}.cpp ../giants_enroll_vapor_functions.cpp
-    ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
-
-  setup_problem(${namel})
+  setup_problem(sedimentation)
 endif()
 
+# 3. Compile jupiter crm
+setup_problem(jupiter_crm)
+
 # 4. Copy input files to run directory
 file(GLOB inputs *.inp *.yaml)
 foreach(input ${inputs})
diff --git a/examples/2019-Li-snap/bryan_pgen.cpp b/examples/2019-Li-snap/bryan.cpp
similarity index 98%
rename from examples/2019-Li-snap/bryan_pgen.cpp
rename to examples/2019-Li-snap/bryan.cpp
index c4fb54ab..5943220d 100644
--- a/examples/2019-Li-snap/bryan_pgen.cpp
+++ b/examples/2019-Li-snap/bryan.cpp
@@ -38,6 +38,9 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
+// special includes
+#include "bryan_vapor_functions.hpp"
+
 int iH2O, iH2Oc;
 Real p0, grav;
 
diff --git a/examples/2019-Li-snap/bryan_mods.cpp b/examples/2019-Li-snap/bryan_vapor_functions.hpp
similarity index 69%
rename from examples/2019-Li-snap/bryan_mods.cpp
rename to examples/2019-Li-snap/bryan_vapor_functions.hpp
index 250eb3a8..cec6c70c 100644
--- a/examples/2019-Li-snap/bryan_mods.cpp
+++ b/examples/2019-Li-snap/bryan_vapor_functions.hpp
@@ -12,17 +12,19 @@ double sat_vapor_p_H2O(double T) {
   return SatVaporPresIdeal(T / tr, pr, betal, gammal);
 }
 
-void Thermodynamics::enrollVaporFunctions() { enrollVaporFunctionH2O(); }
-
 // water svp
-void Thermodynamics::enrollVaporFunctionH2O() {
+void enroll_vapor_Function_H2O(Thermodynamics::SVPFunc1Container &svp_func1) {
   Application::Logger app("snap");
   app->Log("Enrolling H2O vapor pressures");
 
   auto pindex = IndexMap::GetInstance();
   int iH2O = pindex->GetVaporId("H2O");
 
-  svp_func1_[iH2O][0] = [](AirParcel const &qfrac, int, int) {
+  svp_func1[iH2O][0] = [](AirParcel const &qfrac, int, int) {
     return sat_vapor_p_H2O(qfrac.w[IDN]);
   };
 }
+
+void Thermodynamics::enrollVaporFunctions() {
+  enroll_vapor_Function_H2O(svp_func1_);
+}
diff --git a/examples/2019-Li-snap/jupiter2d.inp b/examples/2019-Li-snap/jupiter2d.inp
new file mode 100644
index 00000000..409bf08f
--- /dev/null
+++ b/examples/2019-Li-snap/jupiter2d.inp
@@ -0,0 +1,99 @@
+<comment>
+problem     = test jupiter
+
+<job>
+problem_id  = jupiter2d     # problem ID: basename of output filenames
+
+<output0>
+file_type   = rst
+dt          = 10.E5
+
+<output1>
+file_type   = hst           # History data dump
+dt          = 2.E4          # time increment between outputs
+
+<output2>
+file_type   = netcdf        # Netcdf format output
+variable    = prim          # variables to be output
+dt          = 1.E4          # time increment between outputs
+
+<output3>
+file_type   = netcdf        # Netcdf format output
+variable    = uov           # diagnostic variables
+dt          = 1.E4          # time increment
+
+<time>
+cfl_number  = 0.9
+nlim        = -1            # cycle limit
+tlim        = 2.E5
+xorder      = 5             # horizontal reconstruction order
+integrator  = rk3           # integration method
+
+<mesh>
+nx1         = 80          # Number of zones in X1-direction
+x1min       = -320.E3     # minimum value of X1
+x1max       = 40.E3       # maximum value of X1
+ix1_bc      = reflecting  # Inner-X1 boundary condition flag
+ox1_bc      = reflecting  # Outer-X1 boundary condition flag
+
+nx2         = 320         # Number of zones in X2-direction
+x2min       = 0.          # minimum value of X2
+x2max       = 3200.E3     # maximum value of X2
+ix2_bc      = periodic    # Inner-X2 boundary condition flag
+ox2_bc      = periodic    # Outer-X2 boundary condition flag
+
+nx3         = 1           # Number of zones in X3-direction
+x3min       = -0.5        # minimum value of X3
+x3max       = 0.5         # maximum value of X3
+ix3_bc      = periodic    # Inner-X3 boundary condition flag
+ox3_bc      = periodic    # Outer-X3 boundary condition flag
+
+<meshblock>
+nx1         = 80
+nx2         = 40
+nx3         = 1
+
+<hydro>
+gamma       = 1.42
+grav_acc1   = -24.79
+implicit_flag = 1
+
+#<species>
+#vapor = H2O, NH3
+#cloud = H2O(c), NH3(c), H2O(p), NH3(p)
+
+<species>
+vapor = H2O
+cloud = H2O(c), H2O(p)
+
+<chemistry>
+microphysics_config = water_ammonia.yaml
+
+<thermodynamics>
+Rd          = 3777.   # mu = 2.3175 g/mol
+eps1        = 8.18    8.18    8.18
+beta1       = 0.      24.845  24.845
+rcp1        = 0.15    0.33    0.33
+Ttriple1    = 273.16
+Ptriple1    = 611.7
+
+eps2        = 7.72    7.72    7.72
+rcp2        = 0.078   0.16    0.16
+beta2       = 0.      23.67   23.67
+Ttriple2    = 195.4
+Ptriple2    = 6060.
+
+sa.relax    = 1.0
+sa.max_iter = 4
+sa.ftol     = 1.e-2
+
+<problem>
+qH2O.ppmv   = 25000.
+qNH3.ppmv   = 320.
+
+P0          = 1.E5
+T0          = 169.
+Tmin        = 110.
+hrate       = -100.
+prad        = 2.e5
+init_Ttol   = 2.
diff --git a/examples/2019-Li-snap/jupiter_crm.cpp b/examples/2019-Li-snap/jupiter_crm.cpp
new file mode 100644
index 00000000..71794812
--- /dev/null
+++ b/examples/2019-Li-snap/jupiter_crm.cpp
@@ -0,0 +1,209 @@
+/* -------------------------------------------------------------------------------------
+ * SNAP Example Program
+ *
+ * Contributer:
+ * Cheng Li, University of Michigan
+ *
+ * Year: 2023
+ * Contact: chengcli@umich.edu
+ * Reference: Test Jupiter CRM
+ * -------------------------------------------------------------------------------------
+ */
+
+// athena
+#include <athena/athena.hpp>
+#include <athena/athena_arrays.hpp>
+#include <athena/bvals/bvals.hpp>
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/eos/eos.hpp>
+#include <athena/field/field.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/mesh/mesh.hpp>
+#include <athena/parameter_input.hpp>
+
+// application
+#include <application/application.hpp>
+#include <application/exceptions.hpp>
+
+// canoe
+#include <impl.hpp>
+#include <index_map.hpp>
+
+// climath
+#include <climath/core.h>
+#include <climath/interpolation.h>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// special includes
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
+Real grav, P0, T0, Tmin, prad, hrate;
+int iH2O, iNH3;
+
+void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
+  AllocateUserOutputVariables(4 + NVAPOR);
+  SetUserOutputVariableName(0, "temp");
+  SetUserOutputVariableName(1, "theta");
+  SetUserOutputVariableName(2, "thetav");
+  SetUserOutputVariableName(3, "mse");
+  for (int n = 1; n <= NVAPOR; ++n) {
+    std::string name = "rh" + std::to_string(n);
+    SetUserOutputVariableName(3 + n, name.c_str());
+  }
+}
+
+void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
+  auto pthermo = Thermodynamics::GetInstance();
+
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
+        user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i);
+        // theta_v
+        user_out_var(2, k, j, i) =
+            user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
+        // mse
+        user_out_var(3, k, j, i) =
+            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
+        for (int n = 1; n <= NVAPOR; ++n)
+          user_out_var(3 + n, k, j, i) =
+              pthermo->RelativeHumidity(this, n, k, j, i);
+      }
+}
+
+void Forcing(MeshBlock *pmb, Real const time, Real const dt,
+             AthenaArray<Real> const &w, AthenaArray<Real> const &r,
+             AthenaArray<Real> const &bcc, AthenaArray<Real> &du,
+             AthenaArray<Real> &s) {
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+  auto pthermo = Thermodynamics::GetInstance();
+
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
+
+        Real cv = pthermo->GetCvMass(air, 0);
+
+        if (w(IPR, k, j, i) < prad) {
+          du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv *
+                              (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX));
+        }
+
+        // if (air.w[IDN] < Tmin) {
+        //   u(IEN,k,j,i) += w(IDN,k,j,i)*cv*(Tmin - temp)/sponge_tau*dt;
+        // }
+      }
+}
+
+void Mesh::InitUserMeshData(ParameterInput *pin) {
+  grav = -pin->GetReal("hydro", "grav_acc1");
+
+  P0 = pin->GetReal("problem", "P0");
+  T0 = pin->GetReal("problem", "T0");
+
+  Tmin = pin->GetReal("problem", "Tmin");
+  prad = pin->GetReal("problem", "prad");
+  hrate = pin->GetReal("problem", "hrate") / 86400.;
+
+  // index
+  auto pindex = IndexMap::GetInstance();
+  iH2O = pindex->GetVaporId("H2O");
+  // iNH3 = pindex->GetVaporId("NH3");
+  EnrollUserExplicitSourceFunction(Forcing);
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  srand(Globals::my_rank + time(0));
+
+  Application::Logger app("main");
+  app->Log("ProblemGenerator: jupiter_crm");
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  // mesh limits
+  Real x1min = pmy_mesh->mesh_size.x1min;
+  Real x1max = pmy_mesh->mesh_size.x1max;
+
+  // request temperature and pressure
+  app->Log("request T", T0);
+  app->Log("request P", P0);
+
+  // thermodynamic constants
+  Real gamma = pin->GetReal("hydro", "gamma");
+  Real Rd = pthermo->GetRd();
+  Real cp = gamma / (gamma - 1.) * Rd;
+
+  // set up an adiabatic atmosphere
+  int max_iter = 400, iter = 0;
+  Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01);
+
+  AirParcel air(AirParcel::Type::MoleFrac);
+
+  // estimate surface temperature and pressure
+  Real Ts = T0 - grav / cp * x1min;
+  Real Ps = P0 * pow(Ts / T0, cp / Rd);
+  Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;
+  Real xNH3 = pin->GetReal("problem", "qNH3.ppmv") / 1.E6;
+
+  while (iter++ < max_iter) {
+    // read in vapors
+    air.w[iH2O] = xH2O;
+    // air.w[iNH3] = xNH3;
+    air.w[IPR] = Ps;
+    air.w[IDN] = Ts;
+
+    // stop at just above P0
+    for (int i = is; i <= ie; ++i) {
+      pthermo->Extrapolate(&air, pcoord->dx1f(i),
+                           Thermodynamics::Method::PseudoAdiabat, grav);
+      if (air.w[IPR] < P0) break;
+    }
+
+    // make up for the difference
+    Ts += T0 - air.w[IDN];
+    if (std::abs(T0 - air.w[IDN]) < Ttol) break;
+
+    app->Log("Iteration #", iter);
+    app->Log("T", air.w[IDN]);
+  }
+
+  if (iter > max_iter) {
+    throw RuntimeError("ProblemGenerator", "maximum iteration reached");
+  }
+
+  // construct atmosphere from bottom up
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j) {
+      air.SetZero();
+      air.w[iH2O] = xH2O;
+      // air.w[iNH3] = xNH3;
+      air.w[IPR] = Ps;
+      air.w[IDN] = Ts;
+
+      // half a grid to cell center
+      pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
+                           Thermodynamics::Method::ReversibleAdiabat, grav);
+
+      int i = is;
+      for (; i <= ie; ++i) {
+        if (air.w[IDN] < Tmin) break;
+        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i),
+                             Thermodynamics::Method::PseudoAdiabat, grav,
+                             1.e-5);
+      }
+
+      // Replace adiabatic atmosphere with isothermal atmosphere if temperature
+      // is too low
+      for (; i <= ie; ++i) {
+        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i),
+                             Thermodynamics::Method::Isothermal, grav);
+      }
+    }
+}
diff --git a/examples/2019-Li-snap/sedimentation.cpp b/examples/2019-Li-snap/sedimentation.cpp
index 345e6176..c0e3c697 100644
--- a/examples/2019-Li-snap/sedimentation.cpp
+++ b/examples/2019-Li-snap/sedimentation.cpp
@@ -38,6 +38,9 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
+// special includes
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 Real grav, P0, T0, Tmin;
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
diff --git a/examples/2019-Li-snap/water_ammonia.yaml b/examples/2019-Li-snap/water_ammonia.yaml
index a9a8c1a3..8e6d0308 100644
--- a/examples/2019-Li-snap/water_ammonia.yaml
+++ b/examples/2019-Li-snap/water_ammonia.yaml
@@ -2,14 +2,13 @@ name: giant planet chemistry model
 
 microphysics:
   - water-system
-  - ammonia-system
 
 water-system:
   scheme: Kessler94
   dependent-species: [vapor.H2O, cloud.H2O(c), cloud.H2O(p)]
-  parameters: {autoconversion: 1., accretion: 0.0, evaporation: 1.e-3, sedimentation: -100.}
+  parameters: {autoconversion: 1.e-3, accretion: 0.0, evaporation: 1.e-3, sedimentation: -20.}
 
 ammonia-system:
   scheme: Kessler94
   dependent-species: [vapor.NH3, cloud.NH3(c), cloud.NH3(p)]
-  parameters: {autoconversion: 1., accretion: 0.0, evaporation: 1.e-3, sedimentation: -100.}
+  parameters: {autoconversion: 1.e-3, accretion: 0.0, evaporation: 1.e-3, sedimentation: -20.}
diff --git a/examples/2023-jupiter-mwr-eq/CMakeLists.txt b/examples/2023-jupiter-mwr-eq/CMakeLists.txt
index d782f9df..f415ffc6 100644
--- a/examples/2023-jupiter-mwr-eq/CMakeLists.txt
+++ b/examples/2023-jupiter-mwr-eq/CMakeLists.txt
@@ -1,40 +1,11 @@
-string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
-string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+# =====================================================
+# Juno Microwave Radiometer Forward and Inversion Model
+# =====================================================
 
 # 1. Compile juno_mwr.cpp
+setup_problem(juno_mwr)
 
-set(namel juno_mwr)
-
-add_executable(${namel}.${buildl}
-  ${namel}_pgen.cpp
-  ${namel}_mods.cpp
-  ${CMAKE_BINARY_DIR}/main_${TASK}.cpp
-  )
-
-set_target_properties(
-    ${namel}.${buildl}
-    PROPERTIES
-    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
-    COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
-    )
-
-target_link_libraries(${namel}.${buildl}
-    ${ATHENAPP_LIBRARY_${buildu}}
-    ${CANOE_LIBRARY_${buildu}}
-    ${NETCDF_LIBRARIES}
-    ${CPPDISORT_LIBRARY_${buildu}}
-    )
-
-target_include_directories(${namel}.${buildl}
-    PRIVATE
-    ${CMAKE_BINARY_DIR}
-    ${CANOE_INCLUDE_DIR}
-    ${EIGEN3_INCLUDE_DIR}
-    ${MPI_CXX_INCLUDE_PATH}
-    )
-
-# 1. Copy input file to run directory
-
+# 4. Copy input file to run directory
 file(GLOB inputs *.inp *.yaml)
 foreach(input ${inputs})
     file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin)
diff --git a/examples/2023-jupiter-mwr-eq/juno_mwr_pgen.cpp b/examples/2023-jupiter-mwr-eq/juno_mwr.cpp
similarity index 99%
rename from examples/2023-jupiter-mwr-eq/juno_mwr_pgen.cpp
rename to examples/2023-jupiter-mwr-eq/juno_mwr.cpp
index 9f5a9fc7..b2be09b8 100644
--- a/examples/2023-jupiter-mwr-eq/juno_mwr_pgen.cpp
+++ b/examples/2023-jupiter-mwr-eq/juno_mwr.cpp
@@ -57,6 +57,9 @@
 // inversion
 #include <inversion/profile_inversion.hpp>
 
+// special includes
+#include "juno_mwr_specs.hpp"
+
 Real grav, P0, T0, Tmin, clat;
 Real xHe, xCH4, xH2S, xNa, xKCl, metallicity;
 
diff --git a/examples/2023-jupiter-mwr-eq/juno_mwr_mods.cpp b/examples/2023-jupiter-mwr-eq/juno_mwr_specs.hpp
similarity index 77%
rename from examples/2023-jupiter-mwr-eq/juno_mwr_mods.cpp
rename to examples/2023-jupiter-mwr-eq/juno_mwr_specs.hpp
index a390c2c2..eeea89e1 100644
--- a/examples/2023-jupiter-mwr-eq/juno_mwr_mods.cpp
+++ b/examples/2023-jupiter-mwr-eq/juno_mwr_specs.hpp
@@ -41,20 +41,9 @@
 
 namespace gp = GiantPlanets;
 
-extern Real xHe, xCH4;
-
-void Thermodynamics::enrollVaporFunctions() {
-  Application::Logger app("snap");
-  app->Log("Enrolling vapor functions for juno mwr");
-
-  enrollVaporFunctionH2O();
-  enrollVaporFunctionNH3();
-  enrollVaporFunctionH2S();
-  enrollVaporFunctionNH4SH();
-}
-
 // water svp
-void Thermodynamics::enrollVaporFunctionH2O() {
+void enroll_vapor_function_H2O(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("H2O")) return;
 
@@ -62,20 +51,19 @@ void Thermodynamics::enrollVaporFunctionH2O() {
   app->Log("Enrolling H2O vapor pressures");
 
   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) {
-      svp_func1_[iH2O][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_H2O_BriggsS(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iH2O][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iH2O][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_H2O_BriggsS(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iH2O].size(); ++n) {
+    svp_func1[iH2O][n] = NullSatVaporPres1;
   }
 }
 
 // ammonia svp
-void Thermodynamics::enrollVaporFunctionNH3() {
+void enroll_vapor_function_NH3(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("NH3")) return;
 
@@ -83,20 +71,19 @@ void Thermodynamics::enrollVaporFunctionNH3() {
   app->Log("Enrolling NH3 vapor pressures");
 
   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) {
-      svp_func1_[iNH3][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_NH3_BriggsS(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iNH3][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iNH3][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_NH3_BriggsS(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iNH3].size(); ++n) {
+    svp_func1[iNH3][n] = NullSatVaporPres1;
   }
 }
 
 // hydrogen sulfide svp
-void Thermodynamics::enrollVaporFunctionH2S() {
+void enroll_vapor_function_H2S(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("H2S")) return;
 
@@ -104,20 +91,20 @@ void Thermodynamics::enrollVaporFunctionH2S() {
   app->Log("Enrolling H2S vapor pressures");
 
   int iH2S = pindex->GetVaporId("H2S");
-  for (int n = 0; n < cloud_index_set_[iH2S].size(); ++n) {
-    int j = cloud_index_set_[iH2S][n];
-    if (n == 0) {
-      svp_func1_[iH2S][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_H2S_Antoine(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iH2S][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iH2S][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_H2S_Antoine(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iH2S].size(); ++n) {
+    svp_func1[iH2S][n] = NullSatVaporPres1;
   }
 }
 
 // ammonium hydrosulfide svp
-void Thermodynamics::enrollVaporFunctionNH4SH() {
+void enroll_vapor_function_NH4SH(
+    Thermodynamics::SVPFunc2Container& svp_func2,
+    std::map<IndexPair, ReactionInfo>& cloud_reaction_map) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("NH3") || !pindex->HasVapor("H2S") ||
       !pindex->HasCloud("NH4SH(s)"))
@@ -132,14 +119,26 @@ void Thermodynamics::enrollVaporFunctionNH4SH() {
   int iNH4SH = pindex->GetCloudId("NH4SH(s)");
 
   auto ij = std::minmax(iNH3, iH2S);
-  svp_func2_[ij] = [](AirParcel const& qfrac, int, int, int) {
+  svp_func2[ij] = [](AirParcel const& qfrac, int, int, int) {
     return sat_vapor_p_NH4SH_Lewis(qfrac.w[IDN]);
   };
 
-  cloud_reaction_map_[ij] = std::make_pair<ReactionIndx, ReactionStoi>(
+  cloud_reaction_map[ij] = std::make_pair<ReactionIndx, ReactionStoi>(
       {iNH3, iH2S, iNH4SH}, {1, 1, 1});
 }
 
+void Thermodynamics::enrollVaporFunctions() {
+  Application::Logger app("snap");
+  app->Log("Enrolling vapor functions");
+
+  enroll_vapor_function_H2O(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_NH3(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_H2S(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_NH4SH(svp_func2_, cloud_reaction_map_);
+}
+
+extern Real xHe, xCH4;
+
 // hydrogen heat capacity
 Real Thermodynamics::GetGammad(AirParcel const& qfrac) const {
   Real T = qfrac.w[IDN], cp_h2, cp_he, cp_ch4;
diff --git a/examples/2024-Li-plume/plume_mods.cpp b/examples/2024-Li-plume/plume_mods.cpp
index 84c55396..d346da9f 100644
--- a/examples/2024-Li-plume/plume_mods.cpp
+++ b/examples/2024-Li-plume/plume_mods.cpp
@@ -12,17 +12,19 @@ double sat_vapor_p_H2O_liquid_Ideal(double T) {
   return SatVaporPresIdeal(T / tr, pr, betal, gammal);
 }
 
-void Thermodynamics::enrollVaporFunctions() { enrollVaporFunctionH2O(); }
-
 // water svp
-void Thermodynamics::enrollVaporFunctionH2O() {
+void enroll_vapor_function_H2O(Thermodynamics::SVPFunc1Container &svp_func1) {
   Application::Logger app("snap");
   app->Log("Enrolling H2O vapor pressures");
 
   auto pindex = IndexMap::GetInstance();
   int iH2O = pindex->GetVaporId("H2O");
 
-  svp_func1_[iH2O][0] = [](AirParcel const &qfrac, int, int) {
+  svp_func1[iH2O][0] = [](AirParcel const &qfrac, int, int) {
     return sat_vapor_p_H2O_liquid_Ideal(qfrac.w[IDN]);
   };
 }
+
+void Thermodynamics::enrollVaporFunctions() {
+  enroll_vapor_function_H2O(svp_func1_);
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 81cc335e..f0a6d950 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -51,6 +51,7 @@ add_subdirectory(tasklist)
 add_subdirectory(exo3)
 add_subdirectory(communicator)
 add_subdirectory(nbody)
+add_subdirectory(special)
 
 set(CANOE_LIBRARY_${buildu}
   "canoe_${buildl}"
@@ -79,6 +80,7 @@ set(CANOE_LIBRARY_${buildu}
   ${CFITSIO_LIBRARY}
   ${MPI_CXX_LIBRARIES}
   ${GLOG_LIBRARIES}
+  ${OpenMP_libomp_LIBRARY}
   CACHE INTERNAL "canoe library")
 
 set(CANOE_INCLUDE_DIR
diff --git a/src/microphysics/kessler94.cpp b/src/microphysics/kessler94.cpp
index f18ef5ff..d14a4a81 100644
--- a/src/microphysics/kessler94.cpp
+++ b/src/microphysics/kessler94.cpp
@@ -82,10 +82,9 @@ void Kessler94::EvolveOneStep(AirParcel *air, Real time, Real dt) {
   for (int n = 1; n < Size; ++n) air->w[species_index_[n]] += sol(n);
 }
 
-void Kessler94::SetSedimentationVelocityFromConserved(AthenaArray<Real> vsed[3],
-                                                      Hydro const *phydro,
-                                                      int kl, int ku, int jl,
-                                                      int ju, int il, int iu) {
+void Kessler94::SetVsedFromConserved(AthenaArray<Real> vsed[3],
+                                     Hydro const *phydro, int kl, int ku,
+                                     int jl, int ju, int il, int iu) {
   int ip = species_index_[2] - NHYDRO;
   Real vel = params_.at("sedimentation");
 
diff --git a/src/microphysics/microphysical_schemes.hpp b/src/microphysics/microphysical_schemes.hpp
index d175fd0b..a5f9f57b 100644
--- a/src/microphysics/microphysical_schemes.hpp
+++ b/src/microphysics/microphysical_schemes.hpp
@@ -31,32 +31,46 @@
 class MeshBlock;
 class ParameterInput;
 
+//! \brief virtual base class for all microphysical schemes
 class MicrophysicalSchemeBase : public NamedGroup {
- public:
-  /// constructor and destructor
+ public:  /// constructor and destructor
   MicrophysicalSchemeBase(std::string name) : NamedGroup(name) {}
   virtual ~MicrophysicalSchemeBase() {}
 
-  /// functions
+ public:  /// functions
+  //! \brief Assemble the reaction matrix
+  //!
+  //! \param [in] air air parcel provides data to populate the reaction matrix
+  //! \param [in] time current simulation time
   virtual void AssembleReactionMatrix(AirParcel const &air, Real time) = 0;
-  virtual void EvolveOneStep(AirParcel *Air, Real time, Real dt) = 0;
 
-  /// inbound functions
-  virtual void SetSedimentationVelocityFromConserved(AthenaArray<Real> vsed[3],
-                                                     Hydro const *phydro,
-                                                     int kl, int ku, int jl,
-                                                     int ju, int il,
-                                                     int iu) = 0;
+  //! \brief Evolve the air parcel one time step
+  //!
+  //! \param [in,out] air air parcel to be evolved
+  //! \param [in] time current simulation time
+  //! \param [in] dt time step
+  virtual void EvolveOneStep(AirParcel *air, Real time, Real dt) = 0;
+
+ public:  /// inbound functions
+  //! \brief Set the sedimentation velocity from the conserved variables
+  //!
+  //! \param [in,out] vsed sedimentation velocity
+  //! \param [in] phydro root-level hydrodynamic class
+  virtual void SetVsedFromConserved(AthenaArray<Real> vsed[3],
+                                    Hydro const *phydro, int kl, int ku, int jl,
+                                    int ju, int il, int iu) = 0;
 };
 
 using MicrophysicalSchemePtr = std::shared_ptr<MicrophysicalSchemeBase>;
 using AllMicrophysicalSchemes = std::vector<MicrophysicalSchemePtr>;
 
+//! \brief factory class for contructing microphysical schemes
 class MicrophysicalSchemesFactory {
  public:
   static AllMicrophysicalSchemes Create(MeshBlock *pmb, ParameterInput *pin);
 };
 
+//! \brief base class for all microphysical schemes
 template <int D>
 class MicrophysicalScheme : public MicrophysicalSchemeBase {
  public:
@@ -67,7 +81,7 @@ class MicrophysicalScheme : public MicrophysicalSchemeBase {
   using RealVector = Eigen::Matrix<Real, D, 1>;
   using RealMatrix = Eigen::Matrix<Real, D, D>;
 
-  /// constructor and destructor
+ public:  /// constructor and destructor
   MicrophysicalScheme(std::string name, YAML::Node const &node)
       : MicrophysicalSchemeBase(name) {
     if (node["parameters"]) params_ = ToParameterMap(node["parameters"]);
@@ -85,15 +99,14 @@ class MicrophysicalScheme : public MicrophysicalSchemeBase {
 
   virtual ~MicrophysicalScheme() {}
 
-  /// functions
+ public:  /// functions
   Real const *GetRatePtr() const { return rate_.data(); }
   Real const *GetJacobianPtr() const { return jacb_.data(); }
 
-  /// inbound functions
-  void SetSedimentationVelocityFromConserved(AthenaArray<Real> vsed[3],
-                                             Hydro const *phydro, int kl,
-                                             int ku, int jl, int ju, int il,
-                                             int iu) override {
+ public:  /// inbound functions
+  void SetVsedFromConserved(AthenaArray<Real> vsed[3], Hydro const *phydro,
+                            int kl, int ku, int jl, int ju, int il,
+                            int iu) override {
     for (auto n : species_index_)
       for (int k = kl; k <= ku; ++k)
         for (int j = jl; j <= ju; ++j)
@@ -105,7 +118,7 @@ class MicrophysicalScheme : public MicrophysicalSchemeBase {
   }
 
  protected:
-  //! parameters
+  //! parameters used in microphysics
   ParameterMap params_;
 
   //! indices of species
@@ -116,26 +129,27 @@ class MicrophysicalScheme : public MicrophysicalSchemeBase {
   Eigen::Matrix<Real, D, D> jacb_;
 };
 
+//! \brief Kessler (1994) microphysical scheme
 class Kessler94 : public MicrophysicalScheme<3> {
  public:
   using Base = MicrophysicalScheme<3>;
   enum { Size = Base::Size };
 
-  /// constructor and destructor
+ public:  /// constructor and destructor
   Kessler94(std::string name, YAML::Node const &node);
   ~Kessler94();
 
-  /// functions
+ public:  /// functions
   void AssembleReactionMatrix(AirParcel const &air, Real time) override;
   void EvolveOneStep(AirParcel *air, Real time, Real dt) override;
 
   /// inbound functions
-  void SetSedimentationVelocityFromConserved(AthenaArray<Real> vsed[3],
-                                             Hydro const *phydro, int kl,
-                                             int ku, int jl, int ju, int il,
-                                             int iu) override;
+  void SetVsedFromConserved(AthenaArray<Real> vsed[3], Hydro const *phydro,
+                            int kl, int ku, int jl, int ju, int il,
+                            int iu) override;
 
  protected:
+  //! chemistry solver
   ChemistrySolver<Size> solver_;
 };
 
diff --git a/src/microphysics/microphysics.cpp b/src/microphysics/microphysics.cpp
index fcb52a37..e2706b48 100644
--- a/src/microphysics/microphysics.cpp
+++ b/src/microphysics/microphysics.cpp
@@ -1,3 +1,6 @@
+//! \brief Implementation of the Microphysics class
+//! \file microphysics.cpp
+
 // C/C++
 #include <fstream>
 
@@ -43,10 +46,9 @@ Microphysics::Microphysics(MeshBlock *pmb, ParameterInput *pin)
   vsedf[2].NewAthenaArray(NCLOUD, ncells3 + 1, ncells2, ncells1);
 
   // storage for cloud mass flux at cell boundary
-  vsedf[0].NewAthenaArray(NCLOUD, ncells3, ncells2, ncells1 + 1);
-  mass_flux[0].NewAthenaArray(ncells3, ncells2, ncells1 + 1);
-  mass_flux[1].NewAthenaArray(ncells3, ncells2 + 1, ncells1);
-  mass_flux[2].NewAthenaArray(ncells3 + 1, ncells2, ncells1);
+  mass_flux[0].NewAthenaArray(NCLOUD, ncells3, ncells2, ncells1 + 1);
+  mass_flux[1].NewAthenaArray(NCLOUD, ncells3, ncells2 + 1, ncells1);
+  mass_flux[2].NewAthenaArray(NCLOUD, ncells3 + 1, ncells2, ncells1);
 
   // internal storage for sedimentation velocity at cell center
   vsed_[0].NewAthenaArray(NCLOUD, ncells3, ncells2, ncells1);
@@ -72,23 +74,22 @@ Microphysics::~Microphysics() {
 // void Microphysics::AddFrictionalHeating(
 //     std::vector<AirParcel> &air_column) const {}
 
-void Microphysics::EvolveSystems(AirColumn &air_column, Real time, Real dt) {
+void Microphysics::EvolveSystems(AirColumn &ac, Real time, Real dt) {
   for (auto &system : systems_)
-    for (auto &air : air_column) {
+    for (auto &air : ac) {
       system->AssembleReactionMatrix(air, time);
       system->EvolveOneStep(&air, time, dt);
     }
 }
 
-void Microphysics::SetSedimentationVelocityFromConserved(Hydro const *phydro) {
+void Microphysics::SetVsedFromConserved(Hydro const *phydro) {
   auto pmb = pmy_block_;
 
   int ks = pmb->ks, js = pmb->js, is = pmb->is;
   int ke = pmb->ke, je = pmb->je, ie = pmb->ie;
 
   for (auto &system : systems_) {
-    system->SetSedimentationVelocityFromConserved(vsed_, phydro, ks, ke, js, je,
-                                                  is, ie);
+    system->SetVsedFromConserved(vsed_, phydro, ks, ke, js, je, is, ie);
   }
 
   // interpolation to cell interface
diff --git a/src/microphysics/microphysics.hpp b/src/microphysics/microphysics.hpp
index 421a0a56..9417b05d 100644
--- a/src/microphysics/microphysics.hpp
+++ b/src/microphysics/microphysics.hpp
@@ -19,35 +19,58 @@ class ParameterInput;
 class MicrophysicalSchemeBase;
 class Hydro;
 
+//! \brief root-level management class for microphysics
 class Microphysics {
- public:
-  // tem, v1, v2, v3
-  // enum { NCLOUD_HYDRO = 4 };
+ public:  /// public access members
+  //! \todo(CLI) track cloud temperature and momentum
+  //! tem, v1, v2, v3
+  //! enum { NCLOUD_HYDRO = 4 };
+
+  //! microphysics input key in the input file [microphysics_config]
   static const std::string input_key;
 
-  // access members
-  AthenaArray<Real> w, u, vsedf[3], mass_flux[3];
+  //! primitive variables: mass fraction [kg/kg]
+  AthenaArray<Real> w;
+
+  //! conserved variables: mass concentration [kg/m^3]
+  AthenaArray<Real> u;
+
+  //! sedimentation velocity at cell interface [m/s]
+  AthenaArray<Real> vsedf[3];
 
-  /// constructor and destructor
+  //! mass flux of the dry fluid [kg/m^2/s]
+  AthenaArray<Real> mass_flux[3];
+
+ public:  /// constructor and destructor
   Microphysics(MeshBlock *pmb, ParameterInput *pin);
   ~Microphysics();
 
-  /// functions
+ public:  /// functions
   size_t GetNumSystems() const { return systems_.size(); }
   std::shared_ptr<MicrophysicalSchemeBase> GetSystem(int i) const {
     return systems_[i];
   }
   // void AddFrictionalHeating(std::vector<AirParcel> &air_column) const;
-  void EvolveSystems(AirColumn &air_column, Real time, Real dt);
 
-  /// inbound functions
-  void SetSedimentationVelocityFromConserved(Hydro const *phydro);
+  //! \brief Evolve all microphysical systems
+  //!
+  //! \param [in,out] ac air column to be evolved
+  //! \param [in] time current simulation time
+  //! \param [in] dt time step
+  void EvolveSystems(AirColumn &ac, Real time, Real dt);
+
+ public:  /// inbound functions
+  void SetVsedFromConserved(Hydro const *phydro);
 
  protected:
+  //! sedimentation velocity at cell center [m/s]
   AthenaArray<Real> vsed_[3];
+
+  //! pointers of microphysical systems
   std::vector<std::shared_ptr<MicrophysicalSchemeBase>> systems_;
 
  private:
+  //! meshblock pointer
   MeshBlock const *pmy_block_;
 };
 
diff --git a/src/snap/eos/ideal_moist_hydro.cpp b/src/snap/eos/ideal_moist_hydro.cpp
index 9bb68f4a..f41346fa 100644
--- a/src/snap/eos/ideal_moist_hydro.cpp
+++ b/src/snap/eos/ideal_moist_hydro.cpp
@@ -101,7 +101,7 @@ void EquationOfState::ConservedToPrimitive(
         for (int n = 1; n <= NVAPOR; ++n)
           prim(n, k, j, i) = cons(n, k, j, i) * di;
 
-#ifdef ENABLED_GLOG
+#ifdef ENABLE_GLOG
         LOG_IF(ERROR, std::isnan(w_d) || (w_d < density_floor_))
             << "rank = " << Globals::my_rank << ", (k,j,i) = "
             << "(" << k << "," << j << "," << i << ")"
diff --git a/src/snap/riemann/lmars.cpp b/src/snap/riemann/lmars.cpp
index fb7918ed..3218ca19 100644
--- a/src/snap/riemann/lmars.cpp
+++ b/src/snap/riemann/lmars.cpp
@@ -98,12 +98,18 @@ void Hydro::RiemannSolver(int const k, int const j, int const il, int const iu,
     }
 
     // sedimentation flux
+    // In Athena++ passive tracer module, the tracer mixing ratio is defined as
+    // the ratio of the density of the tracer to the density of the "dry"
+    // species. This ensures conservation of the tracer concentration even with
+    // condensation However, the first component of the primitive variable is
+    // the "total" density To get the denisty of the dry species, the "dry
+    // mixing ratio" (rdl/rdr) is multiplied
     for (int n = 0; n < NCLOUD; ++n) {
-      Real vsed = ubar + pmicro->vsedf[dir](n, k, j, i);
-      if (vsed > 0.) {
-        pmicro->mass_flux[dir](n, k, j, i) = vsed * wli[IDN] * rdl;
+      Real vfld = ubar + pmicro->vsedf[dir](n, k, j, i);
+      if (vfld > 0.) {
+        pmicro->mass_flux[dir](n, k, j, i) = vfld * wli[IDN] * rdl;
       } else {
-        pmicro->mass_flux[dir](n, k, j, i) = vsed * wri[IDN] * rdr;
+        pmicro->mass_flux[dir](n, k, j, i) = vfld * wri[IDN] * rdr;
       }
     }
   }
diff --git a/src/snap/thermodynamics/enroll_vapor_functions.cpp b/src/snap/thermodynamics/enroll_vapor_functions.cpp
index bb507e04..a03658a9 100644
--- a/src/snap/thermodynamics/enroll_vapor_functions.cpp
+++ b/src/snap/thermodynamics/enroll_vapor_functions.cpp
@@ -14,20 +14,10 @@ Real NullSatVaporPres1(AirParcel const& air, int i, int j) {
   return air.w[i] / g * air.w[IPR];
 }
 
-// all vapor functions
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctions() {}
-
-// water svp
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctionH2O() {}
-
-// ammonia svp
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctionNH3() {}
-
-// hydrogen sulfide svp
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctionH2S() {}
-
-// methane svp
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctionCH4() {}
-
-// ammonium hydrosulfide svp
-void __attribute__((weak)) Thermodynamics::enrollVaporFunctionNH4SH() {}
+// enroll null vapor functions
+void __attribute__((weak)) Thermodynamics::enrollVaporFunctions() {
+  for (int i = 0; i <= NVAPOR; ++i)
+    for (int j = 0; j < NPHASE - 1; ++j) {
+      svp_func1_[i][j] = NullSatVaporPres1;
+    }
+}
diff --git a/src/snap/thermodynamics/thermodynamics.cpp b/src/snap/thermodynamics/thermodynamics.cpp
index 8ad398df..b08bd980 100644
--- a/src/snap/thermodynamics/thermodynamics.cpp
+++ b/src/snap/thermodynamics/thermodynamics.cpp
@@ -83,7 +83,7 @@ Thermodynamics* Thermodynamics::fromLegacyInput(ParameterInput* pin) {
   mythermo_->cloud_index_set_.resize(1 + NVAPOR);
   mythermo_->svp_func1_.resize(1 + NVAPOR);
 
-  for (int i = 1; i <= NVAPOR; ++i) {
+  for (int i = 0; i <= NVAPOR; ++i) {
     mythermo_->cloud_index_set_[i].resize(NPHASE - 1);
     mythermo_->svp_func1_[i].resize(NPHASE - 1);
     for (int j = 1; j < NPHASE; ++j) {
diff --git a/src/snap/thermodynamics/thermodynamics.hpp b/src/snap/thermodynamics/thermodynamics.hpp
index cff2bd9e..b16c064e 100644
--- a/src/snap/thermodynamics/thermodynamics.hpp
+++ b/src/snap/thermodynamics/thermodynamics.hpp
@@ -36,6 +36,12 @@ using RealArrayX = std::vector<Real>;
 using SatVaporPresFunc1 = Real (*)(AirParcel const &, int i, int j);
 using SatVaporPresFunc2 = Real (*)(AirParcel const &, int i, int j, int k);
 
+//! \todo(CLI): move to configure.hpp
+enum { MAX_REACTANT = 3 };
+using ReactionIndx = std::array<int, MAX_REACTANT>;
+using ReactionStoi = std::array<int, MAX_REACTANT>;
+using ReactionInfo = std::pair<ReactionIndx, ReactionStoi>;
+
 Real NullSatVaporPres1(AirParcel const &, int, int);
 Real NullSatVaporPres2(AirParcel const &, int, int, int);
 
@@ -66,8 +72,6 @@ class Thermodynamics {
  protected:
   enum { NPHASE = NPHASE_LEGACY };
 
-  enum { MAX_REACTANT = 3 };
-
   //! Constructor for class sets up the initial conditions
   //! Protected ctor access thru static member function Instance
   Thermodynamics() {}
@@ -75,9 +79,8 @@ class Thermodynamics {
   static Thermodynamics *fromYAMLInput(YAML::Node node);
 
  public:
-  using ReactionIndx = std::array<int, MAX_REACTANT>;
-  using ReactionStoi = std::array<int, MAX_REACTANT>;
-  using ReactionInfo = std::pair<ReactionIndx, ReactionStoi>;
+  using SVPFunc1Container = std::vector<std::vector<SatVaporPresFunc1>>;
+  using SVPFunc2Container = std::map<IndexPair, SatVaporPresFunc2>;
 
   enum { Size = 1 + NVAPOR + NCLOUD };
 
@@ -107,34 +110,45 @@ class Thermodynamics {
   //! \return $R_d=\hat{R}/\mu_d$
   Real GetRd() const { return Rd_; }
 
+  //! adiaibatic index of dry air [1]
   Real GetGammad(AirParcel const &var) const;
 
+  //! reference adiabatic index of dry air [1]
   Real GetGammadRef() const { return gammad_ref_; }
 
-  //! Molecular weight
+  //! Molecular weight [kg/mol]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $\mu$
   Real GetMu(int n) const { return mu_[n]; }
 
-  //! Inverse molecular weight
+  //! Inverse molecular weight [mol/kg]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $\mu$
   Real GetInvMu(int n) const { return inv_mu_[n]; }
 
-  //! Ratio of molecular weights
+  //! Ratio of molecular weights [1]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $\epsilon_i=\mu_i/\mu_d$
   Real GetMuRatio(int n) const { return mu_ratio_[n]; }
 
-  //! Ratio of molecular weights
+  //! Ratio of molecular weights [1]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $\epsilon_i^{-1} =\mu_d/\mu_i$
   Real GetInvMuRatio(int n) const { return inv_mu_ratio_[n]; }
 
-  //! Ratio of specific heat capacity [J/(kg K)] at constant volume
+  //! Ratio of specific heat capacity [J/(kg K)] at constant volume [1]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $c_{v,i}/c_{v,d}$
   Real GetCvRatioMass(int n) const { return cv_ratio_mass_[n]; }
 
-  //! Ratio of specific heat capacity [J/(mol K)] at constant volume
+  //! Ratio of specific heat capacity [J/(mol K)] at constant volume [1]
+  //! \param[in] n the index of the thermodynamic species
   //! \return $\hat{c}_{v,i}/\hat{c}_{v,d}$
   Real GetCvRatioMole(int n) const { return cv_ratio_mole_[n]; }
 
+  //! \return the index of the cloud
+  //! \param[in] i the index of the vapor
+  //! \param[in] j the sequential index of the cloud
   int GetCloudIndex(int i, int j) const { return cloud_index_set_[i][j]; }
 
   //! Specific heat capacity [J/(kg K)] at constant volume
@@ -146,6 +160,7 @@ class Thermodynamics {
     return cv_ratio_mass_[n] * cvd;
   }
 
+  //! Reference specific heat capacity [J/(kg K)] at constant volume
   Real GetCvMassRef(int n) const {
     Real cvd = Rd_ / (gammad_ref_ - 1.);
     return cv_ratio_mass_[n] * cvd;
@@ -348,12 +363,6 @@ class Thermodynamics {
 
   // custom functions for vapor (to be overridden by user mods)
   void enrollVaporFunctions();
-  void enrollVaporFunctionH2O();
-  void enrollVaporFunctionNH3();
-  void enrollVaporFunctionH2S();
-  void enrollVaporFunctionCH4();
-  void enrollVaporFunctionCO2();
-  void enrollVaporFunctionNH4SH();
 
  protected:
   //! ideal gas constant of dry air in J/kg
@@ -415,13 +424,13 @@ class Thermodynamics {
   std::array<Real, 1 + NVAPOR> p3_;
 
   //! saturation vapor pressure function: Vapor -> Cloud
-  std::vector<std::vector<SatVaporPresFunc1>> svp_func1_;
+  SVPFunc1Container svp_func1_;
 
   //! cloud index set
   std::vector<IndexSet> cloud_index_set_;
 
   //! saturation vapor pressure function: Vapor + Vapor -> Cloud
-  std::map<IndexPair, SatVaporPresFunc2> svp_func2_;
+  SVPFunc2Container svp_func2_;
 
   //! reaction information map
   std::map<IndexPair, ReactionInfo> cloud_reaction_map_;
diff --git a/src/special/CMakeLists.txt b/src/special/CMakeLists.txt
new file mode 100644
index 00000000..069f8642
--- /dev/null
+++ b/src/special/CMakeLists.txt
@@ -0,0 +1 @@
+# This is a place holder
diff --git a/examples/giants_add_absorber.cpp b/src/special/giants_add_absorber_v1.hpp
similarity index 100%
rename from examples/giants_add_absorber.cpp
rename to src/special/giants_add_absorber_v1.hpp
diff --git a/examples/giants_enroll_vapor_functions.cpp b/src/special/giants_enroll_vapor_functions_v1.hpp
similarity index 52%
rename from examples/giants_enroll_vapor_functions.cpp
rename to src/special/giants_enroll_vapor_functions_v1.hpp
index ab6b5e0f..49336eb4 100644
--- a/examples/giants_enroll_vapor_functions.cpp
+++ b/src/special/giants_enroll_vapor_functions_v1.hpp
@@ -14,18 +14,9 @@
 #include <snap/thermodynamics/vapors/methane_vapors.hpp>
 #include <snap/thermodynamics/vapors/water_vapors.hpp>
 
-void Thermodynamics::enrollVaporFunctions() {
-  Application::Logger app("snap");
-  app->Log("Enrolling vapor functions for juno mwr");
-
-  enrollVaporFunctionH2O();
-  enrollVaporFunctionNH3();
-  enrollVaporFunctionH2S();
-  enrollVaporFunctionNH4SH();
-}
-
 // water svp
-void Thermodynamics::enrollVaporFunctionH2O() {
+void enroll_vapor_function_H2O(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("H2O")) return;
 
@@ -33,20 +24,19 @@ void Thermodynamics::enrollVaporFunctionH2O() {
   app->Log("Enrolling H2O vapor pressures");
 
   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) {
-      svp_func1_[iH2O][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_H2O_BriggsS(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iH2O][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iH2O][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_H2O_BriggsS(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iH2O].size(); ++n) {
+    svp_func1[iH2O][n] = NullSatVaporPres1;
   }
 }
 
 // ammonia svp
-void Thermodynamics::enrollVaporFunctionNH3() {
+void enroll_vapor_function_NH3(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("NH3")) return;
 
@@ -54,20 +44,19 @@ void Thermodynamics::enrollVaporFunctionNH3() {
   app->Log("Enrolling NH3 vapor pressures");
 
   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) {
-      svp_func1_[iNH3][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_NH3_BriggsS(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iNH3][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iNH3][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_NH3_BriggsS(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iNH3].size(); ++n) {
+    svp_func1[iNH3][n] = NullSatVaporPres1;
   }
 }
 
 // hydrogen sulfide svp
-void Thermodynamics::enrollVaporFunctionH2S() {
+void enroll_vapor_function_H2S(Thermodynamics::SVPFunc1Container& svp_func1,
+                               std::vector<IndexSet> const& cloud_index_set) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("H2S")) return;
 
@@ -75,20 +64,20 @@ void Thermodynamics::enrollVaporFunctionH2S() {
   app->Log("Enrolling H2S vapor pressures");
 
   int iH2S = pindex->GetVaporId("H2S");
-  for (int n = 0; n < cloud_index_set_[iH2S].size(); ++n) {
-    int j = cloud_index_set_[iH2S][n];
-    if (n == 0) {
-      svp_func1_[iH2S][n] = [](AirParcel const& qfrac, int, int) {
-        return sat_vapor_p_H2S_Antoine(qfrac.w[IDN]);
-      };
-    } else {
-      svp_func1_[iH2S][n] = NullSatVaporPres1;
-    }
+
+  svp_func1[iH2S][0] = [](AirParcel const& qfrac, int, int) {
+    return sat_vapor_p_H2S_Antoine(qfrac.w[IDN]);
+  };
+
+  for (int n = 1; n < cloud_index_set[iH2S].size(); ++n) {
+    svp_func1[iH2S][n] = NullSatVaporPres1;
   }
 }
 
 // ammonium hydrosulfide svp
-void Thermodynamics::enrollVaporFunctionNH4SH() {
+void enroll_vapor_function_NH4SH(
+    Thermodynamics::SVPFunc2Container& svp_func2,
+    std::map<IndexPair, ReactionInfo>& cloud_reaction_map) {
   auto pindex = IndexMap::GetInstance();
   if (!pindex->HasVapor("NH3") || !pindex->HasVapor("H2S") ||
       !pindex->HasCloud("NH4SH(s)"))
@@ -103,10 +92,21 @@ void Thermodynamics::enrollVaporFunctionNH4SH() {
   int iNH4SH = pindex->GetCloudId("NH4SH(s)");
 
   auto ij = std::minmax(iNH3, iH2S);
-  svp_func2_[ij] = [](AirParcel const& qfrac, int, int, int) {
+
+  svp_func2[ij] = [](AirParcel const& qfrac, int, int, int) {
     return sat_vapor_p_NH4SH_Lewis(qfrac.w[IDN]);
   };
 
-  cloud_reaction_map_[ij] = std::make_pair<ReactionIndx, ReactionStoi>(
+  cloud_reaction_map[ij] = std::make_pair<ReactionIndx, ReactionStoi>(
       {iNH3, iH2S, iNH4SH}, {1, 1, 1});
 }
+
+void Thermodynamics::enrollVaporFunctions() {
+  Application::Logger app("snap");
+  app->Log("Enrolling vapor functions");
+
+  enroll_vapor_function_H2O(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_NH3(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_H2S(svp_func1_, cloud_index_set_);
+  enroll_vapor_function_NH4SH(svp_func2_, cloud_reaction_map_);
+}
diff --git a/examples/jupiter_opacity.yaml b/src/special/jupiter_opacity.yaml
similarity index 100%
rename from examples/jupiter_opacity.yaml
rename to src/special/jupiter_opacity.yaml
diff --git a/examples/mars_absorber.cpp b/src/special/mars_absorber.cpp
similarity index 100%
rename from examples/mars_absorber.cpp
rename to src/special/mars_absorber.cpp
diff --git a/src/tasklist/implicit_hydro_tasks.cpp b/src/tasklist/implicit_hydro_tasks.cpp
index c42ad524..414b00c5 100644
--- a/src/tasklist/implicit_hydro_tasks.cpp
+++ b/src/tasklist/implicit_hydro_tasks.cpp
@@ -296,7 +296,7 @@ TaskStatus ImplicitHydroTasks::ImplicitCorrection(MeshBlock *pmb, int stage) {
 
 TaskStatus ImplicitHydroTasks::UpdateAllConserved(MeshBlock *pmb, int stage) {
   if (stage <= nstages) {
-    pmb->pimpl->pmicro->SetSedimentationVelocityFromConserved(pmb->phydro);
+    pmb->pimpl->pmicro->SetVsedFromConserved(pmb->phydro);
   } else {
     return TaskStatus::fail;
   }
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 73bd7871..7bcdd02e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,114 +1,40 @@
-# ============================
+# ==========================
 # A collection of unit tests
-# ============================
-
-string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
-string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+# ==========================
 
 configure_file(globals.cpp.in globals.cpp @ONLY)
 
-macro(add_test_giants test)
-  get_filename_component(name ${test} NAME_WE)
-
-  add_executable(${name}.${buildl} ${name}.cpp
-	  globals.cpp
-    ${CMAKE_SOURCE_DIR}/examples/giants_add_absorber.cpp
-    ${CMAKE_SOURCE_DIR}/examples/giants_enroll_vapor_functions.cpp
-	  )
-
-  set_target_properties(
-    ${name}.${buildl}
-    PROPERTIES COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
-    )
-
-  target_include_directories(${name}.${buildl}
-    PRIVATE
-    ${CMAKE_BINARY_DIR}
-    ${CANOE_INCLUDE_DIR}
-    ${EIGEN3_INCLUDE_DIR}
-    ${MPI_CXX_INCLUDE_PATH}
-    ${NETCDF_INCLUDES}
-    ${PNETCDF_INCLUDE_DIR}
-    )
-
-  target_link_libraries(${name}.${buildl} gtest_main
-    ${CANOE_LIBRARY_${buildu}}
-    ${NETCDF_LIBRARIES}
-    ${PNETCDF_LIBRARY}
-    ${MPI_CXX_LIBRARIES}
-    ${GLOG_LIBRARIES}
-    )
-
-  add_test(NAME ${name}.${buildl} COMMAND ${name}.${buildl})
-endmacro()
-
-macro(add_test_case test)
-  get_filename_component(name ${test} NAME_WE)
-
-  add_executable(${name}.${buildl} ${name}.cpp globals.cpp)
-
-  set_target_properties(
-    ${name}.${buildl}
-    PROPERTIES COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
-    )
-
-  target_include_directories(${name}.${buildl}
-    PRIVATE
-    $<$<BOOL:${PVFMM}>:${PVFMM_SOURCE_DIR}/SCTL/include>
-    ${CMAKE_BINARY_DIR}
-    ${CANOE_INCLUDE_DIR}
-    ${EIGEN3_INCLUDE_DIR}
-    ${MPI_CXX_INCLUDE_PATH}
-    ${NETCDF_INCLUDES}
-    ${PNETCDF_INCLUDE_DIR}
-    ${OpenMP_CXX_INCLUDE_DIR}
-    ${FFTW_INCLUDE_DIRS}
-    )
-
-  target_link_libraries(${name}.${buildl} gtest_main
-    $<$<BOOL:${PVFMM}>:pvfmmStatic>
-    ${CANOE_LIBRARY_${buildu}}
-    ${NETCDF_LIBRARIES}
-    ${PNETCDF_LIBRARY}
-    #${MPI_CXX_LIBRARIES}
-    ${GLOG_LIBRARIES}
-    ${OpenMP_libomp_LIBRARY}
-    )
-
-  add_test(NAME ${name}.${buildl} COMMAND ${name}.${buildl})
-endmacro()
-
 # Enable testing functionality
 enable_testing()
 
 # Add test cases
-add_test_case(test_weno.cpp)
-add_test_case(test_glog.cpp)
-add_test_giants(test_absorber.cpp)
-add_test_giants(test_yaml_read.cpp)
+setup_test(test_weno)
+setup_test(test_glog)
+setup_test(test_absorber)
+setup_test(test_yaml_read)
 
 if (PVFMM_OPTION STREQUAL "ENABLE_PVFMM")
-  add_test_case(test_nbody.cpp)
+  setup_test(test_nbody)
 endif()
 
 if (${NVAPOR} EQUAL 0)
-  add_test_case(test_mesh.cpp)
+  setup_test(test_mesh)
 endif()
 
 if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 1 AND ${NPHASE_LEGACY} EQUAL 2)
-  add_test_case(test_moist_adiabat.cpp)
+  setup_test(test_moist_adiabat)
 endif()
 
 if (${NVAPOR} EQUAL 2)
   if (${NCLOUD} EQUAL 4)
-    add_test_giants(test_impl.cpp)
-    add_test_giants(test_air_parcel.cpp)
-    add_test_giants(test_radiation.cpp)
-    add_test_giants(test_microwave_opacity.cpp)
-    add_test_giants(test_thermodynamics.cpp)
-    add_test_giants(test_microphysics.cpp)
+    setup_test(test_impl)
+    setup_test(test_air_parcel)
+    setup_test(test_radiation)
+    setup_test(test_microwave_opacity)
+    setup_test(test_thermodynamics)
+    setup_test(test_microphysics)
   elseif (${NCLOUD} EQUAL 5)
-    add_test_giants(test_ammonium_hydrosulfide.cpp)
+    setup_test(test_ammonium_hydrosulfide)
   endif()
 endif()
 
diff --git a/tests/test_air_parcel.cpp b/tests/test_air_parcel.cpp
index 5f562a53..a9d54254 100644
--- a/tests/test_air_parcel.cpp
+++ b/tests/test_air_parcel.cpp
@@ -16,6 +16,10 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestAirParcel : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_ammonium_hydrosulfide.cpp b/tests/test_ammonium_hydrosulfide.cpp
index b489f5e3..6f5b3588 100644
--- a/tests/test_ammonium_hydrosulfide.cpp
+++ b/tests/test_ammonium_hydrosulfide.cpp
@@ -17,6 +17,10 @@
 #include <snap/thermodynamics/vapors/ammonia_vapors.hpp>
 #include <snap/thermodynamics/vapors/water_vapors.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestAmmoniumHydrosulfide : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_impl.cpp b/tests/test_impl.cpp
index 94783438..aa3d530d 100644
--- a/tests/test_impl.cpp
+++ b/tests/test_impl.cpp
@@ -26,6 +26,10 @@
 // tracer
 #include <tracer/tracer.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestImpl : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_microphysics.cpp b/tests/test_microphysics.cpp
index 6d82317d..2e3bc00a 100644
--- a/tests/test_microphysics.cpp
+++ b/tests/test_microphysics.cpp
@@ -19,6 +19,10 @@
 #include <microphysics/microphysical_schemes.hpp>
 #include <microphysics/microphysics.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestMicrophysics : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_microwave_opacity.cpp b/tests/test_microwave_opacity.cpp
index 9bf80607..cce84a7b 100644
--- a/tests/test_microwave_opacity.cpp
+++ b/tests/test_microwave_opacity.cpp
@@ -20,6 +20,10 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestMicrowaveOpacity : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_moist_adiabat.cpp b/tests/test_moist_adiabat.cpp
index 59f213e1..cb422b55 100644
--- a/tests/test_moist_adiabat.cpp
+++ b/tests/test_moist_adiabat.cpp
@@ -20,21 +20,23 @@ double sat_vapor_p_H2O(double T) {
   return SatVaporPresIdeal(T / tr, pr, betal, gammal);
 }
 
-void Thermodynamics::enrollVaporFunctions() { enrollVaporFunctionH2O(); }
-
 // water svp
-void Thermodynamics::enrollVaporFunctionH2O() {
+void enroll_vapor_function_H2O(Thermodynamics::SVPFunc1Container& svp_func1) {
   Application::Logger app("snap");
   app->Log("Enrolling H2O vapor pressures");
 
   auto pindex = IndexMap::GetInstance();
   int iH2O = pindex->GetVaporId("H2O");
 
-  svp_func1_[iH2O][0] = [](AirParcel const& qfrac, int, int) {
+  svp_func1[iH2O][0] = [](AirParcel const& qfrac, int, int) {
     return sat_vapor_p_H2O(qfrac.w[IDN]);
   };
 }
 
+void Thermodynamics::enrollVaporFunctions() {
+  enroll_vapor_function_H2O(svp_func1_);
+}
+
 // expose private members for testing
 class ThermodynamicsTestOnly : public Thermodynamics {
  public:
diff --git a/tests/test_radiation.cpp b/tests/test_radiation.cpp
index 90411065..73934640 100644
--- a/tests/test_radiation.cpp
+++ b/tests/test_radiation.cpp
@@ -20,6 +20,10 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestRadiation : public testing::Test {
  protected:
   ParameterInput *pinput;
diff --git a/tests/test_thermodynamics.cpp b/tests/test_thermodynamics.cpp
index 5a7b8b16..83564a4a 100644
--- a/tests/test_thermodynamics.cpp
+++ b/tests/test_thermodynamics.cpp
@@ -14,6 +14,10 @@
 #include <snap/thermodynamics/vapors/ammonia_vapors.hpp>
 #include <snap/thermodynamics/vapors/water_vapors.hpp>
 
+// special includes
+#include <special/giants_add_absorber_v1.hpp>
+#include <special/giants_enroll_vapor_functions_v1.hpp>
+
 class TestThermodynamics : public testing::Test {
  protected:
   ParameterInput *pinput;