Commit a2f50a16 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Moved vert_interp_lev3d() to Intlevel3d.cc.

parent 75f051b3
......@@ -136,7 +136,7 @@ void partab(FILE *fp, int streamID, int option)
if ( natts > 0 )
{
fprintf(fp, "&parameter\n");
fprintf(fp, " name = _GLOBAL_\n");
fprintf(fp, " name=_GLOBAL_\n");
printHistory(fp, streamID);
printSource(fp, vlistID, 0);
cdo_print_attributes(fp, vlistID, CDI_GLOBAL, 2);
......@@ -160,14 +160,14 @@ void partab(FILE *fp, int streamID, int option)
{
fprintf(fp, "&parameter");
if ( linebreak ) fprintf(fp, "\n");
fprintf(fp, " name = _default_");
fprintf(fp, " name=_default_");
if ( linebreak ) fprintf(fp, "\n");
if ( datatype2str(datatype, pstr) == 0 )
{
fprintf(fp, " datatype = %s", pstr);
fprintf(fp, " datatype=%s", pstr);
if ( linebreak ) fprintf(fp, "\n");
}
fprintf(fp, "/\n");
fprintf(fp, " /\n");
}
}
......@@ -186,59 +186,59 @@ void partab(FILE *fp, int streamID, int option)
vlistInqVarLongname(vlistID, varID, varlongname);
vlistInqVarUnits(vlistID, varID, varunits);
fprintf(fp, " name = %s", varname);
fprintf(fp, " name=%s", varname);
if ( linebreak ) fprintf(fp, "\n");
// if ( code > 0 ) fprintf(fp, " code=%d\n", code);
// if ( tabnum > 0 ) fprintf(fp, " table=%d\n", tabnum);
if ( param >= 0 )
{
cdiParamToString(param, paramstr, sizeof(paramstr));
fprintf(fp, " param = %s", paramstr);
fprintf(fp, " param=%s", paramstr);
if ( linebreak ) fprintf(fp, "\n");
}
if ( strlen(varstdname) )
{
fprintf(fp, " standard_name = %s", varstdname);
fprintf(fp, " standard_name=%s", varstdname);
if ( linebreak ) fprintf(fp, "\n");
}
if ( strlen(varlongname) )
{
fprintf(fp, " long_name = \"%s\"", varlongname);
fprintf(fp, " long_name=\"%s\"", varlongname);
if ( linebreak ) fprintf(fp, "\n");
}
if ( strlen(varunits) )
{
fprintf(fp, " units = \"%s\"", varunits);
fprintf(fp, " units=\"%s\"", varunits);
if ( linebreak ) fprintf(fp, "\n");
}
if ( datatype == -1 )
if ( datatype2str(vlistInqVarDatatype(vlistID, varID), pstr) == 0 )
{
fprintf(fp, " datatype = %s", pstr);
fprintf(fp, " datatype=%s", pstr);
if ( linebreak ) fprintf(fp, "\n");
}
int chunktype = vlistInqVarChunkType(vlistID, varID);
if ( chunktype == CDI_CHUNK_AUTO )
{
fprintf(fp, " chunktype = auto");
fprintf(fp, " chunktype=auto");
if ( linebreak ) fprintf(fp, "\n");
}
else if ( chunktype == CDI_CHUNK_GRID )
{
fprintf(fp, " chunktype = grid");
fprintf(fp, " chunktype=grid");
if ( linebreak ) fprintf(fp, "\n");
}
if ( chunktype == CDI_CHUNK_LINES )
{
fprintf(fp, " chunktype = lines");
fprintf(fp, " chunktype=lines");
if ( linebreak ) fprintf(fp, "\n");
}
if ( option == 2 ) cdo_print_attributes(fp, vlistID, varID, 2);
if ( option == 2 )
fprintf(fp, " missing_value = %g\n", missval);
fprintf(fp, " missing_value=%g\n", missval);
if ( !linebreak ) fprintf(fp, " ");
fprintf(fp, "/\n");
......
......@@ -63,17 +63,11 @@ void vert_interp_lev(int gridsize, double missval, double *vardata1, double *var
}
else if ( IS_EQUAL(w1, 0) )
{
if ( w2 >= 0.5 )
var2[i] = var1L2[i];
else
var2[i] = missval;
var2[i] = (w2 >= 0.5) ? var1L2[i] : missval;
}
else if ( IS_EQUAL(w2, 0) )
{
if ( w1 >= 0.5 )
var2[i] = var1L1[i];
else
var2[i] = missval;
var2[i] = (w1 >= 0.5) ? var1L1[i] : missval;
}
else
{
......@@ -88,7 +82,7 @@ void vert_gen_weights(int expol, int nlev1, double *lev1, int nlev2, double *lev
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i1;
int idx1 = 0, idx2 = 0;
int idx1 = 0, idx2 = 0;
double val1, val2 = 0;
for ( int i2 = 0; i2 < nlev2; ++i2 )
......@@ -119,19 +113,13 @@ void vert_gen_weights(int expol, int nlev1, double *lev1, int nlev2, double *lev
lev_idx1[i2] = 1;
lev_idx2[i2] = 1;
lev_wgt1[i2] = 0;
if ( expol || IS_EQUAL(lev2[i2], val2) )
lev_wgt2[i2] = 1;
else
lev_wgt2[i2] = 0;
lev_wgt2[i2] = (expol || IS_EQUAL(lev2[i2], val2));
}
else if ( i1 == nlev1-1 )
{
lev_idx1[i2] = nlev1-2;
lev_idx2[i2] = nlev1-2;
if ( expol || IS_EQUAL(lev2[i2], val2) )
lev_wgt1[i2] = 1;
else
lev_wgt1[i2] = 0;
lev_wgt1[i2] = (expol || IS_EQUAL(lev2[i2], val2));
lev_wgt2[i2] = 0;
}
else
......@@ -177,14 +165,13 @@ void *Intlevel(void *argument)
int operatorID = cdoOperatorID();
bool expol = false;
if ( operatorID == INTLEVELX ) expol = true;
bool expol = (operatorID == INTLEVELX);
operatorInputArg("levels");
lista_t *flista = lista_new(FLT_LISTA);
int nlev2 = args2flt_lista(operatorArgc(), operatorArgv(), flista);
double *lev2 = (double *) lista_dataptr(flista);
double *lev2 = (double *) lista_dataptr(flista);
if ( cdoVerbose ) for ( i = 0; i < nlev2; ++i ) printf("lev2 %d: %g\n", i, lev2[i]);
......
......@@ -31,6 +31,235 @@
#include "after_vertint.h"
/*
* 3d vertical interpolation routine (see vert_interp_lev() in src/Intlevel.cc)
*/
static
void vert_interp_lev3d(int gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int idx1, idx2;
int offset;
double wgt1, wgt2;
double w1, w2;
double var1L1, var1L2, *var2;
for ( int ilev = 0; ilev < nlev2; ilev++ )
{
offset = ilev*gridsize;
var2 = vardata2 + offset;
for ( int i = 0; i < gridsize; i++ )
{
idx1 = lev_idx1[offset+i];
idx2 = lev_idx2[offset+i];
wgt1 = lev_wgt1[offset+i];
wgt2 = lev_wgt2[offset+i];
/* upper/lower values from input field */
var1L1 = *(vardata1+idx1);
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);
*/
w1 = wgt1;
w2 = wgt2;
if ( DBL_IS_EQUAL(var1L1, missval) ) w1 = 0;
if ( DBL_IS_EQUAL(var1L2, missval) ) w2 = 0;
if ( IS_EQUAL(w1, 0) && IS_EQUAL(w2, 0) )
{
var2[i] = missval;
}
else if ( IS_EQUAL(w1, 0) )
{
var2[i] = (w2 >= 0.5) ? var1L2 : missval;
}
else if ( IS_EQUAL(w2, 0) )
{
var2[i] = (w1 >= 0.5) ? var1L1 : missval;
}
else
{
var2[i] = var1L1*w1 + var1L2*w2;
}
}
}
}
/*
* Create weights for the 3d vertical coordinate
*
* The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
* wrt. the given gridsize. They can directly be used to read values from 3d
* data fields.
*
* 3d version of vert_gen_weights() (src/Intlevel.cc)
*/
static
void vert_gen_weights3d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i,i1, i2;
int idx1 = 0, idx2 = 0;
double val1, val2 = 0;
for ( i = 0; i < gridsize; i++ )
{
for ( i2 = 0; i2 < nlev2; i2++ )
{
/* Because 2 levels were added to the source vertical coordinate (one on
* top, one at the bottom), its loop starts at 1 */
for ( i1 = 1; i1 < nlev1; i1++ )
{
if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
{
idx1 = (i1-1)*gridsize+i;
idx2 = i1*gridsize+i;
}
else
{
idx1 = i1*gridsize+i;
idx2 = (i1-1)*gridsize+i;
}
val1 = lev1[idx1];
val2 = lev1[idx2];
if ( lev2[i2*gridsize+i] > val1 && lev2[i2*gridsize+i] <= val2 ) break;
}
if ( i1 == nlev1 )
{
if ( expol )
cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2*gridsize],i2);
else
cdoAbort("Level %g at index %d not found!");
}
if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
{
lev_idx1[i2*gridsize+i] = gridsize+i;
lev_idx2[i2*gridsize+i] = gridsize+i;
lev_wgt1[i2*gridsize+i] = 0;
lev_wgt2[i2*gridsize+i] = (expol || IS_EQUAL(lev2[i2*gridsize+i], val2));
}
else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
{
lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_wgt1[i2*gridsize+i] = (expol || IS_EQUAL(lev2[i2*gridsize+i], val2));
lev_wgt2[i2*gridsize+i] = 0;
}
else /* target z values has two bounday values in input z field */
{
lev_idx1[i2*gridsize+i] = idx1;
lev_idx2[i2*gridsize+i] = idx2;
lev_wgt1[i2*gridsize+i] = (lev1[idx2] - lev2[i2*gridsize+i]) / (lev1[idx2] - lev1[idx1]);
lev_wgt2[i2*gridsize+i] = (lev2[i2*gridsize+i] - lev1[idx1]) / (lev1[idx2] - lev1[idx1]);
}
/* if (cdoVerbose)
* {
* printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2*gridsize+i]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
* i, i2, i2*gridsize+i, lev2[i2*gridsize+i], idx1, idx2, lev1[idx1], lev1[idx2]);
* printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
* }
*/
/* backshift of the indices because of the two additional levels in input vertical coordinate */
lev_idx1[i2*gridsize+i] -= gridsize;
lev_idx2[i2*gridsize+i] -= gridsize;
}
}
}
/*
* Create weights for the 1d vertical coordinate from a 3d vertical coordinate
*
* The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
* wrt. the given gridsize. They can directly be used to read values from 3d
* data fields.
*
* 3d1d version of vert_gen_weights() (src/Intlevel.cc)
*/
static
void vert_gen_weights3d1d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i,i1, i2;
int idx1 = 0, idx2 = 0;
double val1, val2 = 0;
for ( i = 0; i < gridsize; i++ )
{
for ( i2 = 0; i2 < nlev2; i2++ )
{
/* Because 2 levels were added to the source vertical coordinate (one on
* top, one at the bottom), its loop starts at 1 */
for ( i1 = 1; i1 < nlev1; i1++ )
{
if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
{
idx1 = (i1-1)*gridsize+i;
idx2 = i1*gridsize+i;
}
else
{
idx1 = i1*gridsize+i;
idx2 = (i1-1)*gridsize+i;
}
val1 = lev1[idx1];
val2 = lev1[idx2];
if ( lev2[i2] > val1 && lev2[i2] <= val2 ) break;
}
if ( i1 == nlev1 )
{
if ( expol )
cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2],i2);
else
cdoAbort("Level %g at index %d not found!");
}
if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
{
lev_idx1[i2*gridsize+i] = gridsize+i;
lev_idx2[i2*gridsize+i] = gridsize+i;
lev_wgt1[i2*gridsize+i] = 0;
lev_wgt2[i2*gridsize+i] = (expol || IS_EQUAL(lev2[i2], val2));
}
else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
{
lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_wgt1[i2*gridsize+i] = (expol || IS_EQUAL(lev2[i2], val2));
lev_wgt2[i2*gridsize+i] = 0;
}
else /* target z values has two bounday values in input z field */
{
lev_idx1[i2*gridsize+i] = idx1;
lev_idx2[i2*gridsize+i] = idx2;
lev_wgt1[i2*gridsize+i] = (lev1[idx2] - lev2[i2]) / (lev1[idx2] - lev1[idx1]);
lev_wgt2[i2*gridsize+i] = (lev2[i2] - lev1[idx1]) / (lev1[idx2] - lev1[idx1]);
}
/* if (cdoVerbose)
* {
* printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
* i, i2, i2*gridsize+i, lev2[i2], idx1, idx2, lev1[idx1], lev1[idx2]);
* printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
* }
*/
/* backshift of the indices because of the two additional levels in input vertical coordinate */
lev_idx1[i2*gridsize+i] -= gridsize;
lev_idx2[i2*gridsize+i] -= gridsize;
}
}
}
void *Intlevel3d(void *argument)
{
......@@ -63,11 +292,11 @@ void *Intlevel3d(void *argument)
int INTLEVELX3D = cdoOperatorAdd("intlevelx3d", 0, 0, NULL);
// clang-format on
UNUSED(INTLEVEL3D);
int operatorID = cdoOperatorID();
bool expol = false;
if ( operatorID == INTLEVEL3D ) expol = false;
else if ( operatorID == INTLEVELX3D ) expol = true;
bool expol = (operatorID == INTLEVELX3D);
operatorInputArg("icoordinate");
......
......@@ -377,253 +377,3 @@ void interp_Z(const double *restrict geop, const double *restrict gz, double *pz
}
}
} /* interp_Z */
/*
* 3d vertical interpolation routine (see vert_interp_lev() in src/Intlevel.cc)
*/
void vert_interp_lev3d(int gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i, ilev;
int idx1, idx2;
int offset;
double wgt1, wgt2;
double w1, w2;
double var1L1, var1L2, *var2;
for ( ilev = 0; ilev < nlev2; ilev++ )
{
offset = ilev*gridsize;
var2 = vardata2 + offset;
for ( i = 0; i < gridsize; i++ )
{
idx1 = lev_idx1[offset+i];
idx2 = lev_idx2[offset+i];
wgt1 = lev_wgt1[offset+i];
wgt2 = lev_wgt2[offset+i];
/* upper/lower values from input field */
var1L1 = *(vardata1+idx1);
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);
*/
w1 = wgt1;
w2 = wgt2;
if ( DBL_IS_EQUAL(var1L1, missval) ) w1 = 0;
if ( DBL_IS_EQUAL(var1L2, missval) ) w2 = 0;
if ( IS_EQUAL(w1, 0) && IS_EQUAL(w2, 0) )
{
var2[i] = missval;
}
else if ( IS_EQUAL(w1, 0) )
{
if ( w2 >= 0.5 )
var2[i] = var1L2;
else
var2[i] = missval;
}
else if ( IS_EQUAL(w2, 0) )
{
if ( w1 >= 0.5 )
var2[i] = var1L1;
else
var2[i] = missval;
}
else
{
var2[i] = var1L1*w1 + var1L2*w2;
}
}
}
}
#if defined(CDO)
#include "util.h"
/*
* Create weights for the 3d vertical coordinate
*
* The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
* wrt. the given gridsize. They can directly be used to read values from 3d
* data fields.
*
* 3d version of vert_gen_weights() (src/Intlevel.cc)
*/
void vert_gen_weights3d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i,i1, i2;
int idx1 = 0, idx2 = 0;
double val1, val2 = 0;
for ( i = 0; i < gridsize; i++ )
{
for ( i2 = 0; i2 < nlev2; i2++ )
{
/* Because 2 levels were added to the source vertical coordinate (one on
* top, one at the bottom), its loop starts at 1 */
for ( i1 = 1; i1 < nlev1; i1++ )
{
if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
{
idx1 = (i1-1)*gridsize+i;
idx2 = i1*gridsize+i;
}
else
{
idx1 = i1*gridsize+i;
idx2 = (i1-1)*gridsize+i;
}
val1 = lev1[idx1];
val2 = lev1[idx2];
if ( lev2[i2*gridsize+i] > val1 && lev2[i2*gridsize+i] <= val2 ) break;
}
if ( i1 == nlev1 )
{
if ( expol )
cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2*gridsize],i2);
else
cdoAbort("Level %g at index %d not found!");
}
if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
{
lev_idx1[i2*gridsize+i] = gridsize+i;
lev_idx2[i2*gridsize+i] = gridsize+i;
lev_wgt1[i2*gridsize+i] = 0;
if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
lev_wgt2[i2*gridsize+i] = 1;
else
lev_wgt2[i2*gridsize+i] = 0;
}
else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
{
lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
lev_wgt1[i2*gridsize+i] = 1;
else
lev_wgt1[i2*gridsize+i] = 0;
lev_wgt2[i2*gridsize+i] = 0;
}
else /* target z values has two bounday values in input z field */
{
lev_idx1[i2*gridsize+i] = idx1;
lev_idx2[i2*gridsize+i] = idx2;
lev_wgt1[i2*gridsize+i] = (lev1[idx2] - lev2[i2*gridsize+i]) / (lev1[idx2] - lev1[idx1]);
lev_wgt2[i2*gridsize+i] = (lev2[i2*gridsize+i] - lev1[idx1]) / (lev1[idx2] - lev1[idx1]);
}
/* if (cdoVerbose)
* {
* printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2*gridsize+i]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
* i, i2, i2*gridsize+i, lev2[i2*gridsize+i], idx1, idx2, lev1[idx1], lev1[idx2]);
* printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
* }
*/
/* backshift of the indices because of the two additional levels in input vertical coordinate */
lev_idx1[i2*gridsize+i] -= gridsize;
lev_idx2[i2*gridsize+i] -= gridsize;
}
}
}
/*
* Create weights for the 1d vertical coordinate from a 3d vertical coordinate
*
* The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
* wrt. the given gridsize. They can directly be used to read values from 3d
* data fields.
*
* 3d1d version of vert_gen_weights() (src/Intlevel.cc)
*/
void vert_gen_weights3d1d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i,i1, i2;
int idx1 = 0, idx2 = 0;
double val1, val2 = 0;
for ( i = 0; i < gridsize; i++ )
{
for ( i2 = 0; i2 < nlev2; i2++ )
{
/* Because 2 levels were added to the source vertical coordinate (one on
* top, one at the bottom), its loop starts at 1 */
for ( i1 = 1; i1 < nlev1; i1++ )
{