Canoe Source：Canoe的诊断量输出

It has been 500 days since the last update, the content of the article may be outdated.

Canoe的输出解析

在Juno-MWR-EQ的模拟中，我们设定的输出诊断量如下：

yaml

## /juno.inp
<output2>
file_type   = netcdf      # NetCDF data dump
variable    = prim        # variables to be output
dt          = 1.E-9       # time increment between outputs

<output3>
file_type   = netcdf
variable    = uov
dt          = 1.E-9

<output4>
file_type   = netcdf
variable    = radtoa
dt          = 1.E-9

<output5>
file_type   = netcdf
variable    = radtau
dt          = 1.E-9

dt是输出步长，prim为primitive量，类似的还有cons,见下文。
那么，这些输出设定是如何被模式解析的呢，解析以后通过什么方式被输出的呢？
在粗略研读了源码之后，个人理解：Canoe有三套输出解析系统，分别是：基于Athena++的输出解析模块，Canoe扩展输出模块，以及UOV用户自定义输出模块。

Athena++ Output解析

在athena++模式中，inp文件会被athena/src/outputs/outputs.cpp模块解析。
athena++会按照上述文件的设定的关键字，对输出量进行读取，并扩展到输出对象中：

Required parameters that must be specified in an <output[n]> block are:

variable = cons,prim,D,d,E,e,m,m1,m2,m3,v,v1=vx,v2=vy,v3=vz,p,
bcc,bcc1,bcc2,bcc3,b,b1,b2,b3,phi,uov
file_type = rst,tab,vtk,hst,hdf5

dt = problem time between outputs
EXAMPLE of an <output[n]> block for a VTK dump:

yaml

<output3>
 file_type   = tab       # Tabular data dump
 variable    = prim      # variables to be output
 data_format = %12.5e    # Optional data format string
 dt          = 0.01      # time increment between outputs
 x2_slice    = 0.0       # slice in x2
 x3_slice    = 0.0       # slice in x3

Each <output[n]> block will result in a new node being created in a linked list of OutputType stored in the Outputs class. During a simulation, outputs are made when the simulation time satisfies the criteria implemented in the MakeOutputs() function.

在这里，prim关键字会被以下逻辑进行解析， ContainVariable(output_params.variable, "prim")和或逻辑的使用使得，只要字块中有：

prim关键字的存在，就会扩展以下primitive量到Outputs类中并最终输出：rho,press,vel,phi；
cons关键字的存在，就会扩展以下conserved量到Outputs类中并最终输出：dens,Etot,mom,phi。

cpp

void OutputType::LoadOutputData(MeshBlock *pmb) {
   /* ...省略.... */
  // NEW_OUTPUT_TYPES:
  // (lab-frame) density
  if (ContainVariable(output_params.variable, "D") ||
      ContainVariable(output_params.variable, "cons")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "dens";
    pod->data.InitWithShallowSlice(phyd->u, 4, IDN, 1);
    AppendOutputDataNode(pod);
    num_vars_++;
  }

  // (rest-frame) density
  if (ContainVariable(output_params.variable, "d") ||
      ContainVariable(output_params.variable, "prim")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "rho";
    pod->data.InitWithShallowSlice(phyd->w, 4, IDN, 1);
    AppendOutputDataNode(pod);
    num_vars_++;
  }

  // total energy
  if (NON_BAROTROPIC_EOS) {
    if (ContainVariable(output_params.variable, "E") ||
        ContainVariable(output_params.variable, "cons")) {
      pod = new OutputData;
      pod->type = "SCALARS";
      pod->name = "Etot";
      if (porb->orbital_advection_defined
          && !output_params.orbital_system_output) {
        porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::cons);
        pod->data.InitWithShallowSlice(porb->u_orb, 4, IEN, 1);
      } else {
        pod->data.InitWithShallowSlice(phyd->u, 4, IEN, 1);
      }
      AppendOutputDataNode(pod);
      num_vars_++;
    }

    // pressure
    if (ContainVariable(output_params.variable, "p") ||
        ContainVariable(output_params.variable, "prim")) {
      pod = new OutputData;
      pod->type = "SCALARS";
      pod->name = "press";
      pod->data.InitWithShallowSlice(phyd->w, 4, IPR, 1);
      AppendOutputDataNode(pod);
      num_vars_++;
    }
  }

  // momentum vector
  if (ContainVariable(output_params.variable, "m") ||
      ContainVariable(output_params.variable, "cons")) {
    pod = new OutputData;
    pod->type = "VECTORS";
    pod->name = "mom";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::cons);
      pod->data.InitWithShallowSlice(porb->u_orb, 4, IM1, 3);
    } else {
      pod->data.InitWithShallowSlice(phyd->u, 4, IM1, 3);
    }
    AppendOutputDataNode(pod);
    num_vars_ += 3;
    if (output_params.cartesian_vector) {
      AthenaArray<Real> src;
      if (porb->orbital_advection_defined
          && !output_params.orbital_system_output) {
        src.InitWithShallowSlice(porb->u_orb, 4, IM1, 3);
      } else {
        src.InitWithShallowSlice(phyd->u, 4, IM1, 3);
      }
      pod = new OutputData;
      pod->type = "VECTORS";
      pod->name = "mom_xyz";
      pod->data.NewAthenaArray(3, phyd->u.GetDim3(), phyd->u.GetDim2(),
                               phyd->u.GetDim1());
      CalculateCartesianVector(src,  pod->data,  pmb->pcoord);
      AppendOutputDataNode(pod);
      num_vars_ += 3;
    }
  }

  // each component of momentum
  if (ContainVariable(output_params.variable, "m1")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "mom1";
    pod->data.InitWithShallowSlice(phyd->u, 4, IM1, 1);
    AppendOutputDataNode(pod);
    num_vars_++;
  }
  if (ContainVariable(output_params.variable, "m2")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "mom2";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output
        && porb->orbital_direction == 1) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::cons);
      pod->data.InitWithShallowSlice(porb->u_orb, 4, IM2, 1);
    } else {
      pod->data.InitWithShallowSlice(phyd->u, 4, IM2, 1);
    }
    AppendOutputDataNode(pod);
    num_vars_++;
  }
  if (ContainVariable(output_params.variable, "m3")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "mom3";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output
        && porb->orbital_direction == 2) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::cons);
      pod->data.InitWithShallowSlice(porb->u_orb, 4, IM3, 1);
    } else {
      pod->data.InitWithShallowSlice(phyd->u, 4, IM3, 1);
    }
    AppendOutputDataNode(pod);
    num_vars_++;
  }

  // velocity vector
  if (ContainVariable(output_params.variable, "v") ||
      ContainVariable(output_params.variable, "prim")) {
    pod = new OutputData;
    pod->type = "VECTORS";
    pod->name = "vel";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::prim);
      pod->data.InitWithShallowSlice(porb->w_orb, 4, IVX, 3);
    } else {
      pod->data.InitWithShallowSlice(phyd->w, 4, IVX, 3);
    }
    AppendOutputDataNode(pod);
    num_vars_ += 3;
    if (output_params.cartesian_vector) {
      AthenaArray<Real> src;
      if (porb->orbital_advection_defined
          && !output_params.orbital_system_output) {
        src.InitWithShallowSlice(porb->w_orb, 4, IVX, 3);
      } else {
        src.InitWithShallowSlice(phyd->w, 4, IVX, 3);
      }
      pod = new OutputData;
      pod->type = "VECTORS";
      pod->name = "vel_xyz";
      pod->data.NewAthenaArray(3, phyd->w.GetDim3(), phyd->w.GetDim2(),
                               phyd->w.GetDim1());
      CalculateCartesianVector(src,  pod->data,  pmb->pcoord);
      AppendOutputDataNode(pod);
      num_vars_ += 3;
    }
  }

  // each component of velocity
  if (ContainVariable(output_params.variable, "vx") ||
      ContainVariable(output_params.variable, "v1")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "vel1";
    pod->data.InitWithShallowSlice(phyd->w, 4, IVX, 1);
    AppendOutputDataNode(pod);
    num_vars_++;
  }
  if (ContainVariable(output_params.variable, "vy") ||
      ContainVariable(output_params.variable, "v2")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "vel2";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output
        && porb->orbital_direction == 1) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::prim);
      pod->data.InitWithShallowSlice(porb->w_orb, 4, IVY, 1);
    } else {
      pod->data.InitWithShallowSlice(phyd->w, 4, IVY, 1);
    }
    AppendOutputDataNode(pod);
    num_vars_++;
  }
  if (ContainVariable(output_params.variable, "vz") ||
      ContainVariable(output_params.variable, "v3")) {
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "vel3";
    if (porb->orbital_advection_defined
        && !output_params.orbital_system_output
        && porb->orbital_direction == 2) {
      porb->ConvertOrbitalSystem(phyd->w, phyd->u, OrbitalTransform::prim);
      pod->data.InitWithShallowSlice(porb->w_orb, 4, IVZ, 1);
    } else {
      pod->data.InitWithShallowSlice(phyd->w, 4, IVZ, 1);
    }
    AppendOutputDataNode(pod);
    num_vars_++;
  }

  if (SELF_GRAVITY_ENABLED) {
    if (ContainVariable(output_params.variable, "phi") ||
        ContainVariable(output_params.variable, "prim") ||
        ContainVariable(output_params.variable, "cons")) {
      pod = new OutputData;
      pod->type = "SCALARS";
      pod->name = "phi";
      pod->data.InitWithShallowSlice(pgrav->phi, 4, 0, 1);
      AppendOutputDataNode(pod);
      num_vars_++;
      if (pgrav->output_defect) {
        pod = new OutputData;
        pod->type = "SCALARS";
        pod->name = "defect-phi";
        pod->data.InitWithShallowSlice(pgrav->def, 4, 0, 1);
        AppendOutputDataNode(pod);
        num_vars_++;
      }
    }
  } // endif (SELF_GRAVITY_ENABLED)

那么，上述模块的输出示例：

yaml

netcdf juno_mwr.out2.00000 {
dimensions:
	time = UNLIMITED ; // (1 currently)
	x1 = 1600 ;
	x1f = 1601 ;
	x2 = 5 ;
	x2f = 6 ;
	x3 = 1 ;
	ray = 24 ;
variables:
	float time(time) ;
		time:axis = "T" ;
		time:units = "s" ;
		time:long_name = "time" ;
	float x1(x1) ;
		x1:axis = "Z" ;
		x1:units = "m" ;
		x1:long_name = "Z-coordinate at cell center" ;
	float x1f(x1f) ;
		x1f:units = "m" ;
		x1f:long_name = "Z-coordinate at cell face" ;
	float x2(x2) ;
		x2:axis = "Y" ;
		x2:units = "m" ;
		x2:long_name = "Y-coordinate at cell center" ;
	float x2f(x2f) ;
		x2f:units = "m" ;
		x2f:long_name = "Y-coordinate at cell face" ;
	float x3(x3) ;
		x3:axis = "X" ;
		x3:units = "m" ;
		x3:long_name = "X-coordinate at cell center" ;
	float mu_out(ray) ;
		mu_out:units = "1" ;
		mu_out:long_name = "cosine polar angle" ;
	float phi_out(ray) ;
		phi_out:units = "rad" ;
		phi_out:long_name = "azimuthal angle" ;
	float rho(time, x1, x2, x3) ;
		rho:units = "kg/m^3" ;
		rho:long_name = "density" ;
	float press(time, x1, x2, x3) ;
		press:units = "pa" ;
		press:long_name = "pressure" ;
	float vel1(time, x1, x2, x3) ;
		vel1:units = "m/s" ;
		vel1:long_name = "velocity" ;
	float vel2(time, x1, x2, x3) ;
		vel2:units = "m/s" ;
		vel2:long_name = "velocity" ;
	float vel3(time, x1, x2, x3) ;
		vel3:units = "m/s" ;
		vel3:long_name = "velocity" ;
	float r0(time, x1, x2, x3) ;
	float r1(time, x1, x2, x3) ;
	float r2(time, x1, x2, x3) ;
	float r3(time, x1, x2, x3) ;
	float r4(time, x1, x2, x3) ;
	float r5(time, x1, x2, x3) ;
	float vapor1(time, x1, x2, x3) ;
		vapor1:units = "kg/kg" ;
		vapor1:long_name = "mass mixing ratio of vapor" ;
	float vapor2(time, x1, x2, x3) ;
		vapor2:units = "kg/kg" ;
		vapor2:long_name = "mass mixing ratio of vapor" ;
}

Canoe output扩展

Canoe由于添加了热力学和辐射传输模块，对Athena++的输出模块进行了扩展，这些实现在canoe/src/outputs/load_user_output_data.cpp源码中。

包含了对vapor,rad,radtau,radflux,radtoa等关键字的支持。

cpp

void OutputType::loadUserOutputData(MeshBlock *pmb) {
  OutputData *pod;
  auto phyd = pmb->phydro;
  auto prad = pmb->pimpl->prad;

  // vapor
  if (NVAPOR > 0) {
    if (output_params.variable.compare("prim") == 0 ||
        output_params.variable.compare("vapor") == 0) {
      pod = new OutputData;
      pod->type = "VECTORS";
      pod->name = "vapor";
      pod->data.InitWithShallowSlice(phyd->w, 4, 1, NVAPOR);
      AppendOutputDataNode(pod);
      num_vars_ += NVAPOR;
    }

    if (output_params.variable.compare("cons") == 0) {
      pod = new OutputData;
      pod->type = "VECTORS";
      pod->name = "vapor";
      pod->data.InitWithShallowSlice(phyd->u, 4, 1, NVAPOR);
      AppendOutputDataNode(pod);
      num_vars_ += NVAPOR;
    }
  }

  // radiation
  if (output_params.variable.compare("rad") == 0 ||
      output_params.variable.compare("radtau") == 0) {
    for (int b = 0; b < prad->GetNumBands(); ++b) {
      auto pband = prad->GetBand(b);
      // tau
      pod = new OutputData;
      pod->type = "SCALARS";
      pod->name = pband->GetName() + "-tau";
      pod->data.InitWithShallowSlice(pband->btau, 4, 0, 1);
      AppendOutputDataNode(pod);
      num_vars_ += 1;
    }
  }

  if (output_params.variable.compare("rad") == 0 ||
      output_params.variable.compare("radflux") == 0) {
    // flux up and down
    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "flxup";
    pod->data.InitWithShallowSlice(prad->flxup, 4, 0, 1);
    AppendOutputDataNode(pod);
    num_vars_ += 1;

    pod = new OutputData;
    pod->type = "SCALARS";
    pod->name = "flxdn";
    pod->data.InitWithShallowSlice(prad->flxdn, 4, 0, 1);
    AppendOutputDataNode(pod);
    num_vars_ += 1;
  }

  if (output_params.variable.compare("rad") == 0 ||
      output_params.variable.compare("radtoa") == 0) {
    if (prad->radiance.GetDim3() > 0) {
      pod = new OutputData;
      pod->type = "SCALARS";
      pod->name = "radiance";
      pod->data.InitWithShallowSlice(prad->radiance, 4, 0, 1);
      AppendOutputDataNode(pod);
      num_vars_ += 1;

      // for (auto p : prad->bands) p->writeBinRadiance(&output_params);
    }
  }
  /* ..省略.. */
}

按照和中的设定，输出示例如下：

yaml

netcdf juno_mwr.out4.00000 {
dimensions:
	time = UNLIMITED ; // (1 currently)
	x1 = 1600 ;
	x1f = 1601 ;
	x2 = 5 ;
	x2f = 6 ;
	x3 = 1 ;
	ray = 24 ;
variables:
	float time(time) ;
		time:axis = "T" ;
		time:units = "s" ;
		time:long_name = "time" ;
	float x1(x1) ;
		x1:axis = "Z" ;
		x1:units = "m" ;
		x1:long_name = "Z-coordinate at cell center" ;
	float x1f(x1f) ;
		x1f:units = "m" ;
		x1f:long_name = "Z-coordinate at cell face" ;
	float x2(x2) ;
		x2:axis = "Y" ;
		x2:units = "m" ;
		x2:long_name = "Y-coordinate at cell center" ;
	float x2f(x2f) ;
		x2f:units = "m" ;
		x2f:long_name = "Y-coordinate at cell face" ;
	float x3(x3) ;
		x3:axis = "X" ;
		x3:units = "m" ;
		x3:long_name = "X-coordinate at cell center" ;
	float mu_out(ray) ;
		mu_out:units = "1" ;
		mu_out:long_name = "cosine polar angle" ;
	float phi_out(ray) ;
		phi_out:units = "rad" ;
		phi_out:long_name = "azimuthal angle" ;
	float radiance(time, ray, x2, x3) ;
		radiance:units = "K" ;
		radiance:long_name = "top-of-atmosphere radiance" ;
}

netcdf juno_mwr.out5.00000 {
dimensions:
	time = UNLIMITED ; // (1 currently)
	x1 = 1600 ;
	x1f = 1601 ;
	x2 = 5 ;
	x2f = 6 ;
	x3 = 1 ;
	ray = 24 ;
variables:
	float time(time) ;
		time:axis = "T" ;
		time:units = "s" ;
		time:long_name = "time" ;
	float x1(x1) ;
		x1:axis = "Z" ;
		x1:units = "m" ;
		x1:long_name = "Z-coordinate at cell center" ;
	float x1f(x1f) ;
		x1f:units = "m" ;
		x1f:long_name = "Z-coordinate at cell face" ;
	float x2(x2) ;
		x2:axis = "Y" ;
		x2:units = "m" ;
		x2:long_name = "Y-coordinate at cell center" ;
	float x2f(x2f) ;
		x2f:units = "m" ;
		x2f:long_name = "Y-coordinate at cell face" ;
	float x3(x3) ;
		x3:axis = "X" ;
		x3:units = "m" ;
		x3:long_name = "X-coordinate at cell center" ;
	float mu_out(ray) ;
		mu_out:units = "1" ;
		mu_out:long_name = "cosine polar angle" ;
	float phi_out(ray) ;
		phi_out:units = "rad" ;
		phi_out:long_name = "azimuthal angle" ;
	float CH1tau(time, x1, x2, x3) ;
		CH1tau:units = "1" ;
		CH1tau:long_name = "optical thickness" ;
	float CH2tau(time, x1, x2, x3) ;
		CH2tau:units = "1" ;
		CH2tau:long_name = "optical thickness" ;
	float CH3tau(time, x1, x2, x3) ;
		CH3tau:units = "1" ;
		CH3tau:long_name = "optical thickness" ;
	float CH4tau(time, x1, x2, x3) ;
		CH4tau:units = "1" ;
		CH4tau:long_name = "optical thickness" ;
	float CH5tau(time, x1, x2, x3) ;
		CH5tau:units = "1" ;
		CH5tau:long_name = "optical thickness" ;
	float CH6tau(time, x1, x2, x3) ;
		CH6tau:units = "1" ;
		CH6tau:long_name = "optical thickness" ;
}

这些量的meta内容设定在 canoe/src/outputs/output_utils.cpp中

plaintext

  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"}};
}

User-defined Output Variable (UOV)

在特定的模拟中，我们可能会需要输出我们关心的诊断量，但是往往其不再athena++或者Canoe的输出序列中，那么我们就需要使用UOV功能来注册和装填我们需要的诊断输出量了。
典型的例子有，温度，位温，虚温，湿静能等参数，不属于流体力学的诊断量，或者热力学模块的参数。
<output3>就是对UOV变量的输出设定，使用了uov的关键字来表示我们有一些自定义的诊断量需要输出到output3.nc这个文件，步长为dt。
那么，Canoe(Athena++)通过覆盖MeshBlock类中的成员方法来进行UOV的注册：

cpp

class MeshBlock{
  //! defined in either the prob file or default_pgen.cpp in ../pgen/
  void UserWorkBeforeOutput(ParameterInput *pin); // called in Mesh fn (friend class)
  //! defined in either the prob file or default_pgen.cpp in ../pgen/
  void ProblemGenerator(ParameterInput *pin);
  void InitUserMeshBlockData(ParameterInput *pin);
}

可以理解为：

MeshBlock::InitUserMeshBlockData：可以实现对UOV的注册；
MeshBlock::UserWorkBeforeOutput：计算输出量；
MeshBlock::ProblemGenerator：设置初始值条件。

当然了，这几个函数的功能不仅限与这些。
default_pgen.cpp中没有对这三个函数做具体实现：

cpp

// 4x members of MeshBlock class:

//========================================================================================
//! \fn void MeshBlock::InitUserMeshBlockData(ParameterInput *pin)
//! \brief Function to initialize problem-specific data in MeshBlock class.  Can also be
//! used to initialize variables which are global to other functions in this file.
//! Called in MeshBlock constructor before ProblemGenerator.
//========================================================================================

void __attribute__((weak)) MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  // do nothing
  return;
}

//========================================================================================
//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
//! \brief Should be used to set initial conditions.
//========================================================================================

void __attribute__((weak)) MeshBlock::ProblemGenerator(ParameterInput *pin) {
  // In practice, this function should *always* be replaced by a version
  // that sets the initial conditions for the problem of interest.
  return;
}

//===========================================================================
//! \fn void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin)
//! \brief Function called before generating output files
//========================================================================================

void __attribute__((weak)) MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
  // do nothing
  return;
}

在模拟实际问题时，应该按需求对这三个成员函数进行实现，也就是本专栏中几个example主要实现的代码。

例如，Juno_MWR模拟中，UOV的设置和计算如下：

cpp

// 注册输出量
void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  AllocateUserOutputVariables(4 + NVAPOR);
  SetUserOutputVariableName(0, "temp");
  SetUserOutputVariableName(1, "theta");
  SetUserOutputVariableName(2, "thetav");
  SetUserOutputVariableName(3, "mse");
  for (int n = 1; n <= NVAPOR; ++n) {
    std::string name = "rh" + std::to_string(n);
    SetUserOutputVariableName(3 + n, name.c_str());
  }
}
//计算输出
void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
  auto pthermo = Thermodynamics::GetInstance();
  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j)
      for (int i = is; i <= ie; ++i) {
        user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
        user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i);
        // theta_v
        user_out_var(2, k, j, i) =
            user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
        // mse
        user_out_var(3, k, j, i) =
            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
        for (int n = 1; n <= NVAPOR; ++n)
          user_out_var(3 + n, k, j, i) =
              pthermo->RelativeHumidity(this, n, k, j, i);
      }
}

最终的输出文件meta如下：

yaml

netcdf juno_mwr.out3.00000 {
dimensions:
	time = UNLIMITED ; // (1 currently)
	x1 = 1600 ;
	x1f = 1601 ;
	x2 = 5 ;
	x2f = 6 ;
	x3 = 1 ;
	ray = 24 ;
variables:
	float time(time) ;
		time:axis = "T" ;
		time:units = "s" ;
		time:long_name = "time" ;
	float x1(x1) ;
		x1:axis = "Z" ;
		x1:units = "m" ;
		x1:long_name = "Z-coordinate at cell center" ;
	float x1f(x1f) ;
		x1f:units = "m" ;
		x1f:long_name = "Z-coordinate at cell face" ;
	float x2(x2) ;
		x2:axis = "Y" ;
		x2:units = "m" ;
		x2:long_name = "Y-coordinate at cell center" ;
	float x2f(x2f) ;
		x2f:units = "m" ;
		x2f:long_name = "Y-coordinate at cell face" ;
	float x3(x3) ;
		x3:axis = "X" ;
		x3:units = "m" ;
		x3:long_name = "X-coordinate at cell center" ;
	float mu_out(ray) ;
		mu_out:units = "1" ;
		mu_out:long_name = "cosine polar angle" ;
	float phi_out(ray) ;
		phi_out:units = "rad" ;
		phi_out:long_name = "azimuthal angle" ;
	float temp(time, x1, x2, x3) ;
		temp:units = "K" ;
		temp:long_name = "temperature" ;
	float theta(time, x1, x2, x3) ;
		theta:units = "K" ;
		theta:long_name = "potential temperature" ;
	float thetav(time, x1, x2, x3) ;
		thetav:units = "K" ;
		thetav:long_name = "virtual potential temperature" ;
	float mse(time, x1, x2, x3) ;
		mse:units = "J/kg" ;
		mse:long_name = "moist static energy" ;
	float rh1(time, x1, x2, x3) ;
		rh1:units = "1" ;
		rh1:long_name = "relative humidity 1" ;
	float rh2(time, x1, x2, x3) ;
		rh2:units = "1" ;
		rh2:long_name = "relative humidity 2" ;
}