remap_search_latbins.cc 13.4 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
  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.
*/
17
18
19
20
21

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

22
23
#include <algorithm>

24
#include "cdo_output.h"
25
#include "cdo_options.h"
26
#include <mpim_grid.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
#include "remap.h"
28

29
static void
30
calcBinAddr(GridSearchBins &searchBins)
31
{
32
33
  const size_t ncells = searchBins.ncells;
  const size_t nbins = searchBins.nbins;
34
35
36
37
  size_t *restrict bin_addr = &searchBins.bin_addr[0];
  const float *restrict bin_lats = &searchBins.bin_lats[0];
  const float *restrict cell_bound_box = &searchBins.cell_bound_box[0];

38
  for (size_t n = 0; n < nbins; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
41
      bin_addr[n * 2] = ncells;
      bin_addr[n * 2 + 1] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
    }
43

Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
  for (size_t nele = 0; nele < ncells; ++nele)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
    {
46
47
48
      const size_t nele4 = nele << 2;
      const float cell_bound_box_lat1 = cell_bound_box[nele4];
      const float cell_bound_box_lat2 = cell_bound_box[nele4 + 1];
49
50
      for (size_t n = 0; n < nbins; ++n)
        {
51
          const size_t n2 = n << 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
          if (cell_bound_box_lat1 <= bin_lats[n2 + 1] && cell_bound_box_lat2 >= bin_lats[n2])
53
            {
54
55
              bin_addr[n2] = std::min(nele, bin_addr[n2]);
              bin_addr[n2 + 1] = std::max(nele, bin_addr[n2 + 1]);
56
            }
57
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
59
    }
}
60

61
void
62
calcLatBins(GridSearchBins &searchBins)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
{
64
65
  const size_t nbins = searchBins.nbins;
  const double dlat = PI / nbins;  // lat interval for search bins
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66

67
  if (Options::cdoVerbose) cdoPrint("Using %zu latitude bins to restrict search.", nbins);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68

69
  if (nbins > 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
    {
71
      searchBins.bin_lats.resize(2 * nbins);
72
73
      for (size_t n = 0; n < nbins; ++n)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
          searchBins.bin_lats[n * 2] = (n) *dlat - PIH;
          searchBins.bin_lats[n * 2 + 1] = (n + 1) * dlat - PIH;
76
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77

78
      searchBins.bin_addr.resize(2 * nbins);
79
      calcBinAddr(searchBins);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82

83
size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
get_srch_cells(size_t tgt_cell_addr, GridSearchBins &tgtBins, GridSearchBins &srcBins, float *tgt_cell_bound_box, size_t *srch_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
{
86
87
  const size_t nbins = srcBins.nbins;
  const size_t src_grid_size = srcBins.ncells;
88
89
  const size_t *restrict bin_addr1 = &tgtBins.bin_addr[0];
  const size_t *restrict bin_addr2 = &srcBins.bin_addr[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
  // Restrict searches first using search bins
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

93
94
  size_t min_add = src_grid_size - 1;
  size_t max_add = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

96
  for (size_t n = 0; n < nbins; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
      size_t n2 = n << 1;
99
      if (tgt_cell_addr >= bin_addr1[n2] && tgt_cell_addr <= bin_addr1[n2 + 1])
100
101
102
103
        {
          if (bin_addr2[n2] < min_add) min_add = bin_addr2[n2];
          if (bin_addr2[n2 + 1] > max_add) max_add = bin_addr2[n2 + 1];
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
  // Further restrict searches using bounding boxes
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107

Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
111
  float bound_box_lat1 = tgt_cell_bound_box[0];
  float bound_box_lat2 = tgt_cell_bound_box[1];
  float bound_box_lon1 = tgt_cell_bound_box[2];
  float bound_box_lon2 = tgt_cell_bound_box[3];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
114
  const float *restrict src_cell_bound_box = &srcBins.cell_bound_box[0];
  size_t src_cell_addm4;
115
  size_t num_srch_cells = 0;
116
  for (size_t src_cell_add = min_add; src_cell_add <= max_add; ++src_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
    {
118
      src_cell_addm4 = src_cell_add << 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
      if ((src_cell_bound_box[src_cell_addm4 + 2] <= bound_box_lon2) && (src_cell_bound_box[src_cell_addm4 + 3] >= bound_box_lon1))
120
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
          if ((src_cell_bound_box[src_cell_addm4] <= bound_box_lat2) && (src_cell_bound_box[src_cell_addm4 + 1] >= bound_box_lat1))
122
123
124
125
126
            {
              srch_add[num_srch_cells] = src_cell_add;
              num_srch_cells++;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
  if (bound_box_lon1 < 0.0f || bound_box_lon2 > PI2_f)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
      if (bound_box_lon1 < 0.0f)
132
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
134
          bound_box_lon1 += PI2_f;
          bound_box_lon2 += PI2_f;
135
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
      else
137
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
          bound_box_lon1 -= PI2_f;
          bound_box_lon2 -= PI2_f;
140
141
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
      for (size_t src_cell_add = min_add; src_cell_add <= max_add; ++src_cell_add)
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
        {
          src_cell_addm4 = src_cell_add << 2;
          if ((src_cell_bound_box[src_cell_addm4 + 2] <= bound_box_lon2)
              && (src_cell_bound_box[src_cell_addm4 + 3] >= bound_box_lon1))
            {
              if ((src_cell_bound_box[src_cell_addm4] <= bound_box_lat2)
                  && (src_cell_bound_box[src_cell_addm4 + 1] >= bound_box_lat1))
                {
                  size_t ii;
                  for (ii = 0; ii < num_srch_cells; ++ii)
                    if (srch_add[ii] == src_cell_add) break;

                  if (ii == num_srch_cells)
                    {
                      srch_add[num_srch_cells] = src_cell_add;
                      num_srch_cells++;
                    }
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
164
    }

165
  return num_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167

168
static int
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
gridSearchSquareCurv2dNNScrip(size_t min_add, size_t max_add, size_t *restrict nbr_add, double *restrict nbr_dist, double plat,
                              double plon, const double *restrict src_center_lat, const double *restrict src_center_lon)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
172
{
  int search_result = 0;
173

174
175
176
177
  const double coslat_dst = std::cos(plat);
  const double sinlat_dst = std::sin(plat);
  const double coslon_dst = std::cos(plon);
  const double sinlon_dst = std::sin(plon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178

Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
180
181
  double dist_min = DBL_MAX;
  for (unsigned n = 0; n < 4; ++n) nbr_dist[n] = DBL_MAX;

182
  double distance;
183
  for (size_t srch_add = min_add; srch_add <= max_add; ++srch_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
    {
185
      distance = std::acos(coslat_dst * std::cos(src_center_lat[srch_add])
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
                               * (coslon_dst * std::cos(src_center_lon[srch_add]) + sinlon_dst * std::sin(src_center_lon[srch_add]))
                           + sinlat_dst * std::sin(src_center_lat[srch_add]));
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206

      if (distance < dist_min)
        {
          for (unsigned n = 0; n < 4; ++n)
            {
              if (distance < nbr_dist[n])
                {
                  for (unsigned i = 3; i > n; --i)
                    {
                      nbr_add[i] = nbr_add[i - 1];
                      nbr_dist[i] = nbr_dist[i - 1];
                    }
                  search_result = -1;
                  nbr_add[n] = srch_add;
                  nbr_dist[n] = distance;
                  dist_min = nbr_dist[3];
                  break;
                }
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
209
        }
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
  for (unsigned n = 0; n < 4; ++n) nbr_dist[n] = 1.0 / (nbr_dist[n] + TINY);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
  distance = 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
213
  for (unsigned n = 0; n < 4; ++n) distance += nbr_dist[n];
  for (unsigned n = 0; n < 4; ++n) nbr_dist[n] /= distance;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214

215
  return search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
217
}

218
static unsigned
219
quadCrossProducts(double plon, double plat, double *restrict lons, double *restrict lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220
221
222
223
224
225
226
{
  unsigned n;
  int scross[4], scross_last = 0;
  // Vectors for cross-product check
  double vec1_lat, vec1_lon;
  double vec2_lat, vec2_lon;

227
  // For consistency, we must make sure all lons are in same 2pi interval
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
  vec1_lon = lons[0] - plon;
229
230
231
232
  if (vec1_lon > PI)
    lons[0] -= PI2;
  else if (vec1_lon < -PI)
    lons[0] += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233

234
  for (n = 1; n < 4; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
    {
      vec1_lon = lons[n] - lons[0];
237
238
239
240
      if (vec1_lon > PI)
        lons[n] -= PI2;
      else if (vec1_lon < -PI)
        lons[n] += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
243
    }

  /* corner_loop */
244
  for (n = 0; n < 4; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
    {
246
      const unsigned next_n = (n + 1) % 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
      /*
248
        Here we take the cross product of the vector making up each box side
Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
        with the vector formed by the vertex and search point.
250
        If all the cross products are positive, the point is contained in the box.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
252
253
254
255
256
      */
      vec1_lat = lats[next_n] - lats[n];
      vec1_lon = lons[next_n] - lons[n];
      vec2_lat = plat - lats[n];
      vec2_lon = plon - lons[n];

257
      // Check for 0,2pi crossings
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
      if (vec1_lon > 3.0 * PIH)
259
        vec1_lon -= PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
      else if (vec1_lon < -3.0 * PIH)
261
        vec1_lon += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
262

Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
      if (vec2_lon > 3.0 * PIH)
264
        vec2_lon -= PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
      else if (vec2_lon < -3.0 * PIH)
266
        vec2_lon += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267

268
      double cross_product = vec1_lon * vec2_lat - vec2_lon * vec1_lat;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
270
271
272

      /* If cross product is less than ZERO, this cell doesn't work    */
      /* 2008-10-16 Uwe Schulzweida: bug fix for cross_product eq zero */
      scross[n] = cross_product < 0 ? -1 : cross_product > 0 ? 1 : 0;
273
      if (n == 0) scross_last = scross[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274
      if ((scross[n] < 0 && scross_last > 0) || (scross[n] > 0 && scross_last < 0)) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
276
277
      scross_last = scross[n];
    }

278
  if (n >= 4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
280
    {
      n = 0;
281
282
      if (scross[0] >= 0 && scross[1] >= 0 && scross[2] >= 0 && scross[3] >= 0)
        n = 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
      else if (scross[0] <= 0 && scross[1] <= 0 && scross[2] <= 0 && scross[3] <= 0)
284
        n = 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285
286
287
288
289
    }

  return n;
}

290
bool
Uwe Schulzweida's avatar
Uwe Schulzweida committed
291
292
pointInQuad(bool isCyclic, size_t nx, size_t ny, size_t i, size_t j, size_t adds[4], double lons[4], double lats[4], double plon,
            double plat, const double *restrict centerLon, const double *restrict centerLat)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293
294
{
  bool search_result = false;
295
  size_t ip1 = (i < (nx - 1)) ? i + 1 : isCyclic ? 0 : i;
296
  size_t jp1 = (j < (ny - 1)) ? j + 1 : j;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
297

298
  if (i == ip1 || j == jp1) return search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299

Uwe Schulzweida's avatar
Uwe Schulzweida committed
300
  size_t idx[4];
301
302
303
304
  idx[0] = j * nx + i;
  idx[1] = j * nx + ip1;    // east
  idx[2] = jp1 * nx + ip1;  // north-east
  idx[3] = jp1 * nx + i;    // north
Uwe Schulzweida's avatar
Uwe Schulzweida committed
305

Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
307
  for (unsigned k = 0; k < 4; ++k) lons[k] = centerLon[idx[k]];
  for (unsigned k = 0; k < 4; ++k) lats[k] = centerLat[idx[k]];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308

309
  unsigned n = quadCrossProducts(plon, plat, lons, lats);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310

311
  // If cross products all same sign, we found the location
312
  if (n >= 4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314
      for (unsigned k = 0; k < 4; ++k) adds[k] = idx[k];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
318
319
320
      search_result = true;
    }

  return search_result;
}

321
int
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322
323
gridSearchSquareCurv2dScrip(RemapGrid *src_grid, size_t *restrict src_add, double *restrict src_lats, double *restrict src_lons,
                            double plat, double plon, GridSearchBins &srcBins)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
{
  /*
    Output variables:

    int    src_add[4]              ! address of each corner point enclosing P
    double src_lats[4]             ! latitudes  of the four corner points
    double src_lons[4]             ! longitudes of the four corner points

    Input variables:

    double plat                    ! latitude  of the search point
    double plon                    ! longitude of the search point

    int src_grid_dims[2]           ! size of each src grid dimension

339
    double src_center_lat[]        ! latitude  of each src grid center
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
341
    double src_center_lon[]        ! longitude of each src grid center

Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
    float src_grid_bound_box[][4] ! bound box for source grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343

344
    int src_bin_addr[][2]           ! latitude bins for restricting
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
346
347
  */
  int search_result = 0;

348
  const size_t nbins = srcBins.nbins;
349
350
351
  const size_t *restrict src_bin_addr = &srcBins.bin_addr[0];
  const float *restrict bin_lats = &srcBins.bin_lats[0];
  const float *restrict src_grid_bound_box = &srcBins.cell_bound_box[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
352

353
354
  const float rlat = plat;
  const float rlon = plon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355

356
  // restrict search first using bins
Uwe Schulzweida's avatar
Uwe Schulzweida committed
357

Uwe Schulzweida's avatar
Uwe Schulzweida committed
358
  for (unsigned n = 0; n < 4; ++n) src_add[n] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359

360
  // addresses for restricting search
361
  size_t min_add = src_grid->size - 1;
362
  size_t max_add = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363

364
  for (size_t n = 0; n < nbins; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
    {
366
      const size_t n2 = n << 1;
367
368
      if (rlat >= bin_lats[n2] && rlat <= bin_lats[n2 + 1])
        {
369
370
          if (src_bin_addr[n2] < min_add) min_add = src_bin_addr[n2];
          if (src_bin_addr[n2 + 1] > max_add) max_add = src_bin_addr[n2 + 1];
371
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
    }
373

Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
375
  /* Now perform a more detailed search */

376
377
  const size_t nx = src_grid->dims[0];
  const size_t ny = src_grid->dims[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378

379
  for (size_t srch_add = min_add; srch_add <= max_add; ++srch_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
    {
381
      const size_t srch_add4 = srch_add << 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
382
      /* First check bounding box */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
384
      if (rlon >= src_grid_bound_box[srch_add4 + 2] && rlon <= src_grid_bound_box[srch_add4 + 3]
          && rlat >= src_grid_bound_box[srch_add4] && rlat <= src_grid_bound_box[srch_add4 + 1])
385
386
        {
          /* We are within bounding box so get really serious */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
388

          /* Determine neighbor addresses */
389
390
391
          size_t j = srch_add / nx;
          size_t i = srch_add - j * nx;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
393
          if (pointInQuad(src_grid->is_cyclic, nx, ny, i, j, src_add, src_lons, src_lats, plon, plat, src_grid->cell_center_lon,
                          src_grid->cell_center_lat))
394
395
396
397
            {
              search_result = 1;
              return search_result;
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398

399
          /* Otherwise move on to next cell */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
401

        } /* Bounding box check */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403
404

  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
405
406
407
    If no cell found, point is likely either in a box that straddles either pole or is outside the grid.
    Fall back to a distance-weighted average of the four closest points. Go ahead and compute weights here,
    but store in src_lats and return -add to prevent the parent routine from computing bilinear weights.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408
  */
409
  if (!src_grid->lextrapolate) return search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
411
412
413
414

  /*
    printf("Could not find location for %g %g\n", plat*RAD2DEG, plon*RAD2DEG);
    printf("Using nearest-neighbor average for this point\n");
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
416
  search_result = gridSearchSquareCurv2dNNScrip(min_add, max_add, src_add, src_lats, plat, plon, src_grid->cell_center_lat,
                                                src_grid->cell_center_lon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417

418
  return search_result;
419
}