<p><b>cahrens@lanl.gov</b> 2010-06-17 13:47:05 -0600 (Thu, 17 Jun 2010)</p><p>Version 1.2 of projlatlon: Revised to handle vectors (UReconstructZonal/Meridional)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp        2010-06-17 19:00:55 UTC (rev 358)
+++ branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp        2010-06-17 19:47:05 UTC (rev 359)
@@ -7,20 +7,20 @@
// Assume variable data type is double.
// Assume variable dims are (Time, nCells|nVertices, nVertLevels|nVertLevelsP1)
// Assume no more than 100 vars each for cell and point data
-// Does not deal with tracers.
+// Deals with vectors.
// Does not deal with edge data.
// Only vis up to nVertLevels, not nVertLevelsP1.
// Doubles converted to floats in .vtk files follow these rules:
-// if a NaN, then become -FLT_MAX.
-// if positive infinity, then become FLT_MAX
-// if negative infinity, then become -FLT_MAX
-// if smaller than a float, then become 0.
-// if not zero, but between -FLT_MIN and FLT_MIN make -FLT_MIN or FLT_MIN
-// if outside of -FLT_MAX and FLT_MAX make -FLT_MAX or FLT_MAX
+//        if a NaN, then become -FLT_MAX.
+//        if positive infinity, then become FLT_MAX
+//        if negative infinity, then become -FLT_MAX
+//        if smaller than a float, then become 0.
+//        if not zero, but between -FLT_MIN and FLT_MIN make -FLT_MIN or FLT_MIN
+//        if outside of -FLT_MAX and FLT_MAX make -FLT_MAX or FLT_MAX
//
-// Version 1.0
+// Version 1.2
// Christine Ahrens
-// 4/13/2010
+// 6/17/2010
////////////////////////////////////////////////////////////////////////////////
@@ -38,10 +38,10 @@
using namespace std;
#define CHECK_MALLOC(ptr) \
- if (ptr == NULL) { \
- cerr << "malloc failed!</font>
<font color="red">"; \
- exit(1); \
- }
+        if (ptr == NULL) { \
+                cerr << "malloc failed!</font>
<font color="gray">"; \
+                exit(1); \
+        }
#define MAX_VARS 100
#define DEFAULT_LAYER_THICKNESS 10
@@ -114,6 +114,15 @@
return (float)inputData;
}
+double convertXY2Lon(double x, double y, double lon) {
+ return ((-x * sin(lon)) + (y*cos(lon)));
+}
+
+double convertXYZ2Lat(double x, double y, double z, double lon, double lat) {
+ return ((((x*cos(lon)) + (y*sin(lon)))*sin(lat)) + (z*cos(lat)));
+}
+
+
void printUsage() {
cerr << "USAGE: projlatlon [OPTIONS] infile.nc outfile.vtk" << endl;
cerr << "Create lat/lon projection of MPAS-style data from input file infile.nc " << endl;
@@ -128,16 +137,28 @@
cerr << " -t floatval " << endl;
cerr << " layer_thickness -- use negative value for atmosphere, default is 0" << endl;
cerr << " -c floatval " << endl;
- cerr << " center_longitude in degrees, default is 180" << endl;
+ cerr << " display image with center_longitude in degrees, default is 180" << endl;
+ cerr << " -z " << endl;
+ cerr << " input grid has center longitude at 0, the default is 180 degree or PI radians" << endl;
+ cerr << " -a " << endl;
+ cerr << " atmosphere data (longitude data is in range -PI to PI) " << endl;
+ cerr << " The -a flag also converts the multilayer view to have the " << endl;
+ cerr << " vertical levels stack in the positive Z direction, rather " << endl;
+ cerr << " the negative Z direction, which is the default." << endl;
+ cerr << " -s " << endl;
+ cerr << " timestep to start on (defaults to 0) " << endl;
+ cerr << " -f " << endl;
+ cerr << " timestep to finish on (defaults to last timestep available) " << endl;
cerr << "Variables with time and vertical level are written out." << endl;
cerr << "Tracer vars are named tracer1, tracer2, etc." << endl;
- cerr << "[vtk datafile Version 2.0, projlatlon Version 1.0]" << endl;
+ cerr << "[vtk datafile Version 2.0, projlatlon Version 1.2]" << endl;
exit(1);
}
void getArgs(int argc, char* argv[], int &outputVertLevel,
float &layerThickness, double &centerLon, bool &isMultilayer,
- char* &inputFileName, char* &outputFileName) {
+ bool &isAtmosphere, int &startTimeStep, int &finishTimeStep,
+ bool &isZeroCentered, char* &inputFileName, char* &outputFileName) {
int c;
bool selectedVar = false;
@@ -147,11 +168,15 @@
bool tflag = false;
bool cflag = false;
bool mflag = false;
+ bool aflag = false;
+ bool sflag = false;
+ bool fflag = false;
+ bool zflag = false;
// -l vertical level (if not specified, assume multilayer)
// -t layer thickness (defaults to 1/10 * radius/#layers)
// -c center longitude in degrees for latlon projection
- while ((c=getopt(argc, argv, "l:t:c:m")) != -1) {
+ while ((c=getopt(argc, argv, "l:t:c:mazs:f:")) != -1) {
switch(c) {
case 'l':
lflag = true;
@@ -183,6 +208,22 @@
}
isMultilayer = true;
break;
+ case 'a':
+ aflag = true;
+ isAtmosphere = true;
+ break;
+ case 's':
+ sflag = true;
+ startTimeStep = atoi(optarg);
+ break;
+ case 'f':
+ sflag = true;
+ finishTimeStep = atoi(optarg);
+ break;
+ case 'z':
+ zflag = true;
+ isZeroCentered = true;
+ break;
default:
cerr << "Unrecognized option: -" << c << endl;
printUsage();
@@ -203,11 +244,15 @@
{
bool multilayer = false;
- int outputVertLevel = 0;
+ bool atmosphere = false;
+        int outputVertLevel = 0;
char* inputFileName = NULL;
char* outputFileName = NULL;
- float layerThickness = DEFAULT_LAYER_THICKNESS;
+        float layerThickness = DEFAULT_LAYER_THICKNESS;
double centerLon = 180;
+ int startTimeStep = 0;
+ int finishTimeStep = -1;
+ bool isZeroCentered = false;
if (argc < 3) {
printUsage();
@@ -215,30 +260,30 @@
}
getArgs(argc, argv, outputVertLevel, layerThickness, centerLon, multilayer,
- inputFileName, outputFileName);
-
- NcFile ncFile(inputFileName);
+ atmosphere, startTimeStep, finishTimeStep, isZeroCentered, inputFileName, outputFileName);
+        
+        NcFile ncFile(inputFileName);
- if (!ncFile.is_valid())
- {
- cerr << "Couldn't open file: " << inputFileName << endl;
- exit(1);
- }
+        if (!ncFile.is_valid())
+        {
+                cerr << "Couldn't open file: " << inputFileName << endl;
+                exit(1);
+        }
- NcDim* nCells = ncFile.get_dim("nCells");
- NcDim* nVertices = ncFile.get_dim("nVertices");
- NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
- NcDim* Time = ncFile.get_dim("Time");
- NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
+        NcDim* nCells = ncFile.get_dim("nCells");
+        NcDim* nVertices = ncFile.get_dim("nVertices");
+        NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
+        NcDim* Time = ncFile.get_dim("Time");
+        NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
int maxNVertLevels = nVertLevels->size();
- if ((vertexDegree->size() != 3) && (vertexDegree->size() != 4)) {
- cerr << "This code is only for hexagonal or quad grid projection" << endl;
- exit(1);
- }
+        if ((vertexDegree->size() != 3) && (vertexDegree->size() != 4)) {
+                cerr << "This code is only for hexagonal or quad grid projection" << endl;
+                exit(1);
+        }
- // Can't check for this, b/c if not there it crashes program
- //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
+        // Can't check for this, b/c if not there it crashes program        
+        //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
bool singlelayer = !multilayer;
@@ -267,110 +312,158 @@
maxNVertLevels = 1;
}
- cerr << "centerLon: " << centerLon << " singleLayer: " << singlelayer << " outputVertLevel: " << outputVertLevel << endl;
+ if (atmosphere) {
+ layerThickness = -layerThickness;
+ }
- // figure out what variables to visualize
- NcVar* dualCellVars[MAX_VARS];
- NcVar* dualPointVars[MAX_VARS];
- int dualCellVarIndex = -1;
- int dualPointVarIndex = -1;
- int numDualCells = nVertices->size();
- int numDualPoints = nCells->size()+1;
+ int availTimeSteps = Time->size();
+
+ if (startTimeStep > (availTimeSteps-1)) {
+ cerr << "Start timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
+
+ if (startTimeStep < 0) {
+ cerr << "Start timestep is too small, must be at least 0." << endl;
+ exit(1);
+ }
+
+ if (finishTimeStep > (availTimeSteps-1)) {
+ cerr << "Finish timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
- int numVars = ncFile.num_vars();
+ if (finishTimeStep < 0) {
+ finishTimeStep = availTimeSteps-1;
+ }
- bool tracersExist = false;
+ cerr << "centerLon: " << centerLon << " singleLayer: " << singlelayer << " outputVertLevel: " << outputVertLevel << " start: " << startTimeStep << " finish: " << finishTimeStep << " atmosphere: " << atmosphere << endl;
- for (int i = 0; i < numVars; i++) {
- NcVar* aVar = ncFile.get_var(i);
+        // figure out what variables to visualize
+        NcVar* dualCellVars[MAX_VARS];
+        NcVar* dualPointVars[MAX_VARS];
+ NcVar* vectorVars[3];
+        int dualCellVarIndex = -1;
+        int dualPointVarIndex = -1;
+        int numDualCells = nVertices->size();
+        int numDualPoints = nCells->size()+1;
+        
+        int numVars = ncFile.num_vars();
- // must have 3 dims
- // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
+        bool tracersExist = false;
+ bool isUVector = false;
+ bool isVVector = false;
- int numDims = aVar->num_dims();
- //cout << "Num Dims of var: " << aVar->name() << " is " << numDims << endl;
- if ((numDims != 3) && (strcmp(aVar->name(), "tracers"))) {
- continue; // try the next var
- } else {
- // TODO, check if it is a double
- // assume a double for now
+        for (int i = 0; i < numVars; i++) {
+                NcVar* aVar = ncFile.get_var(i);
- // check for Time dim 0
- NcToken dim0Name = aVar->get_dim(0)->name();
- if (strcmp(dim0Name, "Time"))
- continue;
+                // must have 3 dims
+                // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
- // check for dim 1 being num vertices or cells
- bool isVertexData = false;
- bool isCellData = false;
- NcToken dim1Name = aVar->get_dim(1)->name();
- if (!strcmp(dim1Name, "nVertices"))
- isVertexData = true;
- else if (!strcmp(dim1Name, "nCells"))
- isCellData = true;
- else continue;
+                int numDims = aVar->num_dims();
+                //cout << "Num Dims of var: " << aVar->name() << " is " << numDims << endl;
+                if ((numDims != 3) && (strcmp(aVar->name(), "tracers"))) {
+                        continue; // try the next var
+                } else {
+                        // TODO, check if it is a double
+                        // assume a double for now
- // check if dim 2 is nVertLevels or nVertLevelsP1, too
- NcToken dim2Name = aVar->get_dim(2)->name();
- if ((strcmp(dim2Name, "nVertLevels"))
- && (strcmp(dim2Name, "nVertLevelsP1"))) {
- continue;
- }
+                        // check for Time dim 0
+                        NcToken dim0Name = aVar->get_dim(0)->name();
+                        if (strcmp(dim0Name, "Time"))
+                                continue;
- // Add to cell or point var array
- if (isVertexData) { // means it is dual cell data
- dualCellVarIndex++;
- if (dualCellVarIndex > MAX_VARS-1) {
- cerr << "Exceeded number of cell vars." << endl;
- exit(1);
- }
- dualCellVars[dualCellVarIndex] = aVar;
- //cout << "Adding var " << aVar->name() << " to dualCellVars" << endl;
- } else if (isCellData) { // means it is dual vertex data
- if (strcmp(aVar->name(), "tracers")) {
- dualPointVarIndex++;
- if (dualPointVarIndex > MAX_VARS-1) {
- cerr << "Exceeded number of point vars." << endl;
- exit(1);
+                        // check for dim 1 being num vertices or cells
+                        bool isVertexData = false;
+                        bool isCellData = false;
+                        NcToken dim1Name = aVar->get_dim(1)->name();
+                        if (!strcmp(dim1Name, "nVertices"))
+                                isVertexData = true;
+                        else if (!strcmp(dim1Name, "nCells"))
+                                isCellData = true;
+                        else continue;
+
+                        // check if dim 2 is nVertLevels or nVertLevelsP1, too
+                        NcToken dim2Name = aVar->get_dim(2)->name();
+                        if ((strcmp(dim2Name, "nVertLevels"))
+                                        && (strcmp(dim2Name, "nVertLevelsP1"))) {
+                                continue;
+                        }
+
+                        // Add to cell or point var array
+                        if (isVertexData) { // means it is dual cell data
+                                dualCellVarIndex++;
+                                if (dualCellVarIndex > MAX_VARS-1) {
+                                        cerr << "Exceeded number of cell vars." << endl;
+                                        exit(1);
+                                }
+                                dualCellVars[dualCellVarIndex] = aVar;
+                                //cout << "Adding var " << aVar->name() << " to dualCellVars" << endl;
+                        } else if (isCellData) { // means it is dual vertex data
+ // check for u/v/w vectors
+
+ if (!strncmp(aVar->name(), "uReconstruct", 12)) {
+ if (!strcmp(aVar->name(), "uReconstructZonal")) {
+ vectorVars[0] = aVar;
+ isUVector = true;
+ } else if (!strcmp(aVar->name(), "uReconstructMeridional")) {
+ vectorVars[1] = aVar;
+ isVVector = true;
}
- dualPointVars[dualPointVarIndex] = aVar;
- //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
- } else { // case of tracers, add each as "tracer0", "tracer1", etc.
- tracersExist = true;
- int numTracers = aVar->get_dim(3)->size();
- for (int t = 0; t < numTracers; t++) {
- dualPointVarIndex++;
- if (dualPointVarIndex > MAX_VARS-1) {
- cerr << "Exceeded number of point vars." << endl;
- exit(1);
- }
- dualPointVars[dualPointVarIndex] = aVar;
- //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
- }
- }
- }
- }
- }
+ } else if (strcmp(aVar->name(), "tracers")) {
+                                        dualPointVarIndex++;
+                                        if (dualPointVarIndex > MAX_VARS-1) {
+                                                cerr << "Exceeded number of point vars." << endl;
+                                                exit(1);
+                                        }
+                                        dualPointVars[dualPointVarIndex] = aVar;
+                                        //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
+                                } else { // case of tracers, add each as "tracer0", "tracer1", etc.
+                                        tracersExist = true;
+                                        int numTracers = aVar->get_dim(3)->size();
+                                        for (int t = 0; t < numTracers; t++) {
+                                                dualPointVarIndex++;
+                                                if (dualPointVarIndex > MAX_VARS-1) {
+                                                        cerr << "Exceeded number of point vars." << endl;
+                                                        exit(1);
+                                                }
+                                                dualPointVars[dualPointVarIndex] = aVar;
+                                                //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
+                                        }
+                                }
+                        }
+                }
+        }
- // TODO
- // prompt the user to find out which fields are of interest?
- // for now, assume all are of interest
-
- // get points (centers of primal-mesh cells)
+        // TODO
+         // prompt the user to find out which fields are of interest?
+        // for now, assume all are of interest
+        
+        // get points (centers of primal-mesh cells)
- // TO DO check malloc return vals.
+        // TO DO check malloc return vals.
const float BLOATFACTOR = .5;
- int myCellSize = (int)floor(nCells->size()*(1.0 + BLOATFACTOR));
- double *xCellData = (double*)malloc(myCellSize * sizeof(double));
- CHECK_MALLOC(xCellData);
- NcVar *xCellVar = ncFile.get_var("lonCell");
- xCellVar->get(xCellData+1, nCells->size());
+        int myCellSize = (int)floor(nCells->size()*(1.0 + BLOATFACTOR));
+        double *xCellData = (double*)malloc(myCellSize * sizeof(double));
+        CHECK_MALLOC(xCellData);
+        NcVar *xCellVar = ncFile.get_var("lonCell");
+        xCellVar->get(xCellData+1, nCells->size());
// point 0 is 0.0
*xCellData = 0.0;
const double PI = 3.141592;
+ // if atmospheric data, or zero centered, set center to 180 instead of 0
+ if (atmosphere || isZeroCentered) {
+ for (int j = 1; j <= nCells->size(); j++) {
+ // need to shift over the point so center is at PI
+ if (xCellData[j] < 0) {
+ xCellData[j] += 2*PI;
+ }
+ }
+ }
+
double centerRad = centerLon * PI / 180.0;
if (centerLon != 180) {
@@ -388,51 +481,51 @@
}
}
- double *yCellData = (double*)malloc(myCellSize * sizeof(double));
- CHECK_MALLOC(yCellData);
- NcVar *yCellVar = ncFile.get_var("latCell");
- yCellVar->get(yCellData+1, nCells->size());
+        double *yCellData = (double*)malloc(myCellSize * sizeof(double));
+        CHECK_MALLOC(yCellData);
+        NcVar *yCellVar = ncFile.get_var("latCell");
+        yCellVar->get(yCellData+1, nCells->size());
// point 0 is 0.0
*yCellData = 0.0;
- // get dual-mesh cells
+        // get dual-mesh cells
- int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
- sizeof(int));
- CHECK_MALLOC(cellsOnVertex);
+        int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+                sizeof(int));
+        CHECK_MALLOC(cellsOnVertex);
- int myCOVSize = (int)floor(nVertices->size()*(1.0 + BLOATFACTOR))+1;
- int *myCellsOnVertex = (int *) malloc(myCOVSize * vertexDegree->size() * sizeof(int));
- CHECK_MALLOC(myCellsOnVertex);
+        int myCOVSize = (int)floor(nVertices->size()*(1.0 + BLOATFACTOR))+1;
+        int *myCellsOnVertex = (int *) malloc(myCOVSize * vertexDegree->size() * sizeof(int));
+        CHECK_MALLOC(myCellsOnVertex);
- NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
- //cout << "getting cellsOnVertexVar</font>
<font color="red">";
- cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+        NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
+        //cout << "getting cellsOnVertexVar</font>
<font color="gray">";
+        cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
- // allocate an array to map the extra points and cells to the original
- // so that when obtaining data, we know where to get it
- int *pointMap = (int*)malloc((int)floor(nCells->size()*BLOATFACTOR) * sizeof(int));
- CHECK_MALLOC(pointMap);
- int *cellMap = (int*)malloc((int)floor(nVertices->size()*BLOATFACTOR) * sizeof(int));
- CHECK_MALLOC(cellMap);
+        // allocate an array to map the extra points and cells to the original
+        // so that when obtaining data, we know where to get it
+        int *pointMap = (int*)malloc((int)floor(nCells->size()*BLOATFACTOR) * sizeof(int));
+        CHECK_MALLOC(pointMap);
+        int *cellMap = (int*)malloc((int)floor(nVertices->size()*BLOATFACTOR) * sizeof(int));
+        CHECK_MALLOC(cellMap);
-
- int *myptr = myCellsOnVertex;
- int *myExtraX;
- int *myExtraY;
- int currentExtraPoint = numDualPoints;
- int currentExtraCell = numDualCells;
+        
+        int *myptr = myCellsOnVertex;
+        int *myExtraX;
+        int *myExtraY;
+        int currentExtraPoint = numDualPoints;
+        int currentExtraCell = numDualCells;
// For each cell, examine vertices
- // Add new points and cells where needed to account for wraparound.
+        // Add new points and cells where needed to account for wraparound.
int acell = 775;
int mirrorcell;
int apoint;
int mirrorpoint;
- for (int j = 0; j < numDualCells; j++ ) {
- int *dualCells = cellsOnVertex + (j * vertexDegree->size());
+        for (int j = 0; j < numDualCells; j++ ) {
+                int *dualCells = cellsOnVertex + (j * vertexDegree->size());
// go through and make sure none of the referenced points are
// out of range (<=0 or >nCells->size())
@@ -446,65 +539,89 @@
}
}
- int lastk = vertexDegree->size()-1;
- bool xWrap = false;
- bool yWrap = false;
- for (int k = 0; k < vertexDegree->size(); k++) {
- if (abs(xCellData[dualCells[k]] - xCellData[dualCells[lastk]]) > 5.5) xWrap = true;
- //if (abs(yCellData[dualCells[k]] - yCellData[dualCells[lastk]]) > 2) yWrap = true;
- lastk = k;
- }
+ int lastk;
- if (xWrap || yWrap) {
- //cerr << "Cell wrap: " << j << endl;
+ /////BUG FIX for problem where cells are stretching to a faraway point
+                lastk = vertexDegree->size()-1;
+ //const double thresh = .139626340159546; // 8 degrees
+ const double thresh = .06981317007977; // 8 degrees
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ double xdiff = abs(xCellData[dualCells[k]] - xCellData[dualCells[lastk]]);
+ double ydiff = abs(yCellData[dualCells[k]] - yCellData[dualCells[lastk]]);
+ // Don't look at cells at map border
+                        //if ((xdiff <= 5.5) && (xdiff > (2*.06981317))) {
+                        //if (((xdiff <= 6) && (xdiff > 3)) || (ydiff > .06981317)) {
+                        if (((xdiff > thresh) && (ydiff > thresh))
+ // || ((xdiff <= 6) && (xdiff > 3))
+ || (ydiff > thresh)
+ ) {
+ for (int m = 0; m < vertexDegree->size(); m++) {
+ dualCells[m] = 0;
+ }
+ break;
+ }
}
- if (xWrap) {
- //cerr << "It wrapped in x direction" << endl;
- double anchorX = xCellData[dualCells[0]];
- double anchorY = yCellData[dualCells[0]];
- *myptr = *(dualCells);
- myptr++;
+                lastk = vertexDegree->size()-1;
+                bool xWrap = false;
+                bool yWrap = false;
+                for (int k = 0; k < vertexDegree->size(); k++) {
+                        if (abs(xCellData[dualCells[k]] - xCellData[dualCells[lastk]]) > 5.5) xWrap = true;
+                        //if (abs(yCellData[dualCells[k]] - yCellData[dualCells[lastk]]) > 2) yWrap = true;
+                        lastk = k;
+                }
+
+                if (xWrap || yWrap) {
+                        //cerr << "Cell wrap: " << j << endl;
+                }
+
+                if (xWrap) {
+                        //cerr << "It wrapped in x direction" << endl;
+                        double anchorX = xCellData[dualCells[0]];
+                        double anchorY = yCellData[dualCells[0]];
+                        *myptr = *(dualCells);
+                        myptr++;
+
if (j == acell) {
//cout << "cell " << acell << " anchor x: " << anchorX << " y: " << anchorY << endl;
}
- // modify existing cell, so it doesn't wrap
- // move points to one side
+                        // modify existing cell, so it doesn't wrap
+                        // move points to one side
- // first point is anchor it doesn't move
+                        // first point is anchor it doesn't move
for (int k = 1; k < vertexDegree->size(); k++) {
- double neighX = xCellData[dualCells[k]];
- double neighY = yCellData[dualCells[k]];
+                                double neighX = xCellData[dualCells[k]];
+                                double neighY = yCellData[dualCells[k]];
if (j == acell) {
//cout << "cell " << acell << " k: " << k << " x: " << neighX << " y: " << neighY << endl;
}
- // add a new point, figure out east or west
- if (abs(neighX - anchorX) > 5.5) {
- double neighEastX;
- double neighWestX;
+                                // add a new point, figure out east or west
+                                if (abs(neighX - anchorX) > 5.5) {
+                                        double neighEastX;
+                                        double neighWestX;
- // add on east
- if (neighX < anchorX) {
+                                        // add on east
+                                        if (neighX < anchorX) {
if (j == acell) {
//cout << "add on east" << endl;
}
- neighEastX = neighX + (2*PI);
- xCellData[currentExtraPoint] = neighEastX;
- yCellData[currentExtraPoint] = neighY;
- } else {
- // add on west
+                                                neighEastX = neighX + (2*PI);
+                                                xCellData[currentExtraPoint] = neighEastX;
+                                                yCellData[currentExtraPoint] = neighY;
+                                        } else {
+                                                // add on west
if (j == acell) {
//cout << "add on west" << endl;
}
- neighWestX = neighX - (2*PI);
- xCellData[currentExtraPoint] = neighWestX;
- yCellData[currentExtraPoint] = neighY;
- }
-
+                                                neighWestX = neighX - (2*PI);
+                                                xCellData[currentExtraPoint] = neighWestX;
+                                                yCellData[currentExtraPoint] = neighY;
+                                        }
+                        
if (j == acell) {
//cout << "x: " << xCellData[currentExtraPoint] << " y: " << yCellData[currentExtraPoint] << endl;
}
@@ -513,11 +630,11 @@
//cout << "currentExtraPoint: " << currentExtraPoint << endl;
}
- // add the new point to list of vertices
- *myptr = currentExtraPoint;
+                                        // add the new point to list of vertices
+                                        *myptr = currentExtraPoint;
- // record mapping
- *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+                                        // record mapping
+                                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
if (j == acell) {
mirrorpoint = currentExtraPoint;
@@ -525,28 +642,28 @@
//cout << "mirror point " << mirrorpoint << " has pointMap to: " << dualCells[k] << endl;
}
- myptr++;
- currentExtraPoint++;
- } else {
- // use existing kth point
+                                        myptr++;
+                                        currentExtraPoint++;
+                                } else {
+                                        // use existing kth point
if (j == acell) {
//cout << "use existing point" << endl;
}
- *myptr = dualCells[k];
- myptr++;
- }
- }
+                                        *myptr = dualCells[k];
+                                        myptr++;
+                                }
+                        }
- // add a mirror image cell on other side so there are no
- // gaps in the map
+                        // add a mirror image cell on other side so there are no
+                        // gaps in the map
- // move anchor to other side
- if (anchorX > centerRad) {
- anchorX = anchorX - (2*PI);
- } else {
- anchorX = anchorX + (2*PI);
- }
-
+                        // move anchor to other side
+                        if (anchorX > centerRad) {
+                                anchorX = anchorX - (2*PI);
+                        } else {
+                                anchorX = anchorX + (2*PI);
+                        }
+        
if (j == acell) {
//cout << "add new cell " << currentExtraCell << " to mirror " << acell << " anchorX: " << anchorX << " anchorY: " << anchorY << endl;
mirrorcell = currentExtraCell;
@@ -556,110 +673,110 @@
int* addedCellsPtr = myCellsOnVertex + (currentExtraCell * vertexDegree->size());
// add point coord and add to list of cells
- xCellData[currentExtraPoint] = anchorX;
- yCellData[currentExtraPoint] = anchorY;
- *addedCellsPtr = currentExtraPoint;
+                        xCellData[currentExtraPoint] = anchorX;
+                        yCellData[currentExtraPoint] = anchorY;
+                        *addedCellsPtr = currentExtraPoint;
- // record mapping
- *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[0];
+                        // record mapping
+                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[0];
- addedCellsPtr++;
- currentExtraPoint++;
+                        addedCellsPtr++;
+                        currentExtraPoint++;
- for (int k = 1; k < vertexDegree->size(); k++) {
- double neighX = xCellData[dualCells[k]];
- double neighY = yCellData[dualCells[k]];
+                        for (int k = 1; k < vertexDegree->size(); k++) {
+                                double neighX = xCellData[dualCells[k]];
+                                double neighY = yCellData[dualCells[k]];
- // add a new point, figure out east or west
- if (abs(neighX - anchorX) > 5.5) {
- double neighEastX;
- double neighWestX;
+                                // add a new point, figure out east or west
+                                if (abs(neighX - anchorX) > 5.5) {
+                                        double neighEastX;
+                                        double neighWestX;
- // add on east
- if (neighX < anchorX) {
- neighEastX = neighX + (2*PI);
- xCellData[currentExtraPoint] = neighEastX;
- yCellData[currentExtraPoint] = neighY;
- } else {
- // add on west
- neighWestX = neighX - (2*PI);
- xCellData[currentExtraPoint] = neighWestX;
- yCellData[currentExtraPoint] = neighY;
- }
+                                        // add on east
+                                        if (neighX < anchorX) {
+                                                neighEastX = neighX + (2*PI);
+                                                xCellData[currentExtraPoint] = neighEastX;
+                                                yCellData[currentExtraPoint] = neighY;
+                                        } else {
+                                                // add on west
+                                                neighWestX = neighX - (2*PI);
+                                                xCellData[currentExtraPoint] = neighWestX;
+                                                yCellData[currentExtraPoint] = neighY;
+                                        }
- // add the new point to list of vertices
- *addedCellsPtr = currentExtraPoint;
+                                        // add the new point to list of vertices
+                                        *addedCellsPtr = currentExtraPoint;
- // record mapping
- *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+                                        // record mapping
+                                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
- addedCellsPtr++;
- currentExtraPoint++;
- } else {
- // use existing kth point
- *addedCellsPtr = dualCells[k];
- addedCellsPtr++;
- }
- }
- *(cellMap + (currentExtraCell - numDualCells)) = j;
- currentExtraCell++;
- }
- if (yWrap) {
- //cerr << "It wrapped in y direction" << endl;
- }
+                                        addedCellsPtr++;
+                                        currentExtraPoint++;
+                                } else {
+                                        // use existing kth point
+                                        *addedCellsPtr = dualCells[k];
+                                        addedCellsPtr++;
+                                }
+                        }
+                        *(cellMap + (currentExtraCell - numDualCells)) = j;
+                        currentExtraCell++;
+                }
+                if (yWrap) {
+                        //cerr << "It wrapped in y direction" << endl;
+                }
- // if cell doesn't extend past lat/lon perimeter, then add it to myCellsOnVertex
- if (!xWrap && !yWrap) {
- for (int k=0; k< vertexDegree->size(); k++) {
- *myptr = *(dualCells+k);
- myptr++;
- }
- }
- if (currentExtraCell > myCOVSize) {
- cerr << "Exceeded storage for extra cells!" << endl;
- return 1;
- }
- if (currentExtraPoint > myCellSize) {
- cerr << "Exceeded storage for extra points!" << endl;
- return 1;
- }
- }
-
- // decls for data storage
- double* dualCellVarData;
- double* dualPointVarData;
+                // if cell doesn't extend past lat/lon perimeter, then add it to myCellsOnVertex
+                if (!xWrap && !yWrap) {
+                        for (int k=0; k< vertexDegree->size(); k++) {
+                                *myptr = *(dualCells+k);
+                                myptr++;
+                        }
+                }
+                if (currentExtraCell > myCOVSize) {
+                        cerr << "Exceeded storage for extra cells!" << endl;
+                        return 1;
+                }
+                if (currentExtraPoint > myCellSize) {
+                        cerr << "Exceeded storage for extra points!" << endl;
+                        return 1;
+                }
+        }        
+                                                
+        // decls for data storage
+        double* dualCellVarData;
+        double* dualPointVarData;
- // for each variable, allocate space for variables
+        // for each variable, allocate space for variables
- //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
- int varVertLevels = 0;
-
- // write a file with the geometry.
+        int varVertLevels = 0;
+        
+        // write a file with the geometry.
- ostringstream geoFileName;
+        ostringstream geoFileName;
- geoFileName << "geo_" << outputFileName;
+        geoFileName << "geo_" << outputFileName;         
- ofstream geoFile(geoFileName.str().c_str(), ios::out);
+        ofstream geoFile(geoFileName.str().c_str(), ios::out);
- if (!geoFile)
- {
- cerr << "vtk output file could not be opened" <<endl;
- exit(1);
- }
+        if (!geoFile)
+        {
+                cerr << "vtk output file could not be opened" <<endl;
+                exit(1);
+        }
- // write header
+        // write header
- geoFile << "# vtk DataFile Version 2.0" << endl;
- geoFile << "Project newpop geometry to lat/lon projection from netCDF by Christine Ahrens"
- << endl;
- geoFile << "ASCII" << endl;
- geoFile << "DATASET UNSTRUCTURED_GRID" << endl;
+        geoFile << "# vtk DataFile Version 2.0, projlatlon Version 1.2" << endl;
+        geoFile << "Project newpop geometry to lat/lon projection from netCDF by Christine Ahrens"
+                << endl;
+        geoFile << "ASCII" << endl;
+        geoFile << "DATASET UNSTRUCTURED_GRID" << endl;
- // write points (the points are the primal-grid cell centers)
+        // write points (the points are the primal-grid cell centers)
if (singlelayer) {
geoFile << "POINTS " << currentExtraPoint << " float" << endl;
@@ -667,69 +784,69 @@
geoFile << "POINTS " << currentExtraPoint*(maxNVertLevels+1) << " float" << endl;
}
- // write the point at each vertical level, plus one level for last layer
+        // write the point at each vertical level, plus one level for last layer
- //cout << "Writing points at each level" << endl;
+        //cout << "Writing points at each level" << endl;
- for (int j = 0; j < currentExtraPoint; j++ )
- {
- geoFile.precision(16);
+        for (int j = 0; j < currentExtraPoint; j++ )
+        {
+                geoFile.precision(16);
- if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
- if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
+                if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
+                if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
- float xLon = xCellData[j] * 180.0 / PI;
- float yLon = yCellData[j] * 180.0 / PI;
+                float xLon = xCellData[j] * 180.0 / PI;
+                float yLat = yCellData[j] * 180.0 / PI;
if (singlelayer) {
- geoFile << xLon << "\t" << yLon << "\t" << 0 << endl;
+                        geoFile << xLon << "\t" << yLat << "\t" << 0 << endl;
} else {
for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
- geoFile << xLon << "\t" << yLon << "\t" << -(((float)levelNum)*layerThickness) << endl;
+ geoFile << xLon << "\t" << yLat << "\t" << -(((float)levelNum)*layerThickness) << endl;
}
}
- }
+        }        
- geoFile << endl;
+                geoFile << endl;
- // Write dual-mesh cells
- // Dual-mesh cells are triangles with primal-mesh cell
- // centers as the vertices.
- // The number of dual-mesh cells is the number of vertices in the
- // primal mesh.
+                // Write dual-mesh cells
+                // Dual-mesh cells are triangles with primal-mesh cell
+                // centers as the vertices.
+                // The number of dual-mesh cells is the number of vertices in the
+                // primal mesh.
- int degree = multilayer?vertexDegree->size()*2:vertexDegree->size();
+                int degree = multilayer?vertexDegree->size()*2:vertexDegree->size();
if (singlelayer) {
geoFile << "CELLS " << currentExtraCell << " "
- << currentExtraCell * (degree + 1) << endl;
+ << currentExtraCell * (degree + 1) << endl;        
} else {
geoFile << "CELLS " << currentExtraCell*maxNVertLevels << " "
- << currentExtraCell * (degree + 1) * maxNVertLevels << endl;
+ << currentExtraCell * (degree + 1) * maxNVertLevels << endl;        
}
- // for each dual-mesh cell, write number of points for each
- // and then list the points by number
+                // for each dual-mesh cell, write number of points for each
+                // and then list the points by number
- //cout << "Writing Cells" << endl;
+                //cout << "Writing Cells" << endl;
- for (int j = 0; j < currentExtraCell ; j++) {
+                for (int j = 0; j < currentExtraCell ; j++) {
- // since primal vertex(pt) numbers == dual cell numbers
- // we go through the primal vertices, find the cells around
- // them, and since those primal cell numbers are dual
- // point numbers,
- // we can write the cell numbers for the cellsOnVertex
- // and those will be the numbers of the dual vertices (pts).
-
- int* dualCells = myCellsOnVertex + (j * vertexDegree->size());
+                        // since primal vertex(pt) numbers == dual cell numbers
+                        // we go through the primal vertices, find the cells around
+                        // them, and since those primal cell numbers are dual
+                        // point numbers,
+                        // we can write the cell numbers for the cellsOnVertex
+                        // and those will be the numbers of the dual vertices (pts).
+                        
+                        int* dualCells = myCellsOnVertex + (j * vertexDegree->size());
- // for each level, write the cell
+                        // for each level, write the cell
if (singlelayer) {
- geoFile << degree << "\t" ;
- for (int k = 0; k < vertexDegree->size(); k++)
- geoFile << dualCells[k] << "\t";
+                                geoFile << degree << "\t" ;
+                                for (int k = 0; k < vertexDegree->size(); k++)
+                                        geoFile << dualCells[k] << "\t";
geoFile << endl;
} else {
for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
@@ -737,20 +854,20 @@
geoFile << degree << "\t" ;
for (int k = 0; k < vertexDegree->size(); k++)
- {
+ {                
geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum << "\t";
}
for (int k = 0; k < vertexDegree->size(); k++)
- {
+ {                
geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum+1 << "\t";
}
geoFile << endl;
}
}
- }
+                }
- // write cell types
- int cellType;
+                // write cell types
+                int cellType;
switch (vertexDegree->size()) {
case 3:
if (singlelayer) {
@@ -770,7 +887,7 @@
break;
}
- if (singlelayer) {
+                if (singlelayer) {
geoFile << "CELL_TYPES " << currentExtraCell << endl;
} else {
geoFile << "CELL_TYPES " << currentExtraCell*maxNVertLevels << endl;
@@ -789,45 +906,58 @@
}
}
- // release resources
- geoFile.close();
- free(xCellData);
- free(yCellData);
- free(cellsOnVertex);
- free(myCellsOnVertex);
+        // release resources
+        geoFile.close();
+        //free(xCellData);
+        //free(yCellData);
+        free(cellsOnVertex);
+        free(myCellsOnVertex);
- // For each timestep, write data for each level
+        // For each timestep, write data for each level
- dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
- CHECK_MALLOC(dualCellVarData);
+        dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
+        CHECK_MALLOC(dualCellVarData);
- dualPointVarData = (double*)malloc((sizeof(double)) * currentExtraPoint * maxNVertLevels);
- CHECK_MALLOC(dualPointVarData);
+        dualPointVarData = (double*)malloc((sizeof(double)) * currentExtraPoint * maxNVertLevels);
+        CHECK_MALLOC(dualPointVarData);
- // for each timestep, copy the geometry file, since we don't want to
- // have to recompute the points
- for (int i = 0; i < Time->size(); i++) {
- ostringstream vtkFileName;
- vtkFileName << i << outputFileName;
+ double *vectorData[3] = {NULL, NULL, NULL};
- string geoFileNameString = geoFileName.str();
- string vtkFileNameString = vtkFileName.str();
- int copyRetVal = myCopyFile(&geoFileNameString, &vtkFileNameString);
- //cout << "myCopyFile returned: " << copyRetVal << endl;
-
- ofstream vtkFile(vtkFileName.str().c_str(), ios::out|ios::app);
- if (!vtkFile)
- {
- cerr << "vtk output file could not be opened" <<endl;
- exit(1);
+ if (isUVector) {
+ for (int d = 0; d < 3; d++) {
+ if (singlelayer) {
+ vectorData[d] = (double*)malloc((sizeof(double)) * (nCells->size() + 1));
+ } else {
+ vectorData[d] = (double*)malloc((sizeof(double)) * (nCells->size() + 1) * maxNVertLevels);
+ }
+ CHECK_MALLOC(vectorData[d]);
}
+ }
- vtkFile.precision(16);
+        // for each timestep, copy the geometry file, since we don't want to
+        // have to recompute the points
+        for (int i = startTimeStep; i <= finishTimeStep; i++) {
+                ostringstream vtkFileName;
+                vtkFileName << i << outputFileName;
+                string geoFileNameString = geoFileName.str();
+                string vtkFileNameString = vtkFileName.str();
+                int copyRetVal = myCopyFile(&geoFileNameString, &vtkFileNameString);
+                //cout << "myCopyFile returned: " << copyRetVal << endl;
+
+                ofstream vtkFile(vtkFileName.str().c_str(), ios::out|ios::app);
+                if (!vtkFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+                vtkFile.precision(16);
+
int vertLevelIndex = 0;
if (singlelayer) vertLevelIndex = outputVertLevel;
- // If by point, write out point data
+                // If by point, write out point data
if (singlelayer) {
if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint << endl;
@@ -835,62 +965,62 @@
if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint*(maxNVertLevels+1) << endl;
}
- int printstep = -1;
+                int printstep = -1;
- if (i == printstep) cout << "TIME STEP: " << i << endl;
+                if (i == printstep) cout << "TIME STEP: " << i << endl;
- int tracerNum = 0;
+                int tracerNum = 0;
- for (int v = 0; v <= dualPointVarIndex; v++) {
+                for (int v = 0; v <= dualPointVarIndex; v++) {
- // Read variable number v data for that timestep
+                        // Read variable number v data for that timestep
- varVertLevels = dualPointVars[v]->get_dim(2)->size();
+                        varVertLevels = dualPointVars[v]->get_dim(2)->size();
- bool isTracer = false;
+                        bool isTracer = false;
- if (!strcmp(dualPointVars[v]->name(), "tracers")) {
- isTracer = true;
- // Uncomment if want to exclude tracers.
- // continue;
- }
+                        if (!strcmp(dualPointVars[v]->name(), "tracers")) {
+                                isTracer = true;
+                        // Uncomment if want to exclude tracers.
+                        //        continue;
+                        }
- // Write variable number v data for that timestep
- vtkFile << "SCALARS " << dualPointVars[v]->name();
- if (isTracer) vtkFile << tracerNum+1;
- vtkFile << " float 1" << endl;
- vtkFile << "LOOKUP_TABLE default" << endl;
+                        // Write variable number v data for that timestep
+                        vtkFile << "SCALARS " << dualPointVars[v]->name();
+                        if (isTracer) vtkFile << tracerNum+1;
+                        vtkFile << " float 1" << endl;
+                        vtkFile << "LOOKUP_TABLE default" << endl;
- if (isTracer) {
- dualPointVars[v]->set_cur(i, 0, vertLevelIndex, tracerNum);
+                        if (isTracer) {
+                                dualPointVars[v]->set_cur(i, 0, vertLevelIndex, tracerNum);
if (singlelayer) {
dualPointVars[v]->get(dualPointVarData+1, 1, nCells->size(), 1, 1);
} else {
dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels, 1);
}
- } else {
- dualPointVars[v]->set_cur(i, 0, vertLevelIndex);
+                        } else {
+                                dualPointVars[v]->set_cur(i, 0, vertLevelIndex);
if (singlelayer) {
dualPointVars[v]->get(dualPointVarData+1, 1, nCells->size(), 1);
} else {
dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels);
}
- }
+                        }
- float defaultPointVal = 0.0;
+                        float defaultPointVal = 0.0;
- double *var_target;
- float validData;
+                        double *var_target;
+                        float validData;
- //write dummy
+                        //write dummy
if (singlelayer) {
- var_target = dualPointVarData + 1;
- validData = convertDouble2ValidFloat (*var_target);
+                         var_target = dualPointVarData + 1;
+                                validData = convertDouble2ValidFloat (*var_target);
- // write dummy
- vtkFile << validData << endl;
- } else {
+                                // write dummy
+                                vtkFile << validData << endl;
+                        } else {
//write dummy
var_target = dualPointVarData + maxNVertLevels;
for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
@@ -908,12 +1038,12 @@
var_target = dualPointVarData + maxNVertLevels;
}
- for (int j = 1; j < numDualPoints; j++) {
+                        for (int j = 1; j < numDualPoints; j++) {
if (singlelayer) {
- validData = convertDouble2ValidFloat (*var_target);
- vtkFile << validData << endl;
- var_target++;
+                                        validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+                                        var_target++;
} else {
// write data for one point -- lowest level to highest
@@ -927,20 +1057,20 @@
// Need Mark's input on this one
vtkFile << validData << endl;
}
- }
+                        }
- // put out data for extra points
- for (int j = numDualPoints; j < currentExtraPoint; j++) {
+                        // put out data for extra points
+                        for (int j = numDualPoints; j < currentExtraPoint; j++) {
if (j == mirrorpoint) {
//cout << "data for mirror point " << mirrorpoint << " from point " << *(pointMap + j - numDualPoints) << endl;
}
// use map to find out what point data we are using
- var_target = dualPointVarData + ((*(pointMap + j - numDualPoints))*maxNVertLevels);
+                                var_target = dualPointVarData + ((*(pointMap + j - numDualPoints))*maxNVertLevels);
if (singlelayer) {
- validData = convertDouble2ValidFloat (*var_target);
- vtkFile << validData << endl;
- var_target++;
+                                        validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+                                        var_target++;
} else {
// write data for one point -- lowest level to highest
for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
@@ -953,11 +1083,154 @@
// Need Mark's input on this one
vtkFile << validData << endl;
}
+                        }
+                        if (isTracer) tracerNum++;
+                }
+
+ float validData;
+ int maxDim = 2; // z is always 0 for vectors in lat/lon projection
+
+ if (isUVector) {
+ cout << "isUVector:" << endl;
+
+ for (int d = 0; d < maxDim; d++) {
+ // read data in
+ vectorVars[d]->set_cur(i, 0, vertLevelIndex);
+ if (singlelayer) {
+ vectorVars[d]->get(vectorData[d]+1, 1, nCells->size(), 1);
+ } else {
+ vectorVars[d]->get(vectorData[d]+maxNVertLevels, 1, nCells->size(), maxNVertLevels);
+ }
}
- if (isTracer) tracerNum++;
+
+#ifdef NO_NEED_TO_DO_THIS
+ // Convert all vector data to lat/lon
+ if (singlelayer) {
+ for (int j = 1; j <= nCells->size(); j++) {
+ double x = *(vectorData[0] + j);
+ double y = *(vectorData[1] + j);
+ double z = *(vectorData[2] + j);
+                 float lon = xCellData[j];
+                 float lat = yCellData[j];
+ *(vectorData[0] + j) = convertXY2Lon(x, y, lon);
+ *(vectorData[1] + j) = convertXYZ2Lat(x, y, z, lon, lat);
+ }
+ } else {
+ for (int j = 1; j <= nCells->size(); j++) {
+ float lon = xCellData[j];
+ float lat = yCellData[j];
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ double x = *(vectorData[0] + (j*maxNVertLevels) + levelNum);
+ double y = *(vectorData[1] + (j*maxNVertLevels) + levelNum);
+ double z = *(vectorData[2] + (j*maxNVertLevels) + levelNum);
+ *(vectorData[0] + (j*maxNVertLevels) + levelNum) =
+ convertXY2Lon(x, y, lon);
+ *(vectorData[1] + (j*maxNVertLevels) + levelNum) =
+ convertXYZ2Lat(x, y, z, lon, lat);
+ }
+ }
+ }
+#endif
+
+ vtkFile << "VECTORS uReconstruct double" << endl;
+
+ // dummy first
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + 1));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + maxNVertLevels+ levelNum));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+ }
+ // write highest level dummy point
+ if (multilayer) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + maxNVertLevels + (maxNVertLevels-1)));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+
+ // rest of data
+ for (int j = 1; j < numDualPoints; j++)
+ {
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + j));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + (j * maxNVertLevels) + levelNum));
+ vtkFile << validData << " ";
+ if (j == 1000) {
+ cout << validData << " ";
+ }
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + (j * maxNVertLevels) + maxNVertLevels-1));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+ }
+
+                        // put out data for extra points
+                        for (int j = numDualPoints; j < currentExtraPoint; j++) {
+
+ // use map to find out what point data we are using
+                                int var_target = ((*(pointMap + j - numDualPoints))*maxNVertLevels);
+
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + var_target));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ } else {
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + var_target + levelNum));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ for (int d= 0; d < maxDim; d++) {
+ validData = convertDouble2ValidFloat (*(vectorData[d] + var_target + maxNVertLevels-1));
+ vtkFile << validData << " ";
+ }
+ //vtkFile << endl;
+ vtkFile << 0 << endl;
+ }
+ }
}
- // if by cell, then write out cell data
+                // if by cell, then write out cell data
if (singlelayer) {
if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell << endl;
@@ -965,14 +1238,14 @@
if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell*maxNVertLevels << endl;
}
- for (int v = 0; v <= dualCellVarIndex; v++) {
+                for (int v = 0; v <= dualCellVarIndex; v++) {
- // Write variable number v data for that timestep
- vtkFile << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
- vtkFile << "LOOKUP_TABLE default" << endl;
+                        // Write variable number v data for that timestep
+                        vtkFile << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
+                        vtkFile << "LOOKUP_TABLE default" << endl;
- // Read variable number v data for that timestep and level
- dualCellVars[v]->set_cur(i, 0, vertLevelIndex);
+                        // Read variable number v data for that timestep and level
+                        dualCellVars[v]->set_cur(i, 0, vertLevelIndex);
if (singlelayer) {
dualCellVars[v]->get(dualCellVarData, 1, numDualCells, 1);
} else {
@@ -981,12 +1254,12 @@
double *var_target = dualCellVarData;
- for (int j = 0; j < numDualCells; j++) {
+                        for (int j = 0; j < numDualCells; j++) {
if (singlelayer) {
- float validData = convertDouble2ValidFloat (*var_target);
- vtkFile << validData << endl;
- var_target++;
+                                        float validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+                                        var_target++;
} else {
for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
{
@@ -997,7 +1270,7 @@
}
}
- for (int j = numDualCells; j < currentExtraCell; j++) {
+                        for (int j = numDualCells; j < currentExtraCell; j++) {
if (j == mirrorcell) {
//cout << "data for mirror cell " << mirrorcell << " from cell " << *(cellMap + j - numDualCells) << endl;
}
</font>
</pre>