Skip to content
Snippets Groups Projects
Commit 1bcdc7f9 authored by Oliver Heidmann's avatar Oliver Heidmann
Browse files

class rework Fldstat2.cc

parent 6df167d4
No related branches found
No related tags found
1 merge request!78M300433/class reworks/multiple class reworks 8
......@@ -120,115 +120,142 @@ covariance(const Field &field1, const Field &field2, const Varray<double> &weigh
return covariance(field1.vec_d, field2.vec_d, weight, field1.missval, field2.missval, field1.size);
}
void *
Fldstat2(void *process)
class ModuleFldstat2
{
auto wstatus = false;
auto needWeights = true;
cdo_initialize(process);
// clang-format off
cdo_operator_add("fldcor", FieldFunc_Cor, 0, nullptr);
cdo_operator_add("fldcovar", FieldFunc_Covar, 0, nullptr);
// clang-format on
int operfunc;
const auto operatorID = cdo_operator_id();
const auto operfunc = cdo_operator_f1(operatorID);
CdoStreamID streamID1;
CdoStreamID streamID2;
CdoStreamID streamID3;
const auto streamID1 = cdo_open_read(0);
const auto streamID2 = cdo_open_read(1);
int taxisID1;
int taxisID3;
const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
const auto vlistID3 = vlistDuplicate(vlistID1);
vlist_compare(vlistID1, vlistID2, CMP_ALL);
const auto taxisID1 = vlistInqTaxis(vlistID1);
const auto taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID3);
bool wstatus = false;
bool needWeights = true;
VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
double slon = 0.0, slat = 0.0;
const auto gridID3 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID3, 1);
gridDefYsize(gridID3, 1);
gridDefXvals(gridID3, &slon);
gridDefYvals(gridID3, &slat);
const auto ngrids = vlistNgrids(vlistID1);
for (int index = 0; index < ngrids; ++index) vlistChangeGridIndex(vlistID3, index, gridID3);
Field field1, field2;
const auto streamID3 = cdo_open_write(2);
cdo_def_vlist(streamID3, vlistID3);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> weight;
if (needWeights) weight.resize(gridsizemax);
int lastgridID = -1;
int tsID = 0;
while (true)
{
const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
if (nrecs == 0) break;
const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
if (nrecs2 == 0)
{
cdo_warning("Input streams have different number of time steps!");
break;
}
cdo_taxis_copy_timestep(taxisID3, taxisID1);
cdo_def_timestep(streamID3, tsID);
for (int recID = 0; recID < nrecs; ++recID)
{
int varID, levelID;
cdo_inq_record(streamID1, &varID, &levelID);
field1.init(varList1[varID]);
cdo_inq_record(streamID2, &varID, &levelID);
field2.init(varList2[varID]);
cdo_read_record(streamID1, field1);
cdo_read_record(streamID2, field2);
const auto gridID = varList1[varID].gridID;
if (needWeights && gridID != lastgridID)
{
lastgridID = gridID;
wstatus = (gridcell_weights(gridID, weight) != 0);
}
if (wstatus && tsID == 0 && levelID == 0)
cdo_warning("Using constant grid cell area weights for variable %s!", varList1[varID].name);
double sglval = 0.0;
if (operfunc == FieldFunc_Cor)
sglval = correlation(field1, field2, weight);
else if (operfunc == FieldFunc_Covar)
sglval = covariance(field1, field2, weight);
const auto nmiss3 = DBL_IS_EQUAL(sglval, varList1[varID].missval) ? 1 : 0;
cdo_def_record(streamID3, varID, levelID);
cdo_write_record(streamID3, &sglval, nmiss3);
}
tsID++;
}
cdo_stream_close(streamID3);
cdo_stream_close(streamID2);
cdo_stream_close(streamID1);
cdo_finish();
public:
void
init(void *process)
{
cdo_initialize(process);
// clang-format off
cdo_operator_add("fldcor", FieldFunc_Cor, 0, nullptr);
cdo_operator_add("fldcovar", FieldFunc_Covar, 0, nullptr);
// clang-format on
const auto operatorID = cdo_operator_id();
operfunc = cdo_operator_f1(operatorID);
streamID1 = cdo_open_read(0);
streamID2 = cdo_open_read(1);
const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
const auto vlistID3 = vlistDuplicate(vlistID1);
vlist_compare(vlistID1, vlistID2, CMP_ALL);
taxisID1 = vlistInqTaxis(vlistID1);
taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID3);
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
double slon = 0.0, slat = 0.0;
const auto gridID3 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID3, 1);
gridDefYsize(gridID3, 1);
gridDefXvals(gridID3, &slon);
gridDefYvals(gridID3, &slat);
const auto ngrids = vlistNgrids(vlistID1);
for (int index = 0; index < ngrids; ++index) vlistChangeGridIndex(vlistID3, index, gridID3);
streamID3 = cdo_open_write(2);
cdo_def_vlist(streamID3, vlistID3);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
if (needWeights) weight.resize(gridsizemax);
}
void
run()
{
int lastgridID = -1;
int tsID = 0;
while (true)
{
const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
if (nrecs == 0) break;
const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
if (nrecs2 == 0)
{
cdo_warning("Input streams have different number of time steps!");
break;
}
cdo_taxis_copy_timestep(taxisID3, taxisID1);
cdo_def_timestep(streamID3, tsID);
for (int recID = 0; recID < nrecs; ++recID)
{
int varID, levelID;
cdo_inq_record(streamID1, &varID, &levelID);
field1.init(varList1[varID]);
cdo_inq_record(streamID2, &varID, &levelID);
field2.init(varList2[varID]);
cdo_read_record(streamID1, field1);
cdo_read_record(streamID2, field2);
const auto gridID = varList1[varID].gridID;
if (needWeights && gridID != lastgridID)
{
lastgridID = gridID;
wstatus = (gridcell_weights(gridID, weight) != 0);
}
if (wstatus && tsID == 0 && levelID == 0)
cdo_warning("Using constant grid cell area weights for variable %s!", varList1[varID].name);
double sglval = 0.0;
if (operfunc == FieldFunc_Cor)
sglval = correlation(field1, field2, weight);
else if (operfunc == FieldFunc_Covar)
sglval = covariance(field1, field2, weight);
const auto nmiss3 = DBL_IS_EQUAL(sglval, varList1[varID].missval) ? 1 : 0;
cdo_def_record(streamID3, varID, levelID);
cdo_write_record(streamID3, &sglval, nmiss3);
}
tsID++;
}
}
void
close()
{
cdo_stream_close(streamID3);
cdo_stream_close(streamID2);
cdo_stream_close(streamID1);
cdo_finish();
}
};
void *
Fldstat2(void *process)
{
ModuleFldstat2 fldstat2;
fldstat2.init(process);
fldstat2.run();
fldstat2.close();
return nullptr;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment