xref: /dports/science/cdo/cdo-2.0.0/src/ecautil.cc

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

  Copyright (C) 2006 Brockmann Consult

  Author: Ralf Quast

*/

#include <algorithm>
#include <assert.h>

#include <cdi.h>

#include "process_int.h"
#include "timebase.h"
#include "datetime.h"
#include "ecautil.h"

/**
 * Convert a Gregorian/Julian date to a Julian day number.
 *
 * The Gregorian calendar was adopted midday, October 15, 1582.
 */
static unsigned long
gregdate_to_julday(int year,  /* Gregorian year */
                   int month, /* Gregorian month (1-12) */
                   int day    /* Gregorian day (1-31) */
)
{
#if INT_MAX <= 0X7FFF
  long igreg = 15 + 31 * (10 + (12 * 1582));
  long iy; /* signed, origin 0 year */
  long ja; /* Julian century */
  long jm; /* Julian month */
  long jy; /* Julian year */
#else
  int igreg = 15 + 31 * (10 + (12 * 1582));
  int iy; /* signed, origin 0 year */
  int ja; /* Julian century */
  int jm; /* Julian month */
  int jy; /* Julian year */
#endif
  unsigned long julday; /* returned Julian day number */

  /*
   * Because there is no 0 BC or 0 AD, assume the user wants the start of
   * the common era if they specify year 0.
   */
  if (year == 0) year = 1;

  iy = year;
  if (year < 0) iy++;
  if (month > 2)
    {
      jy = iy;
      jm = month + 1;
    }
  else
    {
      jy = iy - 1;
      jm = month + 13;
    }

  /*
   *  Note: SLIGHTLY STRANGE CONSTRUCTIONS REQUIRED TO AVOID PROBLEMS WITH
   *        OPTIMISATION OR GENERAL ERRORS UNDER VMS!
   */
  julday = day + (int) (30.6001 * jm);
  if (jy >= 0)
    {
      julday += 365 * jy;
      julday += (unsigned long) (0.25 * jy);
    }
  else
    {
      double xi = 365.25 * jy;

      {
        int temp = xi * 100;
        if (temp % 100 != 0) xi -= 1;
      }
      // if (((int) xi) != xi) xi -= 1; // old version, caused warning (25.08.2020 O.Heidmann)
      julday += (int) xi;
    }
  julday += 1720995;

  if (day + (31 * (month + (12 * iy))) >= igreg)
    {
      ja = jy / 100;
      julday -= ja;
      julday += 2;
      julday += ja / 4;
    }

  return julday;
}

/**
 * Computes the day-of-year correspnding a given Gregorian date.
 *
 * @param date a Gregorian date in the form YYYYMMDD
 *
 * @return the day-of-year
 */
static unsigned long
day_of_year(int64_t date)
{
  const int year = date / 10000;
  const int month = (date - year * 10000) / 100;
  const int day = date - year * 10000 - month * 100;

  return gregdate_to_julday(year, month, day) - gregdate_to_julday(year, 1, 1) + 1;
}

/**
 * Counts the number of nonmissing values. The result of the operation
 * is computed according to the following rules:
 *
 * field1  field2  mode  result
 * a       b       0     a + 1
 * a       miss    0     a
 * miss    b       0     1
 * miss    miss    0     0
 *
 * a       b       1     a + 1
 * a       miss    1     0
 * miss    b       1     1
 * miss    miss    1     0
 *
 * a       b       n     b < n ? a : b > n ? a + 1 : a + n
 * a       miss    n     a
 * miss    b       n     b < n ? 0 : b
 * miss    miss    n     0
 *
 * @param field1 the 1st input field, also holds the result
 * @param field2 the 2nd input field
 * @param mode   the counting mode, must be an exact mathematical
 *               integer
 */
static void
count(Field &field1, const Field &field2, double mode)
{
  const auto missval1 = field1.missval;
  const auto missval2 = field2.missval;
  auto &array1 = field1.vec_d;
  const auto &array2 = field2.vec_d;

  const auto len = field1.size;
  if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);

  if (field1.nmiss)
    {
      for (size_t i = 0; i < len; i++)
        {
          if (DBL_IS_EQUAL(array2[i], missval2))
            {
              if (IS_EQUAL(mode, 1.0) || DBL_IS_EQUAL(array1[i], missval1)) array1[i] = 0.0;
              continue;
            }

          if (!DBL_IS_EQUAL(array1[i], missval1))
            {
              if (IS_EQUAL(mode, 0.0) || IS_EQUAL(mode, 1.0))
                array1[i] += 1.0;
              else if (DBL_IS_EQUAL(array2[i], mode) || array2[i] > mode)
                array1[i] += array2[i];
            }
          else
            {
              if (IS_EQUAL(mode, 0.0) || IS_EQUAL(mode, 1.0))
                array1[i] = 1.0;
              else if (array2[i] < mode)
                array1[i] = 0.0;
              else
                array1[i] = array2[i];
            }
        }

      field1.nmiss = field_num_miss(field1);
    }
  else
    {
      for (size_t i = 0; i < len; i++)
        {
          if (DBL_IS_EQUAL(array2[i], missval2))
            {
              if (IS_EQUAL(mode, 1.0)) array1[i] = 0.0;
              continue;
            }

          if (IS_EQUAL(mode, 0.0) || IS_EQUAL(mode, 1.0))
            array1[i] += 1.0;
          else if (DBL_IS_EQUAL(array2[i], mode) || array2[i] > mode)
            array1[i] += array2[i];
        }
    }
}

/**
 * Selects all field elements that compare to the corresponding
 * element of a reference field. The result of the operation is
 * computed according to the following rules:
 *
 * field1  field2  result
 * a       b       comp(a, b) ? a : miss
 * a       miss    miss
 * miss    b       miss
 * miss    miss    miss
 *
 * @param field1  the input field, also holds the result
 * @param field2  the reference field
 * @param compare the comparator
 */
static void
selcomp(Field &field1, const Field &field2, int (*compare)(double, double))
{
  const auto missval1 = field1.missval;
  const auto missval2 = field2.missval;
  auto &array1 = field1.vec_d;
  const auto &array2 = field2.vec_d;

  const auto len = field1.size;
  if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);

  if (field1.nmiss || field2.nmiss)
    {
      for (size_t i = 0; i < len; i++)
        if (DBL_IS_EQUAL(array1[i], missval1) || DBL_IS_EQUAL(array2[i], missval2) || !compare(array1[i], array2[i]))
          array1[i] = missval1;
    }
  else
    {
      for (size_t i = 0; i < len; i++)
        if (!compare(array1[i], array2[i])) array1[i] = missval1;
    }

  field1.nmiss = field_num_miss(field1);
}

/**
 * Selects all field elements that compare to a certain reference
 * value. The result of the operation is computed according to the
 * following rules:
 *
 * field  c      result
 * a      c      comp(a, c) ? a : miss
 * a      miss   miss
 * miss   c      miss
 * miss   miss   miss
 *
 * @param field   the input field, also holds the result
 * @param c       the refence value
 * @param compare the comparator
 */
static void
selcompc(Field &field, double c, int (*compare)(double, double))
{
  const auto missval = field.missval;
  auto &array = field.vec_d;

  const auto len = field.size;

  if (DBL_IS_EQUAL(c, missval))
    {
      for (size_t i = 0; i < len; i++) array[i] = missval;
    }
  else if (field.nmiss)
    {
      for (size_t i = 0; i < len; i++)
        if (DBL_IS_EQUAL(array[i], missval) || !compare(array[i], c)) array[i] = missval;
    }
  else
    {
      for (size_t i = 0; i < len; i++)
        if (!compare(array[i], c)) array[i] = missval;
    }

  field.nmiss = field_num_miss(field);
}

static int
le(double a, double b)
{
  return a <= b;
}

static int
lt(double a, double b)
{
  return a < b;
}

static int
ge(double a, double b)
{
  return a >= b;
}

static int
gt(double a, double b)
{
  return a > b;
}

static int
eq(double a, double b)
{
  return DBL_IS_EQUAL(a, b);
}

static int
ne(double a, double b)
{
  return !DBL_IS_EQUAL(a, b);
}

void
vfarnum(Field &field1, const Field &field2)
{
  count(field1, field2, 0.0);
}

void
vfarnum2(Field &field1, const Field &field2)
{
  count(field1, field2, 1.0);
}

void
vfarnum3(Field &field1, const Field &field2, double n)
{
  count(field1, field2, n);
}

void
vfarsel(Field &field1, const Field &field2)
{
  const auto missval1 = field1.missval;
  const auto missval2 = field2.missval;
  auto &array1 = field1.vec_d;
  const auto &array2 = field2.vec_d;

  const auto len = field1.size;
  if (len != field2.size) cdo_abort("Fields have different gridsize (%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] = missval1;
    }
  else
    {
      for (size_t i = 0; i < len; i++)
        if (IS_EQUAL(array2[i], 0.0)) array1[i] = missval1;
    }

  field1.nmiss = field_num_miss(field1);
}

void
vfarselle(Field &field1, const Field &field2)
{
  selcomp(field1, field2, le);
}

void
vfarsellt(Field &field1, const Field &field2)
{
  selcomp(field1, field2, lt);
}

void
vfarselge(Field &field1, const Field &field2)
{
  selcomp(field1, field2, ge);
}

void
vfarselgt(Field &field1, const Field &field2)
{
  selcomp(field1, field2, gt);
}

void
vfarseleq(Field &field1, const Field &field2)
{
  selcomp(field1, field2, eq);
}

void
vfarselne(Field &field1, const Field &field2)
{
  selcomp(field1, field2, ne);
}

void
vfarsellec(Field &field, double c)
{
  selcompc(field, c, le);
}

void
vfarselltc(Field &field, double c)
{
  selcompc(field, c, lt);
}

void
vfarselgec(Field &field, double c)
{
  selcompc(field, c, ge);
}

void
vfarseleqc(Field &field, double c)
{
  selcompc(field, c, eq);
}

void
vfarselnec(Field &field, double c)
{
  selcompc(field, c, ne);
}

void
vfarselgtc(Field &field, double c)
{
  selcompc(field, c, gt);
}

void
update_hist(FieldVector2D &field, int nlevels, size_t gridsize, const std::vector<double> &yvals, bool onlyNorth)
{
  for (int levelID = 0; levelID < nlevels; levelID++)
    for (size_t i = 0; i < gridsize; i++)
      if (onlyNorth)
        {
          if (yvals[i] >= 0.0) field[1][levelID].vec_d[i] = field[0][levelID].vec_d[i];
        }
      else
        field[1][levelID].vec_d[i] = field[0][levelID].vec_d[i];
}

void
define_mid_of_time(int frequency, int taxisID, int year, int month, int MaxMonths)
{
  int64_t vdate, vdateb, vdatebp1;
  int vtime;
  const auto vtime0 = cdiEncodeTime(0, 0, 0);

  const auto calendar = taxisInqCalendar(taxisID);

  if (frequency == 8)
    {
      vdateb = cdiEncodeDate(year, month, 1);

      int boundmonth = (month + 1 <= MaxMonths) ? month + 1 : 1;
      int boundyear = (boundmonth != 1) ? year : year + 1;
      vdatebp1 = cdiEncodeDate(boundyear, boundmonth, 1);
    }
  else
    {
      vdateb = cdiEncodeDate(year, 1, 1);
      vdatebp1 = cdiEncodeDate(year + 1, 1, 1);
    }

  const auto juldate1 = juldate_encode(calendar, vdateb, vtime0);
  const auto juldate2 = juldate_encode(calendar, vdatebp1, vtime0);

  const auto seconds = juldate_to_seconds(juldate_sub(juldate2, juldate1)) / 2;
  const auto juldatem = juldate_add_seconds(lround(seconds), juldate1);
  juldate_decode(calendar, juldatem, vdate, vtime);

  taxisDefVdate(taxisID, vdate);
  taxisDefVtime(taxisID, vtime);
}

void
adjust_end_date(int nlevels, size_t gridsize, const std::vector<double> &yvals, double missval, int64_t ovdate,
                const FieldVector2D &startDateWithHist, FieldVector2D &endDateWithHist)
{
  const auto ovdateSouth = std::min(cdiEncodeDate(ovdate / 10000, 6, 30), ovdate);

  for (int levelID = 0; levelID < nlevels; levelID++)
    {
      for (size_t i = 0; i < gridsize; i++)
        {
          // start with southern sphere
          if (yvals[i] < 0)
            {
              if (DBL_IS_EQUAL(startDateWithHist[1][levelID].vec_d[i], missval))
                {
                  endDateWithHist[0][levelID].vec_d[i] = missval;
                  continue;
                }
              if (DBL_IS_EQUAL(endDateWithHist[0][levelID].vec_d[i], missval))
                {
                  endDateWithHist[0][levelID].vec_d[i] = ovdateSouth;
                }
            }
          else
            {
              if (DBL_IS_EQUAL(startDateWithHist[0][levelID].vec_d[i], missval))
                {
                  endDateWithHist[0][levelID].vec_d[i] = missval;
                  continue;
                }

              if (DBL_IS_EQUAL(endDateWithHist[0][levelID].vec_d[i], missval))
                {
                  endDateWithHist[0][levelID].vec_d[i] = ovdate;
                }
            }
        }
    }
}

void
compute_gsl(int nlevels, size_t gridsize, std::vector<double> &yvals, double missval, FieldVector2D &startDateWithHist,
            FieldVector2D &endDateWithHist, FieldVector &gslDuration, FieldVector &gslFirstDay, bool useCurrentYear)
{
  double firstDay, duration;

  if (!useCurrentYear)
    {
      for (int levelID = 0; levelID < nlevels; levelID++)
        {
          for (size_t i = 0; i < gridsize; i++)
            {
              // start with southern sphere
              if (yvals[i] < 0.0)
                {
                  duration = (double) (date_to_julday(CALENDAR_PROLEPTIC, (int64_t) endDateWithHist[0][levelID].vec_d[i])
                                       - date_to_julday(CALENDAR_PROLEPTIC, (int64_t) startDateWithHist[1][levelID].vec_d[i]));
                }
              else
                {
                  duration = (double) (date_to_julday(CALENDAR_PROLEPTIC, (int64_t) endDateWithHist[1][levelID].vec_d[i])
                                       - date_to_julday(CALENDAR_PROLEPTIC, (int64_t) startDateWithHist[1][levelID].vec_d[i]));
                }

              if (DBL_IS_EQUAL(startDateWithHist[1][levelID].vec_d[i], missval))
                firstDay = missval;
              else
                firstDay = (double) day_of_year((int64_t) startDateWithHist[1][levelID].vec_d[i]);

              gslDuration[levelID].vec_d[i] = duration;
              gslFirstDay[levelID].vec_d[i] = firstDay;
            }
        }
    }
  else
    {
      // the current year can only have values for the northern hemisphere
      for (int levelID = 0; levelID < nlevels; levelID++)
        {
          for (size_t i = 0; i < gridsize; i++)
            {
              // start with southern sphere
              if (yvals[i] < 0.0)
                {
                  gslDuration[levelID].vec_d[i] = missval;
                  gslFirstDay[levelID].vec_d[i] = missval;
                }
              else
                {
                  duration = (double) (date_to_julday(CALENDAR_PROLEPTIC, (int64_t) endDateWithHist[0][levelID].vec_d[i])
                                       - date_to_julday(CALENDAR_PROLEPTIC, (int64_t) startDateWithHist[0][levelID].vec_d[i]));

                  if (DBL_IS_EQUAL(startDateWithHist[0][levelID].vec_d[i], missval))
                    firstDay = missval;
                  else
                    firstDay = (double) day_of_year((int64_t) startDateWithHist[0][levelID].vec_d[i]);

                  gslDuration[levelID].vec_d[i] = duration;
                  gslFirstDay[levelID].vec_d[i] = firstDay;
                }
            }
        }
    }

  for (int levelID = 0; levelID < nlevels; levelID++)
    {
      gslDuration[levelID].nmiss = field_num_miss(gslDuration[levelID]);
      gslFirstDay[levelID].nmiss = field_num_miss(gslFirstDay[levelID]);
    }
}

void
write_gsl_stream(CdoStreamID ostreamID, int otaxisID, int otsID, int ovarID1, int ovarID2, int ivlistID1, int first_var_id,
                 FieldVector &gslDuration, FieldVector &gslFirstDay, int64_t vdate, int vtime, int nlevels)
{
  (void) ivlistID1;
  (void) first_var_id;

  taxisDefVdate(otaxisID, vdate);
  taxisDefVtime(otaxisID, vtime);
  cdo_def_timestep(ostreamID, otsID);

  for (int levelID = 0; levelID < nlevels; levelID++)
    {
      cdo_def_record(ostreamID, ovarID1, levelID);
      cdo_write_record(ostreamID, gslDuration[levelID].vec_d.data(), gslDuration[levelID].nmiss);
    }

  for (int levelID = 0; levelID < nlevels; levelID++)
    {
      cdo_def_record(ostreamID, ovarID2, levelID);
      cdo_write_record(ostreamID, gslFirstDay[levelID].vec_d.data(), gslFirstDay[levelID].nmiss);
    }
}