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

Merge branch 'develop' of git.mpimet.mpg.de:cdo into develop

parents 4c00a824 49eef9c4
......@@ -47,22 +47,22 @@
static double
adisit_1(double tpot, double sal, double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
const double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
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;
const double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
const double qv = p * (a_b1 - a_d * p);
const double dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
const double dv = a_b2 * p;
const double qnq = -p * (-a_a3 + p * a_c3);
const double qn3 = -p * a_a4;
double tpo = tpot;
double qvs = qv * (sal - 35.) + qc;
double dvs = dv * (sal - 35.) + dc;
const double tpo = tpot;
const double qvs = qv * (sal - 35.) + qc;
const 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);
const double fne = -qvs + t * (dvs + t * (qnq + t * qn3)) - tpo;
const double fst = dvs + t * (2. * qnq + 3. * qn3 * t);
t = t - fne / fst;
return t;
......@@ -73,18 +73,18 @@ adisit_1(double tpot, double sal, double p)
static double
adipot(double t, double s, double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
const double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double s_rel = s - 35.0;
const double 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);
const double aa = (a_a1 + t * (a_a2 - t * (a_a3 - a_a4 * t)));
const double bb = s_rel * (a_b1 - a_b2 * t);
const double cc = (a_c1 + t * (-a_c2 + a_c3 * t));
const double cc1 = a_d * s_rel;
const double dd = (-a_e1 + a_e2 * t);
double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
const double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
return tpot;
}
......@@ -98,9 +98,9 @@ calc_adisit(long gridsize, long nlevel, double *pressure, Field tho, Field sao,
for (long levelID = 0; levelID < nlevel; ++levelID)
{
long offset = gridsize * levelID;
double *thoptr = tho.ptr + offset;
double *saoptr = sao.ptr + offset;
const long offset = gridsize * levelID;
const double *thoptr = tho.ptr + offset;
const double *saoptr = sao.ptr + offset;
double *tisptr = tis.ptr + offset;
for (long i = 0; i < gridsize; ++i)
......@@ -159,24 +159,23 @@ Adisit(void *process)
double *single;
cdoInitialize(process);
int ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
int ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
const int ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
const int ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
UNUSED(ADIPOT);
int operatorID = cdoOperatorID();
const int operatorID = cdoOperatorID();
if (operatorArgc() == 1) pin = parameter2double(operatorArgv()[0]);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int vlistID1 = cdoStreamInqVlist(streamID1);
int nvars = vlistNvars(vlistID1);
const int nvars = vlistNvars(vlistID1);
for (varID = 0; varID < nvars; varID++)
{
int code = vlistInqVarCode(vlistID1, varID);
if (code <= 0)
{
vlistInqVarName(vlistID1, varID, varname);
......@@ -209,13 +208,13 @@ Adisit(void *process)
if (saoID == -1) cdoAbort("Sea water salinity not found!");
if (thoID == -1) cdoAbort("Potential or Insitu temperature not found!");
int gridID = vlistGrid(vlistID1, 0);
size_t gridsize = vlist_check_gridsize(vlistID1);
const int gridID = vlistGrid(vlistID1, 0);
const size_t gridsize = vlist_check_gridsize(vlistID1);
int zaxisID = vlistInqVarZaxis(vlistID1, saoID);
int nlevel1 = zaxisInqSize(zaxisID);
const int nlevel1 = zaxisInqSize(zaxisID);
zaxisID = vlistInqVarZaxis(vlistID1, thoID);
int nlevel2 = zaxisInqSize(zaxisID);
const int nlevel2 = zaxisInqSize(zaxisID);
if (nlevel1 != nlevel2) cdoAbort("temperature and salinity have different number of levels!");
int nlevel = nlevel1;
......@@ -254,9 +253,9 @@ Adisit(void *process)
if (vlistInqVarDatatype(vlistID1, thoID) == CDI_DATATYPE_FLT64 && vlistInqVarDatatype(vlistID1, saoID) == CDI_DATATYPE_FLT64)
datatype = CDI_DATATYPE_FLT64;
int vlistID2 = vlistCreate();
const int vlistID2 = vlistCreate();
int tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
const int tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
if (operatorID == ADISIT)
{
vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
......@@ -284,12 +283,11 @@ Adisit(void *process)
vlistDefVarMissval(vlistID2, saoID2, sao.missval);
vlistDefVarDatatype(vlistID2, saoID2, datatype);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
pstreamDefVlist(streamID2, vlistID2);
int tsID = 0;
......
......@@ -70,9 +70,9 @@ Arith(void *process)
cdoOperatorAdd("setmiss", func_setmiss, 0, nullptr);
// clang-format on
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
bool opercplx = cdoOperatorF2(operatorID);
const int operatorID = cdoOperatorID();
const int operfunc = cdoOperatorF1(operatorID);
const bool opercplx = cdoOperatorF2(operatorID);
operatorCheckArgc(0);
int streamID1 = cdoStreamOpenRead(0);
......@@ -154,9 +154,9 @@ Arith(void *process)
if (filltype == FILL_NONE) vlistCompare(vlistID1, vlistID2, CMP_ALL);
size_t nwpv = (vlistNumber(vlistIDx1) == CDI_COMP && vlistNumber(vlistIDx2) == CDI_COMP) ? 2 : 1;
const size_t nwpv = (vlistNumber(vlistIDx1) == CDI_COMP && vlistNumber(vlistIDx2) == CDI_COMP) ? 2 : 1;
if (nwpv == 2 && !opercplx) cdoAbort("Fields with complex numbers are not supported by this operator!");
size_t gridsizemax = nwpv * vlistGridsizeMax(vlistIDx1);
const size_t gridsizemax = nwpv * vlistGridsizeMax(vlistIDx1);
fieldInit(field1);
fieldInit(field2);
......@@ -192,30 +192,30 @@ Arith(void *process)
if (filltype == FILL_TS)
{
int nvars = vlistNvars(vlistIDx2);
const int nvars = vlistNvars(vlistIDx2);
vardata.resize(nvars);
varnmiss.resize(nvars);
for (varID = 0; varID < nvars; varID++)
{
size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
int nlev = zaxisInqSize(vlistInqVarZaxis(vlistIDx2, varID));
const size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
const int nlev = zaxisInqSize(vlistInqVarZaxis(vlistIDx2, varID));
vardata[varID].resize(nlev * gridsize);
varnmiss[varID].resize(nlev);
}
}
}
int vlistID3 = vlistDuplicate(vlistIDx1);
const int vlistID3 = vlistDuplicate(vlistIDx1);
if (filltype == FILL_TS && vlistIDx1 != vlistID1)
{
int nvars = vlistNvars(vlistID1);
const int nvars = vlistNvars(vlistID1);
for (varID = 0; varID < nvars; varID++) vlistDefVarMissval(vlistID3, varID, vlistInqVarMissval(vlistID1, varID));
}
int taxisID3 = taxisDuplicate(taxisIDx1);
const int taxisID3 = taxisDuplicate(taxisIDx1);
vlistDefTaxis(vlistID3, taxisID3);
int streamID3 = cdoStreamOpenWrite(2);
const int streamID3 = cdoStreamOpenWrite(2);
pstreamDefVlist(streamID3, vlistID3);
int tsID = 0;
......@@ -265,7 +265,7 @@ Arith(void *process)
if (tsID == 0 || filltype == FILL_NONE || filltype == FILL_FILE || filltype == FILL_VARTS)
{
bool lstatus = nlevels2 > 1 ? varID == 0 : recID == 0;
const bool lstatus = nlevels2 > 1 ? varID == 0 : recID == 0;
if (lstatus || (filltype != FILL_VAR && filltype != FILL_VARTS))
{
pstreamInqRecord(streamIDx2, &varID2, &levelID2);
......@@ -277,23 +277,23 @@ Arith(void *process)
if (filltype == FILL_TS)
{
size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
size_t offset = gridsize * levelID;
const size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID));
const size_t offset = gridsize * levelID;
arrayCopy(gridsize, fieldx2->ptr, &vardata[varID][offset]);
varnmiss[varID][levelID] = fieldx2->nmiss;
}
else if (lstatus && (filltype == FILL_VAR || filltype == FILL_VARTS))
{
size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
size_t offset = gridsize * levelID2;
const size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
const size_t offset = gridsize * levelID2;
arrayCopy(gridsize, fieldx2->ptr, &vardata2[offset]);
varnmiss2[levelID2] = fieldx2->nmiss;
}
}
else if (filltype == FILL_TS)
{
size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID2));
size_t offset = gridsize * levelID;
const size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, varID2));
const size_t offset = gridsize * levelID;
arrayCopy(gridsize, &vardata[varID][offset], fieldx2->ptr);
fieldx2->nmiss = varnmiss[varID][levelID];
}
......@@ -304,8 +304,8 @@ Arith(void *process)
if (filltype == FILL_VAR || filltype == FILL_VARTS)
{
levelID2 = (nlevels2 > 1) ? levelID : 0;
size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
size_t offset = gridsize * levelID2;
const size_t gridsize = nwpv * gridInqSize(vlistInqVarGrid(vlistIDx2, 0));
const size_t offset = gridsize * levelID2;
arrayCopy(gridsize, &vardata2[offset], fieldx2->ptr);
fieldx2->nmiss = varnmiss2[levelID2];
fieldx2->grid = vlistInqVarGrid(vlistIDx2, 0);
......
......@@ -84,37 +84,37 @@ Arithc(void *process)
cdoOperatorAdd("mod", func_mod, 0, "divisor");
// clang-format on
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
bool opercplx = cdoOperatorF2(operatorID);
const int operatorID = cdoOperatorID();
const int operfunc = cdoOperatorF1(operatorID);
const bool opercplx = cdoOperatorF2(operatorID);
operatorInputArg(cdoOperatorEnter(operatorID));
if (operatorArgc() < 1) cdoAbort("Too few arguments!");
if (operatorArgc() > 2) cdoAbort("Too many arguments!");
double rconst = parameter2double(operatorArgv()[0]);
const double rconst = parameter2double(operatorArgv()[0]);
double rconstcplx[2];
rconstcplx[0] = rconst;
rconstcplx[1] = 0;
if (operatorArgc() == 2) rconstcplx[1] = parameter2double(operatorArgv()[1]);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = vlistDuplicate(vlistID1);
std::vector<bool> vars;
fill_vars(vlistID1, vars);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
pstreamDefVlist(streamID2, vlistID2);
size_t nwpv = (vlistNumber(vlistID1) == CDI_COMP) ? 2 : 1;
const size_t nwpv = (vlistNumber(vlistID1) == CDI_COMP) ? 2 : 1;
if (nwpv == 2 && !opercplx) cdoAbort("Fields with complex numbers are not supported by this operator!");
size_t gridsizemax = nwpv * vlistGridsizeMax(vlistID1);
const size_t gridsizemax = nwpv * vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
......@@ -134,8 +134,8 @@ Arithc(void *process)
if (vars[varID])
{
int gridID = vlistInqVarGrid(vlistID1, varID);
size_t gridsize = nwpv * gridInqSize(gridID);
const int gridID = vlistInqVarGrid(vlistID1, varID);
const size_t gridsize = nwpv * gridInqSize(gridID);
field.grid = gridID;
field.missval = vlistInqVarMissval(vlistID1, varID);
......
......@@ -35,28 +35,21 @@
static double
dayofyear(int calendar, int64_t vdate, int vtime)
{
int month_360[12] = { 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30 };
int month_365[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
int month_366[12] = { 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
int *dpm;
const int month_360[12] = { 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30 };
const int month_365[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
const int month_366[12] = { 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
int year, month, day;
int hour, minute, second;
double doy = 0;
cdiDecodeDate(vdate, &year, &month, &day);
cdiDecodeTime(vtime, &hour, &minute, &second);
int dpy = days_per_year(calendar, year);
const int dpy = days_per_year(calendar, year);
double doy = 0;
for (int im = 1; im < month; ++im)
{
if (dpy == 360)
dpm = month_360;
else if (dpy == 365)
dpm = month_365;
else
dpm = month_366;
const int *dpm = (dpy == 360) ? month_360 : ((dpy == 365) ? month_365 : month_366);
if (im >= 1 && im <= 12) doy += dpm[im - 1];
}
......@@ -87,39 +80,38 @@ Arithdays(void *process)
int MULDOY = cdoOperatorAdd("muldoy", func_mul, 0, nullptr);
// clang-format on
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
int operfunc2 = cdoOperatorF2(operatorID);
const int operatorID = cdoOperatorID();
const int operfunc = cdoOperatorF1(operatorID);
const int operfunc2 = cdoOperatorF2(operatorID);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = vlistDuplicate(vlistID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int calendar = taxisInqCalendar(taxisID1);
const int calendar = taxisInqCalendar(taxisID1);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
pstreamDefVlist(streamID2, vlistID2);
size_t gridsize = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsize * sizeof(double));
field.ptr = (double *) Malloc(gridsizemax * sizeof(double));
field.weight = nullptr;
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
int64_t vdate = taxisInqVdate(taxisID1);
int vtime = taxisInqVtime(taxisID1);
const int64_t vdate = taxisInqVdate(taxisID1);
const int vtime = taxisInqVtime(taxisID1);
taxisCopyTimestep(taxisID2, taxisID1);
pstreamDefTimestep(streamID2, tsID);
cdiDecodeDate(vdate, &year, &month, &day);
......@@ -130,10 +122,7 @@ Arithdays(void *process)
}
else
{
if (operfunc2 == func_month)
rconst = days_per_month(calendar, year, month);
else
rconst = days_per_year(calendar, year);
rconst = (operfunc2 == func_month) ? days_per_month(calendar, year, month) : days_per_year(calendar, year);
}
if (Options::cdoVerbose) cdoPrint("calendar %d year %d month %d result %g", calendar, year, month, rconst);
......@@ -152,6 +141,7 @@ Arithdays(void *process)
pstreamDefRecord(streamID2, varID, levelID);
pstreamWriteRecord(streamID2, field.ptr, field.nmiss);
}
tsID++;
}
......
......@@ -44,29 +44,28 @@ Arithlat(void *process)
cdoOperatorAdd("mulcoslat", func_mul, 0, nullptr);
cdoOperatorAdd("divcoslat", func_div, 0, nullptr);
int operatorID = cdoOperatorID();
int operfunc = cdoOperatorF1(operatorID);
const int operatorID = cdoOperatorID();
const int operfunc = cdoOperatorF1(operatorID);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = vlistDuplicate(vlistID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
pstreamDefVlist(streamID2, vlistID2);
size_t gridsizemax = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
std::vector<double> array(gridsizemax);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
taxisCopyTimestep(taxisID2, taxisID1);
pstreamDefTimestep(streamID2, tsID);
for (int recID = 0; recID < nrecs; recID++)
......@@ -81,15 +80,15 @@ Arithlat(void *process)
{
gridID0 = gridID;
int gridtype = gridInqType(gridID);
int projtype = (gridtype == GRID_PROJECTION) ? gridInqProjType(gridID) : -1;
const int gridtype = gridInqType(gridID);
const int projtype = (gridtype == GRID_PROJECTION) ? gridInqProjType(gridID) : -1;
if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_LCC)
{
gridID = gridToCurvilinear(gridID, 0);
}
else if (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED)
{
/* No conversion necessary */
// No conversion necessary
}
else if (gridtype == GRID_GME)
{
......@@ -109,10 +108,8 @@ Arithlat(void *process)
scale.resize(gridsize);
gridInqYvals(gridID, &scale[0]);
/* Convert lat/lon units if required */
// Convert lat/lon units if required
gridInqXunits(gridID, units);
grid_to_radian(units, gridsize, &scale[0], "grid latitudes");
if (operfunc == func_mul)
......
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