Commit 2ed91023 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Prepare modules to add operator range.

parent ae9bb879
......@@ -37,6 +37,11 @@
#define MAX_DOY 373
typedef struct {
short varID;
short levelID;
} recinfo_t;
void *Ydaystat(void *argument)
{
......@@ -45,7 +50,6 @@ void *Ydaystat(void *argument)
int nrecs;
int nsets[MAX_DOY];
int nmiss;
int nlevel;
int vdates[MAX_DOY], vtimes[MAX_DOY];
field_type **vars1[MAX_DOY], **vars2[MAX_DOY], **samp1[MAX_DOY];
......@@ -91,17 +95,15 @@ void *Ydaystat(void *argument)
streamDefVlist(streamID2, vlistID2);
int nvars = vlistNvars(vlistID1);
int nrecords = vlistNrecs(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
int *recVarID = (int*) Malloc(nrecords*sizeof(int));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
int gridsize = vlistGridsizeMax(vlistID1);
int gridsizemax = vlistGridsizeMax(vlistID1);
field_type field;
field_init(&field);
field.ptr = (double*) Malloc(gridsize*sizeof(double));
field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
int tsID = 0;
int otsID = 0;
......@@ -137,16 +139,19 @@ void *Ydaystat(void *argument)
if ( tsID == 0 )
{
recVarID[recID] = varID;
recLevelID[recID] = levelID;
recinfo[recID].varID = varID;
recinfo[recID].levelID = levelID;
}
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
field_type *pvars1 = &vars1[dayoy][varID][levelID];
field_type *pvars2 = vars2[dayoy] ? &vars2[dayoy][varID][levelID] : NULL;
int gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
if ( nsets[dayoy] == 0 )
{
streamReadRecord(streamID1, vars1[dayoy][varID][levelID].ptr, &nmiss);
vars1[dayoy][varID][levelID].nmiss = (size_t)nmiss;
streamReadRecord(streamID1, pvars1->ptr, &nmiss);
pvars1->nmiss = (size_t)nmiss;
if ( nmiss > 0 || samp1[dayoy][varID][levelID].ptr )
{
......@@ -154,19 +159,15 @@ void *Ydaystat(void *argument)
samp1[dayoy][varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(vars1[dayoy][varID][levelID].ptr[i],
vars1[dayoy][varID][levelID].missval) )
samp1[dayoy][varID][levelID].ptr[i] = 0;
else
samp1[dayoy][varID][levelID].ptr[i] = 1;
samp1[dayoy][varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
}
}
else
{
streamReadRecord(streamID1, field.ptr, &nmiss);
field.nmiss = (size_t)nmiss;
field.grid = vars1[dayoy][varID][levelID].grid;
field.missval = vars1[dayoy][varID][levelID].missval;
field.grid = pvars1->grid;
field.missval = pvars1->missval;
if ( field.nmiss > 0 || samp1[dayoy][varID][levelID].ptr )
{
......@@ -178,29 +179,33 @@ void *Ydaystat(void *argument)
}
for ( int i = 0; i < gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], vars1[dayoy][varID][levelID].missval) )
if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
samp1[dayoy][varID][levelID].ptr[i]++;
}
if ( lvarstd )
{
farsumq(&vars2[dayoy][varID][levelID], field);
farsum(&vars1[dayoy][varID][levelID], field);
farsumq(pvars2, field);
farsum(pvars1, field);
}
else
{
farfun(&vars1[dayoy][varID][levelID], field, operfunc);
farfun(pvars1, field, operfunc);
}
}
}
if ( nsets[dayoy] == 0 && lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[dayoy][varID][levelID];
field_type *pvars2 = &vars2[dayoy][varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
farmoq(&vars2[dayoy][varID][levelID], vars1[dayoy][varID][levelID]);
farmoq(pvars2, *pvars1);
}
nsets[dayoy]++;
......@@ -227,56 +232,55 @@ void *Ydaystat(void *argument)
if ( nsets[dayoy] )
{
if ( lmean )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[dayoy][varID][levelID].ptr == NULL )
farcdiv(&vars1[dayoy][varID][levelID], (double)nsets[dayoy]);
else
fardiv(&vars1[dayoy][varID][levelID], samp1[dayoy][varID][levelID]);
}
}
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[dayoy][varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[dayoy][varID][levelID].ptr == NULL )
farcdiv(pvars1, (double)nsets[dayoy]);
else
fardiv(pvars1, samp1[dayoy][varID][levelID]);
}
else if ( lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( samp1[dayoy][varID][levelID].ptr == NULL )
{
if ( lstd )
farcstd(&vars1[dayoy][varID][levelID], vars2[dayoy][varID][levelID], nsets[dayoy], divisor);
else
farcvar(&vars1[dayoy][varID][levelID], vars2[dayoy][varID][levelID], nsets[dayoy], divisor);
}
else
{
if ( lstd )
farstd(&vars1[dayoy][varID][levelID], vars2[dayoy][varID][levelID], samp1[dayoy][varID][levelID], divisor);
else
farvar(&vars1[dayoy][varID][levelID], vars2[dayoy][varID][levelID], samp1[dayoy][varID][levelID], divisor);
}
}
}
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[dayoy][varID][levelID];
field_type *pvars2 = &vars2[dayoy][varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
if ( samp1[dayoy][varID][levelID].ptr == NULL )
{
if ( lstd ) farcstd(pvars1, *pvars2, nsets[dayoy], divisor);
else farcvar(pvars1, *pvars2, nsets[dayoy], divisor);
}
else
{
if ( lstd ) farstd(pvars1, *pvars2, samp1[dayoy][varID][levelID], divisor);
else farvar(pvars1, *pvars2, samp1[dayoy][varID][levelID], divisor);
}
}
taxisDefVdate(taxisID2, vdates[dayoy]);
taxisDefVtime(taxisID2, vtimes[dayoy]);
streamDefTimestep(streamID2, otsID);
for ( int recID = 0; recID < nrecords; recID++ )
{
varID = recVarID[recID];
levelID = recLevelID[recID];
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[dayoy][varID][levelID];
if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, vars1[dayoy][varID][levelID].ptr,
(int)vars1[dayoy][varID][levelID].nmiss);
streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
}
otsID++;
......@@ -294,8 +298,7 @@ void *Ydaystat(void *argument)
if ( field.ptr ) Free(field.ptr);
if ( recVarID ) Free(recVarID);
if ( recLevelID ) Free(recLevelID);
Free(recinfo);
streamClose(streamID2);
streamClose(streamID1);
......
......@@ -38,6 +38,11 @@
#define NDAY 373
typedef struct {
short varID;
short levelID;
} recinfo_t;
typedef struct {
int vdate[NDAY];
......@@ -65,8 +70,6 @@ void *Ydrunstat(void *argument)
int tsID;
int inp, its;
int nmiss;
int vdate, vtime;
int dayoy;
cdoInitialize(argument);
......@@ -105,10 +108,9 @@ void *Ydrunstat(void *argument)
streamDefVlist(streamID2, vlistID2);
int nrecords = vlistNrecs(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
int *recVarID = (int*) Malloc(nrecords*sizeof(int));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
cdo_datetime_t *datetime = (cdo_datetime_t*) Malloc((ndates+1)*sizeof(cdo_datetime_t));
......@@ -140,27 +142,30 @@ void *Ydrunstat(void *argument)
if ( tsID == 0 )
{
recVarID[recID] = varID;
recLevelID[recID] = levelID;
recinfo[recID].varID = varID;
recinfo[recID].levelID = levelID;
}
streamReadRecord(streamID1, vars1[tsID][varID][levelID].ptr, &nmiss);
vars1[tsID][varID][levelID].nmiss = nmiss;
field_type *pvars1 = &vars1[tsID][varID][levelID];
field_type *pvars2 = vars2[tsID] ? &vars2[tsID][varID][levelID] : NULL;
streamReadRecord(streamID1, pvars1->ptr, &nmiss);
pvars1->nmiss = nmiss;
if ( lvarstd )
{
farmoq(&vars2[tsID][varID][levelID], vars1[tsID][varID][levelID]);
farmoq(pvars2, *pvars1);
for ( inp = 0; inp < tsID; inp++ )
{
farsumq(&vars2[inp][varID][levelID], vars1[tsID][varID][levelID]);
farsum(&vars1[inp][varID][levelID], vars1[tsID][varID][levelID]);
farsumq(&vars2[inp][varID][levelID], *pvars1);
farsum(&vars1[inp][varID][levelID], *pvars1);
}
}
else
{
for ( inp = 0; inp < tsID; inp++ )
{
farfun(&vars1[inp][varID][levelID], vars1[tsID][varID][levelID], operfunc);
farfun(&vars1[inp][varID][levelID], *pvars1, operfunc);
}
}
}
......@@ -170,8 +175,8 @@ void *Ydrunstat(void *argument)
{
datetime_avg(dpy, ndates, datetime);
vdate = datetime[ndates].date;
vtime = datetime[ndates].time;
int vdate = datetime[ndates].date;
int vtime = datetime[ndates].time;
if ( lvarstd )
ydstatUpdate(stats, vdate, vtime, vars1[0], vars2[0], ndates, operfunc);
......@@ -201,23 +206,26 @@ void *Ydrunstat(void *argument)
{
streamInqRecord(streamID1, &varID, &levelID);
streamReadRecord(streamID1, vars1[ndates-1][varID][levelID].ptr, &nmiss);
vars1[ndates-1][varID][levelID].nmiss = nmiss;
field_type *pvars1 = &vars1[ndates-1][varID][levelID];
field_type *pvars2 = vars2[ndates-1] ? &vars2[ndates-1][varID][levelID] : NULL;
streamReadRecord(streamID1, pvars1->ptr, &nmiss);
pvars1->nmiss = nmiss;
if ( lvarstd )
{
for ( inp = 0; inp < ndates-1; inp++ )
{
farsumq(&vars2[inp][varID][levelID], vars1[ndates-1][varID][levelID]);
farsum(&vars1[inp][varID][levelID], vars1[ndates-1][varID][levelID]);
farsumq(&vars2[inp][varID][levelID], *pvars1);
farsum(&vars1[inp][varID][levelID], *pvars1);
}
farmoq(&vars2[ndates-1][varID][levelID], vars1[ndates-1][varID][levelID]);
farmoq(pvars2, *pvars1);
}
else
{
for ( inp = 0; inp < ndates-1; inp++ )
{
farfun(&vars1[inp][varID][levelID], vars1[ndates-1][varID][levelID], operfunc);
farfun(&vars1[inp][varID][levelID], *pvars1, operfunc);
}
}
}
......@@ -229,13 +237,13 @@ void *Ydrunstat(void *argument)
// set the year to the minimum of years found on output timestep
int outyear = 1e9;
int year, month, day;
for ( dayoy = 0; dayoy < NDAY; dayoy++ )
for ( int dayoy = 0; dayoy < NDAY; dayoy++ )
if ( stats->nsets[dayoy] )
{
cdiDecodeDate(stats->vdate[dayoy], &year, &month, &day);
if ( year < outyear ) outyear = year;
}
for ( dayoy = 0; dayoy < NDAY; dayoy++ )
for ( int dayoy = 0; dayoy < NDAY; dayoy++ )
if ( stats->nsets[dayoy] )
{
cdiDecodeDate(stats->vdate[dayoy], &year, &month, &day);
......@@ -247,23 +255,23 @@ void *Ydrunstat(void *argument)
int otsID = 0;
for ( dayoy = 0; dayoy < NDAY; dayoy++ )
for ( int dayoy = 0; dayoy < NDAY; dayoy++ )
if ( stats->nsets[dayoy] )
{
taxisDefVdate(taxisID2, stats->vdate[dayoy]);
taxisDefVtime(taxisID2, stats->vtime[dayoy]);
streamDefTimestep(streamID2, otsID);
for ( int recID = 0; recID < nrecords; recID++ )
{
varID = recVarID[recID];
levelID = recLevelID[recID];
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &stats->vars1[dayoy][varID][levelID];
if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, stats->vars1[dayoy][varID][levelID].ptr,
stats->vars1[dayoy][varID][levelID].nmiss);
streamWriteRecord(streamID2, pvars1->ptr, pvars1->nmiss);
}
otsID++;
......@@ -281,8 +289,7 @@ void *Ydrunstat(void *argument)
if ( datetime ) Free(datetime);
if ( recVarID ) Free(recVarID);
if ( recLevelID ) Free(recLevelID);
Free(recinfo);
streamClose(streamID2);
streamClose(streamID1);
......@@ -294,12 +301,10 @@ void *Ydrunstat(void *argument)
static
YDAY_STATS *ydstatCreate(int vlistID)
{
int dayoy;
{
YDAY_STATS *stats = (YDAY_STATS*) Malloc(sizeof(YDAY_STATS));
for ( dayoy = 0; dayoy < NDAY; dayoy++ )
for ( int dayoy = 0; dayoy < NDAY; dayoy++ )
{
stats->vdate[dayoy] = 0;
stats->vtime[dayoy] = 0;
......@@ -315,13 +320,13 @@ YDAY_STATS *ydstatCreate(int vlistID)
static
void ydstatDestroy(YDAY_STATS *stats)
{
int dayoy, varID, levelID, nvars, nlevels;
int varID, levelID, nlevels;
if ( stats != NULL )
{
nvars = vlistNvars(stats->vlist);
int nvars = vlistNvars(stats->vlist);
for ( dayoy = 0; dayoy < NDAY; dayoy++ )
for ( int dayoy = 0; dayoy < NDAY; dayoy++ )
{
if ( stats->vars1[dayoy] != NULL )
{
......
......@@ -37,6 +37,12 @@
#define MAX_HOUR 9301 /* 31*12*25 + 1 */
typedef struct {
short varID;
short levelID;
} recinfo_t;
static
int hour_of_year(int vdate, int vtime)
{
......@@ -65,15 +71,11 @@ int hour_of_year(int vdate, int vtime)
void *Yhourstat(void *argument)
{
int i;
int varID;
int vdate, vtime;
int houroy;
int nrecs;
int levelID;
int nsets[MAX_HOUR];
int nmiss;
int nlevel;
int vdates[MAX_HOUR], vtimes[MAX_HOUR];
field_type **vars1[MAX_HOUR], **vars2[MAX_HOUR], **samp1[MAX_HOUR];
......@@ -97,7 +99,7 @@ void *Yhourstat(void *argument)
bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
int divisor = operfunc == func_std1 || operfunc == func_var1;
for ( houroy = 0; houroy < MAX_HOUR; ++houroy )
for ( int houroy = 0; houroy < MAX_HOUR; ++houroy )
{
vars1[houroy] = NULL;
vars2[houroy] = NULL;
......@@ -119,28 +121,26 @@ void *Yhourstat(void *argument)
streamDefVlist(streamID2, vlistID2);
int nvars = vlistNvars(vlistID1);
int nrecords = vlistNrecs(vlistID1);
int maxrecs = vlistNrecs(vlistID1);
int *recVarID = (int*) Malloc(nrecords*sizeof(int));
int *recLevelID = (int*) Malloc(nrecords*sizeof(int));
recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
int gridsize = vlistGridsizeMax(vlistID1);
int gridsizemax = vlistGridsizeMax(vlistID1);
field_type field;
field_init(&field);
field.ptr = (double*) Malloc(gridsize*sizeof(double));
field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
int tsID = 0;
int otsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
vdate = taxisInqVdate(taxisID1);
vtime = taxisInqVtime(taxisID1);
int vdate = taxisInqVdate(taxisID1);
int vtime = taxisInqVtime(taxisID1);
if ( cdoVerbose ) cdoPrint("process timestep: %d %d %d", tsID+1, vdate, vtime);
houroy = hour_of_year(vdate, vtime);
int houroy = hour_of_year(vdate, vtime);
vdates[houroy] = vdate;
vtimes[houroy] = vtime;
......@@ -159,136 +159,138 @@ void *Yhourstat(void *argument)
if ( tsID == 0 )
{
recVarID[recID] = varID;
recLevelID[recID] = levelID;
recinfo[recID].varID = varID;
recinfo[recID].levelID = levelID;
}
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
field_type *pvars1 = &vars1[houroy][varID][levelID];
field_type *pvars2 = vars2[houroy] ? &vars2[houroy][varID][levelID] : NULL;
int gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
if ( nsets[houroy] == 0 )
{
streamReadRecord(streamID1, vars1[houroy][varID][levelID].ptr, &nmiss);
vars1[houroy][varID][levelID].nmiss = (size_t) nmiss;
streamReadRecord(streamID1, pvars1->ptr, &nmiss);
pvars1->nmiss = (size_t) nmiss;
if ( nmiss > 0 || samp1[houroy][varID][levelID].ptr )
{
if ( samp1[houroy][varID][levelID].ptr == NULL )
samp1[houroy][varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( i = 0; i < gridsize; i++ )
if ( DBL_IS_EQUAL(vars1[houroy][varID][levelID].ptr[i],
vars1[houroy][varID][levelID].missval) )
samp1[houroy][varID][levelID].ptr[i] = 0;
else
samp1[houroy][varID][levelID].ptr[i] = 1;
for ( int i = 0; i < gridsize; i++ )
samp1[houroy][varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
}
}
else
{
streamReadRecord(streamID1, field.ptr, &nmiss);
field.nmiss = (size_t) nmiss;
field.grid = vars1[houroy][varID][levelID].grid;
field.missval = vars1[houroy][varID][levelID].missval;
field.grid = pvars1->grid;
field.missval = pvars1->missval;
if ( field.nmiss > 0 || samp1[houroy][varID][levelID].ptr )
{
if ( samp1[houroy][varID][levelID].ptr == NULL )
{
samp1[houroy][varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
for ( i = 0; i < gridsize; i++ )
for ( int i = 0; i < gridsize; i++ )
samp1[houroy][varID][levelID].ptr[i] = nsets[houroy];
}
for ( i = 0; i < gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], vars1[houroy][varID][levelID].missval) )
for ( int i = 0; i < gridsize; i++ )
if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
samp1[houroy][varID][levelID].ptr[i]++;
}
if ( lvarstd )
{
farsumq(&vars2[houroy][varID][levelID], field);
farsum(&vars1[houroy][varID][levelID], field);
farsumq(pvars2, field);
farsum(pvars1, field);
}
else
{
farfun(&vars1[houroy][varID][levelID], field, operfunc);
farfun(pvars1, field, operfunc);
}
}
}
if ( nsets[houroy] == 0 && lvarstd )
for ( varID = 0; varID < nvars; varID++ )
{
for ( int recID = 0; recID < maxrecs; recID++ )
{
int varID = recinfo[recID].varID;
int levelID = recinfo[recID].levelID;
field_type *pvars1 = &vars1[houroy][varID][levelID];
field_type *pvars2 = &vars2[houroy][varID][levelID];
if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
farmoq(&vars2[houroy][varID][levelID], vars1[houroy][varID][levelID]);