/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2006 Brockmann Consult
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

#include <assert.h>
#include <string.h>

#include <cdi.h>

#include "functs.h"
#include "cdo_vlist.h"
#include "cdo_cdi_wrapper.h"
#include <mpim_grid.h>
#include "process_int.h"
#include "ecacore.h"
#include "ecautil.h"
#include "util_files.h"
#include "util_date.h"
#include "datetime.h"

#define FIRST_VAR_ID 0

#define IS_NOT_SET(x) (x == nullptr)
#define IS_SET(x) (x != nullptr)
#define VIS_SET(x) (!x.empty())

void
eca1(const ECA_REQUEST_1 *request)
{
  const auto operatorID = cdoOperatorID();

  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
  int64_t ivdate = 0, ovdate = 0, indate21 = 0;
  int ivtime = 0, ovtime = 0;
  int nrecs;
  int levelID;
  int itsID;
  int otsID;
  long nsets;

  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);

  const auto istreamID = cdoOpenRead(0);

  const auto ivlistID = cdoStreamInqVlist(istreamID);
  const auto ovlistID = vlistCreate();

  const auto gridID = vlistInqVarGrid(ivlistID, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID, FIRST_VAR_ID);
  const auto missval = vlistInqVarMissval(ivlistID, FIRST_VAR_ID);

  auto varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

  vlistDefVarMissval(ovlistID, varID, missval);

  if (IS_SET(request->var1.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var1.name);
  if (IS_SET(request->var1.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var1.longname);
  if (IS_SET(request->var1.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var1.units);

  if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
    {
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

      vlistDefVarMissval(ovlistID, varID, missval);

      if (IS_SET(request->var2.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var2.name);
      if (IS_SET(request->var2.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var2.longname);
      if (IS_SET(request->var2.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var2.units);
    }

  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);

  const auto itaxisID = vlistInqTaxis(ivlistID);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID));
  taxisDefRdate(otaxisID, request->var1.refdate);
  taxisDefRtime(otaxisID, 0);
  vlistDefTaxis(ovlistID, otaxisID);

  const auto ostreamID = cdoOpenWrite(1);
  cdoDefVlist(ostreamID, ovlistID);

  const auto gridsize = gridInqSize(gridID);

  Field field1, field2, field3;
  field1.resize(gridsize);
  field3.resize(gridsize);
  if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3)) field2.resize(gridsize);

  const auto nlevels = zaxisInqSize(zaxisID);

  FieldVector var12(nlevels), samp1(nlevels), samp2(nlevels);
  FieldVector var13, var21, var23;
  if (IS_SET(request->var1.f3)) var13.resize(nlevels);
  if (IS_SET(request->var2.h2)) var21.resize(nlevels);
  if (IS_SET(request->var2.h3)) var23.resize(nlevels);

  for (levelID = 0; levelID < nlevels; levelID++)
    {
      var12[levelID].grid = gridID;
      var12[levelID].missval = missval;
      var12[levelID].resize(gridsize);

      samp1[levelID].grid = gridID;
      samp1[levelID].missval = missval;
      samp1[levelID].resize(gridsize);

      samp2[levelID].grid = gridID;
      samp2[levelID].missval = missval;

      if (IS_SET(request->var1.f3))
        {
          var13[levelID].grid = gridID;
          var13[levelID].missval = missval;
          var13[levelID].resize(gridsize);
        }
      if (IS_SET(request->var2.h2))
        {
          var21[levelID].grid = gridID;
          var21[levelID].missval = missval;
          var21[levelID].resize(gridsize);
        }
      if (IS_SET(request->var2.h3))
        {
          var23[levelID].grid = gridID;
          var23[levelID].missval = missval;
          var23[levelID].resize(gridsize);
        }
    }

  itsID = 0;
  otsID = 0;
  while (true)
    {
      nsets = 0;
      while ((nrecs = cdoStreamInqTimestep(istreamID, itsID)) > 0)
        {
          ivdate = taxisInqVdate(itaxisID);
          ivtime = taxisInqVtime(itaxisID);

          if (nsets == 0)
            {
              SET_DATE(indate2, ivdate, ivtime);
              indate21 = ivdate;
            }
          SET_DATE(indate1, ivdate, ivtime);

          if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;

          for (int recID = 0; recID < nrecs; recID++)
            {
              cdoInqRecord(istreamID, &varID, &levelID);

              if (varID != FIRST_VAR_ID) continue;

              if (nsets == 0)
                {
                  if ( cmplen != 0 && itsID == 0 && request->var1.f2 == &vfarnum2 )
                    {
                      fieldFill(var12[levelID], missval);
                      var12[levelID].nmiss = gridsize;
                    }
                  else if ( request->var1.f2 != &vfarnum2 )
                    {
                      fieldFill(var12[levelID], missval);
                      var12[levelID].nmiss = gridsize;
                    }
                  fieldFill(samp1[levelID], missval);
                  if (!samp2[levelID].empty()) fieldFill(samp2[levelID], 0.0);
                  if (IS_SET(request->var1.f3)) fieldFill(var13[levelID], missval);
                  if (IS_SET(request->var2.h2)) fieldFill(var21[levelID], missval);
                  if (IS_SET(request->var2.h3)) fieldFill(var23[levelID], missval);

                  samp1[levelID].nmiss = gridsize;
                  if (IS_SET(request->var1.f3)) var13[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h2)) var21[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h3)) var23[levelID].nmiss = gridsize;
                }

              cdoReadRecord(istreamID, field1.vec_d.data(), &field1.nmiss);
              field1.grid = var12[levelID].grid;
              field1.missval = var12[levelID].missval;

              vfarnum(samp1[levelID], field1);

              if (IS_SET(request->var2.h2))
                {
                  field2.vec_d = field1.vec_d;
                  field2.nmiss = field1.nmiss;
                  field2.grid = field1.grid;
                  field2.missval = field1.missval;
                }

              if (IS_SET(request->var1.f1)) request->var1.f1(field1, request->var1.f1arg);

              if (field1.nmiss || !samp2[levelID].empty())
                {
                  if (samp2[levelID].empty())
                    {
                      samp2[levelID].resize(gridsize);
                      fieldFill(samp2[levelID], nsets);
                    }
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      if (DBL_IS_EQUAL(field1.vec_d[i], field1.missval)) continue;
                      samp2[levelID].vec_d[i]++;
                    }
                }

              if (IS_NOT_EQUAL(request->var1.mulc, 0.0)) vfarcmul(field1, request->var1.mulc);
              if (IS_NOT_EQUAL(request->var1.addc, 0.0)) vfarcadd(field1, request->var1.addc);

              if (IS_SET(request->var1.f3) && request->var1.f2 == &vfarnum2)
                {
                  varrayCopy(gridsize, var12[levelID].vec_d, field3.vec_d);
                  field3.nmiss = var12[levelID].nmiss;
                  field3.grid = var12[levelID].grid;
                  field3.missval = var12[levelID].missval;
                }
              request->var1.f2(var12[levelID], field1);

              if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
                {
                  // if h2 is null, use the output of f2 as input for h1
                  if (IS_NOT_SET(request->var2.h2))
                    {
                      varrayCopy(gridsize, var12[levelID].vec_d, field2.vec_d);
                      field2.nmiss = var12[levelID].nmiss;
                      field2.grid = var12[levelID].grid;
                      field2.missval = var12[levelID].missval;
                    }

                  if (IS_SET(request->var2.h1)) request->var2.h1(field2, request->var2.h1arg);

                  if (IS_NOT_SET(request->var2.h2))
                    request->var2.h3(var23[levelID], field2);
                  else
                    {
                      request->var2.h2(var21[levelID], field2);
                      if (IS_SET(request->var2.h3)) request->var2.h3(var23[levelID], var21[levelID]);
                    }
                }

              if (IS_SET(request->var1.f3))
                {
                  if (cmplen != 0 && request->var1.f2 == &vfarnum2)
                    {
                      auto &array1 = field3.vec_d;
                      const auto &array2 = var12[levelID].vec_d;
                      const auto missval2 = field2.missval;
                
                      const auto len = field1.size;
                      if (len != field2.size) cdoAbort("Fields have different size (%s)", __func__);
                
                      if (field2.nmiss)
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], missval2) && !DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      else
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      request->var1.f3(var13[levelID], field3);
                    }
                  else
                    request->var1.f3(var13[levelID], var12[levelID]);
                }
            }

          ovdate = ivdate;
          ovtime = ivtime;
          nsets++;
          itsID++;
        }

      if (nrecs == 0 && nsets == 0) break;

      if (request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME)
        for (levelID = 0; levelID < nlevels; levelID++)
          {
            auto &rvar = IS_SET(request->var1.f3) ? var13[levelID] : var12[levelID];

            if (samp2[levelID].empty())
              vfarcdiv(rvar, nsets);
            else
              vfardiv(rvar, samp2[levelID]);

            if (request->var1.epilog == PERCENT_OF_TIME) vfarcmul(rvar, 100.0);
          }
      if ( request->var1.refdate == 19550101 )
        {
          taxisDefVdate(otaxisID, ovdate);
          taxisDefVtime(otaxisID, ovtime);
          cdoDefTimestep(ostreamID, otsID);
        }
      else
        {
          int year, month, day;
          cdiDecodeDate(indate21, &year, &month, &day);
          defineMidOfTime(cdoOperatorF2(operatorID), otaxisID, year, month, 12);
          cdoDefTimestep(ostreamID, otsID);
        }

      if (otsID && vlistInqVarTimetype(ivlistID, FIRST_VAR_ID) == TIME_CONSTANT) continue;

      varID = 0;
      for (levelID = 0; levelID < nlevels; levelID++)
        {
          auto &rvar = IS_SET(request->var1.f3) ? var13[levelID] : var12[levelID];

          vfarsel(rvar, samp1[levelID]);

          cdoDefRecord(ostreamID, varID, levelID);
          cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
        }

      if (IS_SET(request->var2.h2) || IS_SET(request->var2.h3))
        {
          varID = 1;
          for (levelID = 0; levelID < nlevels; levelID++)
            {
              auto &rvar = IS_SET(request->var2.h3) ? var23[levelID] : var21[levelID];

              vfarsel(rvar, samp1[levelID]);

              cdoDefRecord(ostreamID, varID, levelID);
              cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
            }
        }

      if (nrecs == 0) break;
      otsID++;
    }

  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID);
}

void
eca2(const ECA_REQUEST_2 *request)
{
  const auto operatorID = cdoOperatorID();

  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
  int64_t ivdate = 0, ovdate = 0, indate21 = 0;
  int ivtime = 0, ovtime = 0;
  int nrecs;
  int varID, levelID;
  int itsID;
  int otsID;
  long nsets;

  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);

  const auto istreamID1 = cdoOpenRead(0);
  const auto istreamID2 = cdoOpenRead(1);

  const auto ivlistID1 = cdoStreamInqVlist(istreamID1);
  const auto ivlistID2 = cdoStreamInqVlist(istreamID2);
  const auto ovlistID = vlistCreate();

  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);

  const auto gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  const auto missval1 = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);
  const auto missval2 = vlistInqVarMissval(ivlistID2, FIRST_VAR_ID);

  varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

  vlistDefVarMissval(ovlistID, varID, missval1);

  if (IS_SET(request->var1.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var1.name);
  if (IS_SET(request->var1.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var1.longname);
  if (IS_SET(request->var1.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var1.units);

  if (IS_SET(request->var2.h2))
    {
      varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

      vlistDefVarMissval(ovlistID, varID, missval1);

      if (IS_SET(request->var2.name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->var2.name);
      if (IS_SET(request->var2.longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->var2.longname);
      if (IS_SET(request->var2.units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->var2.units);
    }

  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);

  const auto itaxisID1 = vlistInqTaxis(ivlistID1);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
  taxisDefRdate(otaxisID, request->var1.refdate);
  taxisDefRtime(otaxisID, 0);
  vlistDefTaxis(ovlistID, otaxisID);

  const auto ostreamID = cdoOpenWrite(2);
  cdoDefVlist(ostreamID, ovlistID);

  const auto gridsize = gridInqSize(gridID);

  Field field1, field2, field3;
  field1.resize(gridsize);
  field2.resize(gridsize);
  field3.resize(gridsize);
  constexpr int MaxDays = 373;
  FieldVector2D vars2[MaxDays];

  const auto nlevels = zaxisInqSize(zaxisID);

  FieldVector var14(nlevels), samp1(nlevels), samp2(nlevels), samp3(nlevels);
  FieldVector total, var15, var22;
  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) total.resize(nlevels);
  if (IS_SET(request->var1.f5)) var15.resize(nlevels);
  if (IS_SET(request->var2.h2)) var22.resize(nlevels);

  for (levelID = 0; levelID < nlevels; levelID++)
    {
      var14[levelID].grid = gridID;
      var14[levelID].missval = missval1;
      var14[levelID].resize(gridsize);

      samp1[levelID].grid = gridID;
      samp1[levelID].missval = missval1;
      samp1[levelID].resize(gridsize);

      samp2[levelID].grid = gridID;
      samp2[levelID].missval = missval1;
      samp2[levelID].resize(gridsize);

      samp3[levelID].grid = gridID;
      samp3[levelID].missval = missval1;

      if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT)
        {
          total[levelID].grid = gridID;
          total[levelID].missval = missval1;
          total[levelID].resize(gridsize);
        }
      if (IS_SET(request->var1.f5))
        {
          var15[levelID].grid = gridID;
          var15[levelID].missval = missval1;
          var15[levelID].resize(gridsize);
        }
      if (IS_SET(request->var2.h2))
        {
          var22[levelID].grid = gridID;
          var22[levelID].missval = missval1;
          var22[levelID].resize(gridsize);
        }
    }
  itsID = 0;
  while ((nrecs = cdoStreamInqTimestep(istreamID2, itsID)))
    {
      ivdate = taxisInqVdate(vlistInqTaxis(ivlistID2));

      const auto dayoy = decodeDayOfYear(ivdate);
      if (dayoy < 0 || dayoy >= MaxDays) cdoAbort("Day %d out of range!", dayoy);

      if (!vars2[dayoy].size())
        {
          fieldsFromVlist(ivlistID2, vars2[dayoy], FIELD_VEC);
        }

      for (int recID = 0; recID < nrecs; recID++)
        {
          size_t nmiss = 0;
          cdoInqRecord(istreamID2, &varID, &levelID);
          if (varID != FIRST_VAR_ID) continue;
          cdoReadRecord(istreamID2, vars2[dayoy][0][levelID].vec_d.data(), &nmiss);
          vars2[dayoy][0][levelID].nmiss = nmiss;
        }
      itsID++;
    }

  itsID = 0;
  otsID = 0;
  while (true)
    {
      nsets = 0;
      while ((nrecs = cdoStreamInqTimestep(istreamID1, itsID)) > 0)
        {
          ivdate = taxisInqVdate(itaxisID1);
          ivtime = taxisInqVtime(itaxisID1);

          const auto dayoy = decodeDayOfYear(ivdate);
          if (!vars2[dayoy].size())
            cdoAbort("Input streams have different time values!");

          if (nsets == 0) 
            {
              SET_DATE(indate2, ivdate, ivtime);
              indate21 = ivdate;
            }
          SET_DATE(indate1, ivdate, ivtime);

          if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;

          for (int recID = 0; recID < nrecs; recID++)
            {
              cdoInqRecord(istreamID1, &varID, &levelID);

              if (varID != FIRST_VAR_ID) continue;

              if (nsets == 0)
                {
                  if ( cmplen != 0 && itsID == 0 && request->var1.f4 == &vfarnum2 )
                    {
                      fieldFill(var14[levelID], missval1);
                      var14[levelID].nmiss = gridsize;
                    }
                  else if ( request->var1.f4 != &vfarnum2 )
                    {
                      fieldFill(var14[levelID], missval1);
                      var14[levelID].nmiss = gridsize;
                    }
                  fieldFill(samp1[levelID], missval1);
                  fieldFill(samp2[levelID], missval1);
                  if (!samp3[levelID].empty()) fieldFill(samp3[levelID], 0.0);
                  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) fieldFill(total[levelID], 0.0);
                  if (IS_SET(request->var1.f5)) fieldFill(var15[levelID], missval1);
                  if (IS_SET(request->var2.h2)) fieldFill(var22[levelID], missval1);

                  samp1[levelID].nmiss = gridsize;
                  samp2[levelID].nmiss = gridsize;
                  if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) total[levelID].nmiss = gridsize;
                  if (IS_SET(request->var1.f5)) var15[levelID].nmiss = gridsize;
                  if (IS_SET(request->var2.h2)) var22[levelID].nmiss = gridsize;
                }
              cdoReadRecord(istreamID1, field1.vec_d.data(), &field1.nmiss);
              field1.grid = gridID;
              field1.missval = missval1;
              field2.grid = vars2[dayoy][0][levelID].grid;
              field2.nmiss = vars2[dayoy][0][levelID].nmiss;
              field2.missval = missval2;
              vfarcpy(field2, vars2[dayoy][0][levelID]); 


              vfarnum(samp1[levelID], field1);
              vfarnum(samp2[levelID], field2);

              if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT) vfarsum(total[levelID], field1);

              if (IS_SET(request->var1.f1)) request->var1.f1(field1, request->var1.f1arg);
              if (IS_SET(request->var1.f2)) request->var1.f2(field2, request->var1.f2arg);

              if (field1.nmiss || !samp3[levelID].empty())
                {
                  if (samp3[levelID].empty())
                    {
                      samp3[levelID].resize(gridsize);
                      fieldFill(samp3[levelID], nsets);
                    }
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      if (DBL_IS_EQUAL(field1.vec_d[i], field1.missval)) continue;
                      samp3[levelID].vec_d[i]++;
                    }
                }

              if (IS_SET(request->var1.f5) && request->var1.f4 == &vfarnum2)
                {
                  varrayCopy(gridsize, var14[levelID].vec_d, field3.vec_d);
                  field3.nmiss = var14[levelID].nmiss;
                  field3.grid = var14[levelID].grid;
                  field3.missval = var14[levelID].missval;
                }

              request->var1.f3(field1, field2);
              request->var1.f4(var14[levelID], field1);

              if (IS_SET(request->var2.h2))
                {
                  varrayCopy(gridsize, var14[levelID].vec_d, field2.vec_d);
                  field2.nmiss = var14[levelID].nmiss;
                  field2.grid = var14[levelID].grid;
                  field2.missval = var14[levelID].missval;

                  if (IS_SET(request->var2.h1)) request->var2.h1(field2, request->var2.h1arg);

                  request->var2.h2(var22[levelID], field2);
                }

              if (IS_SET(request->var1.f5))
                {
                  if (cmplen != 0 && request->var1.f4 == &vfarnum2)
                    {
                      auto &array1 = field3.vec_d;
                      const auto &array2 = var14[levelID].vec_d;
                      const auto missvaltemp = field1.missval;
                
                      const auto len = field1.size;
                      if (len != field3.size) cdoAbort("Fields have different size (%s)", __func__);
                
                      if (field1.nmiss)
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], missvaltemp) && !DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      else
                        {
                          for (size_t i = 0; i < len; i++)
                            {
                              if (!DBL_IS_EQUAL(array2[i], 0.0))
                                array1[i] = 0.0; 
                            }
                        }
                      request->var1.f5(var15[levelID], field3, request->var1.f5arg);
                    }
                  else
                    request->var1.f5(var15[levelID], var14[levelID], request->var1.f5arg);
                }
            }

          ovdate = ivdate;
          ovtime = ivtime;
          nsets++;
          itsID++;
        }

      if (nrecs == 0 && nsets == 0) break;

      if (request->var1.epilog == MEAN || request->var1.epilog == PERCENT_OF_TIME)
        for (levelID = 0; levelID < nlevels; levelID++)
          {
            auto &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];

            if (samp3[levelID].empty())
              vfarcdiv(rvar, nsets);
            else
              vfardiv(rvar, samp3[levelID]);

            if (request->var1.epilog == PERCENT_OF_TIME) vfarcmul(rvar, 100.0);
          }
      else if (request->var1.epilog == PERCENT_OF_TOTAL_AMOUNT)
        for (levelID = 0; levelID < nlevels; levelID++)
          {
            Field &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];

            vfardiv(rvar, total[levelID]);
            vfarcmul(rvar, 100.0);
          }

      if ( request->var1.refdate == 19550101 )
        {
          taxisDefVdate(otaxisID, ovdate);
          taxisDefVtime(otaxisID, ovtime);
          cdoDefTimestep(ostreamID, otsID);
        }
      else
        {
          int year, month, day;
          cdiDecodeDate(indate21, &year, &month, &day);
          defineMidOfTime(cdoOperatorF2(operatorID), otaxisID, year, month, 12);
          cdoDefTimestep(ostreamID, otsID);
        }

      if (otsID && vlistInqVarTimetype(ivlistID1, FIRST_VAR_ID) == TIME_CONSTANT) continue;

      varID = 0;
      for (levelID = 0; levelID < nlevels; levelID++)
        {
          Field &rvar = IS_SET(request->var1.f5) ? var15[levelID] : var14[levelID];

          vfarsel(rvar, samp1[levelID]);
          vfarsel(rvar, samp2[levelID]);

          cdoDefRecord(ostreamID, varID, levelID);
          cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
        }
      if (IS_SET(request->var2.h2))
        {
          varID = 1;
          for (levelID = 0; levelID < nlevels; levelID++)
            {
              auto &rvar = var22[levelID];

              vfarsel(rvar, samp1[levelID]);
              vfarsel(rvar, samp2[levelID]);

              cdoDefRecord(ostreamID, varID, levelID);
              cdoWriteRecord(ostreamID, rvar.vec_d.data(), rvar.nmiss);
            }
        }

      if (nrecs == 0) break;
      otsID++;
    }

  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID2);
  cdoStreamClose(istreamID1);
}

void
eca3(const ECA_REQUEST_3 *request)
{
  const auto operatorID = cdoOperatorID();

  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
  int64_t ivdate1 = 0, ivdate2 = 0, indate21 = 0;
  int ivtime1 = 0, ivtime2 = 0;
  int64_t ovdate = 0;
  int ovtime = 0;
  int nrecs;
  int varID, levelID;
  int itsID;
  int otsID;
  long nsets;

  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);

  const auto istreamID1 = cdoOpenRead(0);
  const auto istreamID2 = cdoOpenRead(1);

  const auto ivlistID1 = cdoStreamInqVlist(istreamID1);
  const auto ivlistID2 = cdoStreamInqVlist(istreamID2);
  const auto ovlistID = vlistCreate();

  vlistCompare(ivlistID1, ivlistID2, CMP_ALL);

  const auto gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  const auto zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  const auto missval = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);

  varID = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

  vlistDefVarMissval(ovlistID, varID, missval);

  if (IS_SET(request->name)) cdiDefKeyString(ovlistID, varID, CDI_KEY_NAME, request->name);
  if (IS_SET(request->longname)) cdiDefKeyString(ovlistID, varID, CDI_KEY_LONGNAME, request->longname);
  if (IS_SET(request->units)) cdiDefKeyString(ovlistID, varID, CDI_KEY_UNITS, request->units);

  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);

  const auto itaxisID1 = vlistInqTaxis(ivlistID1);
  const auto itaxisID2 = vlistInqTaxis(ivlistID2);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID1));
  taxisDefRdate(otaxisID, request->refdate);
  taxisDefRtime(otaxisID, 0);
  vlistDefTaxis(ovlistID, otaxisID);

  const auto ostreamID = cdoOpenWrite(2);
  cdoDefVlist(ostreamID, ovlistID);

  const auto gridsize = gridInqSize(gridID);

  Field field1, field2;
  field1.resize(gridsize);
  field2.resize(gridsize);

  const int nlevels = zaxisInqSize(zaxisID);

  FieldVector var1(nlevels), var2(nlevels);

  for (levelID = 0; levelID < nlevels; levelID++)
    {
      var1[levelID].grid = gridID;
      var1[levelID].missval = missval;
      var1[levelID].resize(gridsize);

      var2[levelID].grid = gridID;
      var2[levelID].missval = missval;
      var2[levelID].resize(gridsize);
    }

  itsID = 0;
  otsID = 0;
  while (true)
    {
      nsets = 0;
      while ((nrecs = cdoStreamInqTimestep(istreamID1, itsID)) > 0)
        {
          if (!cdoStreamInqTimestep(istreamID2, itsID)) cdoAbort("Input streams have different number of time steps!");

          ivdate1 = taxisInqVdate(itaxisID1);
          ivdate2 = taxisInqVdate(itaxisID2);
          if (ivdate1 != ivdate2) cdoAbort("Input streams have different verification dates at time step %d!", itsID + 1);

          ivtime1 = taxisInqVtime(itaxisID1);
          ivtime2 = taxisInqVtime(itaxisID2);
          if (ivtime1 != ivtime2) cdoAbort("Input streams have different verification times at time step %d!", itsID + 1);

          if (nsets == 0)
            {
              SET_DATE(indate2, ivdate1, ivtime1);
              indate21 = ivdate1;
            }
          SET_DATE(indate1, ivdate1, ivtime1);

          if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;

          for (int recID = 0; recID < nrecs; recID++)
            {
              cdoInqRecord(istreamID1, &varID, &levelID);
              cdoInqRecord(istreamID2, &varID, &levelID);

              if (varID != FIRST_VAR_ID) continue;

              if (nsets == 0)
                {
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      var1[levelID].vec_d[i] = missval;
                      var2[levelID].vec_d[i] = missval;
                    }
                  var1[levelID].nmiss = gridsize;
                  var2[levelID].nmiss = gridsize;
                }

              cdoReadRecord(istreamID1, field1.vec_d.data(), &field1.nmiss);
              field1.grid = var1[levelID].grid;
              field1.missval = var1[levelID].missval;

              cdoReadRecord(istreamID2, field2.vec_d.data(), &field2.nmiss);
              field2.grid = var1[levelID].grid;
              field2.missval = var1[levelID].missval;

              request->f1(var1[levelID], field1);
              request->f2(var2[levelID], field2);
            }

          ovdate = ivdate1;
          ovtime = ivtime1;
          nsets++;
          itsID++;
        }

      if (nrecs == 0 && nsets == 0) break;

      for (levelID = 0; levelID < nlevels; levelID++) request->f3(var1[levelID], var2[levelID]);

      if ( request->refdate == 19550101 )
        {
          taxisDefVdate(otaxisID, ovdate);
          taxisDefVtime(otaxisID, ovtime);
          cdoDefTimestep(ostreamID, otsID);
        }
      else
        {
          int year, month, day;
          cdiDecodeDate(indate21, &year, &month, &day);
          defineMidOfTime(cdoOperatorF2(operatorID), otaxisID, year, month, 12);
          cdoDefTimestep(ostreamID, otsID);
        }

      if (otsID && vlistInqVarTimetype(ivlistID1, FIRST_VAR_ID) == TIME_CONSTANT) continue;

      varID = 0;
      for (levelID = 0; levelID < nlevels; levelID++)
        {
          cdoDefRecord(ostreamID, varID, levelID);
          cdoWriteRecord(ostreamID, var1[levelID].vec_d.data(), var1[levelID].nmiss);
        }

      if (nrecs == 0) break;
      otsID++;
    }

  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID2);
  cdoStreamClose(istreamID1);
}

// check for non missval values
static bool
fldhvs(FieldVector &fieldVector, const size_t nlevels)
{
  for (size_t level = 0; level < nlevels; level++)
    {
      if (fieldVector[level].nmiss != fieldVector[level].size) return true;
    }

  return false;
}

void
eca4(const ECA_REQUEST_4 *request)
{
  const auto operatorID = cdoOperatorID();

  int yearcnt = 0;
  int nrecs;
  int varID, levelID;
  bool resetAtJan = false, resetAtJul = false;
  bool isFirstYear = true;
  int64_t ivdate = 0, ovdate = 0;
  int ivtime = 0, ovtime = 0;
  char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];

  const auto cmplen = DATE_LEN - cdoOperatorF2(operatorID);

  const auto istreamID1 = cdoOpenRead(0);
  const auto istreamID2 = cdoOpenRead(1);

  const auto ivlistID1 = cdoStreamInqVlist(istreamID1);
  const auto ivlistID2 = cdoStreamInqVlist(istreamID2);
  const auto ovlistID = vlistCreate();

  int gridID = vlistInqVarGrid(ivlistID1, FIRST_VAR_ID);
  if (gridInqSize(gridID) != gridInqSize(vlistInqVarGrid(ivlistID2, FIRST_VAR_ID)))
    cdoAbort("Grid sizes of the input fields do not match!");

  const auto zaxisID = vlistInqVarZaxis(ivlistID1, FIRST_VAR_ID);
  const auto missval = vlistInqVarMissval(ivlistID1, FIRST_VAR_ID);

  const auto ovarID1 = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

  vlistDefVarMissval(ovlistID, ovarID1, missval);

  if (IS_SET(request->name)) cdiDefKeyString(ovlistID, ovarID1, CDI_KEY_NAME, request->name);
  if (IS_SET(request->longname)) cdiDefKeyString(ovlistID, ovarID1, CDI_KEY_LONGNAME, request->longname);
  if (IS_SET(request->units)) cdiDefKeyString(ovlistID, ovarID1, CDI_KEY_UNITS, request->units);

  const auto ovarID2 = vlistDefVar(ovlistID, gridID, zaxisID, TIME_VARYING);

  vlistDefVarMissval(ovlistID, ovarID2, missval);

  if (IS_SET(request->name2)) cdiDefKeyString(ovlistID, ovarID2, CDI_KEY_NAME, request->name2);
  if (IS_SET(request->longname2)) cdiDefKeyString(ovlistID, ovarID2, CDI_KEY_LONGNAME, request->longname2);
  if (IS_SET(request->units2)) cdiDefKeyString(ovlistID, ovarID2, CDI_KEY_UNITS, request->units2);

  if (cdoOperatorF2(operatorID) == 16) vlistDefNtsteps(ovlistID, 1);

  const auto itaxisID = vlistInqTaxis(ivlistID1);
  const auto otaxisID = cdoTaxisCreate(TAXIS_RELATIVE);
  taxisDefTunit(otaxisID, TUNIT_DAY);
  //  taxisDefTunit(otaxisID, TUNIT_MINUTE);
  //  taxisDefCalendar(otaxisID, CALENDAR_PROLEPTIC);
  taxisDefCalendar(otaxisID, taxisInqCalendar(itaxisID));
  taxisDefRdate(otaxisID, 19550101);
  taxisDefRtime(otaxisID, 0);
  vlistDefTaxis(ovlistID, otaxisID);

  const auto ostreamID = cdoOpenWrite(2);
  cdoDefVlist(ostreamID, ovlistID);

  bool lyvals = true;
  const auto gridtype = gridInqType(gridID);
  if (gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || gridtype == GRID_PROJECTION)
    {
      gridID = gridToCurvilinear(gridID, 1);
    }
  else if (gridtype == GRID_GME)
    {
      gridID = gridToUnstructured(gridID, 1);
    }
  else
    {
      lyvals = false;
    }

  const auto gridsize = gridInqSize(gridID);
  // for later check on northern\southern hemisphere
  std::vector<double> yvals(gridsize);
  if (lyvals)
    {
      gridInqYvals(gridID, yvals.data());
    }
  else
    {
      for (size_t i = 0; i < gridsize; ++i) yvals[i] = 20;  // Northern hemisphere
    }

  // Two fields are needed because of the definition of gsl for northern and southern hemisphere
  Field fieldGt, fieldLt;
  fieldGt.resize(gridsize);
  fieldLt.resize(gridsize);

  // field for the land-water-distribution
  Field mask;
  mask.resize(gridsize);

  const auto nlevels = zaxisInqSize(zaxisID);

  FieldVector startCount(nlevels), endCount(nlevels);
  FieldVector gslDuration(nlevels), gslFirstDay(nlevels);

  FieldVector2D startDateWithHist(2), endDateWithHist(2);
  /* because of the different definitions for northern and southern hemisphere,
   * the values of the last year have to be present THE LAST YEAR HAS THE INDEX 1 */
  for (int h = 0; h < 2; h++)
    {
      startDateWithHist[h].resize(nlevels);
      endDateWithHist[h].resize(nlevels);
    }

  for (levelID = 0; levelID < nlevels; levelID++)
    {
      startCount[levelID].grid = gridID;
      startCount[levelID].size = gridsize;
      startCount[levelID].missval = missval;
      startCount[levelID].resize(gridsize, 0.0);

      endCount[levelID].grid = gridID;
      endCount[levelID].size = gridsize;
      endCount[levelID].missval = missval;
      endCount[levelID].resize(gridsize, 0.0);

      gslDuration[levelID].grid = gridID;
      gslDuration[levelID].size = gridsize;
      gslDuration[levelID].missval = missval;
      gslDuration[levelID].resize(gridsize);

      gslFirstDay[levelID].grid = gridID;
      gslFirstDay[levelID].size = gridsize;
      gslFirstDay[levelID].missval = missval;
      gslFirstDay[levelID].resize(gridsize);

      for (int h = 0; h < 2; h++)
        {
          startDateWithHist[h][levelID].grid = gridID;
          startDateWithHist[h][levelID].size = gridsize;
          startDateWithHist[h][levelID].missval = missval;
          startDateWithHist[h][levelID].resize(gridsize);

          endDateWithHist[h][levelID].grid = gridID;
          endDateWithHist[h][levelID].size = gridsize;
          endDateWithHist[h][levelID].missval = missval;
          endDateWithHist[h][levelID].resize(gridsize);
        }
    }

  int itsID = 0;
  int otsID = 0;

  if (cdoStreamInqTimestep(istreamID2, itsID))
    {
      cdoInqRecord(istreamID2, &varID, &levelID);
      cdoReadRecord(istreamID2, mask.vec_d.data(), &mask.nmiss);
      mask.grid = gridID;
      mask.missval = vlistInqVarMissval(ivlistID2, 0);

      request->s3(mask, request->s3arg);
    }
  else
    cdoAbort("Could not read land-water mask!");

  while (true)
    {
      long nsets = 0;
      while ((nrecs = cdoStreamInqTimestep(istreamID1, itsID)) > 0)
        {
          ivdate = taxisInqVdate(itaxisID);
          ivtime = taxisInqVtime(itaxisID);

          int month = (ivdate % 10000) / 100;
          if (month < 1 || month > 12) cdoAbort("month %d out of range!", month);

          if (nsets == 0) SET_DATE(indate2, ivdate, ivtime);
          SET_DATE(indate1, ivdate, ivtime);

          if (DATE_IS_NEQ(indate1, indate2, cmplen))
            {
              resetAtJan = false;
              resetAtJul = false;
              break;
            }

          for (int recID = 0; recID < nrecs; recID++)
            {
              cdoInqRecord(istreamID1, &varID, &levelID);

              if (varID != FIRST_VAR_ID) continue;

              if (nsets == 0)
                {
                  fieldFill(gslDuration[levelID], missval);
                  fieldFill(gslFirstDay[levelID], missval);
                  // reinitialize the current year
                  fieldFill(startDateWithHist[0][levelID], missval);
                  fieldFill(endDateWithHist[0][levelID], missval);

                  gslDuration[levelID].nmiss = 0;
                  gslFirstDay[levelID].nmiss = 0;
                  // reinitialize the current year
                  startDateWithHist[0][levelID].nmiss = gridsize;
                  endDateWithHist[0][levelID].nmiss = gridsize;
                }
              // init the history ONCE
              if (0 == itsID)
                {
                  fieldFill(startDateWithHist[1][levelID], missval);
                  fieldFill(endDateWithHist[1][levelID], missval);

                  startDateWithHist[1][levelID].nmiss = gridsize;
                  endDateWithHist[1][levelID].nmiss = gridsize;
                }

              cdoReadRecord(istreamID1, fieldGt.vec_d.data(), &fieldGt.nmiss);
              fieldLt.vec_d = fieldGt.vec_d;
              fieldLt.nmiss = fieldGt.nmiss;
              fieldGt.grid = startCount[levelID].grid;
              fieldGt.missval = startCount[levelID].missval;
              fieldLt.grid = startCount[levelID].grid;
              fieldLt.missval = startCount[levelID].missval;

              // Reinitialization of (start|end)Count variables has to be done different for norther and southern hemisphere
              if (1 == month && !resetAtJan)
                {
                  // reset northern startCount
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      if (yvals[i] >= 0.0)
                        if (!DBL_IS_EQUAL(startCount[levelID].vec_d[i], missval))
                          {
                            startCount[levelID].vec_d[i] = missval;
                            startCount[levelID].nmiss++;
                          }
                    }
                  // reset southern endCount
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      if (yvals[i] < 0.0)
                        if (!DBL_IS_EQUAL(endCount[levelID].vec_d[i], missval))
                          {
                            endCount[levelID].vec_d[i] = missval;
                            endCount[levelID].nmiss++;
                          }
                    }

                  resetAtJan = true;
                }
              if (7 == month && !resetAtJul)
                {
#ifdef _OPENMP
#pragma omp sections
#endif
                  {
#ifdef _OPENMP
#pragma omp section
#endif
                    {
                      // reset northern endCount
                      for (size_t i = 0; i < gridsize; i++)
                        {
                          if (yvals[i] >= 0.0)
                            {
                              if (!DBL_IS_EQUAL(endCount[levelID].vec_d[i], missval))
                                {
                                  endCount[levelID].vec_d[i] = missval;
                                  endCount[levelID].nmiss++;
                                }
                            }
                        }
                    }
#ifdef _OPENMP
#pragma omp section
#endif
                    {
                      // reset southern startCount
                      for (size_t i = 0; i < gridsize; i++)
                        {
                          if (yvals[i] < 0.0)
                            {
                              if (!DBL_IS_EQUAL(startCount[levelID].vec_d[i], missval))
                                {
                                  startCount[levelID].vec_d[i] = missval;
                                  startCount[levelID].nmiss++;
                                }
                            }
                        }
                    }
                  }
                  resetAtJul = true;
                }

                // count the day with temperature larger/smaller than the given limit
#ifdef _OPENMP
#pragma omp sections
#endif
              {
#ifdef _OPENMP
#pragma omp section
#endif
                {
                  vfarsel(fieldGt, mask);
                  request->s1(fieldGt, request->s1arg);
                  vfarnum2(startCount[levelID], fieldGt);
                }
#ifdef _OPENMP
#pragma omp section
#endif
                {
                  vfarsel(fieldLt, mask);
                  request->s2(fieldLt, request->s1arg);
                  vfarnum2(endCount[levelID], fieldLt);
                }
              }

              if (month < 7)
                {
                  for (size_t i = 0; i < gridsize; i++)  // dictinct between northern and southern sphere
                    // start with south
                    if (yvals[i] < 0)
                      {
                        // south: periods can also start in the first half of the year, but this date has already gone into the
                        // history
                        if (DBL_IS_EQUAL(startDateWithHist[1][levelID].vec_d[i], missval)
                            && IS_EQUAL(startCount[levelID].vec_d[i], request->consecutiveDays))
                          {
                            startDateWithHist[1][levelID].vec_d[i] = ivdate;
                            // reset the endCount, because we are only interessted in the end of the eriod, if a start was found
                            endCount[levelID].vec_d[i] = missval;
                            endDateWithHist[0][levelID].vec_d[i] = missval;
                          }
                        if (DBL_IS_EQUAL(endDateWithHist[0][levelID].vec_d[i], missval)
                            && IS_EQUAL(endCount[levelID].vec_d[i], request->consecutiveDays))
                          {
                            endDateWithHist[0][levelID].vec_d[i] = ivdate;
                          }
                      }
                    else
                      {
                        if (DBL_IS_EQUAL(startDateWithHist[0][levelID].vec_d[i], missval)
                            && IS_EQUAL(startCount[levelID].vec_d[i], request->consecutiveDays))
                          {
                            startDateWithHist[0][levelID].vec_d[i] = ivdate;
                          }
                      }
                }
              else
                {
                  for (size_t i = 0; i < gridsize; i++)
                    {
                      if (yvals[i] < 0)
                        {
                          if (DBL_IS_EQUAL(startDateWithHist[0][levelID].vec_d[i], missval)
                              && IS_EQUAL(startCount[levelID].vec_d[i], request->consecutiveDays))
                            {
                              startDateWithHist[0][levelID].vec_d[i] = ivdate;
                            }
                        }
                      else
                        {
                          // north: periods can also start in the second half of the year
                          if (DBL_IS_EQUAL(startDateWithHist[0][levelID].vec_d[i], missval)
                              && IS_EQUAL(startCount[levelID].vec_d[i], request->consecutiveDays))
                            {
                              startDateWithHist[0][levelID].vec_d[i] = ivdate;
                              // reset the endCount, because we are only interessted in the end of the eriod, if a start was found
                              endCount[levelID].vec_d[i] = missval;
                              endDateWithHist[0][levelID].vec_d[i] = missval;
                            }
                          if (DBL_IS_EQUAL(endDateWithHist[0][levelID].vec_d[i], missval)
                              && IS_EQUAL(endCount[levelID].vec_d[i], request->consecutiveDays))
                            {
                              endDateWithHist[0][levelID].vec_d[i] = ivdate;
                            }
                        }
                    }
                }
              // update nmiss for saving data in GRIB
              startCount[levelID].nmiss = fieldNumMiss(startCount[levelID]);
              endCount[levelID].nmiss = fieldNumMiss(endCount[levelID]);
              startDateWithHist[1][levelID].nmiss = fieldNumMiss(startDateWithHist[1][levelID]);
              startDateWithHist[0][levelID].nmiss = fieldNumMiss(startDateWithHist[0][levelID]);
              endDateWithHist[1][levelID].nmiss = fieldNumMiss(endDateWithHist[1][levelID]);
              endDateWithHist[0][levelID].nmiss = fieldNumMiss(endDateWithHist[0][levelID]);
            }

          ovdate = ivdate;
          ovtime = ivtime;
          nsets++;
          itsID++;
        }

      if (nrecs == 0 && nsets == 0) break;

      adjustEndDate(nlevels, gridsize, yvals, missval, ovdate, startDateWithHist, endDateWithHist);

      /*  compute and write GSL for the previous year
       *  AND
       *  write the current start/end dates into the history
       *
       *  this is the default action if more than a year is available */
      if (yearcnt != 0)
        {
          computeGsl(nlevels, gridsize, yvals, missval, startDateWithHist, endDateWithHist, gslDuration, gslFirstDay, false);

          // values of the privous year
          {
            writeGslStream(ostreamID, otaxisID, otsID, ovarID1, ovarID2, ivlistID1, FIRST_VAR_ID, gslDuration, gslFirstDay,
                           cdiEncodeDate(ovdate / 10000 - 1, 12, 31), ovtime, nlevels);
            otsID++;
          }
        }

      /*  if there is a previous year */
      if (ovdate != ivdate)
        {
          /*  if the first year of data was processed, the history has to
           *  be checked befor it get's updated. This is necessary, if a
           *  growing period on the southern hemisphere was found. Otherwise,
           *  it would get overwritten. */
          if (isFirstYear)
            {
              // Check for non missing values, i.e. is there any data for the previous year?
              if (fldhvs(startDateWithHist[1], nlevels))
                {
                  computeGsl(nlevels, gridsize, yvals, missval, startDateWithHist, endDateWithHist, gslDuration, gslFirstDay,
                             false);
                  {
                    writeGslStream(ostreamID, otaxisID, otsID, ovarID1, ovarID2, ivlistID1, FIRST_VAR_ID, gslDuration, gslFirstDay,
                                   cdiEncodeDate(ovdate / 10000 - 1, 12, 31), ovtime, nlevels);
                    otsID++;
                  }
                }
              isFirstYear = false;
            }
#ifdef _OPENMP
#pragma omp sections
#endif
          {
            updateHist(startDateWithHist, nlevels, gridsize, yvals, false);
#ifdef _OPENMP
#pragma omp section
#endif
            updateHist(endDateWithHist, nlevels, gridsize, yvals, true);
          }
        }
      else  // process the current year, this only happens, if the last timestep is reached OR if data for only one year is present
        {
          computeGsl(nlevels, gridsize, yvals, missval, startDateWithHist, endDateWithHist, gslDuration, gslFirstDay, true);
          {
            writeGslStream(ostreamID, otaxisID, otsID, ovarID1, ovarID2, ivlistID1, FIRST_VAR_ID, gslDuration, gslFirstDay, ovdate,
                           ovtime, nlevels);
            otsID++;
          }
        }
      yearcnt++;

      if (nrecs == 0) break;
    }

  cdoStreamClose(ostreamID);
  cdoStreamClose(istreamID2);
  cdoStreamClose(istreamID1);
}