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

Merge branch 'develop'

parents d4eaf3d0 d900de06
2016-09-24 Uwe Schulzweida
* remapnn: optimize sort in kdtree (speedup ~20%)
2016-09-19 Uwe Schulzweida
* New operator shiftx/shifty - Shift fields on rectangular grid in x/y direction
2016-08-29 Uwe Schulzweida
* CDO option -v includes -W
2016-08-18 Uwe Schulzweida
* using CDI library version 1.8.0rc2
......
......@@ -95,7 +95,7 @@ case "${HOSTNAME}" in
${CONFPATH}configure --enable-cxx --prefix=$HOME/local \
$CDOLIBS \
CC=icc CXX=icpc CFLAGS="-g -Wall -O2 -qopt-report=5 -march=native" CXX=icpc
elif test "$COMP" = icpc ; then
elif test "$COMP" = icc ; then
${CONFPATH}configure \
$CDOLIBS \
CC=icc CFLAGS="-g -Wall -Wwrite-strings -O2 -qopt-report=5 -march=native"
......
#! /bin/sh
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.68 for cdo 1.8.0rc2.
# Generated by GNU Autoconf 2.68 for cdo 1.8.0rc3.
#
# Report bugs to <http://mpimet.mpg.de/cdo>.
#
......@@ -570,8 +570,8 @@ MAKEFLAGS=
# Identity of this package.
PACKAGE_NAME='cdo'
PACKAGE_TARNAME='cdo'
PACKAGE_VERSION='1.8.0rc2'
PACKAGE_STRING='cdo 1.8.0rc2'
PACKAGE_VERSION='1.8.0rc3'
PACKAGE_STRING='cdo 1.8.0rc3'
PACKAGE_BUGREPORT='http://mpimet.mpg.de/cdo'
PACKAGE_URL=''
......@@ -1395,7 +1395,7 @@ if test "$ac_init_help" = "long"; then
# Omit some internal or obsolete options to make the list less imposing.
# This message is too long to be a string in the A/UX 3.1 sh.
cat <<_ACEOF
\`configure' configures cdo 1.8.0rc2 to adapt to many kinds of systems.
\`configure' configures cdo 1.8.0rc3 to adapt to many kinds of systems.
Usage: $0 [OPTION]... [VAR=VALUE]...
......@@ -1465,7 +1465,7 @@ fi
if test -n "$ac_init_help"; then
case $ac_init_help in
short | recursive ) echo "Configuration of cdo 1.8.0rc2:";;
short | recursive ) echo "Configuration of cdo 1.8.0rc3:";;
esac
cat <<\_ACEOF
......@@ -1613,7 +1613,7 @@ fi
test -n "$ac_init_help" && exit $ac_status
if $ac_init_version; then
cat <<\_ACEOF
cdo configure 1.8.0rc2
cdo configure 1.8.0rc3
generated by GNU Autoconf 2.68
Copyright (C) 2010 Free Software Foundation, Inc.
......@@ -2206,7 +2206,7 @@ cat >config.log <<_ACEOF
This file contains any messages produced by compilers while
running configure, to aid debugging if configure makes a mistake.
It was created by cdo $as_me 1.8.0rc2, which was
It was created by cdo $as_me 1.8.0rc3, which was
generated by GNU Autoconf 2.68. Invocation command line was
$ $0 $@
......@@ -3155,7 +3155,7 @@ fi
# Define the identity of the package.
PACKAGE='cdo'
VERSION='1.8.0rc2'
VERSION='1.8.0rc3'
cat >>confdefs.h <<_ACEOF
......@@ -21644,7 +21644,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
# report actual input values of CONFIG_FILES etc. instead of their
# values after options handling.
ac_log="
This file was extended by cdo $as_me 1.8.0rc2, which was
This file was extended by cdo $as_me 1.8.0rc3, which was
generated by GNU Autoconf 2.68. Invocation command line was
CONFIG_FILES = $CONFIG_FILES
......@@ -21710,7 +21710,7 @@ _ACEOF
cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
ac_cs_version="\\
cdo config.status 1.8.0rc2
cdo config.status 1.8.0rc3
configured by $0, generated by GNU Autoconf 2.68,
with options \\"\$ac_cs_config\\"
......
......@@ -4,7 +4,7 @@
# autoconf 2.68
# libtool 2.4.2
AC_INIT([cdo], [1.8.0rc2], [http://mpimet.mpg.de/cdo])
AC_INIT([cdo], [1.8.0rc3], [http://mpimet.mpg.de/cdo])
AC_DEFINE_UNQUOTED(CDO, ["$PACKAGE_VERSION"], [CDO version])
......
......@@ -54,26 +54,22 @@ double adisit_1(double tpot, double sal, double p)
a_d = 4.1057E-9,
a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double dc, dv, dvs, fne, fst, qc, qn3, qnq, qv, qvs, t, tpo;
qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
qv = p * (a_b1 - a_d * p);
dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
dv = a_b2 * p;
qnq = -p * (-a_a3 + p * a_c3);
qn3 = -p * a_a4;
{
tpo = tpot;
qvs = qv*(sal - 35.) + qc;
dvs = dv*(sal - 35.) + dc;
t = (tpo + qvs)/dvs;
fne = - qvs + t*(dvs + t*(qnq + t*qn3)) - tpo;
fst = dvs + t*(2.*qnq + 3.*qn3*t);
t = t - fne/fst;
}
return (t);
double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
double qv = p * (a_b1 - a_d * p);
double dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
double dv = a_b2 * p;
double qnq = -p * (-a_a3 + p * a_c3);
double qn3 = -p * a_a4;
double tpo = tpot;
double qvs = qv*(sal - 35.) + qc;
double dvs = dv*(sal - 35.) + dc;
double t = (tpo + qvs)/dvs;
double fne = - qvs + t*(dvs + t*(qnq + t*qn3)) - tpo;
double fst = dvs + t*(2.*qnq + 3.*qn3*t);
t = t - fne/fst;
return t;
}
/* compute potential temperature from insitu temperature */
......@@ -87,19 +83,17 @@ double adipot(double t, double s, double p)
a_d = 4.1057E-9,
a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double aa,bb,cc,cc1,dd,tpot,s_rel;
double s_rel = s - 35.0;
s_rel = s - 35.0;
double aa = (a_a1+ t*(a_a2 - t*(a_a3 - a_a4*t)));
double bb = s_rel*(a_b1 -a_b2*t) ;
double cc = (a_c1 + t*(-a_c2 + a_c3*t));
double cc1 = a_d*s_rel;
double dd = (-a_e1 + a_e2*t);
aa = (a_a1+ t*(a_a2 - t*(a_a3 - a_a4*t)));
bb = s_rel*(a_b1 -a_b2*t) ;
cc = (a_c1 + t*(-a_c2 + a_c3*t));
cc1 = a_d*s_rel;
dd = (-a_e1 + a_e2*t);
double tpot = t-p*(aa + bb + p*(cc - cc1 + p*dd));
tpot=t-p*(aa + bb + p*(cc - cc1 + p*dd));
return (tpot);
return tpot;
}
static
......@@ -109,17 +103,14 @@ void calc_adisit(long gridsize, long nlevel, double *pressure, field_t tho, fiel
/* tho units: Celsius */
/* sao units: psu */
long i, levelID, offset;
double *tisptr, *thoptr, *saoptr;
for ( levelID = 0; levelID < nlevel; ++levelID )
for ( long levelID = 0; levelID < nlevel; ++levelID )
{
offset = gridsize*levelID;
thoptr = tho.ptr + offset;
saoptr = sao.ptr + offset;
tisptr = tis.ptr + offset;
long offset = gridsize*levelID;
double *thoptr = tho.ptr + offset;
double *saoptr = sao.ptr + offset;
double *tisptr = tis.ptr + offset;
for ( i = 0; i < gridsize; ++i )
for ( long i = 0; i < gridsize; ++i )
{
if ( DBL_IS_EQUAL(thoptr[i], tho.missval) ||
DBL_IS_EQUAL(saoptr[i], sao.missval) )
......@@ -141,17 +132,14 @@ void calc_adipot(long gridsize, long nlevel, double *pressure, field_t t, field_
/* t units: Celsius */
/* s units: psu */
long i, levelID, offset;
double *tpotptr, *tptr, *sptr;
for ( levelID = 0; levelID < nlevel; ++levelID )
for ( long levelID = 0; levelID < nlevel; ++levelID )
{
offset = gridsize*levelID;
tptr = t.ptr + offset;
sptr = s.ptr + offset;
tpotptr = tpot.ptr + offset;
long offset = gridsize*levelID;
double *tptr = t.ptr + offset;
double *sptr = s.ptr + offset;
double *tpotptr = tpot.ptr + offset;
for ( i = 0; i < gridsize; ++i )
for ( long i = 0; i < gridsize; ++i )
{
if ( DBL_IS_EQUAL(tptr[i], t.missval) ||
DBL_IS_EQUAL(sptr[i], s.missval) )
......@@ -237,7 +225,7 @@ void *Adisit(void *argument)
int nlevel = nlevel1;
double *pressure = (double*) Malloc(nlevel*sizeof(double));
zaxisInqLevels(zaxisID, pressure);
cdoZaxisInqLevels(zaxisID, pressure);
if ( pin >= 0 )
for ( i = 0; i < nlevel; ++i ) pressure[i] = pin;
......@@ -267,11 +255,15 @@ void *Adisit(void *argument)
sao.missval = vlistInqVarMissval(vlistID1, saoID);
tis.missval = tho.missval;
int datatype = CDI_DATATYPE_FLT32;
if ( vlistInqVarDatatype(vlistID1, thoID) == CDI_DATATYPE_FLT64 &&
vlistInqVarDatatype(vlistID1, saoID) == CDI_DATATYPE_FLT64 )
datatype = CDI_DATATYPE_FLT64;
int vlistID2 = vlistCreate();
int tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARIABLE);
if (operatorID == ADISIT)
if ( operatorID == ADISIT )
{
vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
vlistDefVarName(vlistID2, tisID2, "to");
......@@ -287,6 +279,7 @@ void *Adisit(void *argument)
}
vlistDefVarUnits(vlistID2, tisID2, "K");
vlistDefVarMissval(vlistID2, tisID2, tis.missval);
vlistDefVarDatatype(vlistID2, tisID2, datatype);
int saoID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarParam(vlistID2, saoID2, cdiEncodeParam(5, 255, 255));
......@@ -295,6 +288,7 @@ void *Adisit(void *argument)
vlistDefVarStdname(vlistID2, saoID2, "sea_water_salinity");
vlistDefVarUnits(vlistID2, saoID2, "psu");
vlistDefVarMissval(vlistID2, saoID2, sao.missval);
vlistDefVarDatatype(vlistID2, saoID2, datatype);
int taxisID1 = vlistInqTaxis(vlistID1);
......@@ -331,13 +325,9 @@ void *Adisit(void *argument)
}
if ( operatorID == ADISIT )
{
calc_adisit(gridsize, nlevel, pressure, tho, sao, tis);
}
calc_adisit(gridsize, nlevel, pressure, tho, sao, tis);
else
{
calc_adipot(gridsize, nlevel, pressure, tho, sao, tis);
}
calc_adipot(gridsize, nlevel, pressure, tho, sao, tis);
for ( levelID = 0; levelID < nlevel; ++levelID )
......
......@@ -387,7 +387,7 @@ static
void *after_readTimestep(void *arg)
{
int i;
int recID, varID, gridID, zaxisID, levelID, timeID;
int varID, gridID, zaxisID, levelID, timeID;
int code, leveltype;
int nmiss;
RARG *rarg = (RARG *) arg;
......@@ -402,7 +402,7 @@ void *after_readTimestep(void *arg)
int level = 0;
int levelOffset = 0;
for ( recID = 0; recID < nrecs; recID++ )
for ( int recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(globs->istreamID, &varID, &levelID);
......@@ -688,7 +688,7 @@ void after_control(struct Control *globs, struct Variable *vars)
rtime = after_getTime(globs->StartDate);
}
if ( ofiletype == FILETYPE_NC || ofiletype == FILETYPE_NC2 || ofiletype == FILETYPE_NC4 )
if ( ofiletype == CDI_FILETYPE_NC || ofiletype == CDI_FILETYPE_NC2 || ofiletype == CDI_FILETYPE_NC4 )
{
taxisDefCalendar(globs->taxisID2, CALENDAR_PROLEPTIC);
taxisDefType(globs->taxisID2, TAXIS_RELATIVE);
......@@ -1350,29 +1350,29 @@ void after_parini(struct Control *globs, struct Variable *vars)
#if defined(CDO)
case -1: ofiletype = -1; break;
#else
case -1: ofiletype = FILETYPE_SRV; break;
case -1: ofiletype = CDI_FILETYPE_SRV; break;
#endif
case 0: ofiletype = FILETYPE_SRV; break;
case 1: ofiletype = FILETYPE_GRB; break;
case 2: ofiletype = FILETYPE_NC; break;
case 3: ofiletype = FILETYPE_EXT; break;
case 4: ofiletype = FILETYPE_NC2; break;
case 6: ofiletype = FILETYPE_NC4; break;
case 0: ofiletype = CDI_FILETYPE_SRV; break;
case 1: ofiletype = CDI_FILETYPE_GRB; break;
case 2: ofiletype = CDI_FILETYPE_NC; break;
case 3: ofiletype = CDI_FILETYPE_EXT; break;
case 4: ofiletype = CDI_FILETYPE_NC2; break;
case 6: ofiletype = CDI_FILETYPE_NC4; break;
default: Error( "unknown file format %d", fileFormat);
}
if ( gribFormat ) ofiletype = FILETYPE_GRB;
if ( cdfFormat ) ofiletype = FILETYPE_NC;
if ( gribFormat ) ofiletype = CDI_FILETYPE_GRB;
if ( cdfFormat ) ofiletype = CDI_FILETYPE_NC;
int precision = scan_par(globs->Verbose, namelist, "precision", 0);
if ( precision )
switch ( precision )
{
case 8: DataType = DATATYPE_PACK8; break;
case 16: DataType = DATATYPE_PACK16; break;
case 24: DataType = DATATYPE_PACK24; break;
case 32: DataType = DATATYPE_FLT32; break;
case 64: DataType = DATATYPE_FLT64; break;
case 8: DataType = CDI_DATATYPE_PACK8; break;
case 16: DataType = CDI_DATATYPE_PACK16; break;
case 24: DataType = CDI_DATATYPE_PACK24; break;
case 32: DataType = CDI_DATATYPE_FLT32; break;
case 64: DataType = CDI_DATATYPE_FLT64; break;
default: Error( "unsupported data precision %d", precision);
}
......@@ -2117,7 +2117,7 @@ void after_processing(struct Control *globs, struct Variable *vars)
if ( ! globs->AnalysisData )
for (i = 0; i < globs->NumLevelRequest; i++)
{
if ( (globs->LevelRequest[i] >= 65535) && globs->unitsel && ofiletype == FILETYPE_GRB )
if ( (globs->LevelRequest[i] >= 65535) && globs->unitsel && ofiletype == CDI_FILETYPE_GRB )
{
fprintf(stderr,"\n Level %9.2f out of range (max=65535)!\n", globs->LevelRequest[i]);
exit(1);
......@@ -2149,10 +2149,10 @@ void after_processing(struct Control *globs, struct Variable *vars)
if ( globs->Type == 10 || globs->Type == 40 || globs->Type == 60 )
{
if ( ofiletype == FILETYPE_GRB )
if ( ofiletype == CDI_FILETYPE_GRB )
Error("Can't write fourier coefficients to GRIB!");
else if ( ofiletype == FILETYPE_NC || ofiletype == FILETYPE_NC2 ||
ofiletype == FILETYPE_NC4 )
else if ( ofiletype == CDI_FILETYPE_NC || ofiletype == CDI_FILETYPE_NC2 ||
ofiletype == CDI_FILETYPE_NC4 )
Error("Can't write fourier coefficients to NetCDF!");
}
......
......@@ -27,15 +27,15 @@ const char *filetypestr(int filetype)
{
switch ( filetype )
{
case FILETYPE_GRB: return ("GRIB"); break;
case FILETYPE_GRB2: return ("GRIB2"); break;
case FILETYPE_NC: return ("NetCDF"); break;
case FILETYPE_NC2: return ("NetCDF2"); break;
case FILETYPE_NC4: return ("NetCDF4"); break;
case FILETYPE_NC4C: return ("NetCDF4 classic"); break;
case FILETYPE_SRV: return ("SERVICE"); break;
case FILETYPE_EXT: return ("EXTRA"); break;
case FILETYPE_IEG: return ("IEG"); break;
case CDI_FILETYPE_GRB: return ("GRIB"); break;
case CDI_FILETYPE_GRB2: return ("GRIB2"); break;
case CDI_FILETYPE_NC: return ("NetCDF"); break;
case CDI_FILETYPE_NC2: return ("NetCDF2"); break;
case CDI_FILETYPE_NC4: return ("NetCDF4"); break;
case CDI_FILETYPE_NC4C: return ("NetCDF4 classic"); break;
case CDI_FILETYPE_SRV: return ("SERVICE"); break;
case CDI_FILETYPE_EXT: return ("EXTRA"); break;
case CDI_FILETYPE_IEG: return ("IEG"); break;
default: return ("");
}
}
......@@ -48,18 +48,18 @@ const char *datatypestr(int datatype)
str[0] = 0;
sprintf(str, "%d bit packed", datatype);
if ( datatype == DATATYPE_PACK ) return ("P0");
if ( datatype == CDI_DATATYPE_PACK ) return ("P0");
else if ( datatype > 0 && datatype <= 32 ) return (str);
else if ( datatype == DATATYPE_CPX32 ) return ("C32");
else if ( datatype == DATATYPE_CPX64 ) return ("C64");
else if ( datatype == DATATYPE_FLT32 ) return ("32 bit floats");
else if ( datatype == DATATYPE_FLT64 ) return ("64 bit floats");
else if ( datatype == DATATYPE_INT8 ) return ("I8");
else if ( datatype == DATATYPE_INT16 ) return ("I16");
else if ( datatype == DATATYPE_INT32 ) return ("I32");
else if ( datatype == DATATYPE_UINT8 ) return ("U8");
else if ( datatype == DATATYPE_UINT16 ) return ("U16");
else if ( datatype == DATATYPE_UINT32 ) return ("U32");
else if ( datatype == CDI_DATATYPE_CPX32 ) return ("C32");
else if ( datatype == CDI_DATATYPE_CPX64 ) return ("C64");
else if ( datatype == CDI_DATATYPE_FLT32 ) return ("32 bit floats");
else if ( datatype == CDI_DATATYPE_FLT64 ) return ("64 bit floats");
else if ( datatype == CDI_DATATYPE_INT8 ) return ("I8");
else if ( datatype == CDI_DATATYPE_INT16 ) return ("I16");
else if ( datatype == CDI_DATATYPE_INT32 ) return ("I32");
else if ( datatype == CDI_DATATYPE_UINT8 ) return ("U8");
else if ( datatype == CDI_DATATYPE_UINT16 ) return ("U16");
else if ( datatype == CDI_DATATYPE_UINT32 ) return ("U32");
else return ("");
}
......
......@@ -29,7 +29,7 @@
void *CDItest(void *argument)
{
int nrecs;
int recID, varID, levelID;
int varID, levelID;
int nmiss;
int max_copy = 3;
double s_utime, s_stime;
......@@ -78,7 +78,7 @@ void *CDItest(void *argument)
streamDefTimestep(streamID2, tsID2);
for ( recID = 0; recID < nrecs; recID++ )
for ( int recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID1, &varID, &levelID);
streamDefRecord(streamID2, varID, levelID);
......
......@@ -28,15 +28,15 @@ const char *filetypestr(int filetype)
{
switch ( filetype )
{
case FILETYPE_GRB: return ("GRIB"); break;
case FILETYPE_GRB2: return ("GRIB2"); break;
case FILETYPE_NC: return ("NetCDF"); break;
case FILETYPE_NC2: return ("NetCDF2"); break;
case FILETYPE_NC4: return ("NetCDF4"); break;
case FILETYPE_NC4C: return ("NetCDF4 classic"); break;
case FILETYPE_SRV: return ("SERVICE"); break;
case FILETYPE_EXT: return ("EXTRA"); break;
case FILETYPE_IEG: return ("IEG"); break;
case CDI_FILETYPE_GRB: return ("GRIB"); break;
case CDI_FILETYPE_GRB2: return ("GRIB2"); break;
case CDI_FILETYPE_NC: return ("NetCDF"); break;
case CDI_FILETYPE_NC2: return ("NetCDF2"); break;
case CDI_FILETYPE_NC4: return ("NetCDF4"); break;
case CDI_FILETYPE_NC4C: return ("NetCDF4 classic"); break;
case CDI_FILETYPE_SRV: return ("SERVICE"); break;
case CDI_FILETYPE_EXT: return ("EXTRA"); break;
case CDI_FILETYPE_IEG: return ("IEG"); break;
default: return ("");
}
}
......@@ -49,18 +49,18 @@ const char *datatypestr(int datatype)
str[0] = 0;
sprintf(str, "%d bit packed", datatype);
if ( datatype == DATATYPE_PACK ) return ("P0");
if ( datatype == CDI_DATATYPE_PACK ) return ("P0");
else if ( datatype > 0 && datatype <= 32 ) return (str);
else if ( datatype == DATATYPE_CPX32 ) return ("C32");
else if ( datatype == DATATYPE_CPX64 ) return ("C64");
else if ( datatype == DATATYPE_FLT32 ) return ("32 bit floats");
else if ( datatype == DATATYPE_FLT64 ) return ("64 bit floats");
else if ( datatype == DATATYPE_INT8 ) return ("I8");
else if ( datatype == DATATYPE_INT16 ) return ("I16");
else if ( datatype == DATATYPE_INT32 ) return ("I32");
else if ( datatype == DATATYPE_UINT8 ) return ("U8");
else if ( datatype == DATATYPE_UINT16 ) return ("U16");
else if ( datatype == DATATYPE_UINT32 ) return ("U32");
else if ( datatype == CDI_DATATYPE_CPX32 ) return ("C32");
else if ( datatype == CDI_DATATYPE_CPX64 ) return ("C64");
else if ( datatype == CDI_DATATYPE_FLT32 ) return ("32 bit floats");
else if ( datatype == CDI_DATATYPE_FLT64 ) return ("64 bit floats");
else if ( datatype == CDI_DATATYPE_INT8 ) return ("I8");
else if ( datatype == CDI_DATATYPE_INT16 ) return ("I16");
else if ( datatype == CDI_DATATYPE_INT32 ) return ("I32");
else if ( datatype == CDI_DATATYPE_UINT8 ) return ("U8");
else if ( datatype == CDI_DATATYPE_UINT16 ) return ("U16");
else if ( datatype == CDI_DATATYPE_UINT32 ) return ("U32");
else return ("");
}
......
This diff is collapsed.
......@@ -196,33 +196,36 @@ void *Change(void *argument)
for ( index = 0; index < nzaxis; index++ )
{
zaxisID1 = vlistZaxis(vlistID2, index);
nlevs = zaxisInqSize(zaxisID1);
levels = (double*) Malloc(nlevs*sizeof(double));
newlevels = (double*) Malloc(nlevs*sizeof(double));
zaxisInqLevels(zaxisID1, levels);
for ( k = 0; k < nlevs; k++ ) newlevels[k] = levels[k];
nfound = 0;
for ( i = 0; i < nch; i += 2 )
for ( k = 0; k < nlevs; k++ )
if ( fabs(levels[k] - chlevels[i]) < 0.0001 ) nfound++;
if ( nfound )
{
zaxisID2 = zaxisDuplicate(zaxisID1);
for ( i = 0; i < nch; i += 2 )
for ( k = 0; k < nlevs; k++ )
if ( fabs(levels[k] - chlevels[i]) < 0.001 )
newlevels[k] = chlevels[i+1];
zaxisDefLevels(zaxisID2, newlevels);
vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
}
Free(levels);
Free(newlevels);
}
if ( zaxisInqLevels(zaxisID1, NULL) )
{
nlevs = zaxisInqSize(zaxisID1);
levels = (double*) Malloc(nlevs*sizeof(double));
newlevels = (double*) Malloc(nlevs*sizeof(double));
zaxisInqLevels(zaxisID1, levels);
for ( k = 0; k < nlevs; k++ ) newlevels[k] = levels[k];
nfound = 0;
for ( i = 0; i < nch; i += 2 )
for ( k = 0; k < nlevs; k++ )
if ( fabs(levels[k] - chlevels[i]) < 0.0001 ) nfound++;
if ( nfound )
{
zaxisID2 = zaxisDuplicate(zaxisID1);
for ( i = 0; i < nch; i += 2 )
for ( k = 0; k < nlevs; k++ )
if ( fabs(levels[k] - chlevels[i]) < 0.001 )
newlevels[k] = chlevels[i+1];
zaxisDefLevels(zaxisID2, newlevels);
vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
}
Free(levels);
Free(newlevels);
}
}
}
else if ( operatorID == CHLEVELC || operatorID == CHLEVELV )
{
......@@ -247,27 +250,30 @@ void *Change(void *argument)
}
zaxisID1 = vlistInqVarZaxis(vlistID2, varID);
nlevs = zaxisInqSize(zaxisID1);
levels = (double*) Malloc(nlevs*sizeof(double));
zaxisInqLevels(zaxisID1, levels);
nfound = 0;
for ( k = 0; k < nlevs; k++ )
if ( fabs(levels[k] - chlevels[0]) < 0.0001 ) nfound++;
if ( nfound )