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

Y*pctl: changed to Field vector.

parent 465dac38
......@@ -72,7 +72,7 @@ Runpctl(void *process)
FieldVector3D vars1(ndates + 1);
std::vector<double> array(ndates);
for (int its = 0; its < ndates; its++) fieldsFromVlist(vlistID1, vars1[its], FIELD_PTR);
for (int its = 0; its < ndates; its++) fieldsFromVlist(vlistID1, vars1[its], FIELD_VEC);
int tsID;
for (tsID = 0; tsID < ndates; tsID++)
......@@ -93,7 +93,7 @@ Runpctl(void *process)
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
}
cdoReadRecord(streamID1, vars1[tsID][varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[tsID][varID][levelID].vec.data(), &nmiss);
vars1[tsID][varID][levelID].nmiss = nmiss;
}
}
......@@ -111,27 +111,28 @@ Runpctl(void *process)
for (levelID = 0; levelID < nlevels; levelID++)
{
Field &rvars1_0 = vars1[0][varID][levelID];
nmiss = 0;
for (size_t i = 0; i < gridsize; i++)
{
int j = 0;
for (int inp = 0; inp < ndates; inp++)
{
const double val = vars1[inp][varID][levelID].ptr[i];
const double val = vars1[inp][varID][levelID].vec[i];
if (!DBL_IS_EQUAL(val, missval)) array[j++] = val;
}
if (j > 0)
{
vars1[0][varID][levelID].ptr[i] = percentile(array.data(), j, pn);
rvars1_0.vec[i] = percentile(array.data(), j, pn);
}
else
{
vars1[0][varID][levelID].ptr[i] = missval;
rvars1_0.vec[i] = missval;
nmiss++;
}
}
vars1[0][varID][levelID].nmiss = nmiss;
rvars1_0.nmiss = nmiss;
}
}
......@@ -146,7 +147,7 @@ Runpctl(void *process)
const int levelID = recList[recID].levelID;
cdoDefRecord(streamID2, varID, levelID);
cdoWriteRecord(streamID2, vars1[0][varID][levelID].ptr, vars1[0][varID][levelID].nmiss);
cdoWriteRecord(streamID2, vars1[0][varID][levelID].vec.data(), vars1[0][varID][levelID].nmiss);
}
otsID++;
......@@ -164,15 +165,13 @@ Runpctl(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, vars1[ndates - 1][varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[ndates - 1][varID][levelID].vec.data(), &nmiss);
vars1[ndates - 1][varID][levelID].nmiss = nmiss;
}
tsID++;
}
for (int its = 0; its < ndates; its++) fieldsFree(vlistID1, vars1[its]);
cdoStreamClose(streamID2);
cdoStreamClose(streamID1);
......
......@@ -42,7 +42,6 @@ Seaspctl(void *process)
int nlevels;
int oldmon = 0;
int season_start;
double missval;
cdoInitialize(process);
......@@ -87,32 +86,21 @@ Seaspctl(void *process)
dtlist.setStat(timestat_date);
dtlist.setCalendar(taxisInqCalendar(taxisID1));
size_t gridsize = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsize * sizeof(double));
field.resize(gridsizemax);
Field **vars1 = (Field **) Malloc(nvars * sizeof(Field *));
FieldVector2D vars1;
fieldsFromVlist(vlistID1, vars1, FIELD_VEC);
HISTOGRAM_SET *hset = hsetCreate(nvars);
for (varID = 0; varID < nvars; varID++)
{
gridID = vlistInqVarGrid(vlistID1, varID);
gridsize = gridInqSize(gridID);
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
missval = vlistInqVarMissval(vlistID1, varID);
vars1[varID] = (Field *) Malloc(nlevels * sizeof(Field));
hsetCreateVarLevels(hset, varID, nlevels, gridID);
for (levelID = 0; levelID < nlevels; levelID++)
{
vars1[varID][levelID].grid = gridID;
vars1[varID][levelID].nmiss = 0;
vars1[varID][levelID].missval = missval;
vars1[varID][levelID].ptr = (double *) Malloc(gridsize * sizeof(double));
}
}
int tsID = 0;
......@@ -133,15 +121,14 @@ Seaspctl(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID2, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
}
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID3, &varID, &levelID);
cdoReadRecord(streamID3, field.ptr, &nmiss);
field.nmiss = nmiss;
cdoReadRecord(streamID3, field.vec.data(), &field.nmiss);
field.grid = vars1[varID][levelID].grid;
field.missval = vars1[varID][levelID].missval;
......@@ -183,7 +170,7 @@ Seaspctl(void *process)
recList[recID].levelID = levelID;
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
}
cdoReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hsetAddVarLevelValues(hset, varID, levelID, &vars1[varID][levelID]);
......@@ -214,25 +201,15 @@ Seaspctl(void *process)
const int varID = recList[recID].varID;
const int levelID = recList[recID].levelID;
cdoDefRecord(streamID4, varID, levelID);
cdoWriteRecord(streamID4, vars1[varID][levelID].ptr, vars1[varID][levelID].nmiss);
cdoWriteRecord(streamID4, vars1[varID][levelID].vec.data(), vars1[varID][levelID].nmiss);
}
if (nrecs == 0) break;
otsID++;
}
for (varID = 0; varID < nvars; varID++)
{
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++) Free(vars1[varID][levelID].ptr);
Free(vars1[varID]);
}
Free(vars1);
hsetDestroy(hset);
if (field.ptr) Free(field.ptr);
cdoStreamClose(streamID4);
cdoStreamClose(streamID3);
cdoStreamClose(streamID2);
......
......@@ -29,7 +29,6 @@
#include "cdo_int.h"
#include "param_conversion.h"
#include "percentiles_hist.h"
#include "percentiles.h"
#include "datetime.h"
......@@ -45,54 +44,53 @@ timpctl(int operatorID)
int nlevels;
operatorInputArg("percentile number");
double pn = parameter2double(operatorArgv()[0]);
const double pn = parameter2double(operatorArgv()[0]);
percentile_check_number(pn);
int cmplen = DATE_LEN - cdoOperatorF2(operatorID);
const int cmplen = DATE_LEN - cdoOperatorF2(operatorID);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int streamID2 = cdoStreamOpenRead(cdoStreamName(1));
int streamID3 = cdoStreamOpenRead(cdoStreamName(2));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID2 = cdoStreamOpenRead(cdoStreamName(1));
const int streamID3 = cdoStreamOpenRead(cdoStreamName(2));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = cdoStreamInqVlist(streamID2);
int vlistID3 = cdoStreamInqVlist(streamID3);
int vlistID4 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = cdoStreamInqVlist(streamID2);
const int vlistID3 = cdoStreamInqVlist(streamID3);
const int vlistID4 = vlistDuplicate(vlistID1);
vlistCompare(vlistID1, vlistID2, CMP_ALL);
vlistCompare(vlistID1, vlistID3, CMP_ALL);
if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(vlistID4, 1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = vlistInqTaxis(vlistID2);
int taxisID3 = vlistInqTaxis(vlistID3);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = vlistInqTaxis(vlistID2);
const int taxisID3 = vlistInqTaxis(vlistID3);
/* TODO - check that time axes 2 and 3 are equal */
int taxisID4 = taxisDuplicate(taxisID1);
const int taxisID4 = taxisDuplicate(taxisID1);
taxisWithBounds(taxisID4);
vlistDefTaxis(vlistID4, taxisID4);
int streamID4 = cdoStreamOpenWrite(cdoStreamName(3));
const int streamID4 = cdoStreamOpenWrite(cdoStreamName(3));
cdoDefVlist(streamID4, vlistID4);
int nvars = vlistNvars(vlistID1);
const int nvars = vlistNvars(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
const int maxrecs = vlistNrecs(vlistID1);
std::vector<RecordInfo> recList(maxrecs);
DateTimeList dtlist;
dtlist.setStat(timestat_date);
dtlist.setCalendar(taxisInqCalendar(taxisID1));
size_t gridsize = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsize * sizeof(double));
field.resize(gridsizemax);
FieldVector2D vars1;
fieldsFromVlist(vlistID1, vars1, FIELD_PTR);
fieldsFromVlist(vlistID1, vars1, FIELD_VEC);
HISTOGRAM_SET *hset = hsetCreate(nvars);
for (varID = 0; varID < nvars; varID++)
......@@ -111,25 +109,24 @@ timpctl(int operatorID)
if (nrecs != cdoStreamInqTimestep(streamID3, otsID))
cdoAbort("Number of records at time step %d of %s and %s differ!", otsID + 1, cdoGetStreamName(1), cdoGetStreamName(2));
int64_t vdate2 = taxisInqVdate(taxisID2);
int vtime2 = taxisInqVtime(taxisID2);
int64_t vdate3 = taxisInqVdate(taxisID3);
int vtime3 = taxisInqVtime(taxisID3);
const int64_t vdate2 = taxisInqVdate(taxisID2);
const int vtime2 = taxisInqVtime(taxisID2);
const int64_t vdate3 = taxisInqVdate(taxisID3);
const int vtime3 = taxisInqVtime(taxisID3);
if (vdate2 != vdate3 || vtime2 != vtime3)
cdoAbort("Verification dates at time step %d of %s and %s differ!", otsID + 1, cdoGetStreamName(1), cdoGetStreamName(2));
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID2, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
}
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID3, &varID, &levelID);
cdoReadRecord(streamID3, field.ptr, &nmiss);
field.nmiss = nmiss;
cdoReadRecord(streamID3, field.vec.data(), &field.nmiss);
field.grid = vars1[varID][levelID].grid;
field.missval = vars1[varID][levelID].missval;
......@@ -140,8 +137,8 @@ timpctl(int operatorID)
while (nrecs && (nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
dtlist.taxisInqTimestep(taxisID1, nsets);
int64_t vdate1 = dtlist.getVdate(nsets);
int vtime1 = dtlist.getVtime(nsets);
const int64_t vdate1 = dtlist.getVdate(nsets);
const int vtime1 = dtlist.getVtime(nsets);
if (nsets == 0) SET_DATE(indate2, vdate1, vtime1);
SET_DATE(indate1, vdate1, vtime1);
......@@ -157,7 +154,7 @@ timpctl(int operatorID)
recList[recID].levelID = levelID;
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
}
cdoReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hsetAddVarLevelValues(hset, varID, levelID, &vars1[varID][levelID]);
......@@ -185,21 +182,18 @@ timpctl(int operatorID)
{
if (otsID && recList[recID].lconst) continue;
int varID = recList[recID].varID;
int levelID = recList[recID].levelID;
const int varID = recList[recID].varID;
const int levelID = recList[recID].levelID;
cdoDefRecord(streamID4, varID, levelID);
cdoWriteRecord(streamID4, vars1[varID][levelID].ptr, vars1[varID][levelID].nmiss);
cdoWriteRecord(streamID4, vars1[varID][levelID].vec.data(), vars1[varID][levelID].nmiss);
}
if (nrecs == 0) break;
otsID++;
}
fieldsFree(vlistID1, vars1);
hsetDestroy(hset);
if (field.ptr) Free(field.ptr);
cdoStreamClose(streamID4);
cdoStreamClose(streamID3);
cdoStreamClose(streamID2);
......@@ -209,8 +203,6 @@ timpctl(int operatorID)
void *
Timpctl(void *process)
{
int operatorID;
cdoInitialize(process);
// clang-format off
......@@ -221,7 +213,7 @@ Timpctl(void *process)
cdoOperatorAdd("hourpctl", func_pctl, 4, nullptr);
// clang-format on
operatorID = cdoOperatorID();
const int operatorID = cdoOperatorID();
timpctl(operatorID);
......
......@@ -25,7 +25,6 @@
#include "cdo_int.h"
#include "param_conversion.h"
#include "percentiles_hist.h"
#include "percentiles.h"
#include "datetime.h"
......@@ -46,59 +45,58 @@ Timselpctl(void *process)
operatorInputArg("percentile number, nsets <,noffset <,nskip>>");
int nargc = operatorArgc();
const int nargc = operatorArgc();
if (nargc < 2) cdoAbort("Too few arguments! Need %d found %d.", 2, nargc);
double pn = parameter2double(operatorArgv()[0]);
const double pn = parameter2double(operatorArgv()[0]);
percentile_check_number(pn);
int ndates = parameter2int(operatorArgv()[1]);
const int ndates = parameter2int(operatorArgv()[1]);
int noffset = 0, nskip = 0;
if (nargc > 2) noffset = parameter2int(operatorArgv()[2]);
if (nargc > 3) nskip = parameter2int(operatorArgv()[3]);
if (Options::cdoVerbose) cdoPrint("nsets = %d, noffset = %d, nskip = %d", ndates, noffset, nskip);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int streamID2 = cdoStreamOpenRead(cdoStreamName(1));
int streamID3 = cdoStreamOpenRead(cdoStreamName(2));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID2 = cdoStreamOpenRead(cdoStreamName(1));
const int streamID3 = cdoStreamOpenRead(cdoStreamName(2));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = cdoStreamInqVlist(streamID2);
int vlistID3 = cdoStreamInqVlist(streamID3);
int vlistID4 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = cdoStreamInqVlist(streamID2);
const int vlistID3 = cdoStreamInqVlist(streamID3);
const int vlistID4 = vlistDuplicate(vlistID1);
vlistCompare(vlistID1, vlistID2, CMP_ALL);
vlistCompare(vlistID1, vlistID3, CMP_ALL);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = vlistInqTaxis(vlistID2);
int taxisID3 = vlistInqTaxis(vlistID3);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = vlistInqTaxis(vlistID2);
const int taxisID3 = vlistInqTaxis(vlistID3);
/* TODO - check that time axes 2 and 3 are equal */
int taxisID4 = taxisDuplicate(taxisID1);
const int taxisID4 = taxisDuplicate(taxisID1);
taxisWithBounds(taxisID4);
vlistDefTaxis(vlistID4, taxisID4);
int streamID4 = cdoStreamOpenWrite(cdoStreamName(3));
const int streamID4 = cdoStreamOpenWrite(cdoStreamName(3));
cdoDefVlist(streamID4, vlistID4);
int nvars = vlistNvars(vlistID1);
const int nvars = vlistNvars(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
const int maxrecs = vlistNrecs(vlistID1);
std::vector<RecordInfo> recList(maxrecs);
DateTimeList dtlist;
dtlist.setStat(timestat_date);
dtlist.setCalendar(taxisInqCalendar(taxisID1));
size_t gridsize = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsize * sizeof(double));
field.resize(gridsizemax);
FieldVector2D vars1;
fieldsFromVlist(vlistID1, vars1, FIELD_PTR);
fieldsFromVlist(vlistID1, vars1, FIELD_VEC);
HISTOGRAM_SET *hset = hsetCreate(nvars);
for (varID = 0; varID < nvars; varID++)
......@@ -140,25 +138,24 @@ Timselpctl(void *process)
if (nrecs != cdoStreamInqTimestep(streamID3, otsID))
cdoAbort("Number of records at time step %d of %s and %s differ!", otsID + 1, cdoGetStreamName(1), cdoGetStreamName(2));
int64_t vdate2 = taxisInqVdate(taxisID2);
int vtime2 = taxisInqVtime(taxisID2);
int64_t vdate3 = taxisInqVdate(taxisID3);
int vtime3 = taxisInqVtime(taxisID3);
const int64_t vdate2 = taxisInqVdate(taxisID2);
const int vtime2 = taxisInqVtime(taxisID2);
const int64_t vdate3 = taxisInqVdate(taxisID3);
const int vtime3 = taxisInqVtime(taxisID3);
if (vdate2 != vdate3 || vtime2 != vtime3)
cdoAbort("Verification dates at time step %d of %s and %s differ!", otsID + 1, cdoGetStreamName(1), cdoGetStreamName(2));
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID2, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
}
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID3, &varID, &levelID);
cdoReadRecord(streamID3, field.ptr, &nmiss);
field.nmiss = nmiss;
cdoReadRecord(streamID3, field.vec.data(), &field.nmiss);
field.grid = vars1[varID][levelID].grid;
field.missval = vars1[varID][levelID].missval;
......@@ -185,7 +182,7 @@ Timselpctl(void *process)
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
}
cdoReadRecord(streamID1, vars1[varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hsetAddVarLevelValues(hset, varID, levelID, &vars1[varID][levelID]);
......@@ -212,10 +209,10 @@ Timselpctl(void *process)
{
if (otsID && recList[recID].lconst) continue;
int varID = recList[recID].varID;
int levelID = recList[recID].levelID;
const int varID = recList[recID].varID;
const int levelID = recList[recID].levelID;
cdoDefRecord(streamID4, varID, levelID);
cdoWriteRecord(streamID4, vars1[varID][levelID].ptr, vars1[varID][levelID].nmiss);
cdoWriteRecord(streamID4, vars1[varID][levelID].vec.data(), vars1[varID][levelID].nmiss);
}
if (nrecs == 0) break;
......@@ -233,11 +230,8 @@ Timselpctl(void *process)
LABEL_END:
fieldsFree(vlistID1, vars1);
hsetDestroy(hset);
if (field.ptr) Free(field.ptr);
cdoStreamClose(streamID4);
cdoStreamClose(streamID3);
cdoStreamClose(streamID2);
......
......@@ -94,8 +94,7 @@ Ydaypctl(void *process)
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsizemax * sizeof(double));
field.resize(gridsizemax);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID2, tsID)))
......@@ -120,7 +119,7 @@ Ydaypctl(void *process)
if (!vars1[dayoy].size())
{
fieldsFromVlist(vlistID1, vars1[dayoy], FIELD_PTR);
fieldsFromVlist(vlistID1, vars1[dayoy], FIELD_VEC);
hsets[dayoy] = hsetCreate(nvars);
for (varID = 0; varID < nvars; varID++)
......@@ -135,15 +134,14 @@ Ydaypctl(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, vars1[dayoy][varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID2, vars1[dayoy][varID][levelID].vec.data(), &nmiss);
vars1[dayoy][varID][levelID].nmiss = nmiss;
}
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID3, &varID, &levelID);
cdoReadRecord(streamID3, field.ptr, &nmiss);
field.nmiss = nmiss;
cdoReadRecord(streamID3, field.vec.data(), &field.nmiss);
field.grid = vars1[dayoy][varID][levelID].grid;
field.missval = vars1[dayoy][varID][levelID].missval;
......@@ -182,7 +180,7 @@ Ydaypctl(void *process)
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
}
cdoReadRecord(streamID1, vars1[dayoy][varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID1, vars1[dayoy][varID][levelID].vec.data(), &nmiss);
vars1[dayoy][varID][levelID].nmiss = nmiss;
hsetAddVarLevelValues(hsets[dayoy], varID, levelID, &vars1[dayoy][varID][levelID]);
......@@ -220,7 +218,7 @@ Ydaypctl(void *process)
const int varID = recList[recID].varID;
const int levelID = recList[recID].levelID;
cdoDefRecord(streamID4, varID, levelID);
cdoWriteRecord(streamID4, vars1[dayoy][varID][levelID].ptr, vars1[dayoy][varID][levelID].nmiss);
cdoWriteRecord(streamID4, vars1[dayoy][varID][levelID].vec.data(), vars1[dayoy][varID][levelID].nmiss);
}
otsID++;
......@@ -230,13 +228,10 @@ Ydaypctl(void *process)
{
if (vars1[dayoy].size())
{
fieldsFree(vlistID1, vars1[dayoy]);
hsetDestroy(hsets[dayoy]);
}
}
if (field.ptr) Free(field.ptr);
cdoStreamClose(streamID4);
cdoStreamClose(streamID3);
cdoStreamClose(streamID2);
......
......@@ -104,13 +104,12 @@ Ydrunpctl(void *process)
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
Field field;
fieldInit(field);
field.ptr = (double *) Malloc(gridsizemax * sizeof(double));
field.resize(gridsizemax);
std::vector<CdoDateTime> datetime(ndates + 1);
FieldVector3D vars1(ndates + 1);
for (its = 0; its < ndates; its++) fieldsFromVlist(vlistID1, vars1[its], FIELD_PTR);
for (its = 0; its < ndates; its++) fieldsFromVlist(vlistID1, vars1[its], FIELD_VEC);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID2, tsID)))
......@@ -136,7 +135,7 @@ Ydrunpctl(void *process)
if (!vars2[dayoy].size())
{
fieldsFromVlist(vlistID2, vars2[dayoy], FIELD_PTR);
fieldsFromVlist(vlistID2, vars2[dayoy], FIELD_VEC);
hsets[dayoy] = hsetCreate(nvars);
for (varID = 0; varID < nvars; varID++)
......@@ -151,15 +150,14 @@ Ydrunpctl(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, vars2[dayoy][varID][levelID].ptr, &nmiss);
cdoReadRecord(streamID2, vars2[dayoy][varID][levelID].vec.data(), &nmiss);
vars2[dayoy][varID][levelID].nmiss = nmiss;
}
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID3, &varID, &levelID);
cdoReadRecord(streamID3, field.ptr, &nmiss);
field.nmiss = nmiss;
cdoReadRecord(streamID3, field.vec.data(), &field.nmiss);
field.grid = vars2[dayoy][varID][levelID].grid;
field.missval = vars2[dayoy][varID][levelID].missval;
......@@ -188,7 +186,7 @@ Ydrunpctl(void *process)