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

  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  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 "process_int.h"
#include "cdo_wtime.h"
#include "remap.h"
#include "remap_store_link.h"
#include "cdo_options.h"
#include "progress.h"
#include "cimdOmp.h"
#include <mpim_grid.h>

extern "C"
{
#include "lib/yac/clipping.h"
#include "lib/yac/area.h"
}

struct search_t
{
  enum yac_edge_type *src_edge_type;
  size_t srch_corners;
  size_t max_srch_cells;
  Varray<double> partial_areas;
  Varray<double> partial_weights;
  Varray<grid_cell> src_grid_cells;
  Varray<grid_cell> overlap_buffer;
};

static void
cellSearchRealloc(size_t num_srch_cells, search_t &search)
{
  if (num_srch_cells > search.max_srch_cells)
    {
      search.partial_areas.resize(num_srch_cells);
      search.partial_weights.resize(num_srch_cells);
      search.overlap_buffer.resize(num_srch_cells);
      search.src_grid_cells.resize(num_srch_cells);

      for (size_t n = search.max_srch_cells; n < num_srch_cells; ++n)
        {
          search.overlap_buffer[n].array_size = 0;
          search.overlap_buffer[n].num_corners = 0;
          search.overlap_buffer[n].edge_type = nullptr;
          search.overlap_buffer[n].coordinates_x = nullptr;
          search.overlap_buffer[n].coordinates_y = nullptr;
          search.overlap_buffer[n].coordinates_xyz = nullptr;

          search.src_grid_cells[n].array_size = search.srch_corners;
          search.src_grid_cells[n].num_corners = search.srch_corners;
          search.src_grid_cells[n].edge_type = search.src_edge_type;
          search.src_grid_cells[n].coordinates_x = new double[search.srch_corners];
          search.src_grid_cells[n].coordinates_y = new double[search.srch_corners];
          search.src_grid_cells[n].coordinates_xyz = new double[search.srch_corners][3];
        }

      search.max_srch_cells = num_srch_cells;
    }
}

static void
cellSearchFree(search_t &search)
{
  const auto max_srch_cells = search.max_srch_cells;
  for (size_t n = 0; n < max_srch_cells; n++)
    {
      if (search.overlap_buffer[n].array_size > 0)
        {
          free(search.overlap_buffer[n].coordinates_x);
          free(search.overlap_buffer[n].coordinates_y);
          if (search.overlap_buffer[n].coordinates_xyz) free(search.overlap_buffer[n].coordinates_xyz);
          if (search.overlap_buffer[n].edge_type) free(search.overlap_buffer[n].edge_type);
        }

      delete[] search.src_grid_cells[n].coordinates_x;
      delete[] search.src_grid_cells[n].coordinates_y;
      delete[] search.src_grid_cells[n].coordinates_xyz;
    }

  varrayFree(search.partial_areas);
  varrayFree(search.partial_weights);
  varrayFree(search.overlap_buffer);
  varrayFree(search.src_grid_cells);
}

static void
boundbox_from_corners1r(size_t ic, size_t nc, const Varray<double> &corner_lon, const Varray<double> &corner_lat,
                        float *restrict bound_box)
{
  const auto inc = ic * nc;

  float clat = corner_lat[inc];
  float clon = corner_lon[inc];

  bound_box[0] = clat;
  bound_box[1] = clat;
  bound_box[2] = clon;
  bound_box[3] = clon;

  for (size_t j = 1; j < nc; ++j)
    {
      clat = corner_lat[inc + j];
      clon = corner_lon[inc + j];

      if (clat < bound_box[0]) bound_box[0] = clat;
      if (clat > bound_box[1]) bound_box[1] = clat;
      if (clon < bound_box[2]) bound_box[2] = clon;
      if (clon > bound_box[3]) bound_box[3] = clon;
    }

  if (fabsf(bound_box[3] - bound_box[2]) > PI_f)
    {
      bound_box[2] = 0;
      bound_box[3] = PI2_f;
    }
  /*
  if ( fabsf(bound_box[3] - bound_box[2]) > PI_f )
    {
      if ( bound_box[3] > bound_box[2] && (bound_box[3]-PI2_f) < 0.0f )
        {
          float tmp = bound_box[2];
          bound_box[2] = bound_box[3] - PI2_f;
          bound_box[3] = tmp;
        }
    }
  */
}

static inline double
gridcell_area(const grid_cell &cell)
{
  return yac_huiliers_area(cell);
}

static void
cdo_compute_overlap_areas(size_t n, search_t &search, const grid_cell &target_cell)
{
  auto overlap_cells = search.overlap_buffer.data();

  // Do the clipping and get the cell for the overlapping area
  yac_cell_clipping(n, search.src_grid_cells.data(), target_cell, overlap_cells);

  // Get the partial areas for the overlapping regions
  for (size_t i = 0; i < n; i++) search.partial_areas[i] = gridcell_area(overlap_cells[i]);

#ifdef VERBOSE
  for (size_t i = 0; i < n; i++) printf("overlap area : %lf\n", search.partial_areas[i]);
#endif
}

static constexpr double tol = 1.0e-12;

enum cell_type
{
  UNDEF_CELL,
  LON_LAT_CELL,
  LAT_CELL,
  GREAT_CIRCLE_CELL,
  MIXED_CELL
};

static void
cdo_compute_concave_overlap_areas(size_t N, search_t &search, const grid_cell &target_cell, double target_node_x,
                                  double target_node_y)
{
  auto &partial_areas = search.partial_areas;
  auto overlap_cells = search.overlap_buffer.data();
  auto source_cell = search.src_grid_cells.data();

  double coordinates_x[3] = { -1, -1, -1 };
  double coordinates_y[3] = { -1, -1, -1 };
  double coordinates_xyz[3][3] = { { -1, -1, -1 }, { -1, -1, -1 }, { -1, -1, -1 } };
  enum yac_edge_type edge_types[3] = { GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE };

  grid_cell target_partial_cell;
  target_partial_cell.array_size = 3;
  target_partial_cell.coordinates_x = coordinates_x;
  target_partial_cell.coordinates_y = coordinates_y;
  target_partial_cell.coordinates_xyz = coordinates_xyz;
  target_partial_cell.edge_type = edge_types;
  target_partial_cell.num_corners = 3;

  // Do the clipping and get the cell for the overlapping area

  for (size_t n = 0; n < N; n++) partial_areas[n] = 0.0;

  // common node point to all partial target cells
  target_partial_cell.coordinates_x[0] = target_node_x;
  target_partial_cell.coordinates_y[0] = target_node_y;

  auto partial_cell_xyz = target_partial_cell.coordinates_xyz;
  auto cell_xyz = target_cell.coordinates_xyz;

  gcLLtoXYZ(target_node_x, target_node_y, partial_cell_xyz[0]);


  for (size_t num_corners = 0; num_corners < target_cell.num_corners; ++num_corners)
    {
      const auto corner_a = num_corners;
      const auto corner_b = (num_corners + 1) % target_cell.num_corners;

      // skip clipping and area calculation for degenerated triangles
      //
      // If this is not sufficient, instead we can try something like:
      //
      //     point_list target_list
      //     init_point_list(&target_list);
      //     generate_point_list(&target_list, target_cell);
      //     grid_cell temp_target_cell;
      //     generate_overlap_cell(target_list, temp_target_cell);
      //     free_point_list(&target_list);
      //
      // and use temp_target_cell for triangulation.
      //
      // Compared to the if statement below the alternative seems to be quite costly.

      if (((std::fabs(cell_xyz[corner_a][0] - cell_xyz[corner_b][0]) < tol)
           && (std::fabs(cell_xyz[corner_a][1] - cell_xyz[corner_b][1]) < tol)
           && (std::fabs(cell_xyz[corner_a][2] - cell_xyz[corner_b][2]) < tol))
          || ((std::fabs(cell_xyz[corner_a][0] - partial_cell_xyz[0][0]) < tol)
              && (std::fabs(cell_xyz[corner_a][1] - partial_cell_xyz[0][1]) < tol)
              && (std::fabs(cell_xyz[corner_a][2] - partial_cell_xyz[0][2]) < tol))
          || ((std::fabs(cell_xyz[corner_b][0] - partial_cell_xyz[0][0]) < tol)
              && (std::fabs(cell_xyz[corner_b][1] - partial_cell_xyz[0][1]) < tol)
              && (std::fabs(cell_xyz[corner_b][2] - partial_cell_xyz[0][2]) < tol)))
        continue;

      target_partial_cell.coordinates_x[1] = target_cell.coordinates_x[corner_a];
      target_partial_cell.coordinates_y[1] = target_cell.coordinates_y[corner_a];
      target_partial_cell.coordinates_x[2] = target_cell.coordinates_x[corner_b];
      target_partial_cell.coordinates_y[2] = target_cell.coordinates_y[corner_b];

      partial_cell_xyz[1][0] = cell_xyz[corner_a][0];
      partial_cell_xyz[1][1] = cell_xyz[corner_a][1];
      partial_cell_xyz[1][2] = cell_xyz[corner_a][2];
      partial_cell_xyz[2][0] = cell_xyz[corner_b][0];
      partial_cell_xyz[2][1] = cell_xyz[corner_b][1];
      partial_cell_xyz[2][2] = cell_xyz[corner_b][2];

      yac_cell_clipping(N, source_cell, target_partial_cell, overlap_cells);

      // Get the partial areas for the overlapping regions as sum over the partial target cells.
      for (size_t n = 0; n < N; n++)
        {
          partial_areas[n] += gridcell_area(overlap_cells[n]);
          // we cannot use pole_area because it is rather inaccurate for great
          // circle edges that are nearly circles of longitude
          // partial_areas[n] = pole_area (overlap_cells[n]);
        }
    }

#ifdef VERBOSE
  for (size_t n = 0; n < N; n++) printf("overlap area %zu: %lf\n", n, partial_areas[n]);
#endif
}

static int
get_lonlat_circle_index(RemapGrid *remap_grid)
{
  int lonlat_circle_index = -1;

  if (remap_grid->num_cell_corners == 4)
    {
      if (remap_grid->type == RemapGridType::Reg2D)
        {
          lonlat_circle_index = 1;
        }
      else
        {
          const auto &clon = remap_grid->cell_corner_lon;
          const auto &clat = remap_grid->cell_corner_lat;
          const auto gridsize = remap_grid->size;
          size_t iadd = (gridsize < 100) ? 1 : gridsize / 30 - 1;
          size_t num_i = 0, num_eq0 = 0, num_eq1 = 0;

          for (size_t i = 0; i < gridsize; i += iadd)
            {
              const auto i4 = i * 4;
              num_i++;
              // clang-format off
              if (IS_EQUAL(clon[i4 + 1], clon[i4 + 2]) && IS_EQUAL(clon[i4 + 3], clon[i4 + 0]) &&
                  IS_EQUAL(clat[i4 + 0], clat[i4 + 1]) && IS_EQUAL(clat[i4 + 2], clat[i4 + 3]))
                {
                  num_eq1++;
                }
              else if (IS_EQUAL(clon[i4 + 0], clon[i4 + 1]) && IS_EQUAL(clon[i4 + 2], clon[i4 + 3]) &&
                       IS_EQUAL(clat[i4 + 1], clat[i4 + 2]) && IS_EQUAL(clat[i4 + 3], clat[i4 + 0]))
                {
                  num_eq0++;
                }
              // clang-format on
            }

          if (num_i == num_eq1) lonlat_circle_index = 1;
          if (num_i == num_eq0) lonlat_circle_index = 0;
        }
    }

  // printf("lonlat_circle_index %d\n", lonlat_circle_index);

  return lonlat_circle_index;
}

static void
remapNormalizeWeights(RemapGrid *tgt_grid, RemapVars &rv)
{
  // Include centroids in weights and normalize using destination area if requested
  auto num_links = rv.num_links;
  auto num_wts = rv.num_wts;

  if (rv.normOpt == NormOpt::DESTAREA)
    {
      const auto &cell_area = tgt_grid->cell_area;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(num_wts, num_links, rv, cell_area)
#endif
      for (size_t n = 0; n < num_links; ++n)
        {
          const auto i = rv.tgt_cell_add[n];  // current linear address for target grid cell
          const auto norm_factor = IS_NOT_EQUAL(cell_area[i], 0.0) ? 1.0 / cell_area[i] : 0.0;
          rv.wts[n * num_wts] *= norm_factor;
        }
    }
  else if (rv.normOpt == NormOpt::FRACAREA)
    {
      const auto &cell_frac = tgt_grid->cell_frac;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(num_wts, num_links, rv, cell_frac)
#endif
      for (size_t n = 0; n < num_links; ++n)
        {
          const auto i = rv.tgt_cell_add[n];  // current linear address for target grid cell
          const auto norm_factor = IS_NOT_EQUAL(cell_frac[i], 0.0) ? 1.0 / cell_frac[i] : 0.0;
          rv.wts[n * num_wts] *= norm_factor;
        }
    }
  else if (rv.normOpt == NormOpt::NONE)
    {
    }
}

static void
remapNormalizeWeights(NormOpt normOpt, double cell_area, double cell_frac, size_t num_weights, double *weights)
{
  if (normOpt == NormOpt::DESTAREA)
    {
      const auto norm_factor = IS_NOT_EQUAL(cell_area, 0.0) ? 1.0 / cell_area : 0.0;
      for (size_t i = 0; i < num_weights; ++i) weights[i] *= norm_factor;
    }
  else if (normOpt == NormOpt::FRACAREA)
    {
      const auto norm_factor = IS_NOT_EQUAL(cell_frac, 0.0) ? 1.0 / cell_frac : 0.0;
      for (size_t i = 0; i < num_weights; ++i) weights[i] *= norm_factor;
    }
  else if (normOpt == NormOpt::NONE)
    {
    }
}

static void
setCellCoordinatesYac(RemapGridType remapGridType, const size_t cellAddr, const size_t numCellCorners, RemapGrid *remap_grid,
                      const grid_cell &yac_grid_cell)
{
  auto x = yac_grid_cell.coordinates_x;
  auto y = yac_grid_cell.coordinates_y;

  if (remapGridType == RemapGridType::Reg2D)
    {
      const auto nx = remap_grid->dims[0];
      const auto iy = cellAddr / nx;
      const auto ix = cellAddr - iy * nx;
      const auto reg2d_corner_lon = &remap_grid->reg2d_corner_lon[ix];
      const auto reg2d_corner_lat = &remap_grid->reg2d_corner_lat[iy];
      static const int xi[4] = { 0, 1, 1, 0 };
      static const int yi[4] = { 0, 0, 1, 1 };
      for (int i = 0; i < 4; ++i) x[i] = reg2d_corner_lon[xi[i]];
      for (int i = 0; i < 4; ++i) y[i] = reg2d_corner_lat[yi[i]];
    }
  else
    {
      const auto cell_corner_lon = &remap_grid->cell_corner_lon[cellAddr * numCellCorners];
      const auto cell_corner_lat = &remap_grid->cell_corner_lat[cellAddr * numCellCorners];
      for (size_t ic = 0; ic < numCellCorners; ++ic) x[ic] = cell_corner_lon[ic];
      for (size_t ic = 0; ic < numCellCorners; ++ic) y[ic] = cell_corner_lat[ic];
    }

  double(*restrict xyz)[3] = yac_grid_cell.coordinates_xyz;
  for (size_t ic = 0; ic < numCellCorners; ++ic) gcLLtoXYZ(x[ic], y[ic], xyz[ic]);
}

static void
setCoordinatesYac(const size_t numCells, RemapGridType remapGridType, const std::vector<size_t> &cellIndices,
                  const size_t numCellCorners, RemapGrid *remap_grid, Varray<grid_cell> &yac_grid_cell)
{
  for (size_t n = 0; n < numCells; ++n)
    setCellCoordinatesYac(remapGridType, cellIndices[n], numCellCorners, remap_grid, yac_grid_cell[n]);
}

static void
reg2d_bound_box(RemapGrid *remap_grid, double *grid_bound_box)
{
  const auto nx = remap_grid->dims[0];
  const auto ny = remap_grid->dims[1];
  const auto &reg2d_corner_lon = remap_grid->reg2d_corner_lon;
  const auto &reg2d_corner_lat = remap_grid->reg2d_corner_lat;

  grid_bound_box[0] = reg2d_corner_lat[0];
  grid_bound_box[1] = reg2d_corner_lat[ny];
  if (grid_bound_box[0] > grid_bound_box[1])
    {
      grid_bound_box[0] = reg2d_corner_lat[ny];
      grid_bound_box[1] = reg2d_corner_lat[0];
    }
  grid_bound_box[2] = reg2d_corner_lon[0];
  grid_bound_box[3] = reg2d_corner_lon[nx];
}

static void
greatCircleTypeInit(std::vector<enum yac_edge_type> &great_circle_type, size_t src_num_cell_corners, size_t tgt_num_cell_corners)
{
  auto max_num_cell_corners = src_num_cell_corners;
  if (tgt_num_cell_corners > max_num_cell_corners) max_num_cell_corners = tgt_num_cell_corners;

  great_circle_type.resize(max_num_cell_corners);
  for (size_t i = 0; i < max_num_cell_corners; ++i) great_circle_type[i] = GREAT_CIRCLE;
}

static void
scaleCellFrac(size_t numCells, Varray<double> &cellFrac, const Varray<double> &cellArea)
{
  for (size_t i = 0; i < numCells; ++i)
    if (IS_NOT_EQUAL(cellArea[i], 0)) cellFrac[i] /= cellArea[i];
}

static void
yacGridCellInit(grid_cell &cell, size_t numCellCorners, enum yac_edge_type *edgeType)
{
  cell.array_size = numCellCorners;
  cell.num_corners = numCellCorners;
  cell.edge_type = edgeType;
  cell.coordinates_x = new double[numCellCorners];
  cell.coordinates_y = new double[numCellCorners];
  cell.coordinates_xyz = new double[numCellCorners][3];
}

static void
yacGridCellFree(const grid_cell &cell)
{
  delete[] cell.coordinates_x;
  delete[] cell.coordinates_y;
  delete[] cell.coordinates_xyz;
}

static void
vec_add_weights(Varray<double> &vec, size_t num_weights, const Varray<double> &weights, const std::vector<size_t> &srch_add)
{
  for (size_t i = 0; i < num_weights; ++i)
    {
      const auto partial_weight = weights[i];
      const auto cell_add = srch_add[i];
#ifdef _OPENMP
#pragma omp atomic
#endif
      vec[cell_add] += partial_weight;
    }
}

static size_t
remove_invalid_areas(size_t num_srch_cells, Varray<double> &partial_areas, std::vector<size_t> &srch_add)
{
  size_t n = 0;
  for (size_t i = 0; i < num_srch_cells; ++i)
    {
      if (partial_areas[i] > 0)
        {
          partial_areas[n] = partial_areas[i];
          srch_add[n] = srch_add[i];
          n++;
        }
    }

  return n;
}

static size_t
remove_invalid_weights(size_t grid_size, size_t num_weights, Varray<double> &partial_weights, std::vector<size_t> &srch_add)
{
  size_t n = 0;
  for (size_t i = 0; i < num_weights; ++i)
    {
      const auto cell_add = (partial_weights[i] > 0) ? srch_add[i] : grid_size;
      if (cell_add != grid_size)
        {
          partial_weights[n] = partial_weights[i];
          srch_add[n] = cell_add;
          n++;
        }
    }

  return n;
}

static size_t
remove_unmask_weights(const Varray<int> &grid_mask, size_t num_weights, Varray<double> &partial_weights, std::vector<size_t> &srch_add)
{
  size_t n = 0;
  for (size_t i = 0; i < num_weights; ++i)
    {
      const auto cell_add = srch_add[i];
      /*
        Store the appropriate addresses and weights.
        Also add contributions to cell areas.
        The source grid mask is the master mask.
      */
      if (grid_mask[cell_add])
        {
          partial_weights[n] = partial_weights[i];
          srch_add[n] = cell_add;
          n++;
        }
    }

  return n;
}

void
remapConservWeights(RemapSearch &rsearch, RemapVars &rv)
{
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;

  bool lcheck = true;

  // Variables necessary if segment manages to hit pole
  auto srcGridType = src_grid->type;
  auto tgtGridType = tgt_grid->type;

  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);

  progress::init();

  double start = Options::cdoVerbose ? cdo_get_wtime() : 0;

  auto src_grid_size = src_grid->size;
  auto tgt_grid_size = tgt_grid->size;

  auto src_num_cell_corners = src_grid->num_cell_corners;
  auto tgt_num_cell_corners = tgt_grid->num_cell_corners;

  enum yac_edge_type lonlat_circle_type[] = { LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE };
  std::vector<enum yac_edge_type> great_circle_type;
  greatCircleTypeInit(great_circle_type, src_num_cell_corners, tgt_num_cell_corners);

  auto src_edge_type = great_circle_type.data();
  auto tgt_edge_type = great_circle_type.data();

  enum cell_type target_cell_type = UNDEF_CELL;

  if (src_num_cell_corners == 4)
    {
      const auto lonlat_circle_index = get_lonlat_circle_index(src_grid);
      if (lonlat_circle_index >= 0) src_edge_type = &lonlat_circle_type[lonlat_circle_index];
    }

  if (tgt_num_cell_corners == 4)
    {
      const auto lonlat_circle_index = get_lonlat_circle_index(tgt_grid);
      if (lonlat_circle_index >= 0)
        {
          target_cell_type = LON_LAT_CELL;
          tgt_edge_type = &lonlat_circle_type[lonlat_circle_index];
        }
    }

  if (!(tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL))
    {
      if (tgt_grid->cell_center_lon.empty() || tgt_grid->cell_center_lat.empty())
        cdoAbort("Internal problem (%s): missing target point coordinates!", __func__);
    }

  std::vector<grid_cell> tgt_grid_cell(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i) yacGridCellInit(tgt_grid_cell[i], tgt_num_cell_corners, tgt_edge_type);

  std::vector<search_t> search(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i)
    {
      search[i].srch_corners = src_num_cell_corners;
      search[i].src_edge_type = src_edge_type;
      search[i].max_srch_cells = 0;
    }

  std::vector<std::vector<size_t>> srch_add(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i) srch_add[i].resize(src_grid_size);

  auto srch_corners = src_num_cell_corners;  // num of corners of srch cells

  double src_grid_bound_box[4];
  if (srcGridType == RemapGridType::Reg2D) reg2d_bound_box(src_grid, src_grid_bound_box);

  std::vector<WeightLinks> weightLinks(tgt_grid_size);

  double findex = 0;

  size_t num_srch_cells_stat[3];
  num_srch_cells_stat[0] = 0;
  num_srch_cells_stat[1] = 100000;
  num_srch_cells_stat[2] = 0;

  extern CellSearchMethod cellSearchMethod;
  bool useCellsearch = cellSearchMethod == CellSearchMethod::spherepart || srcGridType == RemapGridType::Reg2D;

  // Loop over destination grid

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic) default(none)                                                                    \
    shared(findex, rsearch, srcGridType, tgtGridType, src_grid_bound_box, rv, Options::cdoVerbose, tgt_num_cell_corners,    \
           target_cell_type, weightLinks, srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, search, srch_add, \
           tgt_grid_cell, num_srch_cells_stat, useCellsearch)
#endif
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
    {
      const auto ompthID = cdo_omp_get_thread_num();

#ifdef _OPENMP
#pragma omp atomic
#endif
      findex++;
      if (ompthID == 0) progress::update(0, 1, findex / tgt_grid_size);

      weightLinks[tgt_cell_add].nlinks = 0;

      if (!tgt_grid->mask[tgt_cell_add]) continue;

      setCellCoordinatesYac(tgtGridType, tgt_cell_add, tgt_num_cell_corners, tgt_grid, tgt_grid_cell[ompthID]);

      // Get search cells
      size_t num_srch_cells;
      if (useCellsearch)
        {
          // num_srch_cells = remapSearchCells(rsearch, target_cell_type == LON_LAT_CELL, tgt_grid_cell[ompthID],
          // srch_add[ompthID].data());
          num_srch_cells
              = remapSearchCells(rsearch, tgtGridType == RemapGridType::Reg2D, tgt_grid_cell[ompthID], srch_add[ompthID].data());
        }
      else
        {
          float tgt_cell_bound_box_r[4];
          boundbox_from_corners1r(tgt_cell_add, tgt_num_cell_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat,
                                  tgt_cell_bound_box_r);

          num_srch_cells
              = get_srch_cells(tgt_cell_add, rsearch.tgtBins, rsearch.srcBins, tgt_cell_bound_box_r, srch_add[ompthID].data());
        }

      if (1 && Options::cdoVerbose)
        {
          num_srch_cells_stat[0] += num_srch_cells;
          if (num_srch_cells < num_srch_cells_stat[1]) num_srch_cells_stat[1] = num_srch_cells;
          if (num_srch_cells > num_srch_cells_stat[2]) num_srch_cells_stat[2] = num_srch_cells;
        }

      if (0 && Options::cdoVerbose) printf("tgt_cell_add %zu  num_srch_cells %zu\n", tgt_cell_add, num_srch_cells);

      if (num_srch_cells == 0) continue;

      // Create search arrays

      cellSearchRealloc(num_srch_cells, search[ompthID]);

      setCoordinatesYac(num_srch_cells, srcGridType, srch_add[ompthID], srch_corners, src_grid, search[ompthID].src_grid_cells);

      if (tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL)
        {
          cdo_compute_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID]);
        }
      else
        {
          cdo_compute_concave_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID],
                                            tgt_grid->cell_center_lon[tgt_cell_add], tgt_grid->cell_center_lat[tgt_cell_add]);
        }

      auto &partial_areas = search[ompthID].partial_areas;
      auto &partial_weights = search[ompthID].partial_weights;

      auto num_weights = remove_invalid_areas(num_srch_cells, partial_areas, srch_add[ompthID]);

      const auto tgt_area = gridcell_area(tgt_grid_cell[ompthID]);
      tgt_grid->cell_area[tgt_cell_add] = tgt_area;

      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] = partial_areas[i] / tgt_area;

      if (rv.normOpt == NormOpt::FRACAREA) yac_correct_weights(num_weights, partial_weights.data());

      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] *= tgt_area;

      num_weights = remove_invalid_weights(src_grid_size, num_weights, partial_weights, srch_add[ompthID]);

      vec_add_weights(src_grid->cell_area, num_weights, partial_weights, srch_add[ompthID]);

      num_weights = remove_unmask_weights(src_grid->mask, num_weights, partial_weights, srch_add[ompthID]);

      vec_add_weights(src_grid->cell_frac, num_weights, partial_weights, srch_add[ompthID]);

      tgt_grid->cell_frac[tgt_cell_add] = varraySum(num_weights, partial_weights);

      store_weightlinks(1, num_weights, srch_add[ompthID].data(), partial_weights.data(), tgt_cell_add, weightLinks);
    }

  progress::update(0, 1, 1);

  if (1 && Options::cdoVerbose)
    {
      cdoPrint("Num search cells min,mean,max :  %zu  %3.1f  %zu", num_srch_cells_stat[1],
               num_srch_cells_stat[0] / (double) tgt_grid_size, num_srch_cells_stat[2]);
    }

  // Finished with all cells: deallocate search arrays
  for (auto ompthID = 0; ompthID < Threading::ompNumThreads; ++ompthID)
    {
      cellSearchFree(search[ompthID]);
      yacGridCellFree(tgt_grid_cell[ompthID]);
    }

  weightLinksToRemapLinks(1, tgt_grid_size, weightLinks, rv);

  // Normalize weights using destination area if requested
  remapNormalizeWeights(tgt_grid, rv);

  if (Options::cdoVerbose) cdoPrint("Total number of links = %zu", rv.num_links);

  scaleCellFrac(src_grid_size, src_grid->cell_frac, src_grid->cell_area);
  scaleCellFrac(tgt_grid_size, tgt_grid->cell_frac, tgt_grid->cell_area);

  // Perform some error checking on final weights
  if (lcheck)
    {
      remapCheckArea(src_grid_size, &src_grid->cell_area[0], "Source");
      remapCheckArea(tgt_grid_size, &tgt_grid->cell_area[0], "Target");

      remapVarsCheckWeights(rv);
    }

  if (Options::cdoVerbose) cdoPrint("Cells search: %.2f seconds", cdo_get_wtime() - start);
}  // remapConservWeights

static double
conservRemap(const double *restrict src_array, size_t num_weights, const Varray<double> &weights, const std::vector<size_t> &src_add)
{
  double tgt_point = 0.;
  for (size_t i = 0; i < num_weights; ++i) tgt_point += src_array[src_add[i]] * weights[i];

  return tgt_point;
}

void
remapConserv(NormOpt normOpt, RemapSearch &rsearch, const double *restrict src_array, double *restrict tgt_array, double missval)
{
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;

  bool lcheck = true;

  // Variables necessary if segment manages to hit pole
  auto srcGridType = src_grid->type;
  auto tgtGridType = tgt_grid->type;

  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);

  progress::init();

  double start = Options::cdoVerbose ? cdo_get_wtime() : 0;

  auto src_grid_size = src_grid->size;
  auto tgt_grid_size = tgt_grid->size;

  Varray<int> src_grid_mask(src_grid_size);
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(src_grid_size, src_array, src_grid_mask, missval)
#endif
  for (size_t i = 0; i < src_grid_size; ++i) src_grid_mask[i] = !DBL_IS_EQUAL(src_array[i], missval);

  auto src_num_cell_corners = src_grid->num_cell_corners;
  auto tgt_num_cell_corners = tgt_grid->num_cell_corners;

  enum yac_edge_type lonlat_circle_type[] = { LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE };
  std::vector<enum yac_edge_type> great_circle_type;
  greatCircleTypeInit(great_circle_type, src_num_cell_corners, tgt_num_cell_corners);

  auto src_edge_type = great_circle_type.data();
  auto tgt_edge_type = great_circle_type.data();

  enum cell_type target_cell_type = UNDEF_CELL;

  if (src_num_cell_corners == 4)
    {
      const auto lonlat_circle_index = get_lonlat_circle_index(src_grid);
      if (lonlat_circle_index >= 0) src_edge_type = &lonlat_circle_type[lonlat_circle_index];
    }

  if (tgt_num_cell_corners == 4)
    {
      const auto lonlat_circle_index = get_lonlat_circle_index(tgt_grid);
      if (lonlat_circle_index >= 0)
        {
          target_cell_type = LON_LAT_CELL;
          tgt_edge_type = &lonlat_circle_type[lonlat_circle_index];
        }
    }

  if (!(tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL))
    {
      if (tgt_grid->cell_center_lon.empty() || tgt_grid->cell_center_lat.empty())
        cdoAbort("Internal problem (%s): missing target point coordinates!", __func__);
    }

  std::vector<grid_cell> tgt_grid_cell(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i) yacGridCellInit(tgt_grid_cell[i], tgt_num_cell_corners, tgt_edge_type);

  std::vector<search_t> search(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i)
    {
      search[i].srch_corners = src_num_cell_corners;
      search[i].src_edge_type = src_edge_type;
      search[i].max_srch_cells = 0;
    }

  std::vector<std::vector<size_t>> srch_add(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i) srch_add[i].resize(src_grid_size);

  auto srch_corners = src_num_cell_corners;  // num of corners of srch cells

  double src_grid_bound_box[4];
  if (srcGridType == RemapGridType::Reg2D) reg2d_bound_box(src_grid, src_grid_bound_box);

  double findex = 0;

  size_t num_srch_cells_stat[3];
  num_srch_cells_stat[0] = 0;
  num_srch_cells_stat[1] = 100000;
  num_srch_cells_stat[2] = 0;

  extern CellSearchMethod cellSearchMethod;
  bool useCellsearch = cellSearchMethod == CellSearchMethod::spherepart || srcGridType == RemapGridType::Reg2D;

  // Loop over destination grid

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic) default(none)                                                                      \
    shared(findex, rsearch, srcGridType, tgtGridType, src_grid_bound_box, Options::cdoVerbose, tgt_num_cell_corners,          \
           target_cell_type, srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, search, srch_add, tgt_grid_cell, \
           num_srch_cells_stat, useCellsearch, src_array, tgt_array, missval, normOpt, src_grid_mask)
#endif
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
    {
      const auto ompthID = cdo_omp_get_thread_num();

#ifdef _OPENMP
#pragma omp atomic
#endif
      findex++;
      if (ompthID == 0) progress::update(0, 1, findex / tgt_grid_size);

      tgt_array[tgt_cell_add] = missval;

      if (!tgt_grid->mask[tgt_cell_add]) continue;

      setCellCoordinatesYac(tgtGridType, tgt_cell_add, tgt_num_cell_corners, tgt_grid, tgt_grid_cell[ompthID]);

      // Get search cells
      size_t num_srch_cells;
      if (useCellsearch)
        {
          // num_srch_cells = remapSearchCells(rsearch, target_cell_type == LON_LAT_CELL, tgt_grid_cell[ompthID],
          // srch_add[ompthID].data());
          num_srch_cells
              = remapSearchCells(rsearch, tgtGridType == RemapGridType::Reg2D, tgt_grid_cell[ompthID], srch_add[ompthID].data());
        }
      else
        {
          float tgt_cell_bound_box_r[4];
          boundbox_from_corners1r(tgt_cell_add, tgt_num_cell_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat,
                                  tgt_cell_bound_box_r);

          num_srch_cells
              = get_srch_cells(tgt_cell_add, rsearch.tgtBins, rsearch.srcBins, tgt_cell_bound_box_r, srch_add[ompthID].data());
        }

      if (1 && Options::cdoVerbose)
        {
          num_srch_cells_stat[0] += num_srch_cells;
          if (num_srch_cells < num_srch_cells_stat[1]) num_srch_cells_stat[1] = num_srch_cells;
          if (num_srch_cells > num_srch_cells_stat[2]) num_srch_cells_stat[2] = num_srch_cells;
        }

      if (0 && Options::cdoVerbose) printf("tgt_cell_add %zu  num_srch_cells %zu\n", tgt_cell_add, num_srch_cells);

      if (num_srch_cells == 0) continue;

      // Create search arrays

      cellSearchRealloc(num_srch_cells, search[ompthID]);

      setCoordinatesYac(num_srch_cells, srcGridType, srch_add[ompthID], srch_corners, src_grid, search[ompthID].src_grid_cells);

      if (tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL)
        {
          cdo_compute_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID]);
        }
      else
        {
          cdo_compute_concave_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID],
                                            tgt_grid->cell_center_lon[tgt_cell_add], tgt_grid->cell_center_lat[tgt_cell_add]);
        }

      auto &partial_areas = search[ompthID].partial_areas;
      auto &partial_weights = search[ompthID].partial_weights;

      auto num_weights = remove_invalid_areas(num_srch_cells, partial_areas, srch_add[ompthID]);

      const auto tgt_area = gridcell_area(tgt_grid_cell[ompthID]);
      tgt_grid->cell_area[tgt_cell_add] = tgt_area;

      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] = partial_areas[i] / tgt_area;

      if (normOpt == NormOpt::FRACAREA) yac_correct_weights(num_weights, partial_weights.data());

      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] *= tgt_area;

      num_weights = remove_invalid_weights(src_grid_size, num_weights, partial_weights, srch_add[ompthID]);

      num_weights = remove_unmask_weights(src_grid_mask, num_weights, partial_weights, srch_add[ompthID]);

      tgt_grid->cell_frac[tgt_cell_add] = varraySum(num_weights, partial_weights);

      if (num_weights)
        {
          sort_weights(num_weights, srch_add[ompthID].data(), partial_weights.data());
          // Normalize weights using destination area if requested
          remapNormalizeWeights(normOpt, tgt_area, tgt_grid->cell_frac[tgt_cell_add], num_weights, partial_weights.data());
          tgt_array[tgt_cell_add] = conservRemap(src_array, num_weights, partial_weights, srch_add[ompthID]);
        }
    }

  progress::update(0, 1, 1);

  if (1 && Options::cdoVerbose)
    {
      cdoPrint("Num search cells min,mean,max :  %zu  %3.1f  %zu", num_srch_cells_stat[1],
               num_srch_cells_stat[0] / (double) tgt_grid_size, num_srch_cells_stat[2]);
    }

  // Finished with all cells: deallocate search arrays

  for (auto ompthID = 0; ompthID < Threading::ompNumThreads; ++ompthID)
    {
      cellSearchFree(search[ompthID]);
      yacGridCellFree(tgt_grid_cell[ompthID]);
    }

  scaleCellFrac(tgt_grid_size, tgt_grid->cell_frac, tgt_grid->cell_area);

  // Perform some error checking on final weights
  if (lcheck) remapCheckArea(tgt_grid_size, &tgt_grid->cell_area[0], "Target");

  if (Options::cdoVerbose) cdoPrint("Cells search: %.2f seconds", cdo_get_wtime() - start);
}  // remapConserv