Skip to content
Snippets Groups Projects
Commit 4074e3fa authored by Oliver Heidmann's avatar Oliver Heidmann Committed by Uwe Schulzweida
Browse files

class rework Varsstat.cc

parent 269189b0
No related branches found
No related tags found
1 merge request!96Develop
Hide whitespace changes
Inline Side-by-side
Showing
with 200 additions and 162 deletions
src/Varsstat.cc
+ 200
− 162
View file @ 4074e3fa
......@@ -80,181 +80,219 @@ add_operators(void)
// clang-format on
}
void *
Varsstat(void *process)
class ModuleVarsstat
{
cdo_initialize(process);
add_operators();
auto operatorID = cdo_operator_id();
auto operfunc = cdo_operator_f1(operatorID);
operator_check_argc(0);
auto lrange = (operfunc == FieldFunc_Range);
auto lmean = (operfunc == FieldFunc_Mean || operfunc == FieldFunc_Avg);
auto lstd = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Std1);
auto lvarstd = (lstd || operfunc == FieldFunc_Var || operfunc == FieldFunc_Var1);
auto lvars2 = (lvarstd || lrange);
int divisor = (operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
auto field2_stdvar_func = lstd ? field2_std : field2_var;
auto fieldc_stdvar_func = lstd ? fieldc_std : fieldc_var;
auto streamID1 = cdo_open_read(0);
auto vlistID1 = cdo_stream_inq_vlist(streamID1);
VarList varList1;
varListInit(varList1, vlistID1);
CdoStreamID streamID1;
int taxisID1;
check_unique_zaxis(vlistID1);
auto zaxisID = vlistZaxis(vlistID1, 0);
auto nlevels = zaxisInqSize(zaxisID);
check_unique_gridsize(vlistID1);
auto gridID = vlistGrid(vlistID1, 0);
auto gridsize = gridInqSize(gridID);
auto timetype = varList1[0].timetype;
auto nvars = vlistNvars(vlistID1);
for (int varID = 1; varID < nvars; ++varID)
{
if (timetype != varList1[varID].timetype) cdo_abort("Number of timesteps differ!");
}
auto vlistID2 = vlistCreate();
vlistDefNtsteps(vlistID2, vlistNtsteps(vlistID1));
CdoStreamID streamID2;
int taxisID2;
int vlistID2;
auto varID2 = vlistDefVar(vlistID2, gridID, zaxisID, timetype);
set_attributes(varList1, vlistID2, varID2, operatorID);
bool lrange;
bool lvarstd;
bool lmean;
bool lstd;
double divisor;
auto taxisID1 = vlistInqTaxis(vlistID1);
auto taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int nlevels;
auto streamID2 = cdo_open_write(1);
cdo_def_vlist(streamID2, vlistID2);
int operfunc;
Field field;
FieldVector vars1, samp1, vars2;
VarList varList1;
FieldVector vars1(nlevels), samp1(nlevels), vars2;
if (lvars2) vars2.resize(nlevels);
for (int levelID = 0; levelID < nlevels; ++levelID)
{
auto missval = varList1[0].missval;
samp1[levelID].grid = gridID;
samp1[levelID].missval = missval;
samp1[levelID].memType = MemType::Double;
vars1[levelID].grid = gridID;
vars1[levelID].missval = missval;
vars1[levelID].memType = MemType::Double;
vars1[levelID].resize(gridsize);
if (lvars2)
{
vars2[levelID].grid = gridID;
vars2[levelID].missval = missval;
vars2[levelID].memType = MemType::Double;
vars2[levelID].resize(gridsize);
}
}
int tsID = 0;
while (true)
{
auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
if (nrecs == 0) break;
cdo_taxis_copy_timestep(taxisID2, taxisID1);
cdo_def_timestep(streamID2, tsID);
for (int recID = 0; recID < nrecs; ++recID)
{
int varID, levelID;
cdo_inq_record(streamID1, &varID, &levelID);
auto &rsamp1 = samp1[levelID];
auto &rvars1 = vars1[levelID];
rvars1.nsamp++;
if (lrange) vars2[levelID].nsamp++;
if (varID == 0)
{
cdo_read_record(streamID1, rvars1);
if (lrange)
{
vars2[levelID].nmiss = rvars1.nmiss;
vars2[levelID].vec_d = rvars1.vec_d;
}
if (lvarstd) field2_moq(vars2[levelID], rvars1);
if (rvars1.nmiss || !rsamp1.empty())
{
if (rsamp1.empty()) rsamp1.resize(rvars1.size);
field2_vinit(rsamp1, rvars1);
}
}
else
{
field.init(varList1[varID]);
cdo_read_record(streamID1, field);
if (field.nmiss || !rsamp1.empty())
{
if (rsamp1.empty()) rsamp1.resize(rvars1.size, rvars1.nsamp);
field2_vincr(rsamp1, field);
}
// clang-format off
public:
void
init(void *process)
{
cdo_initialize(process);
add_operators();
auto operatorID = cdo_operator_id();
operfunc = cdo_operator_f1(operatorID);
operator_check_argc(0);
lrange = (operfunc == FieldFunc_Range);
lmean = (operfunc == FieldFunc_Mean || operfunc == FieldFunc_Avg);
lstd = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Std1);
lvarstd = (lstd || operfunc == FieldFunc_Var || operfunc == FieldFunc_Var1);
auto lvars2 = (lvarstd || lrange);
divisor = (operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
streamID1 = cdo_open_read(0);
auto vlistID1 = cdo_stream_inq_vlist(streamID1);
varListInit(varList1, vlistID1);
check_unique_zaxis(vlistID1);
auto zaxisID = vlistZaxis(vlistID1, 0);
nlevels = zaxisInqSize(zaxisID);
check_unique_gridsize(vlistID1);
auto gridID = vlistGrid(vlistID1, 0);
auto gridsize = gridInqSize(gridID);
auto timetype = varList1[0].timetype;
auto nvars = vlistNvars(vlistID1);
for (int varID = 1; varID < nvars; ++varID)
{
if (timetype != varList1[varID].timetype) cdo_abort("Number of timesteps differ!");
}
vlistID2 = vlistCreate();
vlistDefNtsteps(vlistID2, vlistNtsteps(vlistID1));
auto varID2 = vlistDefVar(vlistID2, gridID, zaxisID, timetype);
set_attributes(varList1, vlistID2, varID2, operatorID);
taxisID1 = vlistInqTaxis(vlistID1);
taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
streamID2 = cdo_open_write(1);
cdo_def_vlist(streamID2, vlistID2);
vars1 = FieldVector(nlevels);
samp1 = FieldVector(nlevels);
if (lvars2) vars2.resize(nlevels);
for (int levelID = 0; levelID < nlevels; ++levelID)
{
auto missval = varList1[0].missval;
samp1[levelID].grid = gridID;
samp1[levelID].missval = missval;
samp1[levelID].memType = MemType::Double;
vars1[levelID].grid = gridID;
vars1[levelID].missval = missval;
vars1[levelID].memType = MemType::Double;
vars1[levelID].resize(gridsize);
if (lvars2)
{
vars2[levelID].grid = gridID;
vars2[levelID].missval = missval;
vars2[levelID].memType = MemType::Double;
vars2[levelID].resize(gridsize);
}
}
}
void
run()
{
auto field2_stdvar_func = lstd ? field2_std : field2_var;
auto fieldc_stdvar_func = lstd ? fieldc_std : fieldc_var;
int tsID = 0;
while (true)
{
auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
if (nrecs == 0) break;
cdo_taxis_copy_timestep(taxisID2, taxisID1);
cdo_def_timestep(streamID2, tsID);
for (int recID = 0; recID < nrecs; ++recID)
{
int varID, levelID;
cdo_inq_record(streamID1, &varID, &levelID);
auto &rsamp1 = samp1[levelID];
auto &rvars1 = vars1[levelID];
rvars1.nsamp++;
if (lrange) vars2[levelID].nsamp++;
if (varID == 0)
{
cdo_read_record(streamID1, rvars1);
if (lrange)
{
vars2[levelID].nmiss = rvars1.nmiss;
vars2[levelID].vec_d = rvars1.vec_d;
}
if (lvarstd) field2_moq(vars2[levelID], rvars1);
if (rvars1.nmiss || !rsamp1.empty())
{
if (rsamp1.empty()) rsamp1.resize(rvars1.size);
field2_vinit(rsamp1, rvars1);
}
}
else
{
field.init(varList1[varID]);
cdo_read_record(streamID1, field);
if (field.nmiss || !rsamp1.empty())
{
if (rsamp1.empty()) rsamp1.resize(rvars1.size, rvars1.nsamp);
field2_vincr(rsamp1, field);
}
// clang-format off
if (lvarstd) field2_sumsumq(rvars1, vars2[levelID], field);
else if (lrange) field2_maxmin(rvars1, vars2[levelID], field);
else field2_function(rvars1, field, operfunc);
// clang-format on
}
}
// clang-format on
}
}
for (int levelID = 0; levelID < nlevels; ++levelID)
{
const auto &rsamp1 = samp1[levelID];
auto &rvars1 = vars1[levelID];
if (rvars1.nsamp)
{
if (lmean)
{
if (!rsamp1.empty())
field2_div(rvars1, rsamp1);
else
fieldc_div(rvars1, (double) rvars1.nsamp);
}
else if (lvarstd)
{
if (!rsamp1.empty())
field2_stdvar_func(rvars1, vars2[levelID], rsamp1, divisor);
else
fieldc_stdvar_func(rvars1, vars2[levelID], rsamp1.nsamp, divisor);
}
else if (lrange) { field2_sub(rvars1, vars2[levelID]); }
cdo_def_record(streamID2, 0, levelID);
cdo_write_record(streamID2, rvars1);
rvars1.nsamp = 0;
}
}
tsID++;
}
}
void
close()
{
cdo_stream_close(streamID2);
cdo_stream_close(streamID1);
vlistDestroy(vlistID2);
cdo_finish();
}
};
for (int levelID = 0; levelID < nlevels; ++levelID)
{
const auto &rsamp1 = samp1[levelID];
auto &rvars1 = vars1[levelID];
if (rvars1.nsamp)
{
if (lmean)
{
if (!rsamp1.empty())
field2_div(rvars1, rsamp1);
else
fieldc_div(rvars1, (double) rvars1.nsamp);
}
else if (lvarstd)
{
if (!rsamp1.empty())
field2_stdvar_func(rvars1, vars2[levelID], rsamp1, divisor);
else
fieldc_stdvar_func(rvars1, vars2[levelID], rsamp1.nsamp, divisor);
}
else if (lrange) { field2_sub(rvars1, vars2[levelID]); }
cdo_def_record(streamID2, 0, levelID);
cdo_write_record(streamID2, rvars1);
rvars1.nsamp = 0;
}
}
tsID++;
}
cdo_stream_close(streamID2);
cdo_stream_close(streamID1);
vlistDestroy(vlistID2);
void *
Varsstat(void *process)
{
ModuleVarsstat varsstat;
cdo_finish();
varsstat.init(process);
varsstat.run();
varsstat.close();
return nullptr;
}
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

Imprint | Privacy Policy