C++ NcVar::put方法代码示例

void PointCounter::output(const string &fileName) const
{
    NcFile file(fileName.c_str(), NcFile::Replace);
    if (!file.is_valid()) {
        string message = "Failed to open file "+fileName+".";
        REPORT_ERROR(message.c_str());
    }

    NcDim *numLonDim = file.add_dim("lon", counters.extent(0));
    NcDim *numLatDim = file.add_dim("lat", counters.extent(1));
    NcDim *numBndsDim = file.add_dim("bnds", 2);
    //NcDim *numLevDim = file.add_dim("lev", counters.extent(2));

    NcVar *lonBndsVar = file.add_var("lon_bnds", ncDouble, numLonDim, numBndsDim);
    NcVar *latBndsVar = file.add_var("lat_bnds", ncDouble, numLatDim, numBndsDim);
    //NcVar *levVar = file.add_var("lev", ncDouble, numLevDim);

    lonBndsVar->add_att("long_name", "longitude bounds");
    lonBndsVar->add_att("units", "degree_east");
    latBndsVar->add_att("long_name", "latitude bounds");
    latBndsVar->add_att("units", "degree_north");

    int numLonBnds = mesh[0].getNumLon()-2, numLatBnds = mesh[0].getNumLat()-1;
    double lonBnds[numLonBnds][2], latBnds[numLatBnds][2];
    for (int i = 0; i < numLonBnds; ++i) {
        lonBnds[i][0] = mesh[0].lon(i)*Rad2Deg;
        lonBnds[i][1] = mesh[0].lon(i+1)*Rad2Deg;
    }
    for (int j = 0; j < numLatBnds; ++j) {
        latBnds[j][0] = mesh[0].lat(j)*Rad2Deg;
        latBnds[j][1] = mesh[0].lat(j+1)*Rad2Deg;
    }

    lonBndsVar->put(&lonBnds[0][0], numLonBnds, 2);
    latBndsVar->put(&latBnds[0][0], numLatBnds, 2);

    NcVar *countersVar = file.add_var("counters", ncInt, numLatDim, numLonDim);

    countersVar->add_att("long_name", "contained point numbers in grid boxes");

    int counters[this->counters.extent(1)][this->counters.extent(0)];
    for (int j = 0; j < this->counters.extent(1); ++j)
        for (int i = 0; i < this->counters.extent(0); ++i)
            counters[j][i] = this->counters(i, j, 0);
    countersVar->put(&counters[0][0], this->counters.extent(1), this->counters.extent(0));
}

void RLLMesh::output(const string &fileName) const
{
    NcFile file(fileName.c_str(), NcFile::Replace);
    if (!file.is_valid()) {
        Message message;
        message << "Failed to create point counter mesh file \"";
        message << fileName << "\"!";
        REPORT_ERROR(message.str());
    }
    NcDim *lonDim = file.add_dim("lon", lon.size());
    NcDim *latDim = file.add_dim("lat", lat.size());
    NcVar *lonVar = file.add_var("lon", ncDouble, lonDim);
    NcVar *latVar = file.add_var("lat", ncDouble, latDim);
    double lonTmp[lon.size()], latTmp[lat.size()];
    for (int i = 0; i < lon.size(); ++i)
        lonTmp[i] = lon(i)*Rad2Deg;
    for (int j = 0; j < lat.size(); ++j)
        latTmp[j] = lat(j)*Rad2Deg;
    lonVar->put(lonTmp, lon.size());
    latVar->put(latTmp, lat.size());
    file.close();
}

int 
main(void)
@{
   // This is the data array we will write. It will just be filled
   // with a progression of numbers for this example.
   int dataOut[NX][NY];
  
   // Create some pretend data. If this wasn't an example program, we
   // would have some real data to write, for example, model output.
   for(int i = 0; i < NX; i++) 
      for(int j = 0; j < NY; j++)
	 dataOut[i][j] = i * NY + j;

   // Create the file. The Replace parameter tells netCDF to overwrite
   // this file, if it already exists.
   NcFile dataFile("simple_xy.nc", NcFile::Replace);

   // You should always check whether a netCDF file creation or open
   // constructor succeeded.
   if (!dataFile.is_valid())
   @{
      cout << "Couldn't open file!\n";
      return NC_ERR;
   @}
  
   // For other method calls, the default behavior of the C++ API is
   // to exit with a message if there is an error.  If that behavior
   // is OK, there is no need to check return values in simple cases
   // like the following.

   // When we create netCDF dimensions, we get back a pointer to an
   // NcDim for each one.
   NcDim* xDim = dataFile.add_dim("x", NX);
   NcDim* yDim = dataFile.add_dim("y", NY);
  
   // Define a netCDF variable. The type of the variable in this case
   // is ncInt (32-bit integer).
   NcVar *data = dataFile.add_var("data", ncInt, xDim, yDim);
     
   // Write the pretend data to the file. Although netCDF supports
   // reading and writing subsets of data, in this case we write all
   // the data in one operation.
   data->put(&dataOut[0][0], NX, NY);

   // The file will be automatically close when the NcFile object goes
   // out of scope. This frees up any internal netCDF resources
   // associated with the file, and flushes any buffers.
   cout << "*** SUCCESS writing example file simple_xy.nc!" << endl;

   return 0;
@}

int TestSuite::testVar()
{


   try
   {

      string FILE_NAME = "tst_vars.nc";
      int NDIMS = 4;
      int NLAT = 6;
      int NLON = 12;

      // Names of things. 
      string LAT_NAME = "latitude";
      string LON_NAME = "longitude";

      int  MAX_ATT_LEN = 80;
      // These are used to construct some example data. 
      float START_LAT = 25.0;
      float START_LON = -125.0;

      string  UNITS = "units";
      string  DEGREES_EAST =  "degrees_east";
      string  DEGREES_NORTH = "degrees_north";

      // For the units attributes. 
      string LAT_UNITS = "degrees_north";
      string LON_UNITS = "degrees_east";

      // Return this code to the OS in case of failure.
#define NC_ERR 2

  
      // We will write latitude and longitude fields. 
      float lats[NLAT],lons[NLON];

      // create some pretend data. If this wasn't an example program, we
      // would have some real data to write for example, model output.
      for (int lat = 0; lat < NLAT; lat++)
	 lats[lat] = START_LAT + 5. * lat;
      for (int lon = 0; lon < NLON; lon++)
	 lons[lon] = START_LON + 5. * lon;

      // Create the file.
      NcFile test(FILE_NAME, NcFile::Replace);

      // Define the dimensions. NetCDF will hand back an ncDim object for
      // each.
      NcDim* latDim = test.addDim(LAT_NAME, NLAT);
      NcDim* lonDim = test.addDim(LON_NAME, NLON);

   
      // Define the coordinate variables.
      NcVar* latVar = test.addVar(LAT_NAME, ncFloat, latDim);
      NcVar* lonVar = test.addVar(LON_NAME, ncFloat, lonDim);
       
      // Define units attributes for coordinate vars. This attaches a
      // text attribute to each of the coordinate variables, containing
      // the units.
      latVar->addAtt(UNITS,ncString, DEGREES_NORTH);
      lonVar->addAtt(UNITS,ncString, DEGREES_EAST);

      // Write the coordinate variable data to the file.
      latVar->put(lats, NLAT);
      lonVar->put(lons, NLON);

      NcValues *latVals = latVar->getValues();
      cout<<"toString returns lats: "<<latVals->toString()<<endl;
      cout<<"toChar returns "<<latVals->toChar(1)<<endl;
      cout<<"toShort returns "<<latVals->toShort(1)<<endl;
      cout<<"toInt returns "<<latVals->toInt(1)<<endl;
      cout<<"toLong returns "<<latVals->toLong(1)<<endl;

      latVals->print(cout);
      
      NcValues *lonVals = lonVar->getValues();
      cout<<"toString returns lats: "<<lonVals->toString()<<endl;
      lonVals->print(cout);
      
	
      cout<<"no segmentation fault thus far"<<endl;
 

      //test varaibles here
   }
   catch(NcException e)
   {
      e.what();
      return 1;
   }
   try
   {
      cout<<"should test adding a variable with more than 5 dimensions here"<<endl;
      // test creating a variable with more than 5 dimensions
   } 
   catch (NcException e)
   {
      e.what();
      return 1;
   }
//.........这里部分代码省略.........

int TestSuite::testFile(string fName, NcFile::FileMode fMode)
{
   string UNITS = "UNITS:";
  
   cout<<"*** Testing tst_file in tst_suite ";
   int NLAT = 6;
   int NLON = 12;
   int NDIMS = 4;
   int NLVL = 2;
   int NREC = 2;
  
   /* These are used to construct some example data. */
   int SAMPLE_PRESSURE = 900;
   float SAMPLE_TEMP = 9.0;
   float START_LAT = 25.0;
   float START_LON = -125.0;
  
  
   /* We will write surface temperature and pressure fields. */
   float pres_out[NLAT][NLON];
   float pres_in[NLAT][NLON];
   float temp_out[NLAT][NLON];
   float temp_in[NLAT][NLON];
   float lats[NLAT], lons[NLON],lats_in[NLAT],lons_in[NLON];
   std::string  chararray []={"I"," hope"," this"," is"," stored "," properly" };
   std::string chararray_in[NLAT];
   int outInts[NLAT],outInts_in[NLAT];
      
   /* It's good practice for each netCDF variable to carry a "units"
    * attribute. */
   char pres_units[] = "hPa";
   char temp_units[] = "celsius";
      
   /* Loop indexes. */
   int lat, lon;
      
   /* Create some pretend data. If this wasn't an example program, we
    * would have some real data to write, for example, model
    * output. */
   for (lat = 0; lat < NLAT; lat++)
   {   
      lats[lat] = START_LAT + 5.*lat;
   }
   for(lat = 0; lat < NLAT; lat++)
   {  
      outInts[lat]= 450;
   }
   for (lon = 0; lon < NLON; lon++)
   {
      lons[lon] = START_LON + 5.*lon;
   }
   for (lat = 0; lat < NLAT; lat++)
   {
      for (lon = 0; lon < NLON; lon++)     
      {
	 pres_out[lat][lon] = SAMPLE_PRESSURE + (lon * NLAT + lat); 
	 temp_out[lat][lon] = SAMPLE_TEMP + .25 * (lon * NLAT + lat);
      }
   }
      
   // nc_set_log_level(3);
   try
   { 
      NcFile f("tst_file.nc",NcFile::Replace);
      NcGroup *root = f.getRootGroup();
      
      NcDim *latDim = root->addDim(string("lat"),NLAT);
      
      NcDim *lonDim = root->addDim(string("lon"),NLON);
      NcVar *latVar = root->addVar(string("latVar"),ncDouble,latDim);
 
      NcVar *lonVar = root->addVar(string("lonVar"),ncFloat,lonDim);
      NcVar *outIntsVar = root->addVar(string("outintsVar"),ncInt,latDim);
      NcVar *charArrVar = root->addVar(string("CharArray"),ncString,latDim);
      
      latVar->addAtt(string(UNITS),ncChar,string("degrees_north"));
      lonVar->addAtt(string(UNITS),ncChar,string("degrees_south"));
      
      outIntsVar->put(&outInts[0],NLAT,0,0,0,0); 
      charArrVar->put(&chararray[0],NLAT,0,0,0,0);
      
      latVar->put(&lats[0],NLAT,0,0,0,0);
      lonVar->put(&lons[0],NLON,0,0,0,0);
      
      NcVar *presVar = root->addVar(string("press"),ncFloat,latDim,lonDim);
      NcVar *tempVar = root->addVar(string("temp"),ncFloat,latDim,lonDim);
      presVar->addAtt(string("UNITS:"),ncChar,string(pres_units));
      tempVar->addAtt(string("UNITS:"),ncChar,string(temp_units));
      
      presVar->put(&pres_out[0][0],NLAT,NLON,0,0,0);
      tempVar->put(&temp_out[0][0],NLAT,NLON,0,0,0);

      //NcValues *ncvalues = presVar->getValues();
      
      {  //another scope for variables 
	 NcGroup::varIterator variableItr;
	 variableItr = root->beginVar();
	 while(variableItr != root->endVar())
	 {
	    variableItr++;
//.........这里部分代码省略.........

void PolygonManager::output(const string &fileName)
{
    // -------------------------------------------------------------------------
    NcFile file(fileName.c_str(), NcFile::Replace);
    if (!file.is_valid()) {
        string message = "Failed to open file "+fileName+".";
        REPORT_ERROR(message.c_str());
    }
    // -------------------------------------------------------------------------
    // time information
    if (TimeManager::onLine()) {
        file.add_att("time", TimeManager::getSeconds());
        file.add_att("time_step", TimeManager::getTimeStep());
        file.add_att("steps", TimeManager::getSteps());
    }
    // -------------------------------------------------------------------------
    // dimensions
    NcDim *numVertexDim = file.add_dim("num_total_vertex", vertices.size());
    NcDim *numEdgeDim = file.add_dim("num_total_edge", edges.size());
    NcDim *numPolygonDim = file.add_dim("num_total_polygon", polygons.size());
    // -------------------------------------------------------------------------
    // vertices part
    NcVar *oldVtxLonVar = file.add_var("old_vertex_lon", ncDouble, numVertexDim);
    NcVar *oldVtxLatVar = file.add_var("old_vertex_lat", ncDouble, numVertexDim);
    NcVar *newVtxLonVar = file.add_var("new_vertex_lon", ncDouble, numVertexDim);
    NcVar *newVtxLatVar = file.add_var("new_vertex_lat", ncDouble, numVertexDim);
    oldVtxLonVar->add_att("long_name", "old vertex longitude");
    oldVtxLonVar->add_att("units", "degree_east");
    oldVtxLatVar->add_att("long_name", "old vertex latitude");
    oldVtxLatVar->add_att("units", "degree_north");
    newVtxLonVar->add_att("long_name", "new vertex longitude");
    newVtxLonVar->add_att("units", "degree_east");
    newVtxLatVar->add_att("long_name", "new vertex latitude");
    newVtxLatVar->add_att("units", "degree_north");
    double *oldVtxLon = new double[vertices.size()];
    double *oldVtxLat = new double[vertices.size()];
    double *newVtxLon = new double[vertices.size()];
    double *newVtxLat = new double[vertices.size()];
    Vertex *vertex = vertices.front();
    for (int i = 0; i < vertices.size(); ++i) {
        oldVtxLon[i] = vertex->getCoordinate(OldTimeLevel).getLon()*Rad2Deg;
        oldVtxLat[i] = vertex->getCoordinate(OldTimeLevel).getLat()*Rad2Deg;
        newVtxLon[i] = vertex->getCoordinate(NewTimeLevel).getLon()*Rad2Deg;
        newVtxLat[i] = vertex->getCoordinate(NewTimeLevel).getLat()*Rad2Deg;
        vertex = vertex->next;
    }
    oldVtxLonVar->put(oldVtxLon, vertices.size());
    oldVtxLatVar->put(oldVtxLat, vertices.size());
    newVtxLonVar->put(newVtxLon, vertices.size());
    newVtxLatVar->put(newVtxLat, vertices.size());
    delete [] oldVtxLon;
    delete [] oldVtxLat;
    delete [] newVtxLon;
    delete [] newVtxLat;
    // -------------------------------------------------------------------------
    // edges part
    NcVar *firstPointVar = file.add_var("first_point_idx", ncInt, numEdgeDim);
    NcVar *secondPointVar = file.add_var("second_point_idx", ncInt, numEdgeDim);
    firstPointVar->add_att("long_name", "first point index of edge");
    secondPointVar->add_att("long_name", "second point index of edge");
    int *idx1 = new int[edges.size()];
    int *idx2 = new int[edges.size()];
    Edge *edge = edges.front();
    for (int i = 0; i < edges.size(); ++i) {
        idx1[i] = edge->getEndPoint(FirstPoint)->getID();
        idx2[i] = edge->getEndPoint(SecondPoint)->getID();
        edge = edge->next;
    }
    firstPointVar->put(idx1, edges.size());
    secondPointVar->put(idx2, edges.size());
    delete [] idx1;
    delete [] idx2;
    // test point
    NcVar *oldTestLonVar = file.add_var("old_testpoint_lon", ncDouble, numEdgeDim);
    NcVar *oldTestLatVar = file.add_var("old_testpoint_lat", ncDouble, numEdgeDim);
    NcVar *newTestLonVar = file.add_var("new_testpoint_lon", ncDouble, numEdgeDim);
    NcVar *newTestLatVar = file.add_var("new_testpoint_lat", ncDouble, numEdgeDim);
    oldTestLonVar->add_att("long_name", "old test point longitude");
    oldTestLonVar->add_att("units", "degree_east");
    oldTestLatVar->add_att("long_name", "old test point latitude");
    oldTestLatVar->add_att("units", "degree_north");
    newTestLonVar->add_att("long_name", "new test point longitude");
    newTestLonVar->add_att("units", "degree_east");
    newTestLatVar->add_att("long_name", "new test point latitude");
    newTestLatVar->add_att("units", "degree_north");
#if defined TTS_ONLINE || PREPROCESS
    double *oldTestLon = new double[edges.size()];
    double *oldTestLat = new double[edges.size()];
    double *newTestLon = new double[edges.size()];
    double *newTestLat = new double[edges.size()];
    edge = edges.front();
    for (int i = 0; i < edges.size(); ++i) {
        Vertex *testPoint = edge->getTestPoint();
        oldTestLon[i] = testPoint->getCoordinate(OldTimeLevel).getLon()*Rad2Deg;
        oldTestLat[i] = testPoint->getCoordinate(OldTimeLevel).getLat()*Rad2Deg;
        newTestLon[i] = testPoint->getCoordinate(NewTimeLevel).getLon()*Rad2Deg;
        newTestLat[i] = testPoint->getCoordinate(NewTimeLevel).getLat()*Rad2Deg;
        edge = edge->next;
    }
    oldTestLonVar->put(oldTestLon, edges.size());
//.........这里部分代码省略.........

int main(int argc, char** argv) {

	NcError error(NcError::verbose_nonfatal);

try {

	// Input file
	std::string strInputFile;

	// Input file list
	std::string strInputFileList;

	// Input file format
	std::string strInputFormat;

	// NetCDF file containing latitude and longitude arrays
	std::string strLatLonFile;

	// Output file (NetCDF)
	std::string strOutputFile;

	// Output variable name
	std::string strOutputVariable;

	// Column in which the longitude index appears
	int iLonIxCol;

	// Column in which the latitude index appears
	int iLatIxCol;

	// Begin latitude
	double dLatBegin;

	// End latitude
	double dLatEnd;

	// Begin longitude
	double dLonBegin;

	// End longitude
	double dLonEnd;

	// Number of latitudes in output
	int nLat;

	// Number of longitudes in output
	int nLon;

	// Parse the command line
	BeginCommandLine()
		CommandLineString(strInputFile, "in", "");
		CommandLineString(strInputFileList, "inlist", "");
		CommandLineStringD(strInputFormat, "in_format", "std", "(std|visit)");
		CommandLineString(strOutputFile, "out", "");
		CommandLineString(strOutputVariable, "outvar", "density");
		CommandLineInt(iLonIxCol, "iloncol", 8);
		CommandLineInt(iLatIxCol, "ilatcol", 9);
		CommandLineDouble(dLatBegin, "lat_begin", -90.0);
		CommandLineDouble(dLatEnd, "lat_end", 90.0);
		CommandLineDouble(dLonBegin, "lon_begin", 0.0);
		CommandLineDouble(dLonEnd, "lon_end", 360.0);
		CommandLineInt(nLat, "nlat", 180);
		CommandLineInt(nLon, "nlon", 360);

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	AnnounceBanner();

	// Check input
	if ((strInputFile == "") && (strInputFileList == "")) {
		_EXCEPTIONT("No input file (--in) or (--inlist) specified");
	}
	if ((strInputFile != "") && (strInputFileList != "")) {
		_EXCEPTIONT("Only one input file (--in) or (--inlist) allowed");
	}
	if (strInputFormat != "std") {
		_EXCEPTIONT("UNIMPLEMENTED:  Only \"--in_format std\" supported");
	}

	// Check output
	if (strOutputFile == "") {
		_EXCEPTIONT("No output file (--out) specified");
	}

	// Check output variable
	if (strOutputVariable == "") {
		_EXCEPTIONT("No output variable name (--outvar) specified");
	}

	// Number of latitudes and longitudes
	if (nLat == 0) {
		_EXCEPTIONT("UNIMPLEMENTED: --nlat must be specified currently");
	}
	if (nLon == 0) {
		_EXCEPTIONT("UNIMPLEMENTED: --nlon must be specified currently");
	}

	// Input file list
	std::vector<std::string> vecInputFiles;
//.........这里部分代码省略.........

void gen(const char* path, NcFile::FileFormat format)		// Generate a netCDF file
{

    NcFile nc(path, NcFile::Replace, NULL, 0, format); // Create, leave in define mode

    // Check if the file was opened successfully
    if (! nc.is_valid()) {
	cerr << "can't create netCDF file " << path << "\n";
	return;
    }

    // Create dimensions
    const int NLATS = 4;
    const int NLONS = 3;
    const int NFRTIMES = 2;
    const int TIMESTRINGLEN = 20;
    NcDim* latd = nc.add_dim("lat", NLATS);
    NcDim* lond = nc.add_dim("lon", NLONS);
    NcDim* frtimed = nc.add_dim("frtime"); // unlimited dimension
    NcDim* timelend = nc.add_dim("timelen", TIMESTRINGLEN); 

    // Create variables and their attributes
    NcVar* P = nc.add_var("P", ncFloat, frtimed, latd, lond);
    P->add_att("long_name", "pressure at maximum wind");
    P->add_att("units", "hectopascals");
    static float range[] = {0., 1500.};
    P->add_att("valid_range", 2, range);
    P->add_att("_FillValue", -9999.0f);

    NcVar* lat = nc.add_var("lat", ncFloat, latd);
    lat->add_att("long_name", "latitude");
    lat->add_att("units", "degrees_north");

    NcVar* lon = nc.add_var("lon", ncFloat, lond);
    lon->add_att("long_name", "longitude");
    lon->add_att("units", "degrees_east");

    NcVar* frtime = nc.add_var("frtime", ncLong, frtimed);
    frtime->add_att("long_name", "forecast time");
    frtime->add_att("units", "hours");

    NcVar* reftime = nc.add_var("reftime",ncChar,timelend);
    reftime->add_att("long_name", "reference time");
    reftime->add_att("units", "text_time");

    NcVar* scalar = nc.add_var("scalarv", ncInt);
    scalar->add_att("scalar_att", 1);

    // Global attributes
    nc.add_att("history", "created by Unidata LDM from NPS broadcast");
    nc.add_att("title", "NMC Global Product Set: Pressure at Maximum Wind");

    // Start writing data, implictly leaves define mode

    static float lats[NLATS] = {-90, -87.5, -85, -82.5};
    lat->put(lats, NLATS);

    static float lons[NLONS] = {-180, -175, -170};
    lon->put(lons, NLONS);

    static int frtimes[NFRTIMES] = {12, 18};
    frtime->put(frtimes, NFRTIMES);

    static const char* s = "1992-3-21 12:00" ;
    reftime->put(s, strlen(s));

    static float P_data[2][4][3] = {
	{{950, 951, 952}, {953, 954, 955}, {956, 957, 958}, {959, 960, 961}},
	{{962, 963, 964}, {965, 966, 967}, {968, 969, 970}, {971, 972, 973}}
      };
    // We could write all P data at once with P->put(&P_data[0][0][0], P->edges()),
    // but instead we write one record at a time, to show use of setcur().
    long rec = 0;                                      // start at zero-th
    const long nrecs = 1;		               // # records to write
    P->put(&P_data[0][0][0], nrecs, NLATS, NLONS);           // write zero-th record
    P->set_cur(++rec);		                       // set to next record
    P->put(&P_data[1][0][0], nrecs, NLATS, NLONS); // write next record

    // close of nc takes place in destructor
}

int main(int argc, char ** argv) {

	MPI_Init(&argc, &argv);

try {

	// Number of latitudes
	int nLat;

	// Number of longitudes
	int nLon;

	// Zonal wavenumber
	int nK;

	// Meridional power
	int nLpow;

	// Output file
	std::string strOutputFile;

	// Parse the command line
	BeginCommandLine()
		CommandLineInt(nLat, "lat", 40);
		CommandLineInt(nLon, "lon", 80);
		CommandLineString(strOutputFile, "out", "topo.nc");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	// Generate longitude and latitude arrays
	AnnounceBanner();
	AnnounceStartBlock("Generating longitude and latitude arrays");

	DataVector<double> dLon;
	dLon.Initialize(nLon);

	Parameters param;
	param.GenerateLatituteArray(nLat);

	std::vector<double> & dLat = param.vecNode;

	double dDeltaLon = 2.0 * M_PI / static_cast<double>(nLon);
	for (int i = 0; i < nLon; i++) {
		dLon[i] = (static_cast<double>(i) + 0.5) * dDeltaLon;
	}

	AnnounceEndBlock("Done");

	// Open NetCDF output file
	AnnounceStartBlock("Writing to file");

	NcFile ncdf_out(strOutputFile.c_str(), NcFile::Replace);

	// Output coordinates
	NcDim * dimLat = ncdf_out.add_dim("lat", nLat);
	NcDim * dimLon = ncdf_out.add_dim("lon", nLon);

	NcVar * varLat = ncdf_out.add_var("lat", ncDouble, dimLat);
	varLat->set_cur((long)0);
	varLat->put(&(param.vecNode[0]), nLat);

	NcVar * varLon = ncdf_out.add_var("lon", ncDouble, dimLon);
	varLon->set_cur((long)0);
	varLon->put(&(dLon[0]), nLon);

	// Generate topography
	DataMatrix<double> dTopo;
	dTopo.Initialize(nLat, nLon);

	double dK = static_cast<double>(nK);
	double dLpow = static_cast<double>(nLpow);

	double dA = 6.37122e6;
	double dX = 500.0;
	double dLatM = 0.0;
	double dLonM = M_PI / 4.0;
	double dD = 5000.0;
	double dH0 = 1.0;
	double dXiM = 4000.0;

	for (int j = 0; j < nLat; j++) {
	for (int i = 0; i < nLon; i++) {

		// Great circle distance
		double dR = dA / dX * acos(sin(dLatM) * sin(dLat[j])
			+ cos(dLatM) * cos(dLat[j]) * cos(dLon[i] - dLonM));

		double dCosXi = 1.0; //cos(M_PI * dR / dXiM);

		dTopo[j][i] = dH0 * exp(- dR * dR / (dD * dD))
			* dCosXi * dCosXi;
	}
	}

	// Write topography
	NcVar * varZs = ncdf_out.add_var("Zs", ncDouble, dimLat, dimLon);
	varZs->set_cur(0, 0);
	varZs->put(&(dTopo[0][0]), nLat, nLon);

//.........这里部分代码省略.........

int main(int argc, char** argv)
{
    if (!cmdline(argc, argv))
    {
        printhelp();
        return EXIT_FAILURE;
    }

    NcFile infile(infilename.c_str(), NcFile::ReadOnly);
    if (!infile.is_valid())
    {
        std::cerr << "Error: invalid input file -- '" << infilename << "'" << std::endl;
        infile.close();
        return EXIT_FAILURE;
    }

    NcFile outfile(outfilename.c_str(), NcFile::Replace);
    if (!outfile.is_valid())
    {
        std::cerr << "Error: cannot open output file -- '" << outfilename << "'" << std::endl;
        outfile.close();
        return EXIT_FAILURE;
    }

    if (varstrings.size() == 0)
    {
        std::cerr << "Warning: no variables specified" << std::endl;
    }

    std::vector<NcVar*> invars;
    for (std::vector<std::string>::const_iterator it = varstrings.begin();
         it != varstrings.end(); ++it)
    {
        NcVar* var = infile.get_var((*it).c_str());
        if (var == NULL)
        {
            std::cerr << "Error: " << *it << ": no such variable" << std::endl;
            infile.close();
            outfile.close();
            return EXIT_FAILURE;
        }
        invars.push_back(var);
    }

    // extract the distinct set of dims
    std::map<std::string, NcDim*> indims;
    for (std::vector<NcVar*>::const_iterator it = invars.begin();
         it != invars.end(); ++it)
    {
        NcVar* var = *it;
        for (int i = 0; i < var->num_dims(); ++i)
        {
            NcDim* dim = var->get_dim(i);
            indims[dim->name()] = dim;
        }
    }

    // add dims to outfile
    std::map<std::string, NcDim*> outdims;
    for (std::map<std::string, NcDim*>::const_iterator it = indims.begin();
         it != indims.end(); ++it)
    {
        NcDim* dim = (*it).second;
        NcDim* outdim = NULL;
        if (dim->is_unlimited())
        {
            outdim = outfile.add_dim(dim->name());
        }
        else
        {
            outdim = outfile.add_dim(dim->name(), dim->size());
        }

        if (outdim != NULL)
        {
            outdims[outdim->name()] = outdim;
        }
    }

    // create variables
    for (std::vector<NcVar*>::const_iterator it = invars.begin();
         it != invars.end(); ++it)
    {
        NcVar* var = *it;
        std::vector<const NcDim*> dims(var->num_dims());
        for (int i = 0; i < var->num_dims(); ++i)
        {
            dims[i] = outdims[var->get_dim(i)->name()];
        }
        NcVar* outvar = outfile.add_var(var->name(), var->type(), var->num_dims(), &dims[0]);

        // identify largest dim, if dim (nearly) exceeds main memory, split along that dim
        int maxdim = -1;
        long maxdimsize = 0;
        long totallen = 1;
        for (int i = 0; i < var->num_dims(); ++i)
        {
            NcDim* dim = var->get_dim(i);
            if (dim->size() > maxdimsize)
            {
//.........这里部分代码省略.........

int main(int argc, char ** argv) {

try {

	// Parameters
	Parameters param;

	// Output filename
	std::string strOutputFile;

	// Horizontal minimum wave number
	int nKmin;

	// Horizontal maximum wave number
	int nKmax;

	// Parse the command line
	BeginCommandLine()
		CommandLineInt(param.nPhiElements, "n", 40);
		CommandLineInt(nKmin, "kmin", 1);
		CommandLineInt(nKmax, "kmax", 20);
		CommandLineDouble(param.dXscale, "X", 1.0);
		CommandLineDouble(param.dT0, "T0", 300.0);
		CommandLineDouble(param.dU0, "U0", 20.0);
		CommandLineDouble(param.dG, "G", 9.80616);
		CommandLineDouble(param.dOmega, "omega", 7.29212e-5);
		CommandLineDouble(param.dGamma, "gamma", 1.4);
		CommandLineString(strOutputFile, "out", "wave.nc");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	AnnounceBanner();

	// Generate latitude values
	param.GenerateLatituteArray(param.nPhiElements);

	// Open NetCDF file
	NcFile ncdf_out(strOutputFile.c_str(), NcFile::Replace);

	NcDim *dimK = ncdf_out.add_dim("k", nKmax - nKmin + 1);
	NcDim *dimLat = ncdf_out.add_dim("lat", param.nPhiElements);
	NcDim *dimEig = ncdf_out.add_dim("eig", param.nPhiElements);

	// Write parameters and latitudes to file
	param.WriteToNcFile(ncdf_out, dimLat, dimLatS);

	// Wave numbers 
	NcVar *varK = ncdf_out.add_var("k", ncInt, dimK);

	DataVector<int> vecK;
	vecK.Initialize(nKmax - nKmin + 1);
	for (int nK = nKmin; nK <= nKmax; nK++) { 
		vecK[nK - nKmin] = nK;
	}

	varK->set_cur((long)0);
	varK->put(vecK, nKmax - nKmin + 1);

	// Eigenvalues
	NcVar *varMR = ncdf_out.add_var("mR", ncDouble, dimK, dimEig);
	NcVar *varMI = ncdf_out.add_var("mI", ncDouble, dimK, dimEig);

	NcVar *varUR = ncdf_out.add_var("uR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varUI = ncdf_out.add_var("uI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varVR = ncdf_out.add_var("vR", ncDouble, dimK, dimEig, dimLatS);
	NcVar *varVI = ncdf_out.add_var("vI", ncDouble, dimK, dimEig, dimLatS);

	NcVar *varPR = ncdf_out.add_var("pR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varPI = ncdf_out.add_var("pI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varWR = ncdf_out.add_var("wR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varWI = ncdf_out.add_var("wI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varRhoR = ncdf_out.add_var("rhoR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varRhoI = ncdf_out.add_var("rhoI", ncDouble, dimK, dimEig, dimLat);

	// Allocate temporary arrays
	DataVector<double> dUR;
	dUR.Initialize(param.nPhiElements);

	DataVector<double> dUI;
	dUI.Initialize(param.nPhiElements);

	DataVector<double> dVR;
	dVR.Initialize(param.nPhiElements-1);

	DataVector<double> dVI;
	dVI.Initialize(param.nPhiElements-1);

	DataVector<double> dPR;
	dPR.Initialize(param.nPhiElements);

	DataVector<double> dPI;
	dPI.Initialize(param.nPhiElements);

	DataVector<double> dWR;
	dWR.Initialize(param.nPhiElements);

//.........这里部分代码省略.........

int main(int argc, char ** argv) {

	MPI_Init(&argc, &argv);

	NcError error(NcError::silent_nonfatal);

try {

	// Input filename
	std::string strInputFile;

	// Output filename
	std::string strOutputFile;

	// Separate topography file
	std::string strTopographyFile;

	// List of variables to extract
	std::string strVariables;

	// Extract geopotential height
	bool fGeopotentialHeight;

	// Pressure levels to extract
	std::string strPressureLevels;

	// Height levels to extract
	std::string strHeightLevels;

	// Extract variables at the surface
	bool fExtractSurface;

	// Extract total energy
	bool fExtractTotalEnergy;

	// Parse the command line
	BeginCommandLine()
		CommandLineString(strInputFile, "in", "");
		CommandLineString(strOutputFile, "out", "");
		CommandLineString(strVariables, "var", "");
		CommandLineBool(fGeopotentialHeight, "output_z");
		CommandLineBool(fExtractTotalEnergy, "output_energy");
		CommandLineString(strPressureLevels, "p", "");
		CommandLineString(strHeightLevels, "z", "");
		CommandLineBool(fExtractSurface, "surf");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	AnnounceBanner();

	// Check command line arguments
	if (strInputFile == "") {
		_EXCEPTIONT("No input file specified");
	}
	if (strOutputFile == "") {
		_EXCEPTIONT("No output file specified");
	}
	if (strVariables == "") {
		_EXCEPTIONT("No variables specified");
	}

	// Parse variable string
	std::vector< std::string > vecVariableStrings;

	ParseVariableList(strVariables, vecVariableStrings);

	// Check variables
	if (vecVariableStrings.size() == 0) {
		_EXCEPTIONT("No variables specified");
	}

	// Parse pressure level string
	std::vector<double> vecPressureLevels;

	ParseLevelArray(strPressureLevels, vecPressureLevels);

	int nPressureLevels = (int)(vecPressureLevels.size());

	for (int k = 0; k < nPressureLevels; k++) {
		if (vecPressureLevels[k] <= 0.0) {
			_EXCEPTIONT("Non-positive pressure values not allowed");
		}
	}

	// Parse height level string
	std::vector<double> vecHeightLevels;

	ParseLevelArray(strHeightLevels, vecHeightLevels);

	int nHeightLevels = (int)(vecHeightLevels.size());

	// Check pressure levels
	if ((nPressureLevels == 0) &&
		(nHeightLevels == 0) &&
		(!fExtractSurface)
	) {
		_EXCEPTIONT("No pressure / height levels to process");
	}

//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Build up output NetCDF file name and open it
  sprintf( NcName, "%s_%4d%02d%02d_%02d%02d.nc", channel,
           tmtime->tm_year + 1900, tmtime->tm_mon + 1, tmtime->tm_mday,
	   tmtime->tm_hour, tmtime->tm_min );
  NcFile ncf ( NcName , NcFile::Replace );
  if (! ncf.is_valid()) return false;

  // Fill arrays on creation
  ncf.set_fill(NcFile::Fill);

  // Add Global Attributes
  if (! ncf.add_att("Satellite", MSG_spacecraft_name(spc).c_str()))
	            return false;
  sprintf(reftime, "%04d-%02d-%02d %02d:%02d:00 UTC",
      tmtime->tm_year + 1900, tmtime->tm_mon + 1, tmtime->tm_mday,
      tmtime->tm_hour, tmtime->tm_min);
  if (! ncf.add_att("Antenna", "Fixed") ) return false;
  if (! ncf.add_att("Receiver", "HIMET") ) return false;
  if (! ncf.add_att("Time", reftime) ) return false;
  if (! ncf.add_att("Area_Name", "SpaceView" ) ) return false;
  sprintf(projname, "GEOS(%3.1f)", sublon);
  if (! ncf.add_att("Projection", projname) ) return false;
  if (! ncf.add_att("Columns", AreaNpix ) ) return false;
  if (! ncf.add_att("Lines", AreaNlin ) ) return false;
  if (! ncf.add_att("SampleX", 1.0 ) ) return false;
  if (! ncf.add_att("SampleY", 1.0 ) ) return false;
  if (! ncf.add_att("AreaStartPix", AreaPixStart ) ) return false;
  if (! ncf.add_att("AreaStartLin", AreaLinStart ) ) return false;
  if (! ncf.add_att("Column_Scale_Factor", cfac) ) return false;
  if (! ncf.add_att("Line_Scale_Factor", lfac) ) return false;
  if (! ncf.add_att("Column_Offset", coff) ) return false;
  if (! ncf.add_att("Line_Offset", loff) ) return false;
  if (! ncf.add_att("Orbit_Radius", sh) ) return false;
  if (! ncf.add_att("Longitude", sublon) ) return false;
  if (! ncf.add_att("NortPolar", 1) ) return false;
  if (! ncf.add_att("NorthSouth", 1) ) return false;
  if (! ncf.add_att("title", TITLE) ) return false;
  if (! ncf.add_att("Institution", INSTITUTION) ) return false;
  if (! ncf.add_att("Type", TYPE) ) return false;
  if (! ncf.add_att("Version", HIMET_VERSION) ) return false;
  if (! ncf.add_att("Conventions", "COARDS") ) return false;
  if (! ncf.add_att("history", "Created from raw data") ) return false;

  // Dimensions
  wd = AreaNpix;
  hg = AreaNlin;
  bpp = header[0].image_structure->number_of_bits_per_pixel;
  ncal = (int) pow(2.0, bpp);

  tdim = ncf.add_dim("time");
  if (!tdim->is_valid()) return false;
  ldim = ncf.add_dim("line", hg);
  if (!ldim->is_valid()) return false;
  cdim = ncf.add_dim("column", wd);
  if (!cdim->is_valid()) return false;
  caldim = ncf.add_dim("calibration", ncal);
  if (!caldim->is_valid()) return false;

  // Get calibration values
  cal = PRO_data->prologue->radiometric_proc.get_calibration((int) chn, bpp);

  // Add Calibration values
  NcVar *cvar = ncf.add_var("calibration", ncFloat, caldim);
  if (!cvar->is_valid()) return false;
  cvar->add_att("long_name", "Calibration coefficients");
  cvar->add_att("variable", channel);
  if (chn > 3 && chn < 12)
    cvar->add_att("units", "K");
  else
    cvar->add_att("units", "mW m^-2 sr^-1 (cm^-1)^-1");
  if (!cvar->put(cal, ncal)) return false;

  tvar = ncf.add_var("time", ncDouble, tdim);
  if (!tvar->is_valid()) return false;
  tvar->add_att("long_name", "Time");
  tvar->add_att("units", "seconds since 2000-01-01 00:00:00 UTC");
  double atime;
  time_t ttime;
  extern long timezone;
  ttime = mktime(tmtime);
  atime = ttime - 946684800 - timezone;
  if (!tvar->put(&atime, 1)) return false;

  ivar = ncf.add_var(channel, ncShort, tdim, ldim, cdim);
  if (!ivar->is_valid()) return false;
  if (!ivar->add_att("add_offset", 0.0)) return false;
  if (!ivar->add_att("scale_factor", 1.0)) return false;
  if (!ivar->add_att("chnum", chn)) return false;

  // Write output values
  if (!ivar->put((const short int *) pixels, 1, hg, wd)) return false;

  // Close NetCDF output
  (void) ncf.close( );

  delete [ ] pixels;
  delete [ ] cal;

  return( true );
}

void TracerManager::output(const string &fileName)
{
    // -------------------------------------------------------------------------
    // output polygon stuffs
    polygonManager.output(fileName);
#ifdef TTS_REMAP
    // -------------------------------------------------------------------------
    NcFile file(fileName.c_str(), NcFile::Write);
    if (!file.is_valid()) {
        Message message;
        message << "Failed to open tracer output file \"";
        message << fileName << "\" for appending meshed density field!";
        REPORT_ERROR(message.str());
    }
    // -------------------------------------------------------------------------
    // output tracer densities on the polygons
    NcDim *numPolygonDim = file.get_dim("num_total_polygon");
    double q0[polygonManager.polygons.size()];
    for (int l = 0; l < tracerNames.size(); ++l) {
        char varName[30];
        sprintf(varName, "q%d", l);
        NcVar *qVar = file.add_var(varName, ncDouble, numPolygonDim);
        Polygon *polygon = polygonManager.polygons.front();
        for (int i = 0; i < polygonManager.polygons.size(); ++i) {
            q0[i] = polygon->tracers[l].getDensity();
            polygon = polygon->next;
        }
        qVar->put(q0, polygonManager.polygons.size());
    }
    // -------------------------------------------------------------------------
    // output tracer densities on the mesh
    const RLLMesh &mesh = tracerDensities[0].getMesh();
    int numLon = mesh.getNumLon()-2;
    int numLat = mesh.getNumLat();
    double lon[numLon], lat[numLat];
    for (int i = 0; i < numLon; ++i)
        lon[i] = mesh.lon(i+1)*Rad2Deg;
    for (int j = 0; j < numLat; ++j)
        lat[j] = mesh.lat(j)*Rad2Deg;
    NcDim *lonDim = file.add_dim("lon", numLon);
    NcDim *latDim = file.add_dim("lat", numLat);
    NcVar *lonVar = file.add_var("lon", ncDouble, lonDim);
    lonVar->add_att("long_name", "longitude");
    lonVar->add_att("units", "degrees_east");
    lonVar->put(lon, numLon);
    NcVar *latVar = file.add_var("lat", ncDouble, latDim);
    latVar->add_att("long_name", "latitude");
    latVar->add_att("units", "degrees_north");
    latVar->put(lat, numLat);
    NcVar *areaVar = file.add_var("area_mesh", ncDouble, latDim, lonDim);
    areaVar->add_att("long_name", "area of fixed mesh cell");
    areaVar->add_att("units", "m2");
    double area[numLat][numLon];
    for (int i = 0; i < numLon; ++i)
        for (int j = 0; j < numLat; ++j)
            area[j][i] = tracerDensities[0].getMesh(Field::Bound).area(i, j);
    areaVar->put(&area[0][0], numLat, numLon);
    double q[numLat][numLon];
    for (int l = 0; l < tracerNames.size(); ++l) {
        char varName[30];
        sprintf(varName, "q%d_mesh", l);
        NcVar *qVar = file.add_var(varName, ncDouble, latDim, lonDim);
        qVar->add_att("long_name", tracerNames[l].c_str());
        for (int i = 0; i < numLon; ++i)
            for (int j = 0; j < numLat; ++j)
                q[j][i] = tracerDensities[l].values(i+1, j, 0).getNew();
        qVar->put(&q[0][0], numLat, numLon);
    }
    // -------------------------------------------------------------------------
    file.close();
#endif
    NOTICE("TracerManager", fileName+" is generated.");
}

bool OutputManagerReference::OpenFile(
	const std::string & strFileName
) {
#ifdef TEMPEST_NETCDF
	// Determine processor rank; only proceed if root node
	int nRank = 0;

#ifdef TEMPEST_MPIOMP
	MPI_Comm_rank(MPI_COMM_WORLD, &nRank);
#endif

	// The active model
	const Model & model = m_grid.GetModel();

	// Open NetCDF file on root process
	if (nRank == 0) {

		// Check for existing NetCDF file
		if (m_pActiveNcOutput != NULL) {
			_EXCEPTIONT("NetCDF file already open");
		}

		// Append .nc extension to file
		std::string strNcFileName = strFileName + ".nc";

		// Open new NetCDF file
		m_pActiveNcOutput = new NcFile(strNcFileName.c_str(), NcFile::Replace);
		if (m_pActiveNcOutput == NULL) {
			_EXCEPTION1("Error opening NetCDF file \"%s\"",
				strNcFileName.c_str());
		}
		if (!m_pActiveNcOutput->is_valid()) {
			_EXCEPTION1("Error opening NetCDF file \"%s\"",
				strNcFileName.c_str());
		}

		// Create nodal time dimension
		NcDim * dimTime = m_pActiveNcOutput->add_dim("time");
		if (dimTime == NULL) {
			_EXCEPTIONT("Error creating \"time\" dimension");
		}

		m_varTime = m_pActiveNcOutput->add_var("time", ncDouble, dimTime);
		if (m_varTime == NULL) {
			_EXCEPTIONT("Error creating \"time\" variable");
		}

		std::string strUnits =
			"days since " + model.GetStartTime().ToDateString();

		std::string strCalendarName = model.GetStartTime().GetCalendarName();

		m_varTime->add_att("long_name", "time");
		m_varTime->add_att("units", strUnits.c_str());
		m_varTime->add_att("calendar", strCalendarName.c_str());
		m_varTime->add_att("bounds", "time_bnds");

		// Create levels dimension
		NcDim * dimLev =
			m_pActiveNcOutput->add_dim("lev", m_dREtaCoord.GetRows());

		// Create interfaces dimension
		NcDim * dimILev =
			m_pActiveNcOutput->add_dim("ilev", m_grid.GetRElements()+1);

		// Create latitude dimension
		NcDim * dimLat =
			m_pActiveNcOutput->add_dim("lat", m_nYReference);

		// Create longitude dimension
		NcDim * dimLon =
			m_pActiveNcOutput->add_dim("lon", m_nXReference);

		// Output physical constants
		const PhysicalConstants & phys = model.GetPhysicalConstants();

		m_pActiveNcOutput->add_att("earth_radius", phys.GetEarthRadius());
		m_pActiveNcOutput->add_att("g", phys.GetG());
		m_pActiveNcOutput->add_att("omega", phys.GetOmega());
		m_pActiveNcOutput->add_att("alpha", phys.GetAlpha());
		m_pActiveNcOutput->add_att("Rd", phys.GetR());
		m_pActiveNcOutput->add_att("Cp", phys.GetCp());
		m_pActiveNcOutput->add_att("T0", phys.GetT0());
		m_pActiveNcOutput->add_att("P0", phys.GetP0());
		m_pActiveNcOutput->add_att("rho_water", phys.GetRhoWater());
		m_pActiveNcOutput->add_att("Rvap", phys.GetRvap());
		m_pActiveNcOutput->add_att("Mvap", phys.GetMvap());
		m_pActiveNcOutput->add_att("Lvap", phys.GetLvap());

		// Output grid parameters
		m_pActiveNcOutput->add_att("Ztop", m_grid.GetZtop());

		// Output equation set
		const EquationSet & eqn = model.GetEquationSet();

		m_pActiveNcOutput->add_att("equation_set", eqn.GetName().c_str());

		// Create variables
		for (int c = 0; c < eqn.GetComponents(); c++) {
			if ((m_fOutputAllVarsOnNodes) ||
//.........这里部分代码省略.........