From fcdba9973120a5ae7d0a6c886f6c4479dd1cc398 Mon Sep 17 00:00:00 2001
From: Cheng Li <69489965+chengcli@users.noreply.github.com>
Date: Sat, 9 Mar 2024 17:34:18 -0500
Subject: [PATCH] Add basic diagnostics (#126)

- add diagnostics that does not depend on communication
---
 .github/workflows/main.yml                    |  48 ++++
 examples/2019-Li-snap/straka.inp              |   7 +
 patches/24.time_integrator.patch              |  34 ++-
 src/CMakeLists.txt                            |   3 +-
 src/astro/celestrial_body.cpp                 |  10 +-
 src/astro/celestrial_body.hpp                 |   3 +-
 src/common.hpp                                |   8 -
 src/diagnostics/CMakeLists.txt                |  19 ++
 src/{snap => }/diagnostics/buoyancy.cpp       |  13 +-
 src/{snap => }/diagnostics/curl.cpp           |   9 +-
 src/diagnostics/diagnostics.cpp               |  46 ++++
 src/{snap => }/diagnostics/diagnostics.hpp    | 129 ++++++-----
 src/diagnostics/diagnostics_factory.cpp       |  54 +++++
 src/{snap => }/diagnostics/divergence.cpp     |  19 +-
 src/diagnostics/hydro_mean.cpp                |  85 +++++++
 .../todo}/convective_heat_flx.cpp             |  39 ++--
 .../todo}/eddy_flux.cpp                       |  14 +-
 .../todo}/eddy_kinetic_energy.cpp             |  12 +-
 .../todo}/horizontal_divergence.cpp           |   7 +-
 .../todo}/latent_heating.cpp                  |  17 +-
 .../todo}/pressure_anomaly.cpp                |   9 +-
 .../todo}/radiative_flux.cpp                  |  14 +-
 .../todo/spherical_angular_momentum.cpp}      |  23 +-
 src/diagnostics/todo/temperature_anomaly.cpp  |  68 ++++++
 .../todo}/tendency.cpp                        |  15 +-
 .../todo}/total_flux.cpp                      |  13 +-
 src/exchanger/exchanger.cpp                   |   3 -
 src/exchanger/linear_exchanger.hpp            |   3 -
 src/harp/radiation.cpp                        |   7 +-
 src/harp/radiation.hpp                        |   3 -
 src/harp/radiation_band.hpp                   |   1 -
 src/harp/spectral_grid.hpp                    |   8 +
 src/impl.cpp                                  |   8 +-
 src/impl.hpp                                  |   2 +
 src/outputs/load_user_output_data.cpp         |  17 ++
 src/outputs/netcdf.cpp                        |  57 +++--
 src/outputs/output_utils.cpp                  |  57 +++--
 src/outputs/pnetcdf.cpp                       |  55 +++--
 src/snap/diagnostics/diagnostics.cpp          | 214 ------------------
 src/snap/diagnostics/hydro_mean.cpp           |  75 ------
 src/snap/diagnostics/temperature_anomaly.cpp  |  40 ----
 src/virtual_groups.hpp                        |   8 +
 42 files changed, 715 insertions(+), 561 deletions(-)
 delete mode 100644 src/common.hpp
 create mode 100644 src/diagnostics/CMakeLists.txt
 rename src/{snap => }/diagnostics/buoyancy.cpp (88%)
 rename src/{snap => }/diagnostics/curl.cpp (97%)
 create mode 100644 src/diagnostics/diagnostics.cpp
 rename src/{snap => }/diagnostics/diagnostics.hpp (67%)
 create mode 100644 src/diagnostics/diagnostics_factory.cpp
 rename src/{snap => }/diagnostics/divergence.cpp (92%)
 create mode 100644 src/diagnostics/hydro_mean.cpp
 rename src/{snap/diagnostics => diagnostics/todo}/convective_heat_flx.cpp (82%)
 rename src/{snap/diagnostics => diagnostics/todo}/eddy_flux.cpp (92%)
 rename src/{snap/diagnostics => diagnostics/todo}/eddy_kinetic_energy.cpp (91%)
 rename src/{snap/diagnostics => diagnostics/todo}/horizontal_divergence.cpp (95%)
 rename src/{snap/diagnostics => diagnostics/todo}/latent_heating.cpp (91%)
 rename src/{snap/diagnostics => diagnostics/todo}/pressure_anomaly.cpp (83%)
 rename src/{snap/diagnostics => diagnostics/todo}/radiative_flux.cpp (88%)
 rename src/{snap/diagnostics/angular_momentum.cpp => diagnostics/todo/spherical_angular_momentum.cpp} (80%)
 create mode 100644 src/diagnostics/todo/temperature_anomaly.cpp
 rename src/{snap/diagnostics => diagnostics/todo}/tendency.cpp (94%)
 rename src/{snap/diagnostics => diagnostics/todo}/total_flux.cpp (91%)
 delete mode 100644 src/snap/diagnostics/diagnostics.cpp
 delete mode 100644 src/snap/diagnostics/hydro_mean.cpp
 delete mode 100644 src/snap/diagnostics/temperature_anomaly.cpp

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index f5442479..5618ec42 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -34,6 +34,54 @@ jobs:
       - name: Run cpplint
         run: ./tools/check_style.sh
 
+  build-cubed-sphere-2d:
+    if: github.event.pull_request.draft == false
+    runs-on: ubuntu-latest
+    needs: style-checker
+
+    steps:
+      - name: install dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y lcov libnetcdf-dev libpnetcdf-dev libeigen3-dev libmpich-dev
+
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: create build directory
+        # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+        # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        run: cmake -B ${{github.workspace}}/build -DTASK=exo2 -DCMAKE_CXX_FLAGS=${{ env.BUILD_FLAG }} -DCMAKE_C_FLAGS=${{ env.BUILD_FLAG }}
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config Debug -- -j2
+
+  build-cubed-sphere-3d:
+    if: github.event.pull_request.draft == false
+    runs-on: ubuntu-latest
+    needs: style-checker
+
+    steps:
+      - name: install dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y lcov libnetcdf-dev libpnetcdf-dev libeigen3-dev libmpich-dev
+
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: create build directory
+        # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+        # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        run: cmake -B ${{github.workspace}}/build -DTASK=exo3 -DCMAKE_CXX_FLAGS=${{ env.BUILD_FLAG }} -DCMAKE_C_FLAGS=${{ env.BUILD_FLAG }}
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config Debug -- -j2
+
   build-non-hydrostatic:
     # The CMake configure and build commands are platform agnostic and should work equally well on Windows or Mac.
     # You can convert this to a matrix build if you need cross-platform coverage.
diff --git a/examples/2019-Li-snap/straka.inp b/examples/2019-Li-snap/straka.inp
index 330371d7..c908b6f8 100644
--- a/examples/2019-Li-snap/straka.inp
+++ b/examples/2019-Li-snap/straka.inp
@@ -21,6 +21,11 @@ file_type  = netcdf
 variable   = uov
 dt         = 300.
 
+<output4>
+file_type  = netcdf
+variable   = diag
+dt         = 300.
+
 <time>
 cfl_number = 0.9          # The Courant, Friedrichs, & Lewy (CFL) Number
 nlim       = -1           # cycle limit
@@ -64,3 +69,5 @@ zr    = 2.E3
 p0    = 1.E5
 Ts    = 300.
 K     = 75.
+
+diagnostics = div,curl
diff --git a/patches/24.time_integrator.patch b/patches/24.time_integrator.patch
index 4af5c1df..329f9851 100644
--- a/patches/24.time_integrator.patch
+++ b/patches/24.time_integrator.patch
@@ -1,19 +1,22 @@
 diff --git a/src/task_list/time_integrator.cpp b/src/task_list/time_integrator.cpp
-index 097afa53..225ab61f 100644
+index 2bfc5178..e0fca29d 100644
 --- a/src/task_list/time_integrator.cpp
 +++ b/src/task_list/time_integrator.cpp
-@@ -36,6 +36,10 @@
+@@ -36,6 +36,13 @@
  #include "../scalars/scalars.hpp"
  #include "task_list.hpp"
  
 +// canoe injection
 +#include <configure.hpp>
 +#include <impl.hpp>
++
++// diagnostics
++#include <diagnostics/diagnostics.hpp>
 +
  //----------------------------------------------------------------------------------------
  //! TimeIntegratorTaskList constructor
  
-@@ -998,6 +1002,14 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -998,6 +1005,14 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
        AddTask(SEND_HYD,src_aterm);
        AddTask(RECV_HYD,NONE);
        AddTask(SETB_HYD,(RECV_HYD|SRC_TERM));
@@ -28,7 +31,7 @@ index 097afa53..225ab61f 100644
      }
  
      if (SHEAR_PERIODIC) {
-@@ -1156,9 +1168,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -1157,9 +1172,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
            }
          } else {
            if (NSCALARS > 0) {
@@ -46,7 +49,7 @@ index 097afa53..225ab61f 100644
            }
          }
        }
-@@ -1224,11 +1244,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1225,11 +1248,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
    //!   SRC_TERM and SourceTerms(), USERWORK, PHY_BVAL, PROLONG, CONS2PRIM,
    //!   ... Although, AMR_FLAG = "flag blocks for AMR" should be FLAG_AMR in VERB_OBJECT
    using namespace HydroIntegratorTaskNames; // NOLINT (build/namespace)
@@ -64,7 +67,7 @@ index 097afa53..225ab61f 100644
    } else if (id == CALC_HYDFLX) {
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
-@@ -1274,7 +1299,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1275,7 +1303,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::AddSourceTerms);
      task_list_[ntasks].lb_time = true;
@@ -73,7 +76,7 @@ index 097afa53..225ab61f 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendHydro);
-@@ -1284,17 +1309,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1285,17 +1313,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendField);
      task_list_[ntasks].lb_time = true;
@@ -93,7 +96,7 @@ index 097afa53..225ab61f 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SetBoundariesHydro);
-@@ -1683,6 +1707,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
+@@ -1688,6 +1715,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
    return TaskStatus::success;
  }
  
@@ -113,7 +116,20 @@ index 097afa53..225ab61f 100644
  //----------------------------------------------------------------------------------------
  // Functions to calculates Hydro fluxes
  
-@@ -2388,6 +2425,7 @@ TaskStatus TimeIntegratorTaskList::IntegrateScalars(MeshBlock *pmb, int stage) {
+@@ -2305,6 +2345,12 @@ TaskStatus TimeIntegratorTaskList::PhysicalBoundary(MeshBlock *pmb, int stage) {
+ TaskStatus TimeIntegratorTaskList::UserWork(MeshBlock *pmb, int stage) {
+   if (stage != nstages) return TaskStatus::success; // only do on last stage
+ 
++  auto all_diags = pmb->pimpl->all_diags;
++
++  for (auto& diag : all_diags) {
++    diag->Progress(pmb->phydro->w);
++  }
++
+   pmb->UserWorkInLoop();
+   return TaskStatus::success;
+ }
+@@ -2395,6 +2441,7 @@ TaskStatus TimeIntegratorTaskList::IntegrateScalars(MeshBlock *pmb, int stage) {
        ave_wghts[2] = stage_wghts[stage-1].gamma_3;
        if (ave_wghts[0] == 0.0 && ave_wghts[1] == 1.0 && ave_wghts[2] == 0.0) {
          ps->s.SwapAthenaArray(ps->s1);
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a9df25d4..5f76b979 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -53,6 +53,7 @@ add_subdirectory(exchanger)
 add_subdirectory(nbody)
 add_subdirectory(special)
 add_subdirectory(single_column)
+add_subdirectory(diagnostics)
 
 set(CANOE_LIBRARY_${buildu}
   "canoe_${buildl}"
@@ -74,8 +75,8 @@ set(CANOE_LIBRARY_${buildu}
   "exchanger_${buildl}"
   "nbody_${buildl}"
   "scm_${buildl}"
+  "diagnostics_${buildl}"
   "yaml-cpp"
-  #${CDISORT_LIBRARY_${buildu}}
   ${CPPDISORT_LIBRARY_${buildu}}
   ${PYTHON_LIBRARY_RELEASE}
   ${NETCDF_LIBRARIES}
diff --git a/src/astro/celestrial_body.cpp b/src/astro/celestrial_body.cpp
index 68101fe9..35b2a588 100644
--- a/src/astro/celestrial_body.cpp
+++ b/src/astro/celestrial_body.cpp
@@ -8,18 +8,19 @@
 
 // athena
 #include <athena/athena.hpp>
+#include <athena/mesh/mesh.hpp>
 #include <athena/parameter_input.hpp>
 
 // climath
 #include <climath/core.h>
 #include <climath/interpolation.h>
 
+// harp
+#include <harp/spectral_grid.hpp>
+
 // utils
 #include <utils/fileio.hpp>
 
-// canoe
-#include <common.hpp>
-
 // astro
 #include "celestrial_body.hpp"
 
@@ -161,7 +162,8 @@ Real CelestrialBody::ParentDistanceInAu(Real time) const {
                                             equinox - M_PI / 2.)));
 }
 
-CelestrialBodyPtr PlanetFactory::Create(ParameterInput *pin) {
+CelestrialBodyPtr PlanetFactory::CreateFrom(MeshBlock *pmb,
+                                            ParameterInput *pin) {
   if (pin->DoesParameterExist("astronomy", "planet")) {
     std::string name = pin->GetString("astronomy", "planet");
     return std::make_shared<CelestrialBody>(pin, name);
diff --git a/src/astro/celestrial_body.hpp b/src/astro/celestrial_body.hpp
index 5ee425ca..1c787d1e 100644
--- a/src/astro/celestrial_body.hpp
+++ b/src/astro/celestrial_body.hpp
@@ -14,6 +14,7 @@
 class ParameterInput;
 struct float_triplet;
 struct Direction;
+class Meshblock;
 
 class CelestrialBody : public NamedGroup {
  public:
@@ -67,7 +68,7 @@ Real GetGravity(char const *name, Real pclat);
 
 class PlanetFactory {
  public:
-  static CelestrialBodyPtr Create(ParameterInput *pin);
+  static CelestrialBodyPtr CreateFrom(MeshBlock *pmb, ParameterInput *pin);
 };
 
 #endif  // SRC_ASTRO_CELESTRIAL_BODY_HPP_
diff --git a/src/common.hpp b/src/common.hpp
deleted file mode 100644
index bad6f15a..00000000
--- a/src/common.hpp
+++ /dev/null
@@ -1,8 +0,0 @@
-#ifndef COMMON_HPP_
-#define COMMON_HPP_
-
-struct Direction {
-  Real mu, phi;
-};
-
-#endif  // COMMON_HPP_
diff --git a/src/diagnostics/CMakeLists.txt b/src/diagnostics/CMakeLists.txt
new file mode 100644
index 00000000..abdb79f2
--- /dev/null
+++ b/src/diagnostics/CMakeLists.txt
@@ -0,0 +1,19 @@
+set(namel diagnostics)
+string(TOUPPER ${namel} nameu)
+
+file(GLOB src_files
+    *.cpp
+    )
+
+string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
+string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+
+add_library(${namel}_${buildl}
+    OBJECT
+    ${src_files}
+    )
+
+set_target_properties(${namel}_${buildl}
+    PROPERTIES
+    COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
+    )
diff --git a/src/snap/diagnostics/buoyancy.cpp b/src/diagnostics/buoyancy.cpp
similarity index 88%
rename from src/snap/diagnostics/buoyancy.cpp
rename to src/diagnostics/buoyancy.cpp
index e9eaaa68..9d08854a 100644
--- a/src/snap/diagnostics/buoyancy.cpp
+++ b/src/diagnostics/buoyancy.cpp
@@ -1,13 +1,14 @@
-#include "../coordinates/coordinates.hpp"
-#include "../hydro/hydro.hpp"
-#include "../hydro/srcterms/hydro_srcterms.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/hydro/srcterms/hydro_srcterms.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// diagnostics
 #include "diagnostics.hpp"
 
 Buoyancy::Buoyancy(MeshBlock* pmb) : Diagnostics(pmb, "b") {
   type = "SCALARS";
-  long_name = "buoyancy";
-  units = "m/s^2";
   data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
   pf_.NewAthenaArray(ncells3_, ncells2_, ncells1_ + 1);
   grav_ = pmb->phydro->hsrc.GetG1();
diff --git a/src/snap/diagnostics/curl.cpp b/src/diagnostics/curl.cpp
similarity index 97%
rename from src/snap/diagnostics/curl.cpp
rename to src/diagnostics/curl.cpp
index 0b6cdedd..4f9ee75a 100644
--- a/src/snap/diagnostics/curl.cpp
+++ b/src/diagnostics/curl.cpp
@@ -1,10 +1,11 @@
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 Curl::Curl(MeshBlock *pmb) : Diagnostics(pmb, "curl") {
-  long_name = "curl";
-  units = "1/s";
   if (pmb->block_size.nx3 > 1) {
     type = "VECTORS";
     data.NewAthenaArray(3, ncells3_, ncells2_, ncells1_);
diff --git a/src/diagnostics/diagnostics.cpp b/src/diagnostics/diagnostics.cpp
new file mode 100644
index 00000000..e88127ab
--- /dev/null
+++ b/src/diagnostics/diagnostics.cpp
@@ -0,0 +1,46 @@
+// C/C++ headers
+#include <cstring>
+#include <sstream>
+#include <stdexcept>
+
+// application
+#include <application/application.hpp>
+
+// Athena++ headers
+#include <athena/mesh/mesh.hpp>
+
+// diagnostics
+#include "diagnostics.hpp"
+
+Diagnostics::Diagnostics(MeshBlock *pmb, std::string name)
+    : NamedGroup(name), pmy_block_(pmb), ncycle_(0) {
+  Application::Logger app("main");
+  app->Log("Initialize Diagnostics");
+
+  ncells1_ = pmb->block_size.nx1 + 2 * (NGHOST);
+  ncells2_ = 1;
+  ncells3_ = 1;
+  if (pmb->pmy_mesh->f2) ncells2_ = pmb->block_size.nx2 + 2 * (NGHOST);
+  if (pmb->pmy_mesh->f3) ncells3_ = pmb->block_size.nx3 + 2 * (NGHOST);
+
+  x1edge_.NewAthenaArray(ncells1_ + 1);
+  x1edge_p1_.NewAthenaArray(ncells1_);
+  x2edge_.NewAthenaArray(ncells1_ + 1);
+  x2edge_p1_.NewAthenaArray(ncells1_);
+  x3edge_.NewAthenaArray(ncells1_ + 1);
+  x3edge_p1_.NewAthenaArray(ncells1_);
+
+  x1area_.NewAthenaArray(ncells1_ + 1);
+  x2area_.NewAthenaArray(ncells1_);
+  x2area_p1_.NewAthenaArray(ncells1_);
+  x3area_.NewAthenaArray(ncells1_);
+  x3area_p1_.NewAthenaArray(ncells1_);
+
+  vol_.NewAthenaArray(ncells1_);
+  total_vol_.NewAthenaArray(ncells1_);
+}
+
+Diagnostics::~Diagnostics() {
+  Application::Logger app("main");
+  app->Log("Destroy Diagnostics");
+}
diff --git a/src/snap/diagnostics/diagnostics.hpp b/src/diagnostics/diagnostics.hpp
similarity index 67%
rename from src/snap/diagnostics/diagnostics.hpp
rename to src/diagnostics/diagnostics.hpp
index c27123b2..5919afa4 100644
--- a/src/snap/diagnostics/diagnostics.hpp
+++ b/src/diagnostics/diagnostics.hpp
@@ -1,75 +1,81 @@
-#ifndef DIAGNOSTICS_HPP_
-#define DIAGNOSTICS_HPP_
+#ifndef SRC_DIAGNOSTICS_HPP_
+#define SRC_DIAGNOSTICS_HPP_
 
-// C++ header
+// C/C++
 #include <string>
 
-// Athena++ headers
-#include "../athena_arrays.hpp"
-#include "../mesh/mesh.hpp"
+// Athena++
+#include <athena/athena_arrays.hpp>
+#include <athena/mesh/mesh.hpp>
+
+// canoe
+#include <configure.hpp>
+#include <virtual_groups.hpp>
+
+// exchanger
+#include <exchanger/exchanger.hpp>
 
 class ParameterInput;
 
-class Diagnostics {
+class Diagnostics : public NamedGroup {
  public:
   // data
-  std::string myname, type, grid;
-  std::string varname, long_name, units;
-  Diagnostics *prev, *next;
+  std::string type;
   AthenaArray<Real> data;
-  int ncycle;
 
   // functions
-  Diagnostics(MeshBlock *pmb, ParameterInput *pin);
   Diagnostics(MeshBlock *pmb, std::string name);
   virtual ~Diagnostics();
 
-  Diagnostics *operator[](std::string name);
-
-  template <typename Dg>
-  Diagnostics *AddDiagnostics(Dg const &d) {
-    Dg *pd = new Dg(d);
-    Diagnostics *p = this;
-    while (p->next != nullptr) p = p->next;
-    p->next = pd;
-    p->next->prev = p;
-    p->next->next = nullptr;
-    return p->next;
-  }
-
+  virtual int GetNumVars() const = 0;
   virtual void Progress(AthenaArray<Real> const &w) {}
   virtual void Finalize(AthenaArray<Real> const &w) {}
 
  protected:
   MeshBlock *pmy_block_;
+
   int ncells1_, ncells2_, ncells3_;
 
+  int ncycle_;
+
   //! mean and eddy component
   AthenaArray<Real> mean_, eddy_;
 
   //! MPI color of each rank
   std::vector<int> color_;
+
   //! rank of the bottom block
   std::vector<int> brank_;
 
   //! scratch geometric arrays
-  AthenaArray<Real> x1area_, x2area_, x2area_p1_, x3area_, x3area_p1_, vol_,
-      total_vol_;
   AthenaArray<Real> x1edge_, x1edge_p1_, x2edge_, x2edge_p1_, x3edge_,
       x3edge_p1_;
 
-  void setColor_(int *color, CoordinateDirection dir);
-  void gatherAllData23_(AthenaArray<Real> &total_vol,
+  AthenaArray<Real> x1area_, x2area_, x2area_p1_, x3area_, x3area_p1_;
+
+  AthenaArray<Real> vol_, total_vol_;
+
+  void GatherVolumnData(AthenaArray<Real> &total_vol,
                         AthenaArray<Real> &total_data);
 };
 
+using DiagnosticsPtr = std::shared_ptr<Diagnostics>;
+using DiagnosticsContainer = std::vector<DiagnosticsPtr>;
+
+class DiagnosticsFactory {
+ public:
+  static DiagnosticsContainer CreateFrom(MeshBlock *pmb, ParameterInput *pin);
+};
+
 // register all diagnostics
 // 1. divergence
 class Divergence : public Diagnostics {
  public:
   Divergence(MeshBlock *pmb);
   virtual ~Divergence() {}
-  void Finalize(AthenaArray<Real> const &w);
+
+  void Finalize(AthenaArray<Real> const &w) override;
+  int GetNumVars() const override { return 1; }
 
  protected:
   AthenaArray<Real> v1f1_, v2f2_, v3f3_;
@@ -80,25 +86,40 @@ class Curl : public Diagnostics {
  public:
   Curl(MeshBlock *pmb);
   virtual ~Curl() {}
-  void Finalize(AthenaArray<Real> const &w);
+
+  void Finalize(AthenaArray<Real> const &w) override;
+  int GetNumVars() const override { return 1; }
 
  protected:
   AthenaArray<Real> v3f2_, v2f3_, v1f3_, v3f1_, v2f1_, v1f2_;
 };
 
-// 3. hydro mean
+// 3. Buoyancy
+class Buoyancy : public Diagnostics {
+ public:
+  Buoyancy(MeshBlock *pmb);
+  virtual ~Buoyancy() {}
+
+  void Finalize(AthenaArray<Real> const &w) override;
+  int GetNumVars() const override { return 1; }
+
+ protected:
+  AthenaArray<Real> pf_;
+  Real grav_;
+};
+
+// 4. hydro mean
 class HydroMean : public Diagnostics {
  public:
   HydroMean(MeshBlock *pmb);
   virtual ~HydroMean();
-  void Progress(AthenaArray<Real> const &w);
-  void Finalize(AthenaArray<Real> const &w);
 
- protected:
-  int last_output_cycle_;
+  void Progress(AthenaArray<Real> const &w) override;
+  void Finalize(AthenaArray<Real> const &w) override;
+  int GetNumVars() const override { return NHYDRO; }
 };
 
-// 4. temperature anomaly
+// 5. temperature anomaly
 class TemperatureAnomaly : public Diagnostics {
  public:
   TemperatureAnomaly(MeshBlock *pmb);
@@ -106,15 +127,18 @@ class TemperatureAnomaly : public Diagnostics {
   void Finalize(AthenaArray<Real> const &w);
 };
 
-// 5. pressure anomaly
+// 6. pressure anomaly
 class PressureAnomaly : public Diagnostics {
  public:
   PressureAnomaly(MeshBlock *pmb);
   virtual ~PressureAnomaly() {}
   void Finalize(AthenaArray<Real> const &w);
+
+ protected:
+  AthenaArray<Real> pf_;
 };
 
-// 6. eddy flux
+/* 6. eddy flux
 class EddyFlux : public Diagnostics {
  public:
   EddyFlux(MeshBlock *pmb);
@@ -143,17 +167,6 @@ class HorizontalDivergence : public Diagnostics {
   AthenaArray<Real> v2f2_, v3f3_;
 };
 
-// 9. Buoyancy
-class Buoyancy : public Diagnostics {
- public:
-  Buoyancy(MeshBlock *pmb);
-  virtual ~Buoyancy() {}
-  void Finalize(AthenaArray<Real> const &w);
-
- protected:
-  AthenaArray<Real> pf_;
-  Real grav_;
-};
 
 // 10. total radiative flux
 class RadiativeFlux : public Diagnostics {
@@ -165,10 +178,10 @@ class RadiativeFlux : public Diagnostics {
 };
 
 // 11. total angular momentum
-class AngularMomentum : public Diagnostics {
+class SphericalAngularMomentum : public Diagnostics {
  public:
-  AngularMomentum(MeshBlock *pmb);
-  virtual ~AngularMomentum() {}
+  SphericalAngularMomentum(MeshBlock *pmb);
+  virtual ~SphericalAngularMomentum() {}
   void Finalize(AthenaArray<Real> const &w);
 };
 
@@ -192,4 +205,12 @@ class Tendency : public Diagnostics {
   AthenaArray<Real> wh_, up_;
 };
 
-#endif
+class ConvectiveHeatFlux : public Diagnostics {
+ public:
+  ConvectiveHeatFlux(MeshBlock *pmb);
+  virtual ~ConvectiveHeatFlux() {}
+  void Progress(AthenaArray<Real> const &w);
+  void Finalize(AthenaArray<Real> const &w);
+};*/
+
+#endif  // SRC_DIAGNOSTICS_HPP_
diff --git a/src/diagnostics/diagnostics_factory.cpp b/src/diagnostics/diagnostics_factory.cpp
new file mode 100644
index 00000000..42b96f56
--- /dev/null
+++ b/src/diagnostics/diagnostics_factory.cpp
@@ -0,0 +1,54 @@
+// C/C++
+#include <cstring>
+
+// athena
+#include <athena/parameter_input.hpp>
+
+// diagnostics
+#include "diagnostics.hpp"
+
+DiagnosticsContainer DiagnosticsFactory::CreateFrom(MeshBlock *pmb,
+                                                    ParameterInput *pin) {
+  DiagnosticsContainer diag;
+
+  char cstr[80];
+  std::string diag_names = pin->GetOrAddString("problem", "diagnostics", "");
+  std::strcpy(cstr, diag_names.c_str());
+  char *p = std::strtok(cstr, " ,");
+
+  while (p != NULL) {
+    std::string name(p);
+    if (name == "div") {  // 1.
+      diag.push_back(std::make_shared<Divergence>(pmb));
+    } else if (name == "curl") {  // 2.
+      diag.push_back(std::make_shared<Curl>(pmb));
+    } else if (name == "b") {  // 3.
+      diag.push_back(std::make_shared<Buoyancy>(pmb));
+    } else if (name == "mean") {  // 4.
+      diag.push_back(std::make_shared<HydroMean>(pmb));
+      /*} else if (name == "tempa") { // 4.
+        diag.push_back(std::make_shared<TemperatureAnomaly>(pmb));
+      } else if (name == "presa") { // 5.
+        diag.push_back(std::make_shared<PressureAnomaly>(pmb));
+      } else if (name == "eddyflux") { // 6.
+        diag.push_back(std::make_shared<EddyFlux>(pmb));
+      } else if (name == "hydroflux") { // 7.
+        diag.push_back(std::make_shared<HydroFlux>(pmb));
+      } else if (name == "div_h") { // 8.
+        diag.push_back(std::make_shared<HorizontalDivergence>(pmb));
+      } else if (name == "radflux") { // 10.
+        diag.push_back(std::make_shared<RadiativeFlux>(pmb));
+      } else if (name == "am") { // 11.
+        diag.push_back(std::make_shared<AngularMomentum>(pmb));
+      } else if (name == "eke") { // 12.
+        diag.push_back(std::make_shared<EddyKineticEnergy>(pmb));
+      } else if (name == "tendency") { // 13.
+        diag.push_back(std::make_shared<Tendency>(pmb)); */
+    } else {
+      throw std::invalid_argument("Invalid diagnostic variable: " + name);
+    }
+    p = std::strtok(NULL, " ,");
+  }
+
+  return diag;
+}
diff --git a/src/snap/diagnostics/divergence.cpp b/src/diagnostics/divergence.cpp
similarity index 92%
rename from src/snap/diagnostics/divergence.cpp
rename to src/diagnostics/divergence.cpp
index 45f37eba..0a5f153e 100644
--- a/src/snap/diagnostics/divergence.cpp
+++ b/src/diagnostics/divergence.cpp
@@ -1,22 +1,29 @@
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 Divergence::Divergence(MeshBlock *pmb) : Diagnostics(pmb, "div") {
   type = "SCALARS";
-  units = "1/s";
-  long_name = "divergence";
+
   data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
   v1f1_.NewAthenaArray(ncells3_, ncells2_, ncells1_ + 1);
-  if (pmb->block_size.nx2 > 1)
+
+  if (pmb->block_size.nx2 > 1) {
     v2f2_.NewAthenaArray(ncells3_, ncells2_ + 1, ncells1_);
-  if (pmb->block_size.nx3 > 1)
+  }
+
+  if (pmb->block_size.nx3 > 1) {
     v3f3_.NewAthenaArray(ncells3_ + 1, ncells2_, ncells1_);
+  }
 }
 
 void Divergence::Finalize(AthenaArray<Real> const &w) {
   MeshBlock *pmb = pmy_block_;
   Coordinates *pcoord = pmb->pcoord;
+
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
 
diff --git a/src/diagnostics/hydro_mean.cpp b/src/diagnostics/hydro_mean.cpp
new file mode 100644
index 00000000..b5c5ba8b
--- /dev/null
+++ b/src/diagnostics/hydro_mean.cpp
@@ -0,0 +1,85 @@
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+#include <athena/stride_iterator.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// canoe
+#include "diagnostics.hpp"
+
+HydroMean::HydroMean(MeshBlock *pmb) : Diagnostics(pmb, "mean") {
+  type = "VECTORS";
+  std::string varname = "rho_bar,";
+
+  for (int n = 1; n <= NVAPOR; ++n)
+    varname += "q" + std::to_string(n) + "_bar,";
+  for (int n = 0; n < 3; ++n)
+    varname += "vel" + std::to_string(n + 1) + "_bar,";
+  varname += "T_bar";
+
+  SetName(varname);
+
+  data.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
+}
+
+HydroMean::~HydroMean() { data.DeleteAthenaArray(); }
+
+void HydroMean::Progress(AthenaArray<Real> const &w) {
+  MeshBlock *pmb = pmy_block_;
+
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  if (ncycle_ == 0) {
+    std::fill(data.data(), data.data() + data.GetSize(), 0.);
+  }
+
+  for (int n = 0; n < IPR; ++n)
+    for (int k = ks; k <= ke; ++k)
+      for (int j = js; j <= je; ++j)
+        for (int i = is; i <= ie; ++i) data(n, k, j, i) += w(n, k, j, i);
+
+  // Temperature field is save in IPR
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i)
+        data(IPR, k, j, i) += pthermo->GetTemp(pmb, k, j, i);
+
+  ncycle_++;
+}
+
+void HydroMean::Finalize(AthenaArray<Real> const &w) {
+  MeshBlock *pmb = pmy_block_;
+
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  // hydro mean
+  if (ncycle_ > 0) {
+    for (int n = 0; n < NHYDRO; ++n)
+      for (int k = ks; k <= ke; ++k)
+        for (int j = js; j <= je; ++j)
+          for (int i = is; i <= ie; ++i) {
+            data(n, k, j, i) /= ncycle_;
+          }
+  } else {
+    for (int n = 0; n < IPR; ++n)
+      for (int k = ks; k <= ke; ++k)
+        for (int j = js; j <= je; ++j)
+          for (int i = is; i <= ie; ++i) data(n, k, j, i) = w(n, k, j, i);
+
+    for (int k = ks; k <= ke; ++k)
+      for (int j = js; j <= je; ++j)
+        for (int i = is; i <= ie; ++i)
+          data(IPR, k, j, i) = pthermo->GetTemp(pmb, k, j, i);
+  }
+
+  // clear cycle
+  ncycle_ = 0;
+}
diff --git a/src/snap/diagnostics/convective_heat_flx.cpp b/src/diagnostics/todo/convective_heat_flx.cpp
similarity index 82%
rename from src/snap/diagnostics/convective_heat_flx.cpp
rename to src/diagnostics/todo/convective_heat_flx.cpp
index d8f8a157..05ccd79e 100644
--- a/src/snap/diagnostics/convective_heat_flx.cpp
+++ b/src/diagnostics/todo/convective_heat_flx.cpp
@@ -4,16 +4,23 @@
 #include <mpi.h>
 #endif
 
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
+// athena
+#include <athena/athena_arrays.hpp>
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+#include <athena/stride_iterator.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
-ConvectiveHeatFlx::ConvectiveHeatFlx(MeshBlock *pmb)
-    : Diagnostics(pmb, "conv_heat_flx") {
+ConvectiveHeatFlux::ConvectiveHeatFlux(MeshBlock *pmb)
+    : Diagnostics(pmb, "conv_heat_flx",
+                  "Z-coordinate convective heat flux (eddy and mean)") {
   type = "VECTORS";
   grid = "--C";
-  long_name = "Z-coordinate convective heat flux (eddy and mean)";
   units = "W/(m^2)";
   eddy_.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
   mean_.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
@@ -21,16 +28,10 @@ ConvectiveHeatFlx::ConvectiveHeatFlx(MeshBlock *pmb)
   data.NewAthenaArray(4, 1, 1, ncells1_);
 }
 
-ConvectiveHeatFlx::~ConvectiveHeatFlx() {
-  eddy_.DeleteAthenaArray();
-  mean_.DeleteAthenaArray();
-  data.DeleteAthenaArray();
-}
-
-void ConvectiveHeatFlx::Progress(AthenaArray<Real> const &w) {
+void ConvectiveHeatFlux::Progress(AthenaArray<Real> const &w) {
   MeshBlock *pmb = pmy_block_;
   Coordinates *pcoord = pmb->pcoord;
-  Thermodynamics *pthermo = pmb->pthermo;
+  auto pthermo = Thermodynamics::GetInstance();
 
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
@@ -72,13 +73,15 @@ void ConvectiveHeatFlx::Progress(AthenaArray<Real> const &w) {
       }
 
   // take horizontal average
-  if (ncycle == 0) std::fill(data.data(), data.data() + data.GetSize(), 0.);
+  if (ncycle_ == 0) {
+    std::fill(data.data(), data.data() + data.GetSize(), 0.);
+  }
 
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
       pcoord->CellVolume(k, j, is, ie, vol_);
       for (int i = is; i <= ie; ++i) {
-        Real Cp = pthermo->GetMeanCp(w.at(k, j, i));
+        Real Cp = pthermo->GetCpMass(w.at(k, j, i));
         data(0, i) += Cp * mean_(IDN, k, j, i) * eddy_(IPR, k, j, i) *
                       eddy_(IVX, k, j, i) * vol_(i);
         data(1, i) += Cp * mean_(IDN, k, j, i) * mean_(IPR, k, j, i) *
@@ -88,11 +91,11 @@ void ConvectiveHeatFlx::Progress(AthenaArray<Real> const &w) {
       }
     }
 
-  ncycle++;
+  ncycle_++;
   delete[] data_sum;
 }
 
-void ConvectiveHeatFlx::Finalize(AthenaArray<Real> const &w) {
+void ConvectiveHeatFlux::Finalize(AthenaArray<Real> const &w) {
   MeshBlock *pmb = pmy_block_;
 
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
diff --git a/src/snap/diagnostics/eddy_flux.cpp b/src/diagnostics/todo/eddy_flux.cpp
similarity index 92%
rename from src/snap/diagnostics/eddy_flux.cpp
rename to src/diagnostics/todo/eddy_flux.cpp
index 8124aead..2a3766d0 100644
--- a/src/snap/diagnostics/eddy_flux.cpp
+++ b/src/diagnostics/todo/eddy_flux.cpp
@@ -4,15 +4,19 @@
 #include <mpi.h>
 #endif
 
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+#include <athena/thermodynamics/thermodynamics.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
-EddyFlux::EddyFlux(MeshBlock *pmb) : Diagnostics(pmb, "eddyflux") {
+EddyFlux::EddyFlux(MeshBlock *pmb)
+    : Diagnostics(pmb, "eddyflux", "Z-coordinate eddy flux");
+{
   type = "VECTORS";
   grid = "--C";
-  long_name = "Z-coordinate eddy flux";
   eddy_.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
   mean_.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
   data.NewAthenaArray(2 * NHYDRO, 1, 1, ncells1_);
diff --git a/src/snap/diagnostics/eddy_kinetic_energy.cpp b/src/diagnostics/todo/eddy_kinetic_energy.cpp
similarity index 91%
rename from src/snap/diagnostics/eddy_kinetic_energy.cpp
rename to src/diagnostics/todo/eddy_kinetic_energy.cpp
index 573e117a..4dd2b50d 100644
--- a/src/snap/diagnostics/eddy_kinetic_energy.cpp
+++ b/src/diagnostics/todo/eddy_kinetic_energy.cpp
@@ -7,9 +7,11 @@
  */
 
 // Athena++ headers
-#include "../coordinates/coordinates.hpp"
-#include "../globals.hpp"
-#include "../utils/utils.hpp"
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/globals.hpp>
+#include <athena/utils/utils.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 // MPI headers
@@ -17,10 +19,10 @@
 #include <mpi.h>
 #endif
 
-EddyKineticEnergy::EddyKineticEnergy(MeshBlock *pmb) : Diagnostics(pmb, "eke") {
+EddyKineticEnergy::EddyKineticEnergy(MeshBlock *pmb)
+    : Diagnostics(pmb, "eke", "eddy kinetic energy, mean kinetic energy") {
   type = "VECTORS";
   grid = "--C";
-  long_name = "eddy kinetic energy,mean kinetic energy";
   units = "J/m^3,J/m^3";
   data.NewAthenaArray(2, 1, 1, ncells1_);
 }
diff --git a/src/snap/diagnostics/horizontal_divergence.cpp b/src/diagnostics/todo/horizontal_divergence.cpp
similarity index 95%
rename from src/snap/diagnostics/horizontal_divergence.cpp
rename to src/diagnostics/todo/horizontal_divergence.cpp
index fff28d44..607ad7d0 100644
--- a/src/snap/diagnostics/horizontal_divergence.cpp
+++ b/src/diagnostics/todo/horizontal_divergence.cpp
@@ -1,5 +1,8 @@
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 HorizontalDivergence::HorizontalDivergence(MeshBlock *pmb)
diff --git a/src/snap/diagnostics/latent_heating.cpp b/src/diagnostics/todo/latent_heating.cpp
similarity index 91%
rename from src/snap/diagnostics/latent_heating.cpp
rename to src/diagnostics/todo/latent_heating.cpp
index 80d2062c..088b07b6 100644
--- a/src/snap/diagnostics/latent_heating.cpp
+++ b/src/diagnostics/todo/latent_heating.cpp
@@ -1,13 +1,18 @@
-#include "../coordinates/coordinates.hpp"
-#include "../hydro/hydro.hpp"
-#include "../hydro/srcterms/hydro_srcterms.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/hydro/srcterms/hydro_srcterms.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
-LatentHeating::LatentHeating(MeshBlock *pmb) : Diagnostics(pmb, "Lvheating") {
+LatentHeating::LatentHeating(MeshBlock *pmb)
+    : Diagnostics(pmb, "Lvheating", "Latent Heating") {
   type = "SCALARS";
   grid = "--C";
-  long_name = "Latent Heating";
   units = "W/m^3";
   data.NewAthenaArray(1, 1, 1, ncells1_);
 }
diff --git a/src/snap/diagnostics/pressure_anomaly.cpp b/src/diagnostics/todo/pressure_anomaly.cpp
similarity index 83%
rename from src/snap/diagnostics/pressure_anomaly.cpp
rename to src/diagnostics/todo/pressure_anomaly.cpp
index 28d56946..79ed2e64 100644
--- a/src/snap/diagnostics/pressure_anomaly.cpp
+++ b/src/diagnostics/todo/pressure_anomaly.cpp
@@ -1,9 +1,12 @@
-#include "../coordinates/coordinates.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
-PressureAnomaly::PressureAnomaly(MeshBlock *pmb) : Diagnostics(pmb, "presa") {
+PressureAnomaly::PressureAnomaly(MeshBlock *pmb)
+    : Diagnostics(pmb, "presa", "Pressure anomaly") {
   type = "SCALARS";
-  long_name = "Pressure anomaly";
   units = "pa";
   mean_.NewAthenaArray(ncells1_);
   data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
diff --git a/src/snap/diagnostics/radiative_flux.cpp b/src/diagnostics/todo/radiative_flux.cpp
similarity index 88%
rename from src/snap/diagnostics/radiative_flux.cpp
rename to src/diagnostics/todo/radiative_flux.cpp
index ad91a4c7..c4be9414 100644
--- a/src/snap/diagnostics/radiative_flux.cpp
+++ b/src/diagnostics/todo/radiative_flux.cpp
@@ -4,14 +4,20 @@
 #include <mpi.h>
 #endif
 
-#include "../coordinates/coordinates.hpp"
-#include "../radiation/radiation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+
+// harp
+#include <harp/radiation.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
-RadiativeFlux::RadiativeFlux(MeshBlock *pmb) : Diagnostics(pmb, "radflux") {
+RadiativeFlux::RadiativeFlux(MeshBlock *pmb)
+    : Diagnostics(pmb, "radflux",
+                  "total upward radiative flux,total downward radiative flux") {
   type = "VECTORS";
   grid = "--F";
-  long_name = "total upward radiative flux,total downward radiative flux";
   units = "w/m^2";
   // 0: upward flux
   // 1: downward flux
diff --git a/src/snap/diagnostics/angular_momentum.cpp b/src/diagnostics/todo/spherical_angular_momentum.cpp
similarity index 80%
rename from src/snap/diagnostics/angular_momentum.cpp
rename to src/diagnostics/todo/spherical_angular_momentum.cpp
index 2633492a..6520e8ae 100644
--- a/src/snap/diagnostics/angular_momentum.cpp
+++ b/src/diagnostics/todo/spherical_angular_momentum.cpp
@@ -12,22 +12,27 @@
 #include <cstring>
 #include <sstream>
 #include <stdexcept>
+
 // MPI headers
 #ifdef MPI_PARALLEL
 #include <mpi.h>
 #endif
 
 // Athena++ header
-#include "../coordinates/coordinates.hpp"
-#include "../math/core.h"
+#include <athena/coordinates/coordinates.hpp>
+
+// climath
+#include <climath/core.h>
+
+// canoe
 #include "diagnostics.hpp"
 
-AngularMomentum::AngularMomentum(MeshBlock *pmb) : Diagnostics(pmb, "am") {
+SphericalAngularMomentum::SphericalAngularMomentum(MeshBlock *pmb)
+    : Diagnostics(pmb, "am",
+                  "mass,moment of inertia relative to a thin spherical "
+                  "shell,mean angular velocity") {
   type = "VECTORS";
   grid = "---";
-  long_name =
-      "mass,moment of inertia relative to a thin spherical shell,mean angular "
-      "velocity";
   units = "kg,1,1/s";
   data.NewAthenaArray(3, 1, 1, 1);
   if (strcmp(COORDINATE_SYSTEM, "spherical_polar") != 0) {
@@ -35,11 +40,11 @@ AngularMomentum::AngularMomentum(MeshBlock *pmb) : Diagnostics(pmb, "am") {
     msg << "### FATAL ERROR in AngularMomentum::AngularMomentum" << std::endl
         << "Diagnostics 'am' can only be used in spherical_polar geometry."
         << std::endl;
-    ATHENA_ERROR(msg);
+    throw std::runtime_error(msg.str().c_str());
   }
 }
 
-void AngularMomentum::Finalize(AthenaArray<Real> const &w) {
+void SphericalAngularMomentum::Finalize(AthenaArray<Real> const &w) {
   MeshBlock *pmb = pmy_block_;
   Coordinates *pcoord = pmb->pcoord;
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
@@ -65,7 +70,7 @@ void AngularMomentum::Finalize(AthenaArray<Real> const &w) {
                           sqr(pcoord->x1v(i) * sin(pcoord->x2v(j)));
         // angular momentum
         mass_moi_am[2] += vol_(i) * w(IDN, k, j, i) * pcoord->x1v(i) *
-                          sin(pcoord->x2v(j)) * w(IV3, k, j, i);
+                          sin(pcoord->x2v(j)) * w(IVZ, k, j, i);
       }
     }
 
diff --git a/src/diagnostics/todo/temperature_anomaly.cpp b/src/diagnostics/todo/temperature_anomaly.cpp
new file mode 100644
index 00000000..31cb742b
--- /dev/null
+++ b/src/diagnostics/todo/temperature_anomaly.cpp
@@ -0,0 +1,68 @@
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/globals.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// canoe
+#include "diagnostics.hpp"
+
+TemperatureAnomaly::TemperatureAnomaly(MeshBlock *pmb)
+    : Diagnostics(pmb, "tempa", "Temperature anomaly") {
+  type = "SCALARS";
+  mean_.NewAthenaArray(ncells1_);
+  data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
+
+  Regroup(X1DIR);
+
+  // calculate total volume
+  total_vol_.ZeroClear();
+  for (int k = pmb->ks; k <= pmb->ke; ++k)
+    for (int j = pmb->js; j <= pmb->je; ++j) {
+      pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_);
+
+      for (int i = pmb->is; i <= pmb->ie; ++i) {
+        total_vol_(i) += vol_(i);
+      }
+    }
+
+#ifdef MPI_PARALLEL
+  MPI_Allreduce(MPI_IN_PLACE, total_vol_.data(), total_vol_.GetSize(),
+                MPI_ATHENA_REAL, MPI_SUM, mpi_comm_);
+#endif
+}
+
+void TemperatureAnomaly::Finalize(AthenaArray<Real> const &w) {
+  MeshBlock *pmb = pmy_block_;
+  Coordinates *pcoord = pmb->pcoord;
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+
+  mean_.ZeroClear();
+
+  // calculate horizontal mean
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j) {
+      pcoord->CellVolume(k, j, is, ie, vol_);
+
+      for (int i = is; i <= ie; ++i)
+        mean_(i) += vol_(i) * pthermo->GetTemp(w.at(k, j, i));
+    }
+
+#ifdef MPI_PARALLEL
+  MPI_Allreduce(MPI_IN_PLACE, mean_.data(), mean_.GetSize(), MPI_ATHENA_REAL,
+                MPI_SUM, mpi_comm_);
+#endif
+
+  // temperature anomaly
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        data(k, j, i) =
+            pthermo->GetTemp(w.at(k, j, i)) - mean_(i) / total_vol_(i);
+      }
+}
diff --git a/src/snap/diagnostics/tendency.cpp b/src/diagnostics/todo/tendency.cpp
similarity index 94%
rename from src/snap/diagnostics/tendency.cpp
rename to src/diagnostics/todo/tendency.cpp
index 56f68857..c602196c 100644
--- a/src/snap/diagnostics/tendency.cpp
+++ b/src/diagnostics/todo/tendency.cpp
@@ -1,6 +1,11 @@
-#include "../coordinates/coordinates.hpp"
-#include "../globals.hpp"
-#include "../particles/particles.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/globals.hpp>
+
+// nbody
+#include <nbody/particles/particles.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 Tendency::Tendency(MeshBlock *pmb) : Diagnostics(pmb, "tendency") {
@@ -27,6 +32,7 @@ Tendency::Tendency(MeshBlock *pmb) : Diagnostics(pmb, "tendency") {
 
   int npart = 0;
   Particles *pp = pmb->ppart;
+
   while (pp != nullptr) {
     npart += pp->u.GetDim4();
     for (int n = 0; n < pp->u.GetDim4(); ++n) {
@@ -38,6 +44,9 @@ Tendency::Tendency(MeshBlock *pmb) : Diagnostics(pmb, "tendency") {
     }
     pp = pp->next;
   }
+
+  SetLongName(long_name);
+
   up_.NewAthenaArray(npart, ncells3_, ncells2_, ncells1_);
   data.NewAthenaArray(NHYDRO + npart, 1, 1, ncells1_);
   last_time_ = 0.;
diff --git a/src/snap/diagnostics/total_flux.cpp b/src/diagnostics/todo/total_flux.cpp
similarity index 91%
rename from src/snap/diagnostics/total_flux.cpp
rename to src/diagnostics/todo/total_flux.cpp
index fb660dfb..53ee221b 100644
--- a/src/snap/diagnostics/total_flux.cpp
+++ b/src/diagnostics/todo/total_flux.cpp
@@ -1,6 +1,11 @@
-#include "../coordinates/coordinates.hpp"
-#include "../globals.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/globals.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 HydroFlux::HydroFlux(MeshBlock *pmb) : Diagnostics(pmb, "hydroflux") {
@@ -24,6 +29,8 @@ HydroFlux::HydroFlux(MeshBlock *pmb) : Diagnostics(pmb, "hydroflux") {
   varname += "v1T";
   long_name += "horizontally averaged heat flux";
 
+  SetLongName(long_name);
+
   data.NewAthenaArray(NHYDRO, 1, 1, ncells1_);
 }
 
diff --git a/src/exchanger/exchanger.cpp b/src/exchanger/exchanger.cpp
index 58a92cdc..529c961a 100644
--- a/src/exchanger/exchanger.cpp
+++ b/src/exchanger/exchanger.cpp
@@ -2,9 +2,6 @@
 #include <athena/bvals/bvals.hpp>
 #include <athena/mesh/mesh.hpp>
 
-// canoe
-#include <common.hpp>
-
 namespace ExchangerHelper {
 
 NeighborBlock const *find_bot_neighbor(MeshBlock const *pmb) {
diff --git a/src/exchanger/linear_exchanger.hpp b/src/exchanger/linear_exchanger.hpp
index 423e3758..397d263c 100644
--- a/src/exchanger/linear_exchanger.hpp
+++ b/src/exchanger/linear_exchanger.hpp
@@ -5,9 +5,6 @@
 #include <athena/globals.hpp>
 #include <athena/mesh/mesh.hpp>
 
-// canoe
-#include "common.hpp"
-
 // exchanger
 #include "exchanger.hpp"
 
diff --git a/src/harp/radiation.cpp b/src/harp/radiation.cpp
index f354a8c2..84048ba9 100644
--- a/src/harp/radiation.cpp
+++ b/src/harp/radiation.cpp
@@ -25,7 +25,6 @@
 #include <utils/vectorize.hpp>
 
 // canoe
-#include <common.hpp>
 #include <impl.hpp>
 
 // harp
@@ -162,7 +161,7 @@ void Radiation::CalRadiance(MeshBlock const *pmb, int k, int j) {
   }
 }
 
-void Radiation::CalTimeStep(MeshBlock const* pmb, int k, int j, int il,
+void Radiation::CalTimeStep(MeshBlock const *pmb, int k, int j, int il,
                             int iu) {
   Real total_flux1 = 0., total_flux2 = 0.;
   auto pcoord = pmb->pcoord;
@@ -178,8 +177,8 @@ void Radiation::CalTimeStep(MeshBlock const* pmb, int k, int j, int il,
     for (size_t b = 0; b < bands_.size(); ++b) {
       total_flux2 += flxup(b, k, j, i + 1) - flxdn(b, k, j, i + 1);
     }
-    rtime(k, j, i) = pcoord->dx1f(i) * phydro->u(IEN, k, j, i) 
-      / (total_flux2 - total_flux1);
+    rtime(k, j, i) =
+        pcoord->dx1f(i) * phydro->u(IEN, k, j, i) / (total_flux2 - total_flux1);
     total_flux1 = total_flux2;
     if (rtime(k, j, i) > 0.) {
       time_step_ = std::min(time_step_, rtime(k, j, i));
diff --git a/src/harp/radiation.hpp b/src/harp/radiation.hpp
index 61e21397..88628480 100644
--- a/src/harp/radiation.hpp
+++ b/src/harp/radiation.hpp
@@ -10,9 +10,6 @@
 // athena
 #include <athena/athena.hpp>
 
-// canoe
-#include <common.hpp>  // Direction
-
 // harp
 #include "radiation_band.hpp"
 
diff --git a/src/harp/radiation_band.hpp b/src/harp/radiation_band.hpp
index 19967591..2f3ac173 100644
--- a/src/harp/radiation_band.hpp
+++ b/src/harp/radiation_band.hpp
@@ -15,7 +15,6 @@
 
 // canoe
 #include <air_parcel.hpp>
-#include <common.hpp>
 #include <virtual_groups.hpp>
 
 // exchanger
diff --git a/src/harp/spectral_grid.hpp b/src/harp/spectral_grid.hpp
index 74c6f195..053251b9 100644
--- a/src/harp/spectral_grid.hpp
+++ b/src/harp/spectral_grid.hpp
@@ -7,6 +7,14 @@
 #include <utility>
 #include <vector>
 
+// external
+#include <yaml-cpp/yaml.h>
+
+//! radiation direction
+struct Direction {
+  Real mu, phi;
+};
+
 //! \brief Smallest unit of a spectral grid
 struct SpectralBin {
   Real rad, wav1, wav2, wght;
diff --git a/src/impl.cpp b/src/impl.cpp
index 312c0391..8dd39fde 100644
--- a/src/impl.cpp
+++ b/src/impl.cpp
@@ -40,6 +40,9 @@
 // single column
 #include "single_column/single_column.hpp"
 
+// diagnostics
+#include "diagnostics/diagnostics.hpp"
+
 // canoe
 #include "impl.hpp"
 #include "index_map.hpp"
@@ -75,6 +78,9 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
   // particle queue
   all_particles = ParticlesFactory::Create(pmb, pin);
 
+  // diagnostics
+  all_diags = DiagnosticsFactory::CreateFrom(pmb, pin);
+
   // cubed sphere
   pexo3 = std::make_shared<CubedSphere>(pmb);
 
@@ -85,7 +91,7 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) {
   scheduler = SchedulerFactory::Create(pmb, pin);
 
   // planet
-  planet = PlanetFactory::Create(pin);
+  planet = PlanetFactory::CreateFrom(pmb, pin);
 }
 
 MeshBlock::Impl::~Impl() {}
diff --git a/src/impl.hpp b/src/impl.hpp
index c70effca..ab2c6146 100644
--- a/src/impl.hpp
+++ b/src/impl.hpp
@@ -40,6 +40,7 @@ class FITSOutputGroup;
 class CelestrialBody;
 class CubedSphere;
 class SingleColumn;
+class Diagnostics;
 
 class ExchangerBase;
 
@@ -65,6 +66,7 @@ class MeshBlock::Impl {
 
   std::vector<std::shared_ptr<Inversion>> all_fits;
   std::vector<std::shared_ptr<ParticleBase>> all_particles;
+  std::vector<std::shared_ptr<Diagnostics>> all_diags;
 
   // scheduler
   std::shared_ptr<Scheduler> scheduler;
diff --git a/src/outputs/load_user_output_data.cpp b/src/outputs/load_user_output_data.cpp
index 0df2cf0e..8513c44c 100644
--- a/src/outputs/load_user_output_data.cpp
+++ b/src/outputs/load_user_output_data.cpp
@@ -12,6 +12,9 @@
 #include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// diagnostics
+#include <diagnostics/diagnostics.hpp>
+
 // outputs
 #include "output_utils.hpp"
 
@@ -19,6 +22,20 @@ void OutputType::loadUserOutputData(MeshBlock *pmb) {
   OutputData *pod;
   auto phyd = pmb->phydro;
   auto prad = pmb->pimpl->prad;
+  auto all_diags = pmb->pimpl->all_diags;
+
+  // diagnostic
+  if (output_params.variable.compare("diag") == 0) {
+    for (auto &diag : all_diags) {
+      diag->Finalize(phyd->w);
+      pod = new OutputData;
+      pod->type = diag->type;
+      pod->name = diag->GetName();
+      pod->data.InitWithShallowSlice(diag->data, 4, 0, diag->data.GetDim4());
+      AppendOutputDataNode(pod);
+      num_vars_ += diag->GetNumVars();
+    }
+  }
 
   // vapor
   if (NVAPOR > 0) {
diff --git a/src/outputs/netcdf.cpp b/src/outputs/netcdf.cpp
index 635c8e21..8c98c63f 100644
--- a/src/outputs/netcdf.cpp
+++ b/src/outputs/netcdf.cpp
@@ -21,6 +21,9 @@
 #include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// utils
+#include <utils/vectorize.hpp>
+
 // output
 #include "output_utils.hpp"
 
@@ -245,34 +248,42 @@ void NetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
 
     pdata = pfirst_data_;
     while (pdata != nullptr) {
-      std::string name, attr;
-      std::vector<std::string> varnames;
-      std::string grid = pmeta->GetGridType(pdata->name);
+      auto names = Vectorize<std::string>(pdata->name.c_str(), ",");
+
+      std::string grid = pmeta->GetGridType(names[0]);
       int nvar = get_num_variables(grid, pdata->data);
 
-      for (int n = 0; n < nvar; ++n) {
-        size_t pos = pdata->name.find('?');
-        if (nvar == 1) {                     // SCALARS
-          if (pos < pdata->name.length()) {  // find '?'
-            varnames.push_back(pdata->name.substr(0, pos) +
-                               pdata->name.substr(pos + 1));
-          } else {
-            varnames.push_back(pdata->name);
-          }
-        } else {  // VECTORS
-          char c[16];
-          snprintf(c, sizeof(c), "%d", n + 1);
-          if (pos < pdata->name.length()) {  // find '?'
-            varnames.push_back(pdata->name.substr(0, pos) + c +
-                               pdata->name.substr(pos + 1));
-          } else {
-            varnames.push_back(pdata->name + c);
+      std::vector<std::string> varnames;
+      if (names.size() >= nvar) {
+        for (int n = 0; n < nvar; ++n) {
+          varnames.push_back(names[n]);
+        }
+      } else {
+        for (int n = 0; n < nvar; ++n) {
+          size_t pos = pdata->name.find('?');
+          if (nvar == 1) {                     // SCALARS
+            if (pos < pdata->name.length()) {  // find '?'
+              varnames.push_back(pdata->name.substr(0, pos) +
+                                 pdata->name.substr(pos + 1));
+            } else {
+              varnames.push_back(pdata->name);
+            }
+          } else {  // VECTORS
+            char c[16];
+            snprintf(c, sizeof(c), "%d", n + 1);
+            if (pos < pdata->name.length()) {  // find '?'
+              varnames.push_back(pdata->name.substr(0, pos) + c +
+                                 pdata->name.substr(pos + 1));
+            } else {
+              varnames.push_back(pdata->name + c);
+            }
           }
         }
       }
 
       for (int n = 0; n < nvar; ++n) {
-        name = varnames[n];
+        auto name = varnames[n];
+
         if (grid == "RCC")  // radiation rays
           nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxisr, ivar);
         else if (grid == "CCF")
@@ -291,13 +302,13 @@ void NetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
           nc_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis, ivar);
 
         // set units
-        attr = pmeta->GetUnits(pdata->name);
+        auto attr = pmeta->GetUnits(name);
         if (attr != "") {
           nc_put_att_text(ifile, *ivar, "units", attr.length(), attr.c_str());
         }
 
         // set long_name
-        attr = pmeta->GetLongName(pdata->name);
+        attr = pmeta->GetLongName(name);
         if (attr != "") {
           nc_put_att_text(ifile, *ivar, "long_name", attr.length(),
                           attr.c_str());
diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp
index fb0ee01a..e6d698ef 100644
--- a/src/outputs/output_utils.cpp
+++ b/src/outputs/output_utils.cpp
@@ -26,27 +26,42 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
   Application::Logger app("outputs");
   app->Log("Initialize MetadataTable");
 
-  table_ = {// short name, long name, units, grid location
-            {"x1", "height at cell center", "m", "--C"},
-            {"x1f", "height at cell boundary", "m", "--F"},
-            {"x2", "distance at cell center", "m", "-C-"},
-            {"x2f", "distance at cell boundary", "m", "-F-"},
-            {"x3", "distance at cell center", "m", "C--"},
-            {"x3f", "distance at cell boundary", "m", "F--"},
-            {"rho", "density", "kg/m^3", "CCC"},
-            {"press", "pressure", "pa", "CCC"},
-            {"vel", "velocity", "m/s", "CCC"},
-            {"vapor", "mass mixing ratio of vapor", "kg/kg", "CCC"},
-            {"temp", "temperature", "K", "CCC"},
-            {"theta", "potential temperature", "K", "CCC"},
-            {"thetav", "virtual potential temperature", "K", "CCC"},
-            {"mse", "moist static energy", "J/kg", "CCC"},
-            {"rh1", "relative humidity 1", "1", "CCC"},
-            {"rh2", "relative humidity 2", "1", "CCC"},
-            {"eps", "turbulent dissipation", "w/kg", "CCC"},
-            {"tke", "turbulent kinetic energy", "J/kg", "CCC"},
-            {"mut", "dynamic turbulent viscosity", "kg/(m.s)", "CCC"},
-            {"radiance", "top-of-atmosphere radiance", "K", "RCC"}};
+  table_ = {
+      // short name, long name, units, grid location
+      {"x1", "height at cell center", "m", "--C"},
+      {"x1f", "height at cell boundary", "m", "--F"},
+      {"x2", "distance at cell center", "m", "-C-"},
+      {"x2f", "distance at cell boundary", "m", "-F-"},
+      {"x3", "distance at cell center", "m", "C--"},
+      {"x3f", "distance at cell boundary", "m", "F--"},
+      {"rho", "density", "kg/m^3", "CCC"},
+      {"press", "pressure", "pa", "CCC"},
+      {"vel", "velocity", "m/s", "CCC"},
+      {"vapor", "mass mixing ratio of vapor", "kg/kg", "CCC"},
+      {"temp", "temperature", "K", "CCC"},
+      {"theta", "potential temperature", "K", "CCC"},
+      {"thetav", "virtual potential temperature", "K", "CCC"},
+      {"mse", "moist static energy", "J/kg", "CCC"},
+      {"rh1", "relative humidity 1", "1", "CCC"},
+      {"rh2", "relative humidity 2", "1", "CCC"},
+      {"eps", "turbulent dissipation", "w/kg", "CCC"},
+      {"tke", "turbulent kinetic energy", "J/kg", "CCC"},
+      {"mut", "dynamic turbulent viscosity", "kg/(m.s)", "CCC"},
+      {"radiance", "top-of-atmosphere radiance", "K", "RCC"},
+      {"curl", "curl", "1/s", "CCC"},
+      {"div", "divergence", "1/s", "CCC"},
+      {"b", "buoyancy", "m/s^2", "CCC"},
+      {"rho_bar", "mean density", "kg/m^3", "CCC"},
+      {"q1_bar", "mean vapor 1 mixing ratio", "kg/kg", "CCC"},
+      {"q2_bar", "mean vapor 2 mixing ratio", "kg/kg", "CCC"},
+      {"q3_bar", "mean vapor 3 mixing ratio", "kg/kg", "CCC"},
+      {"vel1_bar", "mean velocity 1", "m/s", "CCC"},
+      {"vel2_bar", "mean velocity 2", "m/s", "CCC"},
+      {"vel3_bar", "mean velocity 3", "m/s", "CCC"},
+      {"T_bar", "mean temperature", "K", "CCC"},
+      {"tempa", "horizontal temperature anomaly", "K", "CCC"},
+      {"presa", "horizontal pressure anomaly", "pa", "CCC"},
+  };
 }
 
 MetadataTable::~MetadataTable() {
diff --git a/src/outputs/pnetcdf.cpp b/src/outputs/pnetcdf.cpp
index b729c9d5..2d9549b5 100644
--- a/src/outputs/pnetcdf.cpp
+++ b/src/outputs/pnetcdf.cpp
@@ -25,6 +25,9 @@
 #include <harp/radiation.hpp>
 #include <harp/radiation_band.hpp>
 
+// utils
+#include <utils/vectorize.hpp>
+
 // outputs
 #include "output_utils.hpp"
 
@@ -258,34 +261,42 @@ void PnetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
 
   pdata = pfirst_data_;
   while (pdata != nullptr) {
-    std::string grid = pmeta->GetGridType(pdata->name);
-    std::string name, attr;
-    std::vector<std::string> varnames;
+    auto names = Vectorize<std::string>(pdata->name.c_str(), ",");
+
+    std::string grid = pmeta->GetGridType(names[0]);
     int nvar = get_num_variables(grid, pdata->data);
 
-    for (int n = 0; n < nvar; ++n) {
-      size_t pos = pdata->name.find('?');
-      if (nvar == 1) {                     // SCALARS
-        if (pos < pdata->name.length()) {  // find '?'
-          varnames.push_back(pdata->name.substr(0, pos) +
-                             pdata->name.substr(pos + 1));
-        } else {
-          varnames.push_back(pdata->name);
-        }
-      } else {  // VECTORS
-        char c[16];
-        snprintf(c, sizeof(c), "%d", n + 1);
-        if (pos < pdata->name.length()) {  // find '?'
-          varnames.push_back(pdata->name.substr(0, pos) + c +
-                             pdata->name.substr(pos + 1));
-        } else {
-          varnames.push_back(pdata->name + c);
+    std::vector<std::string> varnames;
+    if (names.size() >= nvar) {
+      for (int n = 0; n < nvar; ++n) {
+        varnames.push_back(names[n]);
+      }
+    } else {
+      for (int n = 0; n < nvar; ++n) {
+        size_t pos = pdata->name.find('?');
+        if (nvar == 1) {                     // SCALARS
+          if (pos < pdata->name.length()) {  // find '?'
+            varnames.push_back(pdata->name.substr(0, pos) +
+                               pdata->name.substr(pos + 1));
+          } else {
+            varnames.push_back(pdata->name);
+          }
+        } else {  // VECTORS
+          char c[16];
+          snprintf(c, sizeof(c), "%d", n + 1);
+          if (pos < pdata->name.length()) {  // find '?'
+            varnames.push_back(pdata->name.substr(0, pos) + c +
+                               pdata->name.substr(pos + 1));
+          } else {
+            varnames.push_back(pdata->name + c);
+          }
         }
       }
     }
 
     for (int n = 0; n < nvar; ++n) {
-      name = varnames[n];
+      auto name = varnames[n];
+
       if (grid == "RCC")  // radiation rays
         ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxisr, ivar);
       else if (grid == "CCF")
@@ -304,7 +315,7 @@ void PnetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
         ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis, ivar);
 
       // set units
-      attr = pmeta->GetUnits(pdata->name);
+      auto attr = pmeta->GetUnits(pdata->name);
       if (attr != "") {
         ncmpi_put_att_text(ifile, *ivar, "units", attr.length(), attr.c_str());
       }
diff --git a/src/snap/diagnostics/diagnostics.cpp b/src/snap/diagnostics/diagnostics.cpp
deleted file mode 100644
index 428d102c..00000000
--- a/src/snap/diagnostics/diagnostics.cpp
+++ /dev/null
@@ -1,214 +0,0 @@
-// C/C++ headers
-#include <cstring>
-#include <sstream>
-#include <stdexcept>
-
-// Athena++ headers
-#include "../coordinates/coordinates.hpp"
-#include "../debugger/debugger.hpp"
-#include "../globals.hpp"
-#include "../mesh/mesh.hpp"
-#include "diagnostics.hpp"
-
-Diagnostics::Diagnostics(MeshBlock *pmb, ParameterInput *pin)
-    : myname("HEAD"),
-      type(""),
-      grid(""),
-      varname("HEAD"),
-      long_name(""),
-      units(""),
-      prev(nullptr),
-      next(nullptr),
-      ncycle(0),
-      pmy_block_(pmb) {
-  pmb->pdebug->Enter("Diagnostics List");
-  std::stringstream msg;
-  char cstr[80];
-  std::string diag_names = pin->GetOrAddString("problem", "diagnostics", "");
-  std::strcpy(cstr, diag_names.c_str());
-  char *p = std::strtok(cstr, " ,");
-  while (p != NULL) {
-    std::string name(p);
-    if (name == "div") {  // 1.
-      AddDiagnostics(Divergence(pmb));
-    } else if (name == "curl")  // 2.
-      AddDiagnostics(Curl(pmb));
-    else if (name == "mean")  // 3.
-      AddDiagnostics(HydroMean(pmb));
-    else if (name == "tempa")  // 4.
-      AddDiagnostics(TemperatureAnomaly(pmb));
-    else if (name == "presa")  // 5.
-      AddDiagnostics(PressureAnomaly(pmb));
-    else if (name == "eddyflux")  // 6.
-      AddDiagnostics(EddyFlux(pmb));
-    else if (name == "hydroflux")  // 7.
-      AddDiagnostics(HydroFlux(pmb));
-    else if (name == "div_h")  // 8.
-      AddDiagnostics(HorizontalDivergence(pmb));
-    else if (name == "b")  // 9.
-      AddDiagnostics(Buoyancy(pmb));
-    else if (name == "radflux")  // 10.
-      AddDiagnostics(RadiativeFlux(pmb));
-    else if (name == "am")  // 11.
-      AddDiagnostics(AngularMomentum(pmb));
-    else if (name == "eke")  // 12.
-      AddDiagnostics(EddyKineticEnergy(pmb));
-    else if (name == "tendency")  // 13.
-      AddDiagnostics(Tendency(pmb));
-    else {
-      msg << "### FATAL ERROR in function Diagnostics::Diagnostics" << std::endl
-          << "Diagnostic variable " << name << " not defined";
-      ATHENA_ERROR(msg);
-    }
-    // msg << "- add diagnostics " + name << std::endl;
-    p = std::strtok(NULL, " ,");
-  }
-
-#if MAIN_TASKLIST != InversionTaskList
-  if (NGHOST < 2) {
-    msg << "### FATAL ERROR in function Diagnostics::Diagnostics" << std::endl
-        << "Most diagnostic variables require at least 2 ghost cells";
-    ATHENA_ERROR(msg);
-  }
-#endif
-
-  // pmb->pdebug->WriteMessage(msg.str());
-  pmb->pdebug->Leave();
-}
-
-Diagnostics::Diagnostics(MeshBlock *pmb, std::string name)
-    : myname(name),
-      varname(name),
-      prev(nullptr),
-      next(nullptr),
-      ncycle(0),
-      pmy_block_(pmb) {
-  pmb->pdebug->Enter("Diagnostics-" + name);
-  std::stringstream msg;
-  ncells1_ = pmb->block_size.nx1 + 2 * (NGHOST);
-  ncells2_ = 1;
-  ncells3_ = 1;
-  if (pmb->pmy_mesh->f2) ncells2_ = pmb->block_size.nx2 + 2 * (NGHOST);
-  if (pmb->pmy_mesh->f3) ncells3_ = pmb->block_size.nx3 + 2 * (NGHOST);
-
-  x1edge_.NewAthenaArray(ncells1_ + 1);
-  x1edge_p1_.NewAthenaArray(ncells1_);
-  x2edge_.NewAthenaArray(ncells1_ + 1);
-  x2edge_p1_.NewAthenaArray(ncells1_);
-  x3edge_.NewAthenaArray(ncells1_ + 1);
-  x3edge_p1_.NewAthenaArray(ncells1_);
-
-  x1area_.NewAthenaArray(ncells1_ + 1);
-  x2area_.NewAthenaArray(ncells1_);
-  x2area_p1_.NewAthenaArray(ncells1_);
-  x3area_.NewAthenaArray(ncells1_);
-  x3area_p1_.NewAthenaArray(ncells1_);
-
-  vol_.NewAthenaArray(ncells1_);
-  total_vol_.NewAthenaArray(ncells1_);
-  brank_.resize(Globals::nranks);
-  color_.resize(Globals::nranks);
-
-  for (int i = 0; i < Globals::nranks; ++i) {
-    brank_[i] = -1;
-    color_[i] = -1;
-  }
-
-  pmb->pdebug->Leave();
-}
-
-Diagnostics::~Diagnostics() {
-  if (prev != nullptr) prev->next = next;
-  if (next != nullptr) next->prev = prev;
-}
-
-Diagnostics *Diagnostics::operator[](std::string name) {
-  Diagnostics *p = this;
-  while (p != nullptr) {
-    if (p->myname == name) return p;
-    p = p->next;
-  }
-  return p;
-}
-
-void Diagnostics::setColor_(int *color, CoordinateDirection dir) {
-  MeshBlock *pmb = pmy_block_;
-  NeighborBlock bblock, tblock;
-  pmb->FindNeighbors(dir, bblock, tblock);
-
-  if (dir == X1DIR) {
-    if (pmb->block_size.x1min <= pmb->pmy_mesh->mesh_size.x1min) {
-      bblock.snb.gid = -1;
-      bblock.snb.rank = -1;
-    }
-    if (pmb->block_size.x1max >= pmb->pmy_mesh->mesh_size.x1max) {
-      tblock.snb.gid = -1;
-      tblock.snb.rank = -1;
-    }
-  } else if (dir == X2DIR) {
-    if (pmb->block_size.x2min <= pmb->pmy_mesh->mesh_size.x2min) {
-      bblock.snb.gid = -1;
-      bblock.snb.rank = -1;
-    }
-    if (pmb->block_size.x2max >= pmb->pmy_mesh->mesh_size.x2max) {
-      tblock.snb.gid = -1;
-      tblock.snb.rank = -1;
-    }
-  } else {  // X3DIR
-    if (pmb->block_size.x3min <= pmb->pmy_mesh->mesh_size.x3min) {
-      bblock.snb.gid = -1;
-      bblock.snb.rank = -1;
-    }
-    if (pmb->block_size.x3max >= pmb->pmy_mesh->mesh_size.x3max) {
-      tblock.snb.gid = -1;
-      tblock.snb.rank = -1;
-    }
-  }
-
-#ifdef MPI_PARALLEL
-  MPI_Allgather(&bblock.snb.rank, 1, MPI_INT, brank_.data(), 1, MPI_INT,
-                MPI_COMM_WORLD);
-#else
-  brank_[0] = -1;
-#endif
-
-  int c = 0;
-  for (int i = 0; i < Globals::nranks; ++i) {
-    // color[i] = brank_[i] == -1 ? color[i] : color[brank_[i]];
-    if (brank_[i] == -1) {
-      if (color[i] == -1) color[i] = c++;
-    } else
-      color[i] = color[brank_[i]];
-  }
-}
-
-void Diagnostics::gatherAllData23_(AthenaArray<Real> &total_vol,
-                                   AthenaArray<Real> &total_data) {
-  MeshBlock *pmb = pmy_block_;
-
-  // calculate total volume
-  total_vol.ZeroClear();
-  for (int k = pmb->ks; k <= pmb->ke; ++k)
-    for (int j = pmb->js; j <= pmb->je; ++j) {
-      pmb->pcoord->CellVolume(k, j, pmb->is, pmb->ie, vol_);
-      for (int i = pmb->is; i <= pmb->ie; ++i) total_vol(i) += vol_(i);
-    }
-
-    // sum over all ranks
-#ifdef MPI_PARALLEL
-  MPI_Comm comm;
-  std::fill(color_.data(), color_.data() + Globals::nranks, -1);
-  setColor_(color_.data(), X2DIR);
-  setColor_(color_.data(), X3DIR);
-  MPI_Comm_split(MPI_COMM_WORLD, color_[Globals::my_rank], Globals::my_rank,
-                 &comm);
-  // int size;
-  // MPI_Comm_size(comm, &size);
-  // std::cout << size << std::endl;
-  MPI_Allreduce(MPI_IN_PLACE, total_vol.data(), total_vol.GetSize(),
-                MPI_ATHENA_REAL, MPI_SUM, comm);
-  MPI_Allreduce(MPI_IN_PLACE, total_data.data(), total_data.GetSize(),
-                MPI_ATHENA_REAL, MPI_SUM, comm);
-  MPI_Comm_free(&comm);
-#endif
-}
diff --git a/src/snap/diagnostics/hydro_mean.cpp b/src/snap/diagnostics/hydro_mean.cpp
deleted file mode 100644
index b8d03c0e..00000000
--- a/src/snap/diagnostics/hydro_mean.cpp
+++ /dev/null
@@ -1,75 +0,0 @@
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
-#include "diagnostics.hpp"
-
-HydroMean::HydroMean(MeshBlock *pmb) : Diagnostics(pmb, "mean") {
-  type = "VECTORS";
-  varname = "rho_bar,";
-  long_name = "mean density,";
-  units = "kg/m^3,";
-
-  for (int n = 1; n <= NVAPOR; ++n) {
-    varname += "q" + std::to_string(n) + "_bar,";
-    units += "kg/kg,";
-    long_name += "mean vapor,";
-  }
-  for (int n = 0; n < 3; ++n) {
-    units += "m/s,";
-    varname += "vel" + std::to_string(n + 1) + "_bar,";
-    long_name += "mean velocity,";
-  }
-  units += "K";
-  varname += "T_bar";
-  long_name += "mean temperature";
-
-  data.NewAthenaArray(NHYDRO, ncells3_, ncells2_, ncells1_);
-}
-
-HydroMean::~HydroMean() { data.DeleteAthenaArray(); }
-
-void HydroMean::Progress(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
-
-  if (ncycle == 0) std::fill(data.data(), data.data() + data.GetSize(), 0.);
-
-  for (int n = 0; n < IPR; ++n)
-    for (int k = pmb->ks; k <= pmb->ke; ++k)
-      for (int j = pmb->js; j <= pmb->je; ++j)
-        for (int i = pmb->is; i <= pmb->ie; ++i)
-          data(n, k, j, i) += w(n, k, j, i);
-
-  // Temperature field is save in IPR
-  for (int k = pmb->ks; k <= pmb->ke; ++k)
-    for (int j = pmb->js; j <= pmb->je; ++j)
-      for (int i = pmb->is; i <= pmb->ie; ++i)
-        data(IPR, k, j, i) += pmb->pthermo->GetTemp(w.at(k, j, i));
-
-  ncycle++;
-}
-
-void HydroMean::Finalize(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
-
-  // hydro mean
-  if (ncycle > 0) {
-    for (int n = 0; n < NHYDRO; ++n)
-      for (int k = pmb->ks; k <= pmb->ke; ++k)
-        for (int j = pmb->js; j <= pmb->je; ++j)
-          for (int i = pmb->is; i <= pmb->ie; ++i) data(n, k, j, i) /= ncycle;
-  } else {
-    for (int n = 0; n < IPR; ++n)
-      for (int k = pmb->ks; k <= pmb->ke; ++k)
-        for (int j = pmb->js; j <= pmb->je; ++j)
-          for (int i = pmb->is; i <= pmb->ie; ++i)
-            data(n, k, j, i) = w(n, k, j, i);
-
-    for (int k = pmb->ks; k <= pmb->ke; ++k)
-      for (int j = pmb->js; j <= pmb->je; ++j)
-        for (int i = pmb->is; i <= pmb->ie; ++i)
-          data(IPR, k, j, i) = pmb->pthermo->GetTemp(w.at(k, j, i));
-  }
-
-  // clear cycle
-  ncycle = 0;
-}
diff --git a/src/snap/diagnostics/temperature_anomaly.cpp b/src/snap/diagnostics/temperature_anomaly.cpp
deleted file mode 100644
index a042bf46..00000000
--- a/src/snap/diagnostics/temperature_anomaly.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#include "../coordinates/coordinates.hpp"
-#include "../globals.hpp"
-#include "../thermodynamics/thermodynamics.hpp"
-#include "diagnostics.hpp"
-
-TemperatureAnomaly::TemperatureAnomaly(MeshBlock *pmb)
-    : Diagnostics(pmb, "tempa") {
-  type = "SCALARS";
-  long_name = "Temperature anomaly";
-  units = "K";
-  mean_.NewAthenaArray(ncells1_);
-  data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
-}
-
-void TemperatureAnomaly::Finalize(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
-  Coordinates *pcoord = pmb->pcoord;
-  Thermodynamics *pthermo = pmb->pthermo;
-  int is = pmb->is, js = pmb->js, ks = pmb->ks;
-  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
-
-  mean_.ZeroClear();
-
-  // calculate horizontal mean
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j) {
-      pcoord->CellVolume(k, j, is, ie, vol_);
-      for (int i = is; i <= ie; ++i)
-        mean_(i) += vol_(i) * pthermo->GetTemp(w.at(k, j, i));
-    }
-
-  gatherAllData23_(total_vol_, mean_);
-
-  // temperature anomaly
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j)
-      for (int i = is; i <= ie; ++i)
-        data(k, j, i) =
-            pthermo->GetTemp(w.at(k, j, i)) - mean_(i) / total_vol_(i);
-}
diff --git a/src/virtual_groups.hpp b/src/virtual_groups.hpp
index c44e40fd..c83bd963 100644
--- a/src/virtual_groups.hpp
+++ b/src/virtual_groups.hpp
@@ -24,11 +24,19 @@ class OutputParameters;
 class NamedGroup {
  public:
   explicit NamedGroup(std::string name) : myname_(name) {}
+  NamedGroup(std::string name, std::string long_name)
+      : myname_(name), long_name_(long_name) {}
+
   virtual ~NamedGroup() {}
+
   std::string GetName() const { return myname_; }
+  std::string GetLongName() const { return long_name_; }
+  void SetName(std::string name) { myname_ = name; }
+  void SetLongName(std::string long_name) { long_name_ = long_name; }
 
  private:
   std::string myname_;
+  std::string long_name_;
 };
 
 class FlagGroup {