From 6bb06f43defd7e8d8df06b9fc0d0e81e2d786205 Mon Sep 17 00:00:00 2001
From: stormy/cli <chengcli@umich.edu>
Date: Wed, 10 Apr 2024 23:51:05 -0400
Subject: [PATCH] wip

---
 examples/2024-YZhou-jupiter/jupiter_dry.inp | 26 +++---
 src/diagnostics/anomaly.cpp                 | 91 +++++++++++++++++++++
 src/diagnostics/diagnostics.hpp             | 53 ++++++------
 src/diagnostics/diagnostics_factory.cpp     |  8 +-
 src/diagnostics/hydro_flux.cpp              |  3 +
 src/diagnostics/pressure_anomaly.cpp        | 67 ---------------
 src/diagnostics/radiative_flux.cpp          |  3 +
 src/diagnostics/temperature_anomaly.cpp     | 73 -----------------
 src/diagnostics/v1_moments.cpp              | 73 +++++++++++++++++
 src/outputs/output_utils.cpp                | 19 ++++-
 src/outputs/pnetcdf.cpp                     |  4 +-
 11 files changed, 231 insertions(+), 189 deletions(-)
 create mode 100644 src/diagnostics/anomaly.cpp
 delete mode 100644 src/diagnostics/pressure_anomaly.cpp
 delete mode 100644 src/diagnostics/temperature_anomaly.cpp
 create mode 100644 src/diagnostics/v1_moments.cpp

diff --git a/examples/2024-YZhou-jupiter/jupiter_dry.inp b/examples/2024-YZhou-jupiter/jupiter_dry.inp
index 997bc962..b4215e24 100644
--- a/examples/2024-YZhou-jupiter/jupiter_dry.inp
+++ b/examples/2024-YZhou-jupiter/jupiter_dry.inp
@@ -6,7 +6,7 @@ problem_id  = jupiter2d     # problem ID: basename of output filenames
 
 <output0>
 file_type   = rst
-dt          = 100.E5
+dt          = 10.E5
 
 <output1>
 file_type   = hst           # History data dump
@@ -35,33 +35,33 @@ dt          = 4.E4          # time increment
 <time>
 cfl_number  = 0.9
 nlim        = -1            # cycle limit
-tlim        = 200.E5
+tlim        = 1000.E5
 xorder      = 5             # horizontal reconstruction order
 integrator  = rk3           # integration method
 
 <mesh>
-nx1         = 80          # Number of zones in X1-direction
+nx1         = 40          # 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
+nx2         = 80         # Number of zones in X2-direction
 x2min       = 0.          # minimum value of X2
-x2max       = 3200.E3     # maximum value of X2
+x2max       = 1600.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
+nx3         = 80         # Number of zones in X3-direction
+x3min       = 0.        	# minimum value of X3
+x3max       = 1600.E3    	# 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
+nx1         = 40
+nx2         = 20
+nx3         = 20
 
 <hydro>
 gamma       = 1.42
@@ -81,7 +81,7 @@ Rd          = 3777.   # mu = 2.3175 g/mol
 P0          = 1.E5
 T0          = 169.
 Tmin        = 110.
-hrate       = -100.
+hrate       = -20.
 prad        = 2.e5
 init_Ttol   = 2.
-diagnostics = div, curl, b, mean, div_h, tempa, presa, radflux, hydroflux
+diagnostics = div, curl, b, mean, div_h, anomaly, radflux, hydroflux
diff --git a/src/diagnostics/anomaly.cpp b/src/diagnostics/anomaly.cpp
new file mode 100644
index 00000000..54de0ee3
--- /dev/null
+++ b/src/diagnostics/anomaly.cpp
@@ -0,0 +1,91 @@
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// exchanger
+#include <exchanger/exchanger.hpp>
+
+// diagnostics
+#include "diagnostics.hpp"
+
+Anomaly::Anomaly(MeshBlock *pmb) : Diagnostics(pmb, "rhoa,tempa,v1a,presa") {
+  type = "VECTORS";
+  data.NewAthenaArray(4, ncells3_, ncells2_, ncells1_);
+  mean_.resize(4 * ncells1_);
+
+  pexh = std::make_shared<PlanarExchanger<Real, 2>>("diag/anomaly");
+
+  // volumn
+  pexh->send_buffer[0].resize(total_vol_.size());
+  pexh->recv_buffer[0].resize(total_vol_.size());
+
+  // value
+  pexh->send_buffer[1].resize(4 * ncells1_);
+  pexh->recv_buffer[1].resize(4 * ncells1_);
+}
+
+void Anomaly::Finalize(MeshBlock *pmb) {
+  Coordinates *pcoord = pmb->pcoord;
+  auto const &w = pmb->phydro->w;
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+
+  std::fill(total_vol_.begin(), total_vol_.end(), 0.);
+  std::fill(mean_.begin(), mean_.end(), 0.);
+
+  // 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) {
+        total_vol_[i] += vol_(i);
+        // density anomaly
+        mean_[i] += vol_(i) * w(IDN, k, j, i);
+        // temperature anomaly
+        mean_[ncells1_ + i] += vol_(i) * pthermo->GetTemp(pmb, k, j, i);
+        // vertical velocity anomaly
+        mean_[2 * ncells1_ + i] += vol_(i) * w(IVX, k, j, i);
+        // pressure anomaly
+        mean_[3 * ncells1_ + i] += vol_(i) * w(IPR, k, j, i);
+      }
+    }
+
+  pexh->Regroup(pmb, X1DIR);
+  packData(pmb);
+  pexh->ReduceAll(pmb, MPI_SUM);
+  unpackData(pmb);
+
+  // pressure anomaly
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        // density anomaly
+        data(0, k, j, i) = w(IDN, k, j, i) - mean_[i] / total_vol_[i];
+        // temperature anomaly
+        data(1, k, j, i) = pthermo->GetTemp(pmb, k, j, i) -
+                           mean_[ncells1_ + i] / total_vol_[i];
+        // vertical velocity anomaly
+        data(2, k, j, i) =
+            w(IVX, k, j, i) - mean_[2 * ncells1_ + i] / total_vol_[i];
+        // pressure anomaly
+        data(3, k, j, i) =
+            w(IPR, k, j, i) - mean_[3 * ncells1_ + i] / total_vol_[i];
+      }
+}
+
+void Anomaly::packData(MeshBlock const *pmb) {
+  pexh->send_buffer[0].swap(total_vol_);
+  pexh->send_buffer[1].swap(mean_);
+}
+
+void Anomaly::unpackData(MeshBlock const *pmb) {
+  total_vol_.swap(pexh->recv_buffer[0]);
+  mean_.swap(pexh->recv_buffer[1]);
+}
diff --git a/src/diagnostics/diagnostics.hpp b/src/diagnostics/diagnostics.hpp
index 994fb8f0..dfdc8955 100644
--- a/src/diagnostics/diagnostics.hpp
+++ b/src/diagnostics/diagnostics.hpp
@@ -12,12 +12,12 @@
 #include <configure.hpp>
 #include <virtual_groups.hpp>
 
-// exchanger
-#include <exchanger/exchanger.hpp>
-
 class ParameterInput;
 class Meshlock;
 
+template <typename T, int N>
+class PlanarExchanger;
+
 class Diagnostics : public NamedGroup {
  public:
   // data
@@ -128,37 +128,17 @@ class HorizontalDivergence : public Diagnostics {
   AthenaArray<Real> v2f2_, v3f3_;
 };
 
-// 6. temperature anomaly
-class TemperatureAnomaly : public Diagnostics {
+// 6. anomaly
+class Anomaly : public Diagnostics {
  public:
   std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
 
  public:
-  TemperatureAnomaly(MeshBlock *pmb);
-  virtual ~TemperatureAnomaly() {}
+  Anomaly(MeshBlock *pmb);
+  virtual ~Anomaly() {}
 
   void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
-
- protected:
-  void packData(MeshBlock const *pmb);
-  void unpackData(MeshBlock const *pmb);
-
-  //! mean component
-  std::vector<Real> mean_;
-};
-
-// 7. pressure anomaly
-class PressureAnomaly : public Diagnostics {
- public:
-  std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
-
- public:
-  PressureAnomaly(MeshBlock *pmb);
-  virtual ~PressureAnomaly() {}
-
-  void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
+  int GetNumVars() const override { return 4; }
 
  protected:
   void packData(MeshBlock const *pmb);
@@ -208,6 +188,23 @@ class HydroFlux : public Diagnostics {
   int ncycle_;
 };
 
+// 10. w_avg
+class V1Moments : public Diagnostics {
+ public:
+  std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
+
+ public:
+  V1Moments(MeshBlock *pmb);
+  virtual ~V1Moments() {}
+
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return 3; }
+
+ protected:
+  void packData(MeshBlock const *pmb);
+  void unpackData(MeshBlock const *pmb);
+};
+
 /* 6. eddy flux
 class EddyFlux : public Diagnostics {
  public:
diff --git a/src/diagnostics/diagnostics_factory.cpp b/src/diagnostics/diagnostics_factory.cpp
index 2b44d87d..c4de7c88 100644
--- a/src/diagnostics/diagnostics_factory.cpp
+++ b/src/diagnostics/diagnostics_factory.cpp
@@ -28,14 +28,14 @@ DiagnosticsContainer DiagnosticsFactory::CreateFrom(MeshBlock *pmb,
       diag.push_back(std::make_shared<HydroMean>(pmb));
     } else if (name == "div_h") {  // 5.
       diag.push_back(std::make_shared<HorizontalDivergence>(pmb));
-    } else if (name == "tempa") {  // 6.
-      diag.push_back(std::make_shared<TemperatureAnomaly>(pmb));
-    } else if (name == "presa") {  // 7.
-      diag.push_back(std::make_shared<PressureAnomaly>(pmb));
+    } else if (name == "anomaly") {  // 6.
+      diag.push_back(std::make_shared<Anomaly>(pmb));
     } else if (name == "radflux") {  // 8.
       diag.push_back(std::make_shared<RadiativeFlux>(pmb));
     } else if (name == "hydroflux") {  // 9.
       diag.push_back(std::make_shared<HydroFlux>(pmb));
+    } else if (name == "w_avg") {  // 10.
+      diag.push_back(std::make_shared<V1Moments>(pmb));
       /*} else if (name == "eddyflux") { // 6.
         diag.push_back(std::make_shared<EddyFlux>(pmb));
       } else if (name == "am") { // 11.
diff --git a/src/diagnostics/hydro_flux.cpp b/src/diagnostics/hydro_flux.cpp
index 581f2c6f..e33424d8 100644
--- a/src/diagnostics/hydro_flux.cpp
+++ b/src/diagnostics/hydro_flux.cpp
@@ -2,6 +2,9 @@
 #include <athena/coordinates/coordinates.hpp>
 #include <athena/hydro/hydro.hpp>
 
+// exchanger
+#include <exchanger/exchanger.hpp>
+
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
diff --git a/src/diagnostics/pressure_anomaly.cpp b/src/diagnostics/pressure_anomaly.cpp
deleted file mode 100644
index 6e97c267..00000000
--- a/src/diagnostics/pressure_anomaly.cpp
+++ /dev/null
@@ -1,67 +0,0 @@
-// athena
-#include <athena/coordinates/coordinates.hpp>
-#include <athena/hydro/hydro.hpp>
-
-// diagnostics
-#include "diagnostics.hpp"
-
-PressureAnomaly::PressureAnomaly(MeshBlock *pmb) : Diagnostics(pmb, "presa") {
-  type = "SCALARS";
-  data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
-  mean_.resize(ncells1_);
-
-  pexh = std::make_shared<PlanarExchanger<Real, 2>>("diag/presa");
-
-  // volumn
-  pexh->send_buffer[0].resize(total_vol_.size());
-  pexh->recv_buffer[0].resize(total_vol_.size());
-
-  // value
-  pexh->send_buffer[1].resize(ncells1_);
-  pexh->recv_buffer[1].resize(ncells1_);
-}
-
-void PressureAnomaly::Finalize(MeshBlock *pmb) {
-  Coordinates *pcoord = pmb->pcoord;
-  auto const &w = pmb->phydro->w;
-
-  int is = pmb->is, js = pmb->js, ks = pmb->ks;
-  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
-
-  std::fill(total_vol_.begin(), total_vol_.end(), 0.);
-  std::fill(mean_.begin(), mean_.end(), 0.);
-
-  // 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) {
-        total_vol_[i] += vol_(i);
-        mean_[i] += vol_(i) * w(IPR, k, j, i);
-      }
-    }
-
-  pexh->Regroup(pmb, X1DIR);
-
-  packData(pmb);
-  pexh->ReduceAll(pmb, MPI_SUM);
-  unpackData(pmb);
-
-  // pressure 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) = w(IPR, k, j, i) - mean_[i] / total_vol_[i];
-      }
-}
-
-void PressureAnomaly::packData(MeshBlock const *pmb) {
-  pexh->send_buffer[0].swap(total_vol_);
-  pexh->send_buffer[1].swap(mean_);
-}
-
-void PressureAnomaly::unpackData(MeshBlock const *pmb) {
-  total_vol_.swap(pexh->recv_buffer[0]);
-  mean_.swap(pexh->recv_buffer[1]);
-}
diff --git a/src/diagnostics/radiative_flux.cpp b/src/diagnostics/radiative_flux.cpp
index 045a6cbb..741fe9c5 100644
--- a/src/diagnostics/radiative_flux.cpp
+++ b/src/diagnostics/radiative_flux.cpp
@@ -1,6 +1,9 @@
 // athena
 #include <athena/coordinates/coordinates.hpp>
 
+// exchanger
+#include <exchanger/exchanger.hpp>
+
 // harp
 #include <harp/radiation.hpp>
 
diff --git a/src/diagnostics/temperature_anomaly.cpp b/src/diagnostics/temperature_anomaly.cpp
deleted file mode 100644
index d06566b7..00000000
--- a/src/diagnostics/temperature_anomaly.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-// athena
-#include <athena/coordinates/coordinates.hpp>
-#include <athena/hydro/hydro.hpp>
-
-// snap
-#include <snap/thermodynamics/thermodynamics.hpp>
-
-// diagnostics
-#include "diagnostics.hpp"
-
-TemperatureAnomaly::TemperatureAnomaly(MeshBlock *pmb)
-    : Diagnostics(pmb, "tempa") {
-  type = "SCALARS";
-  data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
-  mean_.resize(ncells1_);
-
-  pexh = std::make_shared<PlanarExchanger<Real, 2>>("diag/tempa");
-
-  // volumn
-  pexh->send_buffer[0].resize(total_vol_.size());
-  pexh->recv_buffer[0].resize(total_vol_.size());
-
-  // value
-  pexh->send_buffer[1].resize(ncells1_);
-  pexh->recv_buffer[1].resize(ncells1_);
-}
-
-void TemperatureAnomaly::Finalize(MeshBlock *pmb) {
-  Coordinates *pcoord = pmb->pcoord;
-  auto const &w = pmb->phydro->w;
-
-  auto pthermo = Thermodynamics::GetInstance();
-
-  int is = pmb->is, js = pmb->js, ks = pmb->ks;
-  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
-
-  std::fill(total_vol_.begin(), total_vol_.end(), 0.);
-  std::fill(mean_.begin(), mean_.end(), 0.);
-
-  // 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) {
-        total_vol_[i] += vol_(i);
-        mean_[i] += vol_(i) * pthermo->GetTemp(pmb, k, j, i);
-      }
-    }
-
-  pexh->Regroup(pmb, X1DIR);
-  packData(pmb);
-  pexh->ReduceAll(pmb, MPI_SUM);
-  unpackData(pmb);
-
-  // pressure 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(pmb, k, j, i) - mean_[i] / total_vol_[i];
-      }
-}
-
-void TemperatureAnomaly::packData(MeshBlock const *pmb) {
-  pexh->send_buffer[0].swap(total_vol_);
-  pexh->send_buffer[1].swap(mean_);
-}
-
-void TemperatureAnomaly::unpackData(MeshBlock const *pmb) {
-  total_vol_.swap(pexh->recv_buffer[0]);
-  mean_.swap(pexh->recv_buffer[1]);
-}
diff --git a/src/diagnostics/v1_moments.cpp b/src/diagnostics/v1_moments.cpp
new file mode 100644
index 00000000..3a79ccd6
--- /dev/null
+++ b/src/diagnostics/v1_moments.cpp
@@ -0,0 +1,73 @@
+// athena++
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/mesh/mesh.hpp>
+
+// exchanger
+#include <exchanger/exchanger.hpp>
+
+// canoe
+#include "diagnostics.hpp"
+
+V1Moments::V1Moments(MeshBlock *pmb) : Diagnostics(pmb, "w_avg,w2_avg,w3_avg") {
+  type = "VECTORS";
+
+  data.NewAthenaArray(3, ncells1_);
+
+  pexh = std::make_shared<PlanarExchanger<Real, 2>>("diag/v1_moments");
+
+  // volumn
+  pexh->send_buffer[0].resize(total_vol_.size());
+  pexh->recv_buffer[0].resize(total_vol_.size());
+  // value
+  pexh->send_buffer[1].resize(3 * ncells1_);
+  pexh->recv_buffer[1].resize(3 * ncells1_);
+}
+
+void V1Moments::Finalize(MeshBlock *pmb) {
+  std::fill(total_vol_.begin(), total_vol_.end(), 0.);
+  std::fill(data.data(), data.data() + data.GetSize(), 0.);
+
+  // calculate total volume
+  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);
+        Real v1 = pmb->phydro->w(IVX, k, j, i);
+        data(0, i) += vol_(i) * v1;
+        data(1, i) += vol_(i) * v1 * v1;
+        data(2, i) += vol_(i) * v1 * v1 * v1;
+      }
+    }
+
+  pexh->Regroup(pmb, X1DIR);
+  packData(pmb);
+  pexh->ReduceAll(pmb, MPI_SUM);
+  unpackData(pmb);
+
+  for (int i = pmb->is; i <= pmb->ie; ++i)
+    for (int n = 0; n < 3; ++n) {
+      data(n, i) /= total_vol_[i];
+    }
+}
+
+void V1Moments::packData(MeshBlock const *pmb) {
+  pexh->send_buffer[0].swap(total_vol_);
+
+  for (int n = 0; n < 3; ++n)
+    for (int i = pmb->is; i <= pmb->ie; ++i) {
+      pexh->send_buffer[1][n * ncells1_ + i] = data(n, i);
+    }
+}
+
+void V1Moments::unpackData(MeshBlock const *pmb) {
+  total_vol_.swap(pexh->recv_buffer[0]);
+
+  for (int n = 0; n < 3; ++n) {
+    for (int i = pmb->is; i <= pmb->ie; ++i) {
+      data(n, i) = pexh->recv_buffer[1][n * ncells1_ + i];
+    }
+  }
+}
diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp
index 1dcfbcd8..2a6a1512 100644
--- a/src/outputs/output_utils.cpp
+++ b/src/outputs/output_utils.cpp
@@ -42,15 +42,22 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
       {"theta", "potential temperature", "K", "CCC"},
       {"thetav", "virtual potential temperature", "K", "CCC"},
       {"mse", "moist static energy", "J/kg", "CCC"},
+      // relative humidity
       {"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"},
+      // radiation
       {"radiance", "top-of-atmosphere radiance", "K", "RCC"},
+      // curl
       {"curl", "curl", "1/s", "CCC"},
+      // divergence
       {"div", "divergence", "1/s", "CCC"},
+      {"div_h", "horizontal divergence", "1/s", "CCC"},
+      // buoyancy
       {"b", "buoyancy", "m/s^2", "CCC"},
+      // hydro mean
       {"rho_bar", "mean density", "kg/m^3", "CCC"},
       {"q1_bar", "mean vapor1 mixing ratio", "kg/kg", "CCC"},
       {"q2_bar", "mean vapor2 mixing ratio", "kg/kg", "CCC"},
@@ -59,12 +66,16 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
       {"vel2_bar", "mean velocity 2", "m/s", "CCC"},
       {"vel3_bar", "mean velocity 3", "m/s", "CCC"},
       {"T_bar", "mean temperature", "K", "CCC"},
-      {"div_h", "horizontal divergence", "1/s", "CCC"},
+      // anomalies
+      {"rhoa", "horizontal density anomaly", "kg/m^3", "CCC"},
       {"tempa", "horizontal temperature anomaly", "K", "CCC"},
+      {"v1a", "horizontal vertical velocity anomaly", "m/s", "CCC"},
       {"presa", "horizontal pressure anomaly", "pa", "CCC"},
+      // radiative flux
       {"rflx_up", "total upward radiative flux", "w/m^2", "--F"},
       {"rflx_dn", "total downward radiative flux", "w/m^2", "--F"},
-      {"v1rho", "total upward mass flux", "kg/(m^2.s)", "--F"},
+      // hydro flux
+      {"v1rho", "total upward mass flux", "kg/(m^2.s)", "--C"},
       {"v1q1", "total upward vapor1 flux", "1/(m^2.s)", "--C"},
       {"v1q2", "total upward vapor2 flux", "1/(m^2.s)", "--C"},
       {"v1q3", "total upward vapor3 flux", "1/(m^2.s)", "--C"},
@@ -72,6 +83,10 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
       {"v1v2", "total upward velocity 2 flux", "(m/s)/(m^2.s)", "--C"},
       {"v1v3", "total upward velocity 3 flux", "(m/s)/(m^2.s)", "--C"},
       {"v1T", "total upward temperature flux", "K/(m^2.s)", "--C"},
+      // v1 moments
+      {"w_avg", "mean vertical velocity", "m/s", "--C"},
+      {"w2_avg", "mean squard vertical velocity", "(m/s)^2", "--C"},
+      {"w3_avg", "mean cubed vertical velocity", "(m/s)^3", "--C"},
       // CMPI6 1D variables description
       {"ta_avg", "domain avg. air temperature profile", "K", "--C"},
       {"ua_avg", "domain avg. zonal wind profile", "m/s", "--C"},
diff --git a/src/outputs/pnetcdf.cpp b/src/outputs/pnetcdf.cpp
index 1d3ec8c5..010d558d 100644
--- a/src/outputs/pnetcdf.cpp
+++ b/src/outputs/pnetcdf.cpp
@@ -310,13 +310,13 @@ void PnetcdfOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, bool flag) {
         ncmpi_def_var(ifile, name.c_str(), NC_FLOAT, 4, iaxis, ivar);
 
       // set units
-      auto attr = pmeta->GetUnits(pdata->name);
+      auto attr = pmeta->GetUnits(name);
       if (attr != "") {
         ncmpi_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 != "") {
         ncmpi_put_att_text(ifile, *ivar, "long_name", attr.length(),
                            attr.c_str());