Commit 024ab16b authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

intlevel: optionally use 3d z-var from input file.

parent f57d8469
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.3
* Version 1.9.3 release
2018-01-27 Uwe Schulzweida
* intlevel: optionally use 3d z-var from input file
2018-01-25 Uwe Schulzweida
* Replaced isnan() by std::isnan()
......
......@@ -84,6 +84,36 @@ void vert_interp_lev(size_t gridsize, double missval, double *vardata1, double *
}
}
/*
* 3d vertical interpolation routine (see vert_interp_lev() in src/Intlevel.cc)
*/
void vert_interp_lev3d(size_t gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
for ( int ilev = 0; ilev < nlev2; ilev++ )
{
size_t offset = ilev*gridsize;
double *var2 = vardata2 + offset;
for ( size_t i = 0; i < gridsize; i++ )
{
int idx1 = lev_idx1[offset+i];
int idx2 = lev_idx2[offset+i];
double wgt1 = lev_wgt1[offset+i];
double wgt2 = lev_wgt2[offset+i];
/* upper/lower values from input field */
double var1L1 = *(vardata1+idx1);
double var1L2 = *(vardata1+idx2);
/* if (cdoVerbose) printf("i:%d level %d: idx1=%d idx2=%d (offset+i:%d) wgt1=%g wgt2=%g var1L1:%g var1L2:%g ",
* i, ilev, idx1, idx2, offset+i, wgt1, wgt2, var1L1, var1L2);
*/
var2[i] = vert_interp_lev_kernel(wgt1, wgt2, var1L1, var1L2, missval);
}
}
}
void vert_gen_weights(int expol, int nlev1, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
......@@ -147,6 +177,29 @@ void vert_gen_weights(int expol, int nlev1, double *lev1, int nlev2, double *lev
}
}
static
void vert_gen_weights3d1d(bool expol, int nlev1, size_t gridsize, double *xlev1, int nlev2, double *lev2,
int *xlev_idx1, int *xlev_idx2, double *xlev_wgt1, double *xlev_wgt2)
{
std::vector<double> lev1(nlev1);
std::vector<int> lev_idx1(nlev2);
std::vector<int> lev_idx2(nlev2);
std::vector<double> lev_wgt1(nlev2);
std::vector<double> lev_wgt2(nlev2);
for ( size_t i = 0; i < gridsize; i++ )
{
for ( int k = 0; k < nlev1; ++k ) lev1[k] = xlev1[k*gridsize+i];
vert_gen_weights(expol, nlev1, &lev1[0], nlev2, &lev2[0], &lev_idx1[0], &lev_idx2[0], &lev_wgt1[0], &lev_wgt2[0]);
for ( int k = 0; k < nlev2; ++k ) xlev_idx1[k*gridsize+i] = lev_idx1[k]*gridsize+i;
for ( int k = 0; k < nlev2; ++k ) xlev_idx2[k*gridsize+i] = lev_idx2[k]*gridsize+i;
for ( int k = 0; k < nlev2; ++k ) xlev_wgt1[k*gridsize+i] = lev_wgt1[k];
for ( int k = 0; k < nlev2; ++k ) xlev_wgt2[k*gridsize+i] = lev_wgt2[k];
}
}
bool levelDirUp(int nlev, double *lev)
{
......@@ -167,6 +220,40 @@ bool levelDirDown(int nlev, double *lev)
return ldown;
}
void vert_init_level_0_and_N(int nlev, size_t gridsize, double *zlevels)
{
/*
* Check monotony of vertical levels
*/
std::vector<double> zlev(nlev);
for ( int i = 0; i < nlev; ++i ) zlev[i] = zlevels[(i+1)*gridsize];
bool lup = levelDirUp(nlev, &zlev[0]);
bool ldown = levelDirDown(nlev, &zlev[0]);
/* Add artificial values for indication of extrapolation areas (lowermost + upmost levels) */
if ( lup )
{
for ( size_t i = 0; i < gridsize ;i++ )
{
zlevels[i] = -1.e33;
zlevels[(nlev+1)*gridsize + i] = 1.e33;
}
}
else if ( ldown )
{
for ( size_t i = 0; i < gridsize ;i++ )
{
zlevels[i] = 1.e33;
zlevels[(nlev+1)*gridsize + i] = -1.e33;
}
}
else
cdoWarning("Non monotonic zaxis!");
if ( cdoVerbose ) for ( int i = 0; i < nlev+2; ++i ) cdoPrint("lev1 %d: %g", i, zlevels[i*gridsize]);
}
void *Intlevel(void *argument)
{
int nrecs;
......@@ -214,77 +301,16 @@ void *Intlevel(void *argument)
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int i;
int nzaxis = vlistNzaxis(vlistID1);
for ( i = 0; i < nzaxis; i++ )
{
int zaxisID = vlistZaxis(vlistID1, i);
nlevel = zaxisInqSize(zaxisID);
if ( zaxisInqType(zaxisID) != ZAXIS_HYBRID && zaxisInqType(zaxisID) != ZAXIS_HYBRID_HALF )
if ( nlevel > 1 )
{
zaxisID1 = zaxisID;
break;
}
}
if ( i == nzaxis ) cdoAbort("No processable variable found!");
int nlev1 = nlevel;
double *lev1 = (double*) Malloc((nlev1+2)*sizeof(double));
cdoZaxisInqLevels(zaxisID1, lev1+1);
bool lup = levelDirUp(nlev1, lev1+1);
bool ldown = levelDirDown(nlev1, lev1+1);
if ( lup )
{
lev1[0] = -1.e33;
lev1[nlev1+1] = 1.e33;
}
else if ( ldown )
{
lev1[0] = 1.e33;
lev1[nlev1+1] = -1.e33;
}
else
cdoWarning("Non monotonic zaxis!");
if ( cdoVerbose ) for ( int i = 0; i < nlev1+2; ++i ) cdoPrint("lev1 %d: %g", i, lev1[i]);
int *lev_idx1 = (int*) Malloc(nlev2*sizeof(int));
int *lev_idx2 = (int*) Malloc(nlev2*sizeof(int));
double *lev_wgt1 = (double*) Malloc(nlev2*sizeof(double));
double *lev_wgt2 = (double*) Malloc(nlev2*sizeof(double));
vert_gen_weights(expol, nlev1+2, lev1, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
int zaxisID2 = zaxisCreate(zaxisInqType(zaxisID1), nlev2);
zaxisDefLevels(zaxisID2, lev2);
{
char str[CDI_MAX_NAME];
zaxisInqName(zaxisID1, str);
zaxisDefName(zaxisID2, str);
str[0] = 0;
zaxisInqLongname(zaxisID1, str);
if ( str[0] ) zaxisDefLongname(zaxisID2, str);
str[0] = 0;
zaxisInqUnits(zaxisID1, str);
if ( str[0] ) zaxisDefUnits(zaxisID2, str);
zaxisDefDatatype(zaxisID2, zaxisInqDatatype(zaxisID1));
}
for ( int i = 0; i < nzaxis; i++ )
if ( zaxisID1 == vlistZaxis(vlistID1, i) )
vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
pstreamDefVlist(streamID2, vlistID2);
int nvars = vlistNvars(vlistID1);
// Find z-variable
int zaxisID2 = CDI_UNDEFID;
int nlev1 = 0;
double *lev1 = NULL;
int zvarID = CDI_UNDEFID;
int zvarIsConstant = true;
bool zvarIsVarying = false;
size_t zvarGridsize = 0;
size_t wisize = 0;
if ( zvarname )
{
for ( varID = 0; varID < nvars; varID++ )
......@@ -299,9 +325,97 @@ void *Intlevel(void *argument)
}
if ( zvarID == CDI_UNDEFID ) cdoAbort("Variable %s not found!", zvarname);
zvarIsConstant = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
zvarIsVarying = vlistInqVarTimetype(vlistID1, zvarID) == TIME_VARYING;
zvarGridsize = gridInqSize(vlistInqVarGrid(vlistID1, zvarID));
zaxisID1 = vlistInqVarZaxis(vlistID1, zvarID);
nlev1 = zaxisInqSize(zaxisID1);
lev1 = (double*) Malloc(zvarGridsize*(nlev1+2)*sizeof(double));
zaxisID2 = zaxisCreate(ZAXIS_GENERIC, nlev2);
char str[CDI_MAX_NAME];
strcpy(str, "zlev");
zaxisDefName(zaxisID2, str);
str[0] = 0;
vlistInqVarLongname(vlistID1, zvarID, str);
if ( str[0] ) zaxisDefLongname(zaxisID2, str);
str[0] = 0;
vlistInqVarUnits(vlistID1, zvarID, str);
if ( str[0] ) zaxisDefUnits(zaxisID2, str);
wisize = zvarGridsize*nlev2;
}
else
{
int i;
int nzaxis = vlistNzaxis(vlistID1);
for ( i = 0; i < nzaxis; i++ )
{
int zaxisID = vlistZaxis(vlistID1, i);
nlevel = zaxisInqSize(zaxisID);
if ( zaxisInqType(zaxisID) != ZAXIS_HYBRID && zaxisInqType(zaxisID) != ZAXIS_HYBRID_HALF )
if ( nlevel > 1 )
{
zaxisID1 = zaxisID;
break;
}
}
if ( i == nzaxis ) cdoAbort("No processable variable found!");
zaxisID2 = zaxisCreate(zaxisInqType(zaxisID1), nlev2);
char str[CDI_MAX_NAME];
zaxisInqName(zaxisID1, str);
zaxisDefName(zaxisID2, str);
str[0] = 0;
zaxisInqLongname(zaxisID1, str);
if ( str[0] ) zaxisDefLongname(zaxisID2, str);
str[0] = 0;
zaxisInqUnits(zaxisID1, str);
if ( str[0] ) zaxisDefUnits(zaxisID2, str);
zaxisDefDatatype(zaxisID2, zaxisInqDatatype(zaxisID1));
nlev1 = nlevel;
lev1 = (double*) Malloc((nlev1+2)*sizeof(double));
cdoZaxisInqLevels(zaxisID1, lev1+1);
bool lup = levelDirUp(nlev1, lev1+1);
bool ldown = levelDirDown(nlev1, lev1+1);
if ( lup )
{
lev1[0] = -1.e33;
lev1[nlev1+1] = 1.e33;
}
else if ( ldown )
{
lev1[0] = 1.e33;
lev1[nlev1+1] = -1.e33;
}
else
cdoWarning("Non monotonic zaxis!");
if ( cdoVerbose ) for ( int i = 0; i < nlev1+2; ++i ) cdoPrint("lev1 %d: %g", i, lev1[i]);
wisize = nlev2;
}
zaxisDefLevels(zaxisID2, lev2);
int nzaxis = vlistNzaxis(vlistID1);
for ( int i = 0; i < nzaxis; i++ )
if ( zaxisID1 == vlistZaxis(vlistID1, i) )
vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
int *lev_idx1 = (int*) Malloc(wisize*sizeof(int));
int *lev_idx2 = (int*) Malloc(wisize*sizeof(int));
double *lev_wgt1 = (double*) Malloc(wisize*sizeof(double));
double *lev_wgt2 = (double*) Malloc(wisize*sizeof(double));
int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
pstreamDefVlist(streamID2, vlistID2);
bool *vars = (bool*) Malloc(nvars*sizeof(bool));
bool *varinterp = (bool*) Malloc(nvars*sizeof(bool));
size_t **varnmiss = (size_t**) Malloc(nvars*sizeof(size_t*));
......@@ -334,7 +448,6 @@ void *Intlevel(void *argument)
}
}
int tsID = 0;
while ( (nrecs = pstreamInqTimestep(streamID1, tsID)) )
{
......@@ -350,11 +463,29 @@ void *Intlevel(void *argument)
size_t gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
size_t offset = gridsize*levelID;
double *single1 = vardata1[varID] + offset;
pstreamReadRecord(streamID1, single1, &varnmiss[varID][levelID]);
vars[varID] = true;
}
if ( tsID == 0 || zvarIsVarying )
{
if ( zvarname )
{
for ( levelID = 0; levelID < nlev1; ++levelID )
{
size_t offset1 = zvarGridsize*levelID;
size_t offset2 = zvarGridsize*(levelID+1);
memcpy(lev1+offset2, vardata1[zvarID]+offset1, zvarGridsize*sizeof(double));
}
vert_init_level_0_and_N(nlev1, zvarGridsize, lev1);
vert_gen_weights3d1d(expol, nlev1+2, zvarGridsize, lev1, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
}
else
vert_gen_weights(expol, nlev1+2, lev1, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
}
for ( varID = 0; varID < nvars; varID++ )
{
if ( vars[varID] && varinterp[varID] )
......@@ -363,8 +494,12 @@ void *Intlevel(void *argument)
double missval = vlistInqVarMissval(vlistID1, varID);
size_t gridsize = gridInqSize(gridID);
vert_interp_lev(gridsize, missval, vardata1[varID], vardata2[varID],
nlev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
if ( zvarname )
vert_interp_lev3d(gridsize, missval, vardata1[varID], vardata2[varID],
nlev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
else
vert_interp_lev(gridsize, missval, vardata1[varID], vardata2[varID],
nlev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
for ( levelID = 0; levelID < nlev2; levelID++ )
{
......
......@@ -33,40 +33,11 @@
bool levelDirUp(int nlev, double *lev);
bool levelDirDown(int nlev, double *lev);
double vert_interp_lev_kernel(double w1, double w2, double var1L1, double var1L2, double missval);
void vert_interp_lev3d(size_t gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2);
void vert_gen_weights(int expol, int nlev1, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2);
/*
* 3d vertical interpolation routine (see vert_interp_lev() in src/Intlevel.cc)
*/
static
void vert_interp_lev3d(size_t gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
for ( int ilev = 0; ilev < nlev2; ilev++ )
{
size_t offset = ilev*gridsize;
double *var2 = vardata2 + offset;
for ( size_t i = 0; i < gridsize; i++ )
{
int idx1 = lev_idx1[offset+i];
int idx2 = lev_idx2[offset+i];
double wgt1 = lev_wgt1[offset+i];
double wgt2 = lev_wgt2[offset+i];
/* upper/lower values from input field */
double var1L1 = *(vardata1+idx1);
double var1L2 = *(vardata1+idx2);
/* if (cdoVerbose) printf("i:%d level %d: idx1=%d idx2=%d (offset+i:%d) wgt1=%g wgt2=%g var1L1:%g var1L2:%g ",
* i, ilev, idx1, idx2, offset+i, wgt1, wgt2, var1L1, var1L2);
*/
var2[i] = vert_interp_lev_kernel(wgt1, wgt2, var1L1, var1L2, missval);
}
}
}
void vert_init_level_0_and_N(int nlev, size_t gridsize, double *zlevels);
/*
* Create weights for the 3d vertical coordinate
......@@ -262,33 +233,7 @@ void *Intlevel3d(void *argument)
}
if ( i == nzaxis ) cdoAbort("No processable variable found (grid coordinate differ)!");
/*
* Check monotony of vertical levels
*/
std::vector<double> lev1(nlevi);
for ( int i = 0; i < nlevi; ++i ) lev1[i] = zlevels_in[(i+1)*gridSize];
bool lup = levelDirUp(nlevi, &lev1[0]);
bool ldown = levelDirDown(nlevi, &lev1[0]);
/* Add artificial values for intication of extrapolation areas (lowermost + upmost levels) */
if ( lup )
{
for ( size_t i = 0; i < gridSize ;i++ )
{
zlevels_in[i] = -1.e33;
zlevels_in[(nlevi+1)*gridSize + i] = 1.e33;
}
}
else if ( ldown )
{
for ( size_t i = 0; i < gridSize ;i++ )
{
zlevels_in[i] = 1.e33;
zlevels_in[(nlevi+1)*gridSize + i] = -1.e33;
}
}
else
cdoWarning("Non monotonic zaxis!");
vert_init_level_0_and_N(nlevi, gridSize, zlevels_in);
/*
* Create weights for later interpolation - assumption: input vertical correct is constant in time
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment