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

Added operator dTdivdt (experimental).

parent c56ebe39
......@@ -25,6 +25,8 @@
#include "cdo_int.h"
#include "pstream_int.h"
#include "datetime.h"
void *
Deltat(void *process)
......@@ -34,27 +36,39 @@ Deltat(void *process)
cdoInitialize(process);
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
// clang-format off
cdoOperatorAdd("deltat", 0, 0, nullptr);
cdoOperatorAdd("dTdivdt", 0, 1, nullptr);
// clang-format on
int operatorID = cdoOperatorID();
bool ldivdt = cdoOperatorF2(operatorID);
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = vlistDuplicate(vlistID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
const int calendar = taxisInqCalendar(taxisID1);
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1));
pstreamDefVlist(streamID2, vlistID2);
FieldVector2D vars;
fieldsFromVlist(vlistID1, vars, FIELD_PTR);
size_t gridsizemax = vlistGridsizeMax(vlistID1);
std::vector<double> array1(gridsizemax);
std::vector<double> array2(gridsizemax);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
std::vector<double> array1(gridsizemax), array2(gridsizemax);
int tsID = 0;
int nrecs = cdoStreamInqTimestep(streamID1, tsID);
juldate_t juldate0 = juldate_encode(calendar, taxisInqVdate(taxisID1), taxisInqVtime(taxisID1));
for (int recID = 0; recID < nrecs; ++recID)
{
pstreamInqRecord(streamID1, &varID, &levelID);
......@@ -66,6 +80,10 @@ Deltat(void *process)
int tsID2 = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
juldate_t juldate1 = juldate_encode(calendar, taxisInqVdate(taxisID1), taxisInqVtime(taxisID1));
double idt_in_sec = ldivdt ? 1./juldate_to_seconds(juldate_sub(juldate1, juldate0)) : 1;
juldate0 = juldate1;
taxisCopyTimestep(taxisID2, taxisID1);
pstreamDefTimestep(streamID2, tsID2);
......@@ -74,9 +92,9 @@ Deltat(void *process)
pstreamInqRecord(streamID1, &varID, &levelID);
pstreamReadRecord(streamID1, array1.data(), &nmiss);
double missval = vars[varID][levelID].missval;
const double missval = vars[varID][levelID].missval;
const size_t gridsize = vars[varID][levelID].size;
double *array0 = vars[varID][levelID].ptr;
size_t gridsize = vars[varID][levelID].size;
if (nmiss || vars[varID][levelID].nmiss)
{
for (size_t i = 0; i < gridsize; ++i)
......@@ -84,14 +102,14 @@ Deltat(void *process)
if (DBL_IS_EQUAL(array0[i], missval) || DBL_IS_EQUAL(array1[i], missval))
array2[i] = missval;
else
array2[i] = array1[i] - array0[i];
array2[i] = (array1[i] - array0[i]) * idt_in_sec;
}
nmiss = arrayNumMV(gridsize, array2.data(), missval);
}
else
{
for (size_t i = 0; i < gridsize; ++i) array2[i] = array1[i] - array0[i];
for (size_t i = 0; i < gridsize; ++i) array2[i] = (array1[i] - array0[i]) * idt_in_sec;
}
arrayCopy(gridsize, array1.data(), array0);
......
......@@ -124,10 +124,10 @@ Sorttimestamp(void *process)
int nts = xtsID;
timeinfo_t *timeinfo = (timeinfo_t *) Malloc(nts * sizeof(timeinfo_t));
const int calendar = taxisInqCalendar(taxisID2);
for (tsID = 0; tsID < nts; tsID++)
{
const int calendar = taxisInqCalendar(taxisID2);
const int64_t julday = date_to_julday(calendar, vdate[tsID]);
const int secofday = time_to_sec(vtime[tsID]);
const double vdatetime = julday + secofday / 86400.;
......
......@@ -24,11 +24,10 @@ juldate_t
juldate_encode(int calendar, int64_t date, int time)
{
int year, month, day, hour, minute, second;
juldate_t juldate;
cdiDecodeDate(date, &year, &month, &day);
cdiDecodeTime(time, &hour, &minute, &second);
juldate_t juldate;
encode_caldaysec(calendar, year, month, day, hour, minute, second, &juldate.julday, &juldate.secofday);
return juldate;
......@@ -38,7 +37,6 @@ void
juldate_decode(int calendar, juldate_t juldate, int64_t *date, int *time)
{
int year, month, day, hour, minute, second;
decode_caldaysec(calendar, juldate.julday, juldate.secofday, &year, &month, &day, &hour, &minute, &second);
*date = cdiEncodeDate(year, month, day);
......@@ -69,5 +67,4 @@ double
juldate_to_seconds(juldate_t juldate)
{
return juldate.julday * 86400. + juldate.secofday;
;
}
......@@ -44,7 +44,7 @@
#define CondcOperators {"ifthenc", "ifnotthenc"}
#define ConsecstatOperators {"consects", "consecsum"}
#define CopyOperators {"copy", "selall", "szip"}
#define DeltatOperators {"deltat"}
#define DeltatOperators {"deltat", "dTdivdt"}
#define DeltimeOperators {"delday", "del29feb"}
#define DeriveparOperators {"gheight", "sealevelpressure"}
#define DetrendOperators {"detrend"}
......
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