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

void NetCdfConfigureDialog::getDimEdges(int dimId, unsigned &size, double &firstValue, double &lastValue)
{
    if ((_currentFile->get_var(_currentVar->get_dim(dimId)->name())) != NULL)
    {
        NcVar *tmpVarOfDim = _currentFile->get_var(_currentVar->get_dim(dimId)->name());
        if ((tmpVarOfDim->num_dims()) == 1)
        {
            int sizeOfDim = tmpVarOfDim->get_dim(0)->size();
            size = sizeOfDim;
            double arrayOfDimStart[1] = {0};
            size_t edgeOfArray[1] = {1};
            long edgeOrigin[1] = {0};
            tmpVarOfDim->set_cur(edgeOrigin);
            tmpVarOfDim->get(arrayOfDimStart,edgeOfArray);
            firstValue = arrayOfDimStart[0];
            double arrayOfDimEnd[1] = {0};
            edgeOrigin[0] = sizeOfDim - 1;
            tmpVarOfDim->set_cur(edgeOrigin);
            tmpVarOfDim->get(arrayOfDimEnd,edgeOfArray);
            lastValue = arrayOfDimEnd[0];
        }
    } else {
        size = 0;
        firstValue = 0;
        lastValue = 0;
    }
}

FieldPtr FileArome::getFieldCore(std::string iVariable, int iTime) const {
   // Not cached, retrieve data
   NcVar* var = getVar(iVariable);
   int nLat  = mNLat;
   int nLon  = mNLon;

   int numDims = var->num_dims();

   long* count;
   long totalCount = nLat*nLon;
   if(numDims == 4) {
      // Variable has a surface dimension
      count = new long[4];
      count[0] = 1;
      count[1] = 1;
      count[2] = nLat;
      count[3] = nLon;
      var->set_cur(iTime, 0, 0, 0);
   }
   else if(numDims == 3) {
      count = new long[3];
      count[0] = 1;
      count[1] = nLat;
      count[2] = nLon;
      var->set_cur(iTime, 0, 0);
   }
   else {
      std::stringstream ss;
      ss << "Cannot read variable '" << iVariable << "' from '" << getFilename() << "'";
      Util::error(ss.str());
   }
   float* values = new float[nLat*nLon];
   var->get(values, count);
   float MV = getMissingValue(var);

   float offset = getOffset(var);
   float scale = getScale(var);
   int index = 0;
   FieldPtr field = getEmptyField();
   for(int lat = 0; lat < nLat; lat++) {
      for(int lon = 0; lon < nLon; lon++) {
         float value = values[index];
         assert(index < totalCount);
         if(value == MV) {
            // Field has missing value indicator and the value is missing
            // Save values using our own internal missing value indicator
            value = Util::MV;
         }
         else {
            value = scale*values[index] + offset;
         }
         (*field)(lat,lon,0) = value;
         index++;
      }
   }
   delete[] values;
   delete[] count;
   return field;
}

template<class T> static bool load_nc_array(const NcFile& ncf, const string& name, vector<T>& dest, bool required = true, int offset = 0, int count = -1)
{
	NcVar *v = load_nc_variable(ncf, name.c_str(), required);
	if (v)
	{
		vector<long> offsets = list_of(offset).repeat(v->num_dims()-1, 0);
		v->set_cur(&offsets.front());
		vector<long> counts (v->num_dims());
		long* shape = v->edges();
		transform(shape, shape + v->num_dims(), offsets.begin(), counts.begin(), minus<long>());
		delete shape;
		if (count > 0)
		{
			counts[0] = count;
		}
		dest.resize(product(counts));
		bool success = v->get(&dest.front(), &counts.front());
		if (!success)
		{
			dest.resize(0);
			check(!required, string("NetcdfDataset::load_nc_array<") + typeid(T).name() + "> " + name + '\n' + "failed with offset " + str(offsets) + ", counts " + str(counts));
		}
		return success;
	}
	return false;
}

//YUAN: recid - the order (from ZERO) in the .nc file, chtid - the cohort id
int SiteinInputer::getRecID(const int &siteid){
	NcError err(NcError::silent_nonfatal);

	NcFile siteFile(siteinfname.c_str(), NcFile::ReadOnly);
 	NcVar* siteidV = siteFile.get_var("CHTID");

 	int id = -1;
	for (int i=0; i<(int)siteidV->num_vals(); i++){
		siteidV->set_cur(i);
		siteidV->get(&id, 1);
		if(id==siteid) return i;
	}
	return -1;
}

void Regioner::createCohorList4Run(){
	// read in a list of cohorts to run

	//netcdf error
	NcError err(NcError::silent_nonfatal);

	//open file and check if valid
	string filename = md.runchtfile;
	NcFile runFile(filename.c_str(), NcFile::ReadOnly);
 	if(!runFile.is_valid()){
 		string msg = filename+" is not valid";
 		char* msgc = const_cast< char* > ( msg.c_str());
 		throw Exception(msgc, I_NCFILE_NOT_EXIST);
 	}
 	
 	NcDim* chtD = runFile.get_dim("CHTID");
 	if(!chtD->is_valid()){
 		string msg="CHT Dimension is not valid in createCohortList4Run";
 		char* msgc = const_cast<char*> (msg.c_str());
 		throw Exception(msgc, I_NCDIM_NOT_EXIST);
 	}
 	
 	NcVar* chtV = runFile.get_var("CHTID");
 	if(chtV==NULL){
 		string msg="Cannot get CHTID in createCohortList4Run ";
 		char* msgc = const_cast<char*> (msg.c_str());
 		throw Exception(msgc, I_NCVAR_NOT_EXIST);
 	}

 	int numcht = chtD->size();
 	
	int chtid  = -1;
	int chtid0 = -1;
	int chtidx = -1;
	for (int i=0; i<numcht; i++){
		chtV->set_cur(i);
   		chtV->get(&chtid, 1);
   		runchtlist.push_back(chtid);
	   	
	   	if (i==0) chtid0=chtid;
	   	if (i==numcht-1) chtidx=chtid;
   	}

	cout <<md.casename << ": " <<numcht <<"  cohorts to be run @" <<md.runstages<< "\n";
	cout <<"   from:  " <<chtid0<<"  to:  " <<chtidx <<"\n";
   
};

static string get_nc_string(const NcFile& ncf, const string& name, int offset = 0, bool required = true)
{
	static array<long, 2> offsets = {{0, 0}};
	static array<long, 2> counts = {{1, 0}};
	NcVar *v = load_nc_variable(ncf, name.c_str(), required);
	if (v)
	{
		long* shape = v->edges();
		offsets.front() = offset;
		counts.back() = shape[1];
		v->set_cur(&offsets.front());
		char* temp = new char [shape[1]];
		delete shape;
		bool success = v->get(temp, &counts.front());
		if(!success)
		{
			check(!required, " index " + str(offset) + " out of bounds for " + name + " in netcdf file");
		}
		string s(temp);
		delete [] temp;
		return s;
	}
	return "";
}

RadarData_t ReadRadarFile(const string &filename,
                          const float latmin, const float latmax,
                          const float lonmin, const float lonmax)
{
    RadarData_t inputData;

    // Sets the internal netcdf error handling to nonfatal for this scope
    NcError error_handler(NcError::verbose_nonfatal);

    NcFile radarFile(filename.c_str());

    if (!radarFile.is_valid())
    {
        cerr << "ERROR: Could not open radar file: " << filename << " for reading.\n";

        // Error is indicated by the lack of initialization of
        // the filename member of the struct.
        return(inputData);
    }

    NcVar* latVar = radarFile.get_var("lat");

    if (NULL == latVar)
    {
        cerr << "ERROR: invalid data file.  No variable called 'lat'!\n";
        radarFile.close();
        return(inputData);
    }

    long latCnt = latVar->num_vals();
    double* latVals = new double[latCnt];
    latVar->get(latVals, latCnt);
    inputData.latUnits = GrabAttribute(latVar, 0);
    inputData.latSpacing = strtod(GrabAttribute(latVar, 1).c_str(), NULL);

    const long minLatIndex = lower_bound(latVals, latmin, latCnt);
    const long maxLatIndex = upper_bound(latVals, latmax, latCnt);

    delete latVals;
    latCnt = (maxLatIndex - minLatIndex) + 1;
    latVar->set_cur(minLatIndex);
    inputData.latVals = new double[latCnt];
    latVar->get(inputData.latVals, latCnt);


    NcVar* lonVar = radarFile.get_var("lon");

    if (NULL == lonVar)
    {
        cerr << "ERROR: invalid data file.  No variable called 'lon'!\n";
        radarFile.close();
        return(inputData);
    }

    long lonCnt = lonVar->num_vals();
    double* lonVals = new double[lonCnt];
    lonVar->get(lonVals, lonCnt);
    inputData.lonUnits = GrabAttribute(lonVar, 0);
    inputData.lonSpacing = strtod(GrabAttribute(lonVar, 1).c_str(), NULL);

    const long minLonIndex = lower_bound(lonVals, lonmin, lonCnt);
    const long maxLonIndex = upper_bound(lonVals, lonmax, lonCnt);

    delete lonVals;
    lonCnt = (maxLonIndex - minLonIndex) + 1;
    lonVar->set_cur(minLonIndex);
    inputData.lonVals = new double[lonCnt];
    lonVar->get(inputData.lonVals, lonCnt);

    NcVar* reflectVar = NULL;
    reflectVar = radarFile.get_var("value");

    if ( reflectVar == NULL )
    {
        // Try this variable name
        reflectVar = radarFile.get_var("Reflectivity");
    }

    if (reflectVar == NULL)
    {
        cerr << "ERROR: invalid data file.  No variable called 'value'!\n";
        radarFile.close();
        return(inputData);
    }


    inputData.dataEdges = reflectVar->edges();	// [0] - time, [1] - lat, [2] - lon
    inputData.dataEdges[1] = latCnt;
    inputData.dataEdges[2] = lonCnt;
    inputData.dataVals = new double[inputData.dataEdges[0] * inputData.dataEdges[1] * inputData.dataEdges[2]];
    reflectVar->set_cur(0, minLatIndex, minLonIndex);
    reflectVar->get(inputData.dataVals, inputData.dataEdges);

    inputData.var_LongName = GrabAttribute(reflectVar, 0);
    inputData.var_Units = "dBZ";//GrabAttribute(reflectVar, 1);


    NcVar* timeVar = radarFile.get_var("time");

    if (NULL == timeVar)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      GetInputFileList(inList,vecFiles);
    }

    //open up first file
    NcFile readin(vecFiles[0].c_str());
    if (!readin.is_valid()){
      _EXCEPTION1("Unable to open file %s for reading",\
        vecFiles[0].c_str());
    }
    int tLen,latLen,lonLen;

    NcDim * time = readin.get_dim("time");
    tLen = time->size();
    NcVar * timeVar = readin.get_var("time");

    NcDim * lat = readin.get_dim("lat");
    latLen = lat->size();
    NcVar * latVar = readin.get_var("lat");

    NcDim * lon = readin.get_dim("lon");
    lonLen = lon->size();
    NcVar * lonVar = readin.get_var("lon");

    //read input variable
    NcVar * inVar = readin.get_var(varName.c_str());
    //Create output matrix
    DataMatrix<double> outMat(latLen,lonLen);
    densCalc(inVar,outMat);
    //Option for calculating the yearly standard deviation 
 /*   if (calcStdDev){
      for (int a=0; a<latLen; a++){
        for (int b=0; b<lonLen; b++){
          storeMat[0][a][b] = outMat[a][b];
        }
      }
    }
*/
    //If multiple files, add these values to the output
    if (vecFiles.size()>1){
      DataMatrix<double> addMat(latLen,lonLen);
      std::cout<<"There are "<<vecFiles.size()<<" files."<<std::endl;
      for (int v=1; v<vecFiles.size(); v++){
        NcFile addread(vecFiles[v].c_str());
        NcVar * inVar = addread.get_var(varName.c_str());
        densCalc(inVar,addMat);
        for (int a=0; a<latLen; a++){
          for (int b=0; b<lonLen; b++){
            outMat[a][b]+=addMat[a][b];
          }
        }
/*        if (calcStdDev){
          for (int a=0; a<latLen; a++){
            for (int b=0; b<lonLen; b++){
              storeMat[v][a][b] = addMat[a][b];
            }
          }
        }*/
        addread.close(); 
      }
      //Divide output by number of files
      double div = 1./((double) vecFiles.size());
      for (int a=0; a<latLen; a++){
        for (int b=0; b<lonLen; b++){
          outMat[a][b]*=div;
        }
      }
    }
 
    NcFile readout(outFile.c_str(),NcFile::Replace, NULL,0,NcFile::Offset64Bits);
    NcDim * outLat = readout.add_dim("lat", latLen);
    NcDim * outLon = readout.add_dim("lon", lonLen);
    NcVar * outLatVar = readout.add_var("lat",ncDouble,outLat);
    NcVar * outLonVar = readout.add_var("lon",ncDouble,outLon);
    std::cout<<"Copying dimension attributes."<<std::endl;
    copy_dim_var(latVar,outLatVar);
    copy_dim_var(lonVar,outLonVar);

    std::cout<<"Creating density variable."<<std::endl;

    NcVar * densVar = readout.add_var("dens",ncDouble,outLat,outLon);
    densVar->set_cur(0,0);
    densVar->put((&outMat[0][0]),latLen,lonLen);

/*    if (calcStdDev){
      NcVar * stdDevVar = readout.add_var("stddev", ncDouble,outLat,outLon);
      DataMatrix<double> stdDevMat(latLen,lonLen);
      yearlyStdDev(storeMat,vecFiles.size(),latLen,lonLen,stdDevMat);
      stdDevVar->set_cur(0,0);
      stdDevVar->put(&(stdDevMat[0][0]),latLen,lonLen);
      std::cout<<" created sd variable"<<std::endl;

    }
*/
    readout.close();
    readin.close();
  } 
  catch (Exception &e){
    std::cout<<e.ToString()<<std::endl;
  }
}

void ReadCFTimeDataFromNcFile(
	NcFile * ncfile,
	const std::string & strFilename,
	std::vector<Time> & vecTimes,
	bool fWarnOnMissingCalendar
) {
	// Empty existing Time vector
	vecTimes.clear();

	// Get time dimension
	NcDim * dimTime = ncfile->get_dim("time");
	if (dimTime == NULL) {
		_EXCEPTION1("Dimension \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}

	// Get time variable
	NcVar * varTime = ncfile->get_var("time");
	if (varTime == NULL) {
		_EXCEPTION1("Variable \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}
	if (varTime->num_dims() != 1) {
		_EXCEPTION1("Variable \"time\" has more than one dimension in file \"%s\"",
			strFilename.c_str());
	}
	if (strcmp(varTime->get_dim(0)->name(), "time") != 0) {
		_EXCEPTION1("Variable \"time\" does not have dimension \"time\" in file \"%s\"",
			strFilename.c_str());
	}

	// Calendar attribute
	NcAtt * attTimeCal = varTime->get_att("calendar");
	std::string strCalendar;
	if (attTimeCal == NULL) {
		if (fWarnOnMissingCalendar) {
			Announce("WARNING: Variable \"time\" is missing \"calendar\" attribute; assuming \"standard\"");
		}
		strCalendar = "standard";
	} else {
		strCalendar = attTimeCal->as_string(0);
	}
	Time::CalendarType eCalendarType =
		Time::CalendarTypeFromString(strCalendar);

	// Units attribute
	NcAtt * attTimeUnits = varTime->get_att("units");
	if (attTimeUnits == NULL) {
		_EXCEPTION1("Variable \"time\" is missing \"units\" attribute in file \"%s\"",
			strFilename.c_str());
	}
	std::string strTimeUnits = attTimeUnits->as_string(0);

	// Load in time data
	DataVector<int> vecTimeInt;
	DataVector<float> vecTimeFloat;
	DataVector<double> vecTimeDouble;
	DataVector<ncint64> vecTimeInt64;

	if (varTime->type() == ncInt) {
		vecTimeInt.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt[0]), dimTime->size());

	} else if (varTime->type() == ncFloat) {
		vecTimeFloat.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeFloat[0]), dimTime->size());

	} else if (varTime->type() == ncDouble) {
		vecTimeDouble.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeDouble[0]), dimTime->size());

	} else if (varTime->type() == ncInt64) {
		vecTimeInt64.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt64[0]), dimTime->size());

	} else {
		_EXCEPTION1("Variable \"time\" has invalid type "
			"(expected \"int\", \"int64\", \"float\" or \"double\")"
			" in file \"%s\"", strFilename.c_str());
	}

	for (int t = 0; t < dimTime->size(); t++) {
		Time time(eCalendarType);
		if (varTime->type() == ncInt) {
			time.FromCFCompliantUnitsOffsetInt(
				strTimeUnits,
				vecTimeInt[t]);

		} else if (varTime->type() == ncFloat) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
				static_cast<double>(vecTimeFloat[t]));

		} else if (varTime->type() == ncDouble) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            std::vector<long> counts;
            for (int i = 0; i < var->num_dims(); ++i)
            {
                NcDim* dim = var->get_dim(i);
                long current = 0;
                long count = dim->size();
                if (i == maxdim)
                {
                    current = pagesizedim * p;
                    count = pagesizedim;
                    if (p == pages -1)
                    {
                        if (dim->size() % pagesizedim != 0)
                        {
                            count = dim->size() % pagesizedim;
                        }
                    }
                }
                cur.push_back(current);
                counts.push_back(count);
                len *= count;
            }
            curvec.push_back(cur);
            countsvec.push_back(counts);
            lengths.push_back(len);
        }

        std::vector< std::vector<long> >::const_iterator it1;
        std::vector< std::vector<long> >::const_iterator it2;
        std::vector<long>::const_iterator it3;

        for (it1 = curvec.begin(), it2 = countsvec.begin(), it3 = lengths.begin();
             it1 != curvec.end() && it2 != countsvec.end() && it3 != lengths.end(); ++it1, ++it2, ++it3)
        {
            std::vector<long> cur = *it1;
            std::vector<long> counts = *it2;
            long len = *it3;

            var->set_cur(&cur[0]);
            outvar->set_cur(&cur[0]);
            switch (outvar->type())
            {
            case ncByte:
            {
                ncbyte* barr = new ncbyte[len];
                var->get(barr, &counts[0]);
                outvar->put(barr, &counts[0]);
                delete[] barr;
                break;
            }
            case ncChar:
            {
                char* carr = new char[len];
                var->get(carr, &counts[0]);
                outvar->put(carr, &counts[0]);
                delete[] carr;
                break;
            }
            case ncShort:
            {
                short* sarr = new short[len];
                var->get(sarr, &counts[0]);
                outvar->put(sarr, &counts[0]);
                delete[] sarr;
                break;
            }
            case ncInt:
            {
                long* larr = new long[len];
                var->get(larr, &counts[0]);
                outvar->put(larr, &counts[0]);
                delete[] larr;
                break;
            }
            case ncFloat:
            {
                float* farr = new float[len];
                var->get(farr, &counts[0]);
                outvar->put(farr, &counts[0]);
                delete[] farr;
                break;
            }
            case ncDouble:
            {
                double* darr = new double[len];
                var->get(darr, &counts[0]);
                outvar->put(darr, &counts[0]);
                delete[] darr;
                break;
            }
            default:
                break;
            }
        }
    }

    infile.close();
    outfile.close();
    return 0;
}

//.........这里部分代码省略.........
	// 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");
	}

	// Open input file
	AnnounceStartBlock("Loading input file");
	NcFile ncdf_in(strInputFile.c_str(), NcFile::ReadOnly);
	if (!ncdf_in.is_valid()) {
		_EXCEPTION1("Unable to open file \"%s\" for reading",
			strInputFile.c_str());
	}

	// Load time array
	Announce("Time");
	NcVar * varTime = ncdf_in.get_var("time");
	if (varTime == NULL) {
		_EXCEPTION1("File \"%s\" does not contain variable \"time\"",
			strInputFile.c_str());
	}
	int nTime = varTime->get_dim(0)->size();

	DataArray1D<double> dTime(nTime);
	varTime->set_cur((long)0);
	varTime->get(&(dTime[0]), nTime);

	// Load latitude array
	Announce("Latitude");
	NcVar * varLat = ncdf_in.get_var("lat");
	if (varLat == NULL) {
		_EXCEPTION1("File \"%s\" does not contain variable \"lat\"",
			strInputFile.c_str());
	}
	int nLat = varLat->get_dim(0)->size();

	DataArray1D<double> dLat(nLat);
	varLat->set_cur((long)0);
	varLat->get(&(dLat[0]), nLat);

	// Load longitude array
	Announce("Longitude");
	NcVar * varLon = ncdf_in.get_var("lon");
	if (varLon == NULL) {
		_EXCEPTION1("File \"%s\" does not contain variable \"lon\"",
			strInputFile.c_str());
	}
	int nLon = varLon->get_dim(0)->size();

	DataArray1D<double> dLon(nLon);
	varLon->set_cur((long)0);
	varLon->get(&(dLon[0]), nLon);

	// Load level array
	Announce("Level");
	NcVar * varLev = ncdf_in.get_var("lev");
	if (varLev == NULL) {

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) {

	NcError error(NcError::silent_nonfatal);

try {
	// Input filename
	std::string strInputFile;

	// Output mesh filename
	std::string strOutputFile;

	// Polynomial degree per element
	int nP = 2;

	// Parse the command line
	BeginCommandLine()
		CommandLineString(strInputFile, "in", "");
		CommandLineString(strOutputFile, "out", "");
		//CommandLineInt(nP, "np", 2);
		//CommandLineBool(fCGLL, "cgll");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	// Check file names
	if (strInputFile == "") {
		std::cout << "ERROR: No input file specified" << std::endl;
		return (-1);
	}
	if (strOutputFile == "") {
		std::cout << "ERROR: No output file specified" << std::endl;
		return (-1);
	}
	if (nP < 1) {
		std::cout << "ERROR: --np must be >= 2" << std::endl;
		return (-1);
	}

	AnnounceBanner();

	// Load input mesh
	AnnounceStartBlock("Loading input mesh");

	Mesh meshIn(strInputFile);
	meshIn.RemoveZeroEdges();

	AnnounceEndBlock("Done");

	// Construct edge map
	AnnounceStartBlock("Constructing edge map");

	meshIn.ConstructEdgeMap();

	AnnounceEndBlock("Done");

	// Build connectivity vector using edge map
	AnnounceStartBlock("Constructing connectivity");

    std::vector< std::set<int> > vecConnectivity;
    int err = GenerateConnectivityData(meshIn, vecConnectivity);
    if (err) return err;

	AnnounceEndBlock("Done");

	// Open output file
	AnnounceStartBlock("Writing connectivity file");

	NcFile ncmesh(strInputFile.c_str(), NcFile::ReadOnly);

	NcVar * varLat = ncmesh.get_var("grid_center_lat");
	NcVar * varLon = ncmesh.get_var("grid_center_lon");

	// Check if center latitudes and longitudes are already available
	DataArray1D<double> dAllLats;
	DataArray1D<double> dAllLons;

	bool fConvertLatToDegrees = true;
	bool fConvertLonToDegrees = true;

	if ((varLat == NULL) || (varLon == NULL)) {
		Announce("grid_center_lat not found, recalculating face centers");
	} else {
		Announce("grid_center_lat found in file, loading values");

		if (varLat->get_dim(0)->size() != vecConnectivity.size()) {
			_EXCEPTIONT("grid_center_lat dimension mismatch");
		}
		if (varLon->get_dim(0)->size() != vecConnectivity.size()) {
			_EXCEPTIONT("grid_center_lon dimension mismatch");
		}

		dAllLats.Allocate(vecConnectivity.size());
		varLat->set_cur((long)0);
		varLat->get(dAllLats, vecConnectivity.size());

		NcAtt * attLatUnits = varLat->get_att("units");
		std::string strLatUnits = attLatUnits->as_string(0);
		if (strLatUnits == "degrees") {
			fConvertLatToDegrees = false;
		}
//.........这里部分代码省略.........

int SiteinInputer::getSiteinData(SiteIn* sid, const int & recid){
	if (recid==-1) return -1;     //CID not exists in sitein.nc, so back to calling
 
	NcFile siteFile(siteinfname.c_str(), NcFile::ReadOnly);
 		
	NcVar* ysfV = siteFile.get_var("YSF");
	if(ysfV==NULL){
	  cout <<"cannot get ysf\n"; 
	}
	ysfV->set_cur(recid);
	ysfV->get(&sid->ysf, 1);
	 
	NcVar* envlaiV = siteFile.get_var("LAI");
	if(envlaiV==NULL){
	  cout <<"cannot get lai \n";
	}
	envlaiV->set_cur(recid);
	envlaiV->get(sid->lai,1, 12);

	NcVar* vegcV = siteFile.get_var("VEGC");
	if(vegcV==NULL){
	  cout <<"cannot get vegc \n";
	}
	vegcV->set_cur(recid);
	vegcV->get(&sid->vegc, 1);

	NcVar* vegnV = siteFile.get_var("VEGN");
	if(vegcV==NULL){
	  cout <<"cannot get vegn \n"; 
	}
	vegnV->set_cur(recid);
	vegnV->get(&sid->vegn, 1);

	NcVar* mossthickV = siteFile.get_var("MOSSTHICK");
	if(mossthickV==NULL){
	  cout <<"cannot get mossthick \n";	 
	}
	mossthickV->set_cur(recid);
	mossthickV->get(&sid->mossthick, 1);
	 
	NcVar* fibthickV = siteFile.get_var("FIBTHICK");
	if(fibthickV==NULL){
	  cout <<"cannot get fibthick \n";  
	}
	fibthickV->set_cur(recid);
	fibthickV->get(&sid->fibthick, 1);
	 
	NcVar* humthickV = siteFile.get_var("HUMTHICK");
	if(humthickV==NULL){
	  cout <<"cannot get humthick \n";  
	}
	humthickV->set_cur(recid);
	humthickV->get(&sid->humthick, 1);

	NcVar* soilcV = siteFile.get_var("SOILC");
	if(soilcV==NULL){
	  cout <<"cannot get soilc \n";
	}
	soilcV->set_cur(recid);
	soilcV->get(&sid->soilc, 1);

	NcVar* fibcV = siteFile.get_var("FIBC");
	if(fibcV==NULL){
	  cout <<"cannot get fibc \n";
	}
	fibcV->set_cur(recid);
	fibcV->get(&sid->fibc, 1);

	NcVar* humcV = siteFile.get_var("HUMC");
	if(humcV==NULL){
	  cout <<"cannot get humc \n";
	}
	humcV->set_cur(recid);
	humcV->get(&sid->humc, 1);

	NcVar* mincV = siteFile.get_var("MINC");
	if(mincV==NULL){
	  cout <<"cannot get minc \n";
	}
	mincV->set_cur(recid);
	mincV->get(&sid->minc, 1);

	NcVar* orgnV = siteFile.get_var("ORGN");
	if(orgnV==NULL){
	  cout <<"cannot get orgn \n";
	}
	orgnV->set_cur(recid);
	orgnV->get(&sid->orgn, 1);

	NcVar* avlnV = siteFile.get_var("AVLN");
	if(avlnV==NULL){
	  cout <<"cannot get avln \n";
	}
	avlnV->set_cur(recid);
	avlnV->get(&sid->avln, 1);

	NcVar* mintypeV = siteFile.get_var("MINTYPE");
	if(mintypeV==NULL){
	  cout <<"cannot get mineral type \n"; 
	}
//.........这里部分代码省略.........

int main(int argc, const char* argv[]) {
    string fileName;
    int maxDataSavePts, boundHandlingMethod;
    double stoppingTol;
    
    if (argc == 8) {
        for (int i = 0; i < argc; i++) {            
            if (argv[i][0] == '-') {
                switch(argv[i][1]) {
                    case 'n':
                        maxDataSavePts = atoi(argv[i+1]);
                        i++;
                        break;
                    case 'e':
                        stoppingTol = atof(argv[i+1]);
                        i++;
                        break;
                    case 'b':
                        boundHandlingMethod = atof(argv[i+1]);
                        i++;
                        break;
                }
            } else {                
                fileName = argv[i];
            }
        }
    } else {
        fprintf(stderr, "\nUsage: ./FilePrep <fileName> -n <maxDataSavePts> -e <stoppingTol> -b <boundHandlingMethod>\n");
        fprintf(stderr, "                  <fileName> Name of the NetCDF file to be prepared for simulation.\n");
        fprintf(stderr, "    -n      <maxDataSavePts> Maximum number of initial states to be held in the file.\n");
        fprintf(stderr, "    -e         <stoppingTol> Maximum allowable L2 distance between successive probability distributions.\n\n");
        fprintf(stderr, "    -b <boundHandlingMethod> 1 for deletion. 2 for resampling.");
        abort();
    }
        
    NcFile file(fileName.c_str(), NcFile::Write);
    if (!file.is_valid()) {
        fprintf(stderr, "Error: %s could not be opened.\n", fileName.c_str());
        abort();
    }
        
    NcDim* dim = file.add_dim("maxDataSavePts", maxDataSavePts);
    
    NcVar* stoppingTolVar = file.add_var("stoppingTol", ncDouble);
    NcVar* numDataSavePtsVar = file.add_var("numDataSavePts", ncInt, file.get_dim("numMdls")); 
    NcVar* boundHandlingMethodVar = file.add_var("boundHandlingMethod", ncInt, file.get_dim("numMdls"));
    file.add_var("numBoundedSpeciesStates", ncInt, file.get_dim("numMdls"));
    file.add_var("speciesStateBounded", ncInt, file.get_dim("numMdls"), file.get_dim("maxSpecies"));
    file.add_var("speciesStateLowerBounds", ncDouble, file.get_dim("numMdls"), file.get_dim("maxSpecies"));
    file.add_var("speciesStateUpperBounds", ncDouble, file.get_dim("numMdls"), file.get_dim("maxSpecies"));
    file.add_var("dataSavePts", ncInt, file.get_dim("numMdls"), dim);
    file.add_var("speciesStateChanges", ncInt, file.get_dim("numMdls"), file.get_dim("maxSpecies"));
    
    file.get_var("state")->rename("initFwdData");    
    file.add_var("initRevData", ncDouble, file.get_dim("numMdls"), file.get_dim("maxTrials"), file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    file.add_var("fwdData", ncDouble, file.get_dim("numMdls"), dim, file.get_dim("maxTrials"), file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    file.add_var("revData", ncDouble, file.get_dim("numMdls"), dim, file.get_dim("maxTrials"), file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    
    file.add_var("initFwdAbsCurr", ncInt, file.get_dim("numMdls"), file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    file.add_var("initRevAbsCurr", ncInt, file.get_dim("numMdls"), file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    file.add_var("fwdAbsCurr", ncInt, file.get_dim("numMdls"), dim, file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    file.add_var("revAbsCurr", ncInt, file.get_dim("numMdls"), dim, file.get_dim("maxTimePts"), file.get_dim("maxSaveSpecies"));
    
    stoppingTolVar->put(&stoppingTol, 1);
    
    int numMdls = (int) file.get_dim("numMdls")->size();    
    for (int i = 0; i < numMdls; i++) {
        numDataSavePtsVar->set_cur(i);
        numDataSavePtsVar->put(&maxDataSavePts, 1);
    }
    
    boundHandlingMethodVar->put(&boundHandlingMethod, 1);
}