1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <assert.h>
6 #include <string.h>
7
8 #include "dmemory.h"
9 #include "cdi.h"
10 #include "cdi_cksum.h"
11 #include "cdi_int.h"
12 #include "cdi_uuid.h"
13 #include "grid.h"
14 #include "resource_handle.h"
15 #include "resource_unpack.h"
16 #include "namespace.h"
17 #include "serialize.h"
18 #include "vlist.h"
19
20 int (*proj_lonlat_to_lcc_func)(struct CDI_GridProjParams gpp, size_t, double*, double*) = NULL;
21 int (*proj_lcc_to_lonlat_func)(struct CDI_GridProjParams *gpp, double, double, size_t, double*, double*) = NULL;
22 int (*proj_lonlat_to_stere_func)(struct CDI_GridProjParams gpp, size_t, double*, double*) = NULL;
23 int (*proj_stere_to_lonlat_func)(struct CDI_GridProjParams *gpp, double, double, size_t, double*, double*) = NULL;
24
25 /* the value in the second pair of brackets must match the length of
26 * the longest string (including terminating NUL) */
27 static const char Grids[][17] = {
28 /* 0 */ "undefined",
29 /* 1 */ "generic",
30 /* 2 */ "gaussian",
31 /* 3 */ "gaussian_reduced",
32 /* 4 */ "lonlat",
33 /* 5 */ "spectral",
34 /* 6 */ "fourier",
35 /* 7 */ "gme",
36 /* 8 */ "trajectory",
37 /* 9 */ "unstructured",
38 /* 10 */ "curvilinear",
39 /* 11 */ "lcc",
40 /* 12 */ "projection",
41 /* 13 */ "characterXY",
42 };
43
44 /* must match table below */
45 enum xystdname_idx {
46 grid_xystdname_grid_latlon,
47 grid_xystdname_latlon,
48 grid_xystdname_projection,
49 grid_xystdname_char,
50 };
51 static const char xystdname_tab[][2][24] = {
52 [grid_xystdname_grid_latlon] = { "grid_longitude",
53 "grid_latitude" },
54 [grid_xystdname_latlon] = { "longitude",
55 "latitude" },
56 [grid_xystdname_projection] = { "projection_x_coordinate",
57 "projection_y_coordinate" },
58 [grid_xystdname_char] = { "region",
59 "region" },
60 };
61
62
63
64 static int gridCompareP ( void * gridptr1, void * gridptr2 );
65 static void gridDestroyP ( void * gridptr );
66 static void gridPrintP ( void * gridptr, FILE * fp );
67 static int gridGetPackSize ( void * gridptr, void *context);
68 static void gridPack ( void * gridptr, void * buff, int size, int *position, void *context);
69 static int gridTxCode ( void );
70
71 static const resOps gridOps = {
72 gridCompareP,
73 gridDestroyP,
74 gridPrintP,
75 gridGetPackSize,
76 gridPack,
77 gridTxCode
78 };
79
80 static int GRID_Debug = 0; /* If set to 1, debugging */
81
82
grid_to_pointer(int gridID)83 grid_t *grid_to_pointer(int gridID)
84 {
85 return (grid_t *)reshGetVal(gridID, &gridOps);
86 }
87
88 #define gridMark4Update(gridID) reshSetStatus(gridID, &gridOps, RESH_DESYNC_IN_USE)
89
90 static
cdiInqAttConvertedToFloat(int gridID,int atttype,const char * attname,int attlen,double * attflt)91 bool cdiInqAttConvertedToFloat(int gridID, int atttype, const char *attname, int attlen, double *attflt)
92 {
93 bool status = true;
94
95 if (atttype == CDI_DATATYPE_INT32)
96 {
97 int attint;
98 int *pattint = attlen > 1 ? (int*) malloc(attlen*sizeof(int)) : &attint;
99 cdiInqAttInt(gridID, CDI_GLOBAL, attname, attlen, pattint);
100 for ( int i = 0; i < attlen; ++i ) attflt[i] = (double)pattint[i];
101 if (attlen > 1) free(pattint);
102 }
103 else if (atttype == CDI_DATATYPE_FLT32 || atttype == CDI_DATATYPE_FLT64)
104 {
105 cdiInqAttFlt(gridID, CDI_GLOBAL, attname, attlen, attflt);
106 }
107 else
108 {
109 status = false;
110 }
111
112 return status;
113 }
114
115 static
grid_axis_init(struct gridaxis_t * axisptr)116 void grid_axis_init(struct gridaxis_t *axisptr)
117 {
118 axisptr->size = 0;
119 axisptr->vals = NULL;
120 axisptr->bounds = NULL;
121 axisptr->flag = 0;
122 axisptr->first = 0.0;
123 axisptr->last = 0.0;
124 axisptr->inc = 0.0;
125 #ifndef USE_MPI
126 axisptr->clength = 0;
127 axisptr->cvals = NULL;
128 #endif
129 cdiInitKeys(&axisptr->keys);
130 }
131
grid_init(grid_t * gridptr)132 void grid_init(grid_t *gridptr)
133 {
134 gridptr->self = CDI_UNDEFID;
135 gridptr->type = CDI_UNDEFID;
136 gridptr->datatype = CDI_UNDEFID;
137 gridptr->proj = CDI_UNDEFID;
138 gridptr->projtype = CDI_UNDEFID;
139 gridptr->mask = NULL;
140 gridptr->mask_gme = NULL;
141 gridptr->size = 0;
142
143 grid_axis_init(&gridptr->x);
144 grid_axis_init(&gridptr->y);
145
146 gridptr->area = NULL;
147 gridptr->reducedPoints = NULL;
148 gridptr->reducedPointsSize = 0;
149
150 gridptr->gme.nd = 0;
151 gridptr->gme.ni = 0;
152 gridptr->gme.ni2 = 0;
153 gridptr->gme.ni3 = 0;
154
155 gridptr->trunc = 0;
156 gridptr->nvertex = 0;
157 gridptr->np = 0;
158 gridptr->isCyclic = CDI_UNDEFID;
159
160 gridptr->lcomplex = false;
161 gridptr->hasdims = true;
162 gridptr->name = NULL;
163 gridptr->vtable = &cdiGridVtable;
164
165 cdiInitKeys(&gridptr->keys);
166 gridptr->atts.nalloc = MAX_ATTRIBUTES;
167 gridptr->atts.nelems = 0;
168
169 cdiDefVarKeyInt(&gridptr->keys, CDI_KEY_DATATYPE, CDI_DATATYPE_FLT64);
170 cdiDefVarKeyInt(&gridptr->keys, CDI_KEY_SCANNINGMODE, 64);
171 }
172
173
174 static
grid_free_components(grid_t * gridptr)175 void grid_free_components(grid_t *gridptr)
176 {
177 void *p2free[] = { gridptr->mask, gridptr->mask_gme,
178 gridptr->x.vals, gridptr->y.vals,
179 #ifndef USE_MPI
180 gridptr->x.cvals, gridptr->y.cvals,
181 #endif
182 gridptr->x.bounds, gridptr->y.bounds,
183 gridptr->reducedPoints, gridptr->area,
184 gridptr->name};
185
186 for ( size_t i = 0; i < sizeof(p2free)/sizeof(p2free[0]); ++i )
187 if ( p2free[i] ) Free(p2free[i]);
188
189 /*
190 int gridID = gridptr->self;
191 cdiDeleteKeys(gridID, CDI_XAXIS);
192 cdiDeleteKeys(gridID, CDI_YAXIS);
193 cdiDeleteKeys(gridID, CDI_GLOBAL);
194 cdiDeleteAtts(gridID, CDI_GLOBAL);
195 */
196 }
197
grid_free(grid_t * gridptr)198 void grid_free(grid_t *gridptr)
199 {
200 grid_free_components(gridptr);
201 grid_init(gridptr);
202 }
203
204 static
gridNewEntry(cdiResH resH)205 grid_t *gridNewEntry(cdiResH resH)
206 {
207 grid_t *gridptr = (grid_t*) Malloc(sizeof(grid_t));
208 grid_init(gridptr);
209
210 if ( resH == CDI_UNDEFID )
211 gridptr->self = reshPut(gridptr, &gridOps);
212 else
213 {
214 gridptr->self = resH;
215 reshReplace(resH, gridptr, &gridOps);
216 }
217
218 return gridptr;
219 }
220
221 static
gridInit(void)222 void gridInit(void)
223 {
224 static bool gridInitialized = false;
225 if ( gridInitialized ) return;
226 gridInitialized = true;
227
228 const char *env = getenv("GRID_DEBUG");
229 if ( env ) GRID_Debug = atoi(env);
230 }
231
232 static
grid_copy_base_scalar_fields(grid_t * gridptrOrig,grid_t * gridptrDup)233 void grid_copy_base_scalar_fields(grid_t *gridptrOrig, grid_t *gridptrDup)
234 {
235 memcpy(gridptrDup, gridptrOrig, sizeof(grid_t));
236 gridptrDup->self = CDI_UNDEFID;
237 cdiInitKeys(&gridptrDup->keys);
238 cdiCopyVarKeys(&gridptrOrig->keys, &gridptrDup->keys);
239 cdiInitKeys(&gridptrDup->x.keys);
240 cdiCopyVarKeys(&gridptrOrig->x.keys, &gridptrDup->x.keys);
241 cdiInitKeys(&gridptrDup->y.keys);
242 cdiCopyVarKeys(&gridptrOrig->y.keys, &gridptrDup->y.keys);
243 }
244
245
246 static
grid_copy_base(grid_t * gridptrOrig)247 grid_t *grid_copy_base(grid_t *gridptrOrig)
248 {
249 grid_t *gridptrDup = (grid_t *)Malloc(sizeof (*gridptrDup));
250 gridptrOrig->vtable->copyScalarFields(gridptrOrig, gridptrDup);
251 gridptrOrig->vtable->copyArrayFields(gridptrOrig, gridptrDup);
252 return gridptrDup;
253 }
254
255
cdiGridCount(void)256 unsigned cdiGridCount(void)
257 {
258 return reshCountType(&gridOps);
259 }
260
261 static inline
gridaxisSetKey(struct gridaxis_t * axisptr,int key,const char * name)262 void gridaxisSetKey(struct gridaxis_t *axisptr, int key, const char *name)
263 {
264 if (find_key(&axisptr->keys, key) == NULL)
265 cdiDefVarKeyBytes(&axisptr->keys, key, (const unsigned char*)name, (int)strlen(name)+1);
266 }
267
cdiGridTypeInit(grid_t * gridptr,int gridtype,size_t size)268 void cdiGridTypeInit(grid_t *gridptr, int gridtype, size_t size)
269 {
270 gridptr->type = gridtype;
271 gridptr->size = size;
272
273 if ( gridtype == GRID_LONLAT ) gridptr->nvertex = 2;
274 else if ( gridtype == GRID_GAUSSIAN ) gridptr->nvertex = 2;
275 else if ( gridtype == GRID_GAUSSIAN_REDUCED ) gridptr->nvertex = 2;
276 else if ( gridtype == GRID_CURVILINEAR ) gridptr->nvertex = 4;
277 else if ( gridtype == GRID_UNSTRUCTURED ) gridptr->x.size = size;
278
279 switch (gridtype)
280 {
281 case GRID_LONLAT:
282 case GRID_GAUSSIAN:
283 case GRID_GAUSSIAN_REDUCED:
284 case GRID_TRAJECTORY:
285 case GRID_CURVILINEAR:
286 case GRID_UNSTRUCTURED:
287 case GRID_GME:
288 {
289 if ( gridtype == GRID_TRAJECTORY )
290 {
291 gridaxisSetKey(&gridptr->x, CDI_KEY_NAME, "tlon");
292 gridaxisSetKey(&gridptr->y, CDI_KEY_NAME, "tlat");
293 }
294 else
295 {
296 gridaxisSetKey(&gridptr->x, CDI_KEY_NAME, "lon");
297 gridaxisSetKey(&gridptr->y, CDI_KEY_NAME, "lat");
298 }
299
300 gridaxisSetKey(&gridptr->x, CDI_KEY_LONGNAME, "longitude");
301 gridaxisSetKey(&gridptr->y, CDI_KEY_LONGNAME, "latitude");
302
303 gridaxisSetKey(&gridptr->x, CDI_KEY_UNITS, "degrees_east");
304 gridaxisSetKey(&gridptr->y, CDI_KEY_UNITS, "degrees_north");
305
306 gridaxisSetKey(&gridptr->x, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_latlon][0]);
307 gridaxisSetKey(&gridptr->y, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_latlon][1]);
308
309 break;
310 }
311 #ifndef USE_MPI
312 case GRID_CHARXY:
313 {
314 if ( gridptr->x.cvals ) gridaxisSetKey(&gridptr->x, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_char][0]);
315 if ( gridptr->y.cvals ) gridaxisSetKey(&gridptr->y, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_char][1]);
316
317 break;
318 }
319 #endif
320 case GRID_GENERIC:
321 case GRID_PROJECTION:
322 {
323 gridaxisSetKey(&gridptr->x, CDI_KEY_NAME, "x");
324 gridaxisSetKey(&gridptr->y, CDI_KEY_NAME, "y");
325 if (gridtype == GRID_PROJECTION)
326 {
327 gridaxisSetKey(&gridptr->x, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_projection][0]);
328 gridaxisSetKey(&gridptr->y, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_projection][1]);
329 gridaxisSetKey(&gridptr->x, CDI_KEY_UNITS, "m");
330 gridaxisSetKey(&gridptr->y, CDI_KEY_UNITS, "m");
331 }
332 break;
333 }
334 }
335 }
336
337
338 // used also in CDO
gridGenXvals(int xsize,double xfirst,double xlast,double xinc,double * restrict xvals)339 void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *restrict xvals)
340 {
341 if (fabs(xinc) <= 0 && xsize > 1)
342 {
343 if (xfirst >= xlast)
344 {
345 while ( xfirst >= xlast ) xlast += 360;
346 xinc = (xlast-xfirst)/(xsize);
347 }
348 else
349 {
350 xinc = (xlast-xfirst)/(xsize-1);
351 }
352 }
353
354 for ( int i = 0; i < xsize; ++i )
355 xvals[i] = xfirst + i*xinc;
356 }
357
358 static
calc_gaussgrid(double * restrict yvals,int ysize,double yfirst,double ylast)359 void calc_gaussgrid(double *restrict yvals, int ysize, double yfirst, double ylast)
360 {
361 double *restrict yw = (double *) malloc((size_t)ysize * sizeof(double));
362 gaussianLatitudes((size_t)ysize, yvals, yw);
363 free(yw);
364 for (int i = 0; i < ysize; i++ )
365 yvals[i] = asin(yvals[i])/M_PI*180.0;
366
367 if (yfirst < ylast && yfirst > -90.0 && ylast < 90.0)
368 {
369 const int yhsize = ysize/2;
370 for (int i = 0; i < yhsize; i++ )
371 {
372 const double ytmp = yvals[i];
373 yvals[i] = yvals[ysize-i-1];
374 yvals[ysize-i-1] = ytmp;
375 }
376 }
377 }
378
379 static
gridGenYvalsGaussian(int ysize,double yfirst,double ylast,double * restrict yvals)380 void gridGenYvalsGaussian(int ysize, double yfirst, double ylast, double *restrict yvals)
381 {
382 const double deleps = 0.002;
383
384 calc_gaussgrid(yvals, ysize, yfirst, ylast);
385
386 if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
387 if ( fabs(yvals[0] - yfirst) > deleps || fabs(yvals[ysize-1] - ylast) > deleps )
388 {
389 bool lfound = false;
390 int ny = (int) (180./(fabs(ylast-yfirst)/(ysize-1)) + 0.5);
391 ny -= ny%2;
392 if ( ny > ysize && ny < 4096 )
393 {
394 double *ytmp = (double *) Malloc((size_t)ny * sizeof (double));
395 calc_gaussgrid(ytmp, ny, yfirst, ylast);
396
397 int i;
398 for ( i = 0; i < (ny-ysize); i++ )
399 if ( fabs(ytmp[i] - yfirst) < deleps ) break;
400 int nstart = i;
401
402 lfound = (nstart+ysize-1) < ny && fabs(ytmp[nstart+ysize-1] - ylast) < deleps;
403 if ( lfound )
404 {
405 for (i = 0; i < ysize; i++) yvals[i] = ytmp[i+nstart];
406 }
407
408 if ( ytmp ) Free(ytmp);
409 }
410
411 if ( !lfound )
412 {
413 Warning("Cannot calculate gaussian latitudes for lat1 = %g latn = %g!", yfirst, ylast);
414 for (int i = 0; i < ysize; i++ ) yvals[i] = 0;
415 yvals[0] = yfirst;
416 yvals[ysize-1] = ylast;
417 }
418 }
419 }
420
421 static
gridGenYvalsRegular(int ysize,double yfirst,double ylast,double yinc,double * restrict yvals)422 void gridGenYvalsRegular(int ysize, double yfirst, double ylast, double yinc, double *restrict yvals)
423 {
424 if (fabs(yinc) <= 0 && ysize > 1)
425 {
426 if ( IS_EQUAL(yfirst, ylast) && IS_NOT_EQUAL(yfirst, 0) ) ylast *= -1;
427
428 if ( yfirst > ylast )
429 yinc = (yfirst-ylast)/(ysize-1);
430 else if ( yfirst < ylast )
431 yinc = (ylast-yfirst)/(ysize-1);
432 else
433 {
434 if ( ysize%2 != 0 )
435 {
436 yinc = 180.0/(ysize-1);
437 yfirst = -90;
438 }
439 else
440 {
441 yinc = 180.0/ysize;
442 yfirst = -90 + yinc/2;
443 }
444 }
445 }
446
447 if ( yfirst > ylast && yinc > 0 ) yinc = -yinc;
448
449 for (int i = 0; i < ysize; i++ )
450 yvals[i] = yfirst + i*yinc;
451 }
452
453 // used also in CDO
gridGenYvals(int gridtype,int ysize,double yfirst,double ylast,double yinc,double * restrict yvals)454 void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *restrict yvals)
455 {
456 if (gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED)
457 {
458 if (ysize > 2)
459 {
460 gridGenYvalsGaussian(ysize, yfirst, ylast, yvals);
461 }
462 else
463 {
464 yvals[0] = yfirst;
465 yvals[ysize-1] = ylast;
466 }
467 }
468 // else if (gridtype == GRID_LONLAT || gridtype == GRID_GENERIC)
469 else
470 {
471 gridGenYvalsRegular(ysize, yfirst, ylast, yinc, yvals);
472 }
473 /*
474 else
475 Error("unable to calculate values for %s grid!", gridNamePtr(gridtype));
476 */
477 }
478
479 /*
480 @Function gridCreate
481 @Title Create a horizontal Grid
482
483 @Prototype int gridCreate(int gridtype, SizeType size)
484 @Parameter
485 @Item gridtype The type of the grid, one of the set of predefined CDI grid types.
486 The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_LONLAT},
487 @func{GRID_GAUSSIAN}, @func{GRID_PROJECTION}, @func{GRID_SPECTRAL},
488 @func{GRID_GME}, @func{GRID_CURVILINEAR} and @func{GRID_UNSTRUCTURED}.
489 @Item size Number of gridpoints.
490
491 @Description
492 The function @func{gridCreate} creates a horizontal Grid.
493
494 @Result
495 @func{gridCreate} returns an identifier to the Grid.
496
497 @Example
498 Here is an example using @func{gridCreate} to create a regular lon/lat Grid:
499
500 @Source
501 #include "cdi.h"
502 ...
503 #define nlon 12
504 #define nlat 6
505 ...
506 double lons[nlon] = {0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330};
507 double lats[nlat] = {-75, -45, -15, 15, 45, 75};
508 int gridID;
509 ...
510 gridID = gridCreate(GRID_LONLAT, nlon*nlat);
511 gridDefXsize(gridID, nlon);
512 gridDefYsize(gridID, nlat);
513 gridDefXvals(gridID, lons);
514 gridDefYvals(gridID, lats);
515 ...
516 @EndSource
517 @EndFunction
518 */
gridCreate(int gridtype,SizeType size)519 int gridCreate(int gridtype, SizeType size)
520 {
521 if (CDI_Debug) Message("gridtype=%s size=%zu", gridNamePtr(gridtype), size);
522
523 xassert(size);
524 gridInit();
525
526 grid_t *gridptr = gridNewEntry(CDI_UNDEFID);
527 if (! gridptr) Error("No memory");
528
529 const int gridID = gridptr->self;
530
531 if (CDI_Debug) Message("gridID: %d", gridID);
532
533 cdiGridTypeInit(gridptr, gridtype, (size_t)size);
534
535 return gridID;
536 }
537
538 static
gridDestroyKernel(grid_t * gridptr)539 void gridDestroyKernel(grid_t * gridptr)
540 {
541 xassert(gridptr);
542
543 const int id = gridptr->self;
544
545 grid_free_components(gridptr);
546 Free(gridptr);
547
548 reshRemove(id, &gridOps);
549 }
550
551 /*
552 @Function gridDestroy
553 @Title Destroy a horizontal Grid
554
555 @Prototype void gridDestroy(int gridID)
556 @Parameter
557 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
558
559 @EndFunction
560 */
gridDestroy(int gridID)561 void gridDestroy(int gridID)
562 {
563 //cdiDeleteKeys(gridID, CDI_GLOBAL);
564 //cdiDeleteAtts(gridID, CDI_GLOBAL);
565
566 grid_t *gridptr = grid_to_pointer(gridID);
567 gridptr->vtable->destroy(gridptr);
568 }
569
570 static
gridDestroyP(void * gridptr)571 void gridDestroyP(void * gridptr)
572 {
573 ((grid_t *)gridptr)->vtable->destroy((grid_t *)gridptr);
574 }
575
576
gridNamePtr(int gridtype)577 const char *gridNamePtr(int gridtype)
578 {
579 const int size = (int) (sizeof(Grids)/sizeof(Grids[0]));
580
581 const char *name = (gridtype >= 0 && gridtype < size) ? Grids[gridtype] : Grids[GRID_GENERIC];
582
583 return name;
584 }
585
586
gridName(int gridtype,char * gridname)587 void gridName(int gridtype, char *gridname)
588 {
589 strcpy(gridname, gridNamePtr(gridtype));
590 }
591
592 /*
593 @Function gridDefXname
594 @Title Define the name of a X-axis
595
596 @Prototype void gridDefXname(int gridID, const char *name)
597 @Parameter
598 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
599 @Item name Name of the X-axis.
600
601 @Description
602 The function @func{gridDefXname} defines the name of a X-axis.
603
604 @EndFunction
605 */
gridDefXname(int gridID,const char * name)606 void gridDefXname(int gridID, const char *name)
607 {
608 (void)cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, name);
609 }
610
611 /*
612 @Function gridDefXlongname
613 @Title Define the longname of a X-axis
614
615 @Prototype void gridDefXlongname(int gridID, const char *longname)
616 @Parameter
617 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
618 @Item longname Longname of the X-axis.
619
620 @Description
621 The function @func{gridDefXlongname} defines the longname of a X-axis.
622
623 @EndFunction
624 */
gridDefXlongname(int gridID,const char * longname)625 void gridDefXlongname(int gridID, const char *longname)
626 {
627 (void)cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_LONGNAME, longname);
628 }
629
630 /*
631 @Function gridDefXunits
632 @Title Define the units of a X-axis
633
634 @Prototype void gridDefXunits(int gridID, const char *units)
635 @Parameter
636 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
637 @Item units Units of the X-axis.
638
639 @Description
640 The function @func{gridDefXunits} defines the units of a X-axis.
641
642 @EndFunction
643 */
gridDefXunits(int gridID,const char * units)644 void gridDefXunits(int gridID, const char *units)
645 {
646 (void)cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, units);
647 }
648
649 /*
650 @Function gridDefYname
651 @Title Define the name of a Y-axis
652
653 @Prototype void gridDefYname(int gridID, const char *name)
654 @Parameter
655 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
656 @Item name Name of the Y-axis.
657
658 @Description
659 The function @func{gridDefYname} defines the name of a Y-axis.
660
661 @EndFunction
662 */
gridDefYname(int gridID,const char * name)663 void gridDefYname(int gridID, const char *name)
664 {
665 (void)cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, name);
666 }
667
668 /*
669 @Function gridDefYlongname
670 @Title Define the longname of a Y-axis
671
672 @Prototype void gridDefYlongname(int gridID, const char *longname)
673 @Parameter
674 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
675 @Item longname Longname of the Y-axis.
676
677 @Description
678 The function @func{gridDefYlongname} defines the longname of a Y-axis.
679
680 @EndFunction
681 */
gridDefYlongname(int gridID,const char * longname)682 void gridDefYlongname(int gridID, const char *longname)
683 {
684 (void)cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_LONGNAME, longname);
685 }
686
687 /*
688 @Function gridDefYunits
689 @Title Define the units of a Y-axis
690
691 @Prototype void gridDefYunits(int gridID, const char *units)
692 @Parameter
693 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
694 @Item units Units of the Y-axis.
695
696 @Description
697 The function @func{gridDefYunits} defines the units of a Y-axis.
698
699 @EndFunction
700 */
gridDefYunits(int gridID,const char * units)701 void gridDefYunits(int gridID, const char *units)
702 {
703 (void)cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, units);
704 }
705
706 /*
707 @Function gridInqXname
708 @Title Get the name of a X-axis
709
710 @Prototype void gridInqXname(int gridID, char *name)
711 @Parameter
712 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
713 @Item name Name of the X-axis. The caller must allocate space for the
714 returned string. The maximum possible length, in characters, of
715 the string is given by the predefined constant @func{CDI_MAX_NAME}.
716
717 @Description
718 The function @func{gridInqXname} returns the name of a X-axis.
719
720 @Result
721 @func{gridInqXname} returns the name of the X-axis to the parameter name.
722
723 @EndFunction
724 */
gridInqXname(int gridID,char * name)725 void gridInqXname(int gridID, char *name)
726 {
727 int length = CDI_MAX_NAME;
728 (void)cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, name, &length);
729 }
730
731 /*
732 @Function gridInqXlongname
733 @Title Get the longname of a X-axis
734
735 @Prototype void gridInqXlongname(int gridID, char *longname)
736 @Parameter
737 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
738 @Item longname Longname of the X-axis. The caller must allocate space for the
739 returned string. The maximum possible length, in characters, of
740 the string is given by the predefined constant @func{CDI_MAX_NAME}.
741
742 @Description
743 The function @func{gridInqXlongname} returns the longname of a X-axis.
744
745 @Result
746 @func{gridInqXlongname} returns the longname of the X-axis to the parameter longname.
747
748 @EndFunction
749 */
gridInqXlongname(int gridID,char * longname)750 void gridInqXlongname(int gridID, char *longname)
751 {
752 int length = CDI_MAX_NAME;
753 (void)cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_LONGNAME, longname, &length);
754 }
755
756 /*
757 @Function gridInqXunits
758 @Title Get the units of a X-axis
759
760 @Prototype void gridInqXunits(int gridID, char *units)
761 @Parameter
762 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
763 @Item units Units of the X-axis. The caller must allocate space for the
764 returned string. The maximum possible length, in characters, of
765 the string is given by the predefined constant @func{CDI_MAX_NAME}.
766
767 @Description
768 The function @func{gridInqXunits} returns the units of a X-axis.
769
770 @Result
771 @func{gridInqXunits} returns the units of the X-axis to the parameter units.
772
773 @EndFunction
774 */
gridInqXunits(int gridID,char * units)775 void gridInqXunits(int gridID, char *units)
776 {
777 int length = CDI_MAX_NAME;
778 (void)cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, units, &length);
779 }
780
781 /*
782 @Function gridInqYname
783 @Title Get the name of a Y-axis
784
785 @Prototype void gridInqYname(int gridID, char *name)
786 @Parameter
787 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
788 @Item name Name of the Y-axis. The caller must allocate space for the
789 returned string. The maximum possible length, in characters, of
790 the string is given by the predefined constant @func{CDI_MAX_NAME}.
791
792 @Description
793 The function @func{gridInqYname} returns the name of a Y-axis.
794
795 @Result
796 @func{gridInqYname} returns the name of the Y-axis to the parameter name.
797
798 @EndFunction
799 */
gridInqYname(int gridID,char * name)800 void gridInqYname(int gridID, char *name)
801 {
802 int length = CDI_MAX_NAME;
803 (void)cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, name, &length);
804 }
805
806 /*
807 @Function gridInqYlongname
808 @Title Get the longname of a Y-axis
809
810 @Prototype void gridInqYlongname(int gridID, char *longname)
811 @Parameter
812 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
813 @Item longname Longname of the Y-axis. The caller must allocate space for the
814 returned string. The maximum possible length, in characters, of
815 the string is given by the predefined constant @func{CDI_MAX_NAME}.
816
817 @Description
818 The function @func{gridInqYlongname} returns the longname of a Y-axis.
819
820 @Result
821 @func{gridInqYlongname} returns the longname of the Y-axis to the parameter longname.
822
823 @EndFunction
824 */
gridInqYlongname(int gridID,char * longname)825 void gridInqYlongname(int gridID, char *longname)
826 {
827 int length = CDI_MAX_NAME;
828 (void)cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_LONGNAME, longname, &length);
829 }
830
831 /*
832 @Function gridInqYunits
833 @Title Get the units of a Y-axis
834
835 @Prototype void gridInqYunits(int gridID, char *units)
836 @Parameter
837 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
838 @Item units Units of the Y-axis. The caller must allocate space for the
839 returned string. The maximum possible length, in characters, of
840 the string is given by the predefined constant @func{CDI_MAX_NAME}.
841
842 @Description
843 The function @func{gridInqYunits} returns the units of a Y-axis.
844
845 @Result
846 @func{gridInqYunits} returns the units of the Y-axis to the parameter units.
847
848 @EndFunction
849 */
gridInqYunits(int gridID,char * units)850 void gridInqYunits(int gridID, char *units)
851 {
852 int length = CDI_MAX_NAME;
853 (void)cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, units, &length);
854 }
855
856
gridDefProj(int gridID,int projID)857 void gridDefProj(int gridID, int projID)
858 {
859 grid_t *gridptr = grid_to_pointer(gridID);
860 gridptr->proj = projID;
861
862 if (gridptr->type == GRID_CURVILINEAR)
863 {
864 grid_t *projptr = grid_to_pointer(projID);
865 const char *xdimname = cdiInqVarKeyStringPtr(&gridptr->x.keys, CDI_KEY_DIMNAME);
866 const char *ydimname = cdiInqVarKeyStringPtr(&gridptr->y.keys, CDI_KEY_DIMNAME);
867 if (xdimname && find_key(&projptr->x.keys, CDI_KEY_NAME)) cdiDefKeyString(projID, CDI_XAXIS, CDI_KEY_NAME, xdimname);
868 if (ydimname && find_key(&projptr->y.keys, CDI_KEY_NAME)) cdiDefKeyString(projID, CDI_YAXIS, CDI_KEY_NAME, ydimname);
869 }
870 }
871
872
gridInqProj(int gridID)873 int gridInqProj(int gridID)
874 {
875 grid_t *gridptr = grid_to_pointer(gridID);
876 return gridptr->proj;
877 }
878
879
gridInqProjType(int gridID)880 int gridInqProjType(int gridID)
881 {
882 grid_t *gridptr = grid_to_pointer(gridID);
883
884 int projtype = gridptr->projtype;
885 if (projtype == -1)
886 {
887 char gmapname[CDI_MAX_NAME];
888 int length = CDI_MAX_NAME;
889 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname, &length);
890 if (gmapname[0])
891 {
892 // clang-format off
893 if (strIsEqual(gmapname, "rotated_latitude_longitude")) projtype = CDI_PROJ_RLL;
894 else if (strIsEqual(gmapname, "lambert_azimuthal_equal_area")) projtype = CDI_PROJ_LAEA;
895 else if (strIsEqual(gmapname, "lambert_conformal_conic")) projtype = CDI_PROJ_LCC;
896 else if (strIsEqual(gmapname, "sinusoidal")) projtype = CDI_PROJ_SINU;
897 else if (strIsEqual(gmapname, "polar_stereographic")) projtype = CDI_PROJ_STERE;
898 // clang-format on
899 gridptr->projtype = projtype;
900 }
901 }
902
903 return projtype;
904 }
905
906
gridVerifyProj(int gridID)907 void gridVerifyProj(int gridID)
908 {
909 grid_t *gridptr = grid_to_pointer(gridID);
910
911 const int projtype = gridInqProjType(gridID);
912 if (projtype == CDI_PROJ_RLL)
913 {
914 gridaxisSetKey(&gridptr->x, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_grid_latlon][0]);
915 gridaxisSetKey(&gridptr->y, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_grid_latlon][1]);
916 gridaxisSetKey(&gridptr->x, CDI_KEY_UNITS, "degrees");
917 gridaxisSetKey(&gridptr->y, CDI_KEY_UNITS, "degrees");
918 }
919 else if (projtype == CDI_PROJ_LCC)
920 {
921 gridaxisSetKey(&gridptr->x, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_projection][0]);
922 gridaxisSetKey(&gridptr->y, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_projection][1]);
923 gridaxisSetKey(&gridptr->x, CDI_KEY_UNITS, "m");
924 gridaxisSetKey(&gridptr->y, CDI_KEY_UNITS, "m");
925 }
926 }
927
928 /*
929 @Function gridInqType
930 @Title Get the type of a Grid
931
932 @Prototype int gridInqType(int gridID)
933 @Parameter
934 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
935
936 @Description
937 The function @func{gridInqType} returns the type of a Grid.
938
939 @Result
940 @func{gridInqType} returns the type of the grid,
941 one of the set of predefined CDI grid types.
942 The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_LONLAT},
943 @func{GRID_GAUSSIAN}, @func{GRID_PROJECTION}, @func{GRID_SPECTRAL}, @func{GRID_GME},
944 @func{GRID_CURVILINEAR} and @func{GRID_UNSTRUCTURED}.
945
946 @EndFunction
947 */
gridInqType(int gridID)948 int gridInqType(int gridID)
949 {
950 grid_t *gridptr = grid_to_pointer(gridID);
951 return gridptr->type;
952 }
953
954
955 /*
956 @Function gridInqSize
957 @Title Get the size of a Grid
958
959 @Prototype SizeType gridInqSize(int gridID)
960 @Parameter
961 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
962
963 @Description
964 The function @func{gridInqSize} returns the size of a Grid.
965
966 @Result
967 @func{gridInqSize} returns the number of grid points of a Grid.
968
969 @EndFunction
970 */
gridInqSize(int gridID)971 SizeType gridInqSize(int gridID)
972 {
973 grid_t *gridptr = grid_to_pointer(gridID);
974
975 size_t size = gridptr->size;
976 if ( size == 0 )
977 {
978 size_t xsize = gridptr->x.size;
979 size_t ysize = gridptr->y.size;
980
981 size = ysize ? xsize * ysize : xsize;
982
983 gridptr->size = size;
984 }
985
986 return (SizeType)size;
987 }
988
989 static
nsp2trunc(int nsp)990 int nsp2trunc(int nsp)
991 {
992 /* nsp = (trunc+1)*(trunc+1) */
993 /* => trunc^2 + 3*trunc - (x-2) = 0 */
994 /* */
995 /* with: y^2 + p*y + q = 0 */
996 /* y = -p/2 +- sqrt((p/2)^2 - q) */
997 /* p = 3 and q = - (x-2) */
998 int trunc = (int) (sqrt(nsp*4 + 1.) - 3) / 2;
999 return trunc;
1000 }
1001
1002
gridInqTrunc(int gridID)1003 int gridInqTrunc(int gridID)
1004 {
1005 grid_t *gridptr = grid_to_pointer(gridID);
1006
1007 if ( gridptr->trunc == 0 )
1008 {
1009 if ( gridptr->type == GRID_SPECTRAL )
1010 gridptr->trunc = nsp2trunc(gridptr->size);
1011 /*
1012 else if ( gridptr->type == GRID_GAUSSIAN )
1013 gridptr->trunc = nlat2trunc(gridptr->y.size);
1014 */
1015 }
1016
1017 return gridptr->trunc;
1018 }
1019
1020
gridDefTrunc(int gridID,int trunc)1021 void gridDefTrunc(int gridID, int trunc)
1022 {
1023 grid_t *gridptr = grid_to_pointer(gridID);
1024
1025 if ( gridptr->trunc != trunc )
1026 {
1027 gridMark4Update(gridID);
1028 gridptr->trunc = trunc;
1029 }
1030 }
1031
1032 /*
1033 @Function gridDefXsize
1034 @Title Define the number of values of a X-axis
1035
1036 @Prototype void gridDefXsize(int gridID, SizeType xsize)
1037 @Parameter
1038 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
1039 @Item xsize Number of values of a X-axis.
1040
1041 @Description
1042 The function @func{gridDefXsize} defines the number of values of a X-axis.
1043
1044 @EndFunction
1045 */
gridDefXsize(int gridID,SizeType xsize)1046 void gridDefXsize(int gridID, SizeType xsize)
1047 {
1048 grid_t *gridptr = grid_to_pointer(gridID);
1049
1050 const size_t gridSize = gridInqSize(gridID);
1051 if ((size_t)xsize > gridSize)
1052 Error("xsize %zu is greater then gridsize %zu", (size_t)xsize, gridSize);
1053
1054 const int gridType = gridInqType(gridID);
1055 if (gridType == GRID_UNSTRUCTURED && (size_t)xsize != gridSize)
1056 Error("xsize %zu must be equal to gridsize %zu for gridtype: %s", (size_t)xsize, gridSize, gridNamePtr(gridType));
1057 if (gridType == GRID_GAUSSIAN_REDUCED && xsize != 2 && (size_t)xsize != gridSize)
1058 Error("xsize %zu must be equal to gridsize %zu for gridtype: %s", (size_t)xsize, gridSize, gridNamePtr(gridType));
1059
1060 if (gridptr->x.size != (size_t)xsize)
1061 {
1062 gridMark4Update(gridID);
1063 gridptr->x.size = (size_t)xsize;
1064 }
1065
1066 if (gridType != GRID_UNSTRUCTURED && gridType != GRID_GAUSSIAN_REDUCED && gridType != GRID_PROJECTION)
1067 {
1068 const size_t axisproduct = gridptr->x.size*gridptr->y.size;
1069 if (axisproduct > 0 && axisproduct != gridSize)
1070 Error("Inconsistent grid declaration! (xsize=%zu ysize=%zu gridsize=%zu)",
1071 gridptr->x.size, gridptr->y.size, gridSize);
1072 }
1073 }
1074
gridDefDatatype(int gridID,int datatype)1075 void gridDefDatatype(int gridID, int datatype)
1076 {
1077 cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
1078 }
1079
gridInqDatatype(int gridID)1080 int gridInqDatatype(int gridID)
1081 {
1082 int datatype = CDI_UNDEFID;
1083 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
1084 return datatype;
1085 }
1086
1087 /*
1088 @Function gridInqXsize
1089 @Title Get the number of values of a X-axis
1090
1091 @Prototype SizeType gridInqXsize(int gridID)
1092 @Parameter
1093 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
1094
1095 @Description
1096 The function @func{gridInqXsize} returns the number of values of a X-axis.
1097
1098 @Result
1099 @func{gridInqXsize} returns the number of values of a X-axis.
1100
1101 @EndFunction
1102 */
gridInqXsize(int gridID)1103 SizeType gridInqXsize(int gridID)
1104 {
1105 grid_t *gridptr = grid_to_pointer(gridID);
1106 return (SizeType)gridptr->x.size;
1107 }
1108
1109 /*
1110 @Function gridDefYsize
1111 @Title Define the number of values of a Y-axis
1112
1113 @Prototype void gridDefYsize(int gridID, SizeType ysize)
1114 @Parameter
1115 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
1116 @Item ysize Number of values of a Y-axis.
1117
1118 @Description
1119 The function @func{gridDefYsize} defines the number of values of a Y-axis.
1120
1121 @EndFunction
1122 */
gridDefYsize(int gridID,SizeType ysize)1123 void gridDefYsize(int gridID, SizeType ysize)
1124 {
1125 grid_t *gridptr = grid_to_pointer(gridID);
1126
1127 const size_t gridSize = gridInqSize(gridID);
1128
1129 if ((size_t)ysize > gridSize)
1130 Error("ysize %zu is greater then gridsize %zu", (size_t)ysize, gridSize);
1131
1132 const int gridType = gridInqType(gridID);
1133 if (gridType == GRID_UNSTRUCTURED && (size_t)ysize != gridSize)
1134 Error("ysize %zu must be equal gridsize %zu for gridtype: %s", gridNamePtr(gridType), (size_t)ysize, gridSize);
1135
1136 if (gridptr->y.size != (size_t)ysize)
1137 {
1138 gridMark4Update(gridID);
1139 gridptr->y.size = (size_t)ysize;
1140 }
1141
1142 if (gridType != GRID_UNSTRUCTURED && gridType != GRID_GAUSSIAN_REDUCED && gridType != GRID_PROJECTION)
1143 {
1144 const size_t axisproduct = gridptr->x.size*gridptr->y.size;
1145 if (axisproduct > 0 && axisproduct != gridSize)
1146 Error("Inconsistent grid declaration! (xsize=%zu ysize=%zu gridsize=%zu)",
1147 gridptr->x.size, gridptr->y.size, gridSize);
1148 }
1149 }
1150
1151 /*
1152 @Function gridInqYsize
1153 @Title Get the number of values of a Y-axis
1154
1155 @Prototype SizeType gridInqYsize(int gridID)
1156 @Parameter
1157 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
1158
1159 @Description
1160 The function @func{gridInqYsize} returns the number of values of a Y-axis.
1161
1162 @Result
1163 @func{gridInqYsize} returns the number of values of a Y-axis.
1164
1165 @EndFunction
1166 */
gridInqYsize(int gridID)1167 SizeType gridInqYsize(int gridID)
1168 {
1169 grid_t *gridptr = grid_to_pointer(gridID);
1170 return (SizeType)gridptr->y.size;
1171 }
1172
1173 /*
1174 @Function gridDefNP
1175 @Title Define the number of parallels between a pole and the equator
1176
1177 @Prototype void gridDefNP(int gridID, int np)
1178 @Parameter
1179 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
1180 @Item np Number of parallels between a pole and the equator.
1181
1182 @Description
1183 The function @func{gridDefNP} defines the number of parallels between a pole and the equator
1184 of a Gaussian grid.
1185
1186 @EndFunction
1187 */
gridDefNP(int gridID,int np)1188 void gridDefNP(int gridID, int np)
1189 {
1190 grid_t *gridptr = grid_to_pointer(gridID);
1191
1192 if ( gridptr->np != np )
1193 {
1194 gridMark4Update(gridID);
1195 gridptr->np = np;
1196 }
1197 }
1198
1199 /*
1200 @Function gridInqNP
1201 @Title Get the number of parallels between a pole and the equator
1202
1203 @Prototype int gridInqNP(int gridID)
1204 @Parameter
1205 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
1206
1207 @Description
1208 The function @func{gridInqNP} returns the number of parallels between a pole and the equator
1209 of a Gaussian grid.
1210
1211 @Result
1212 @func{gridInqNP} returns the number of parallels between a pole and the equator.
1213
1214 @EndFunction
1215 */
gridInqNP(int gridID)1216 int gridInqNP(int gridID)
1217 {
1218 grid_t *gridptr = grid_to_pointer(gridID);
1219 return gridptr->np;
1220 }
1221
1222 /*
1223 @Function
1224 @Title
1225
1226 @Prototype
1227 @Parameter
1228 @Item Grid identifier
1229
1230 @EndFunction
1231 */
gridDefReducedPoints(int gridID,int reducedPointsSize,const int reducedPoints[])1232 void gridDefReducedPoints(int gridID, int reducedPointsSize, const int reducedPoints[])
1233 {
1234 grid_t *gridptr = grid_to_pointer(gridID);
1235
1236 gridptr->reducedPoints = (int *) Malloc((size_t)reducedPointsSize * sizeof(int));
1237 gridptr->reducedPointsSize = reducedPointsSize;
1238 memcpy(gridptr->reducedPoints, reducedPoints, (size_t)reducedPointsSize * sizeof(int));
1239 gridMark4Update(gridID);
1240 }
1241
1242 /*
1243 @Function
1244 @Title
1245
1246 @Prototype
1247 @Parameter
1248 @Item Grid identifier
1249
1250 @EndFunction
1251 */
gridInqReducedPoints(int gridID,int * reducedPoints)1252 void gridInqReducedPoints(int gridID, int *reducedPoints)
1253 {
1254 grid_t *gridptr = grid_to_pointer(gridID);
1255
1256 if ( gridptr->reducedPoints == 0 ) Error("undefined pointer!");
1257
1258 memcpy(reducedPoints, gridptr->reducedPoints, (size_t)gridptr->reducedPointsSize * sizeof(int));
1259 }
1260
1261 static size_t
gridInqMaskSerialGeneric(grid_t * gridptr,mask_t ** internalMask,int * restrict mask)1262 gridInqMaskSerialGeneric(grid_t *gridptr, mask_t **internalMask, int *restrict mask)
1263 {
1264 size_t size = gridptr->size;
1265
1266 if (CDI_Debug && size == 0)
1267 Warning("Size undefined for gridID = %d", gridptr->self);
1268
1269 const mask_t *restrict mask_src = *internalMask;
1270 if (mask_src)
1271 {
1272 if (mask && size > 0)
1273 for (size_t i = 0; i < size; ++i)
1274 mask[i] = (int)mask_src[i];
1275 }
1276 else
1277 size = 0;
1278
1279 return size;
1280 }
1281
1282 static SizeType
gridInqMaskSerial(grid_t * gridptr,int * mask)1283 gridInqMaskSerial(grid_t *gridptr, int *mask)
1284 {
1285 return (SizeType)gridInqMaskSerialGeneric(gridptr, &gridptr->mask, mask);
1286 }
1287
1288
gridInqMask(int gridID,int * mask)1289 int gridInqMask(int gridID, int *mask)
1290 {
1291 grid_t *gridptr = grid_to_pointer(gridID);
1292 return gridptr->vtable->inqMask(gridptr, mask);
1293 }
1294
1295 static void
gridDefMaskSerial(grid_t * gridptr,const int * mask)1296 gridDefMaskSerial(grid_t *gridptr, const int *mask)
1297 {
1298 const size_t size = gridptr->size;
1299 if (size == 0) Error("Size undefined for gridID = %d", gridptr->self);
1300
1301 if (mask == NULL)
1302 {
1303 if (gridptr->mask)
1304 {
1305 Free(gridptr->mask);
1306 gridptr->mask = NULL;
1307 }
1308 }
1309 else
1310 {
1311 if (gridptr->mask == NULL)
1312 gridptr->mask = (mask_t *) Malloc(size*sizeof(mask_t));
1313 else if (CDI_Debug)
1314 Warning("grid mask already defined!");
1315
1316 for (size_t i = 0; i < size; ++i )
1317 gridptr->mask[i] = (mask_t)(mask[i] != 0);
1318 }
1319 }
1320
gridDefMask(int gridID,const int * mask)1321 void gridDefMask(int gridID, const int *mask)
1322 {
1323 grid_t *gridptr = grid_to_pointer(gridID);
1324 gridptr->vtable->defMask(gridptr, mask);
1325 gridMark4Update(gridID);
1326 }
1327
1328 static int
gridInqMaskGMESerial(grid_t * gridptr,int * mask_gme)1329 gridInqMaskGMESerial(grid_t *gridptr, int *mask_gme)
1330 {
1331 return gridInqMaskSerialGeneric(gridptr, &gridptr->mask_gme, mask_gme);
1332 }
1333
gridInqMaskGME(int gridID,int * mask)1334 int gridInqMaskGME(int gridID, int *mask)
1335 {
1336 grid_t *gridptr = grid_to_pointer(gridID);
1337 return gridptr->vtable->inqMaskGME(gridptr, mask);
1338 }
1339
1340 static void
gridDefMaskGMESerial(grid_t * gridptr,const int * mask)1341 gridDefMaskGMESerial(grid_t *gridptr, const int *mask)
1342 {
1343 const size_t size = gridptr->size;
1344 if (size == 0) Error("Size undefined for gridID = %d", gridptr->self);
1345
1346 if (gridptr->mask_gme == NULL)
1347 gridptr->mask_gme = (mask_t *) Malloc(size * sizeof (mask_t));
1348 else if ( CDI_Debug )
1349 Warning("mask already defined!");
1350
1351 for (size_t i = 0; i < size; ++i)
1352 gridptr->mask_gme[i] = (mask_t)(mask[i] != 0);
1353 }
1354
gridDefMaskGME(int gridID,const int * mask)1355 void gridDefMaskGME(int gridID, const int *mask)
1356 {
1357 grid_t *gridptr = grid_to_pointer(gridID);
1358 gridptr->vtable->defMaskGME(gridptr, mask);
1359 gridMark4Update(gridID);
1360 }
1361
1362 static
gridInqXValsSerial(grid_t * gridptr,double * xvals)1363 SizeType gridInqXValsSerial(grid_t *gridptr, double *xvals)
1364 {
1365 const int gridtype = gridptr->type;
1366 size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->x.size;
1367
1368 if (CDI_Debug && size == 0) Warning("size undefined for gridID = %d", gridptr->self);
1369
1370 if (gridptr->x.vals)
1371 {
1372 if (size && xvals)
1373 {
1374 const double *gridptr_xvals = gridptr->vtable->inqXValsPtr(gridptr);
1375 memcpy(xvals, gridptr_xvals, size * sizeof (double));
1376 }
1377 }
1378 else
1379 size = 0;
1380
1381 return (SizeType)size;
1382 }
1383
1384 static
gridInqXValsPartSerial(grid_t * gridptr,int start,SizeType length,double * xvals)1385 SizeType gridInqXValsPartSerial(grid_t *gridptr, int start, SizeType length, double *xvals)
1386 {
1387 const int gridtype = gridptr->type;
1388 size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->x.size;
1389
1390 if (CDI_Debug && size == 0) Warning("size undefined for gridID = %d", gridptr->self);
1391
1392 if (gridptr->x.vals)
1393 {
1394 if (size && xvals && (size_t)length <= size)
1395 {
1396 const double *gridptr_xvals = gridptr->vtable->inqXValsPtr(gridptr);
1397 memcpy(xvals, gridptr_xvals+start, (size_t)length * sizeof (double));
1398 }
1399 }
1400 else
1401 size = 0;
1402
1403 return (SizeType)size;
1404 }
1405
1406 #ifndef USE_MPI
1407 static
gridInqXCvalsSerial(grid_t * gridptr,char ** xcvals)1408 SizeType gridInqXCvalsSerial(grid_t *gridptr, char **xcvals)
1409 {
1410 if (gridptr->type != GRID_CHARXY)
1411 Error("Function only valid for grid type 'GRID_CHARXY'.");
1412
1413 const size_t size = gridptr->x.size;
1414 size_t maxclength = 0;
1415
1416 const char **gridptr_xcvals = gridptr->vtable->inqXCvalsPtr(gridptr);
1417 if (gridptr_xcvals && size && xcvals)
1418 {
1419 maxclength = gridptr->x.clength;
1420 for ( size_t i = 0; i < size; i++ )
1421 memcpy(xcvals[i], gridptr_xcvals[i], maxclength*sizeof(char));
1422 }
1423
1424 return (SizeType)maxclength;
1425 }
1426
1427 static
gridInqXIscSerial(grid_t * gridptr)1428 int gridInqXIscSerial(grid_t *gridptr)
1429 {
1430 /*
1431 if ( gridptr->type != GRID_CHARXY )
1432 Error("Axis type is 'char' but grid is not type 'GRID_CHARXY'.");
1433 */
1434 return gridptr->x.clength;
1435 }
1436 #endif
1437
1438 /*
1439 @Function gridInqXvals
1440 @Title Get all values of a X-axis
1441
1442 @Prototype SizeType gridInqXvals(int gridID, double *xvals)
1443 @Parameter
1444 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
1445 @Item xvals Pointer to the location into which the X-values are read.
1446 The caller must allocate space for the returned values.
1447
1448 @Description
1449 The function @func{gridInqXvals} returns all values of the X-axis.
1450
1451 @Result
1452 Upon successful completion @func{gridInqXvals} returns the number of values and
1453 the values are stored in @func{xvals}.
1454 Otherwise, 0 is returned and @func{xvals} is empty.
1455
1456 @EndFunction
1457 */
gridInqXvals(int gridID,double * xvals)1458 SizeType gridInqXvals(int gridID, double *xvals)
1459 {
1460 grid_t *gridptr = grid_to_pointer(gridID);
1461 return gridptr->vtable->inqXVals(gridptr, xvals);
1462 }
1463
1464
gridInqXvalsPart(int gridID,int start,SizeType length,double * xvals)1465 SizeType gridInqXvalsPart(int gridID, int start, SizeType length, double *xvals)
1466 {
1467 grid_t *gridptr = grid_to_pointer(gridID);
1468 return gridptr->vtable->inqXValsPart(gridptr, start, length, xvals);
1469 }
1470
1471
gridInqXCvals(int gridID,char ** xcvals)1472 SizeType gridInqXCvals(int gridID, char **xcvals)
1473 {
1474 grid_t *gridptr = grid_to_pointer(gridID);
1475 #ifndef USE_MPI
1476 return gridptr->vtable->inqXCvals(gridptr, xcvals);
1477 #else
1478 return 0;
1479 #endif
1480 }
1481
1482
gridInqXIsc(int gridID)1483 int gridInqXIsc(int gridID)
1484 {
1485 grid_t *gridptr = grid_to_pointer(gridID);
1486 #ifndef USE_MPI
1487 return gridptr->vtable->inqXIsc(gridptr);
1488 #else
1489 return 0;
1490 #endif
1491 }
1492
1493 static
gridDefXValsSerial(grid_t * gridptr,const double * xvals)1494 void gridDefXValsSerial(grid_t *gridptr, const double *xvals)
1495 {
1496 const int gridtype = gridptr->type;
1497 const size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->x.size;
1498
1499 if (size == 0) Error("Size undefined for gridID = %d", gridptr->self);
1500
1501 if (gridptr->x.vals && CDI_Debug) Warning("values already defined!");
1502 gridptr->x.vals = (double *)Realloc(gridptr->x.vals, size * sizeof(double));
1503 memcpy(gridptr->x.vals, xvals, size * sizeof (double));
1504 }
1505
1506 #ifndef USE_MPI
1507 static
gridInqYCvalsSerial(grid_t * gridptr,char ** ycvals)1508 SizeType gridInqYCvalsSerial(grid_t *gridptr, char **ycvals)
1509 {
1510 if ( gridptr->type != GRID_CHARXY )
1511 Error("Function only valid for grid type 'GRID_CHARXY'.");
1512
1513 const size_t size = gridptr->y.size;
1514 size_t maxclength = 0;
1515
1516 const char **gridptr_ycvals = gridptr->vtable->inqYCvalsPtr(gridptr);
1517 if (gridptr_ycvals && size && ycvals)
1518 {
1519 maxclength = gridptr->y.clength;
1520 for (size_t i = 0; i < size; i++)
1521 memcpy(ycvals[i], gridptr_ycvals[i], maxclength*sizeof(char));
1522 }
1523
1524 return (SizeType)maxclength;
1525 }
1526
1527 static
gridInqYIscSerial(grid_t * gridptr)1528 int gridInqYIscSerial(grid_t *gridptr)
1529 {
1530 // if ( gridptr->type != GRID_CHARXY ) Error("Axis type is 'char' but grid is not type 'GRID_CHARXY'.");
1531 return gridptr->y.clength;
1532 }
1533 #endif
1534
1535 /*
1536 @Function gridDefXvals
1537 @Title Define the values of a X-axis
1538
1539 @Prototype void gridDefXvals(int gridID, const double *xvals)
1540 @Parameter
1541 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
1542 @Item xvals X-values of the grid.
1543
1544 @Description
1545 The function @func{gridDefXvals} defines all values of the X-axis.
1546
1547 @EndFunction
1548 */
gridDefXvals(int gridID,const double * xvals)1549 void gridDefXvals(int gridID, const double *xvals)
1550 {
1551 grid_t *gridptr = grid_to_pointer(gridID);
1552 gridptr->vtable->defXVals(gridptr, xvals);
1553 gridMark4Update(gridID);
1554 }
1555
1556 static
gridInqYValsSerial(grid_t * gridptr,double * yvals)1557 SizeType gridInqYValsSerial(grid_t *gridptr, double *yvals)
1558 {
1559 const int gridtype = gridptr->type;
1560 size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->y.size;
1561
1562 if (CDI_Debug && size == 0) Warning("size undefined for gridID = %d!", gridptr->self);
1563
1564 if (gridptr->y.vals)
1565 {
1566 if (size && yvals)
1567 {
1568 const double *gridptr_yvals = gridptr->vtable->inqYValsPtr(gridptr);
1569 memcpy(yvals, gridptr_yvals, size * sizeof (double));
1570 }
1571 }
1572 else
1573 size = 0;
1574
1575 return (SizeType)size;
1576 }
1577
1578 static
gridInqYValsPartSerial(grid_t * gridptr,int start,SizeType length,double * yvals)1579 SizeType gridInqYValsPartSerial(grid_t *gridptr, int start, SizeType length, double *yvals)
1580 {
1581 const int gridtype = gridptr->type;
1582 size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->y.size;
1583
1584 if (CDI_Debug && size == 0) Warning("size undefined for gridID = %d!", gridptr->self);
1585
1586 if (gridptr->y.vals)
1587 {
1588 if (size && yvals && (size_t)length <= size)
1589 {
1590 const double *gridptr_yvals = gridptr->vtable->inqYValsPtr(gridptr);
1591 memcpy(yvals, gridptr_yvals+start, (size_t)length * sizeof (double));
1592 }
1593 }
1594 else
1595 size = 0;
1596
1597 return (SizeType)size;
1598 }
1599
1600 /*
1601 @Function gridInqYvals
1602 @Title Get all values of a Y-axis
1603
1604 @Prototype SizeType gridInqYvals(int gridID, double *yvals)
1605 @Parameter
1606 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
1607 @Item yvals Pointer to the location into which the Y-values are read.
1608 The caller must allocate space for the returned values.
1609
1610 @Description
1611 The function @func{gridInqYvals} returns all values of the Y-axis.
1612
1613 @Result
1614 Upon successful completion @func{gridInqYvals} returns the number of values and
1615 the values are stored in @func{yvals}.
1616 Otherwise, 0 is returned and @func{yvals} is empty.
1617
1618 @EndFunction
1619 */
gridInqYvals(int gridID,double * yvals)1620 SizeType gridInqYvals(int gridID, double *yvals)
1621 {
1622 grid_t *gridptr = grid_to_pointer(gridID);
1623 return gridptr->vtable->inqYVals(gridptr, yvals);
1624 }
1625
1626
gridInqYvalsPart(int gridID,int start,SizeType size,double * yvals)1627 SizeType gridInqYvalsPart(int gridID, int start, SizeType size, double *yvals)
1628 {
1629 grid_t *gridptr = grid_to_pointer(gridID);
1630 return gridptr->vtable->inqYValsPart(gridptr, start, size, yvals);
1631 }
1632
1633
gridInqYCvals(int gridID,char ** ycvals)1634 SizeType gridInqYCvals(int gridID, char **ycvals)
1635 {
1636 grid_t *gridptr = grid_to_pointer(gridID);
1637 #ifndef USE_MPI
1638 return gridptr->vtable->inqYCvals(gridptr, ycvals);
1639 #else
1640 return 0;
1641 #endif
1642 }
1643
1644
gridInqYIsc(int gridID)1645 int gridInqYIsc(int gridID)
1646 {
1647 grid_t *gridptr = grid_to_pointer(gridID);
1648 #ifndef USE_MPI
1649 return gridptr->vtable->inqYIsc(gridptr);
1650 #else
1651 return 0;
1652 #endif
1653 }
1654
1655 static
gridDefYValsSerial(grid_t * gridptr,const double * yvals)1656 void gridDefYValsSerial(grid_t *gridptr, const double *yvals)
1657 {
1658 const int gridtype = gridptr->type;
1659 const size_t size = (gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED) ? gridptr->size : gridptr->y.size;
1660
1661 if (size == 0)
1662 Error("Size undefined for gridID = %d!", gridptr->self);
1663
1664 if (gridptr->y.vals && CDI_Debug)
1665 Warning("Values already defined!");
1666
1667 gridptr->y.vals = (double *)Realloc(gridptr->y.vals, size * sizeof(double));
1668 memcpy(gridptr->y.vals, yvals, size * sizeof(double));
1669 }
1670
1671
1672 /*
1673 @Function gridDefYvals
1674 @Title Define the values of a Y-axis
1675
1676 @Prototype void gridDefYvals(int gridID, const double *yvals)
1677 @Parameter
1678 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
1679 @Item yvals Y-values of the grid.
1680
1681 @Description
1682 The function @func{gridDefYvals} defines all values of the Y-axis.
1683
1684 @EndFunction
1685 */
gridDefYvals(int gridID,const double * yvals)1686 void gridDefYvals(int gridID, const double *yvals)
1687 {
1688 grid_t *gridptr = grid_to_pointer(gridID);
1689 gridptr->vtable->defYVals(gridptr, yvals);
1690 gridMark4Update(gridID);
1691 }
1692
1693 static double
gridInqXValSerial(grid_t * gridptr,SizeType index)1694 gridInqXValSerial(grid_t *gridptr, SizeType index)
1695 {
1696 const double xval = gridptr->x.vals ? gridptr->x.vals[index] : 0;
1697 return xval;
1698 }
1699
1700
gridInqXval(int gridID,SizeType index)1701 double gridInqXval(int gridID, SizeType index)
1702 {
1703 grid_t *gridptr = grid_to_pointer(gridID);
1704 return gridptr->vtable->inqXVal(gridptr, index);
1705 }
1706
1707 static double
gridInqYValSerial(grid_t * gridptr,SizeType index)1708 gridInqYValSerial(grid_t *gridptr, SizeType index)
1709 {
1710 const double yval = gridptr->y.vals ? gridptr->y.vals[index] : 0;
1711 return yval;
1712 }
1713
1714 /*
1715 @Function
1716 @Title
1717
1718 @Prototype
1719 @Parameter
1720 @Item Grid identifier
1721
1722 @EndFunction
1723 */
gridInqYval(int gridID,SizeType index)1724 double gridInqYval(int gridID, SizeType index)
1725 {
1726 grid_t *gridptr = grid_to_pointer(gridID);
1727 return gridptr->vtable->inqYVal(gridptr, index);
1728 }
1729
1730 static
grid_calc_increment(size_t size,const double * vals)1731 double grid_calc_increment(size_t size, const double *vals)
1732 {
1733 if (size > 1)
1734 {
1735 double inc = (vals[size-1] - vals[0])/(size-1);
1736 const double abs_inc = fabs(inc);
1737 for (size_t i = 1; i < size; ++i)
1738 if (fabs(fabs(vals[i-1] - vals[i]) - abs_inc) > 0.01*abs_inc)
1739 {
1740 inc = 0.0;
1741 break;
1742 }
1743
1744 return inc;
1745 }
1746
1747 return 0.0;
1748 }
1749
1750 static
grid_calc_increment_in_meter(SizeType size,const double * vals)1751 double grid_calc_increment_in_meter(SizeType size, const double *vals)
1752 {
1753 if (size > 1)
1754 {
1755 const double inc = (vals[size-1] - vals[0])/(size-1);
1756 return round(fabs(inc));
1757 }
1758
1759 return 0.0;
1760 }
1761
1762 /*
1763 @Function
1764 @Title
1765
1766 @Prototype
1767 @Parameter
1768 @Item Grid identifier
1769
1770 @EndFunction
1771 */
gridInqXinc(int gridID)1772 double gridInqXinc(int gridID)
1773 {
1774 grid_t *gridptr = grid_to_pointer(gridID);
1775 const double *restrict xvals = gridptr->vtable->inqXValsPtr(gridptr);
1776
1777 if (fabs(gridptr->x.inc) <= 0 && xvals)
1778 {
1779 const size_t xsize = gridptr->x.size;
1780 if (xsize > 1) gridptr->x.inc = grid_calc_increment(xsize, xvals);
1781 }
1782
1783 return gridptr->x.inc;
1784 }
1785
gridInqXincInMeter(int gridID)1786 double gridInqXincInMeter(int gridID)
1787 {
1788 grid_t *gridptr = grid_to_pointer(gridID);
1789 const double *restrict xvals = gridptr->vtable->inqXValsPtr(gridptr);
1790
1791 if (fabs(gridptr->x.inc) <= 0 && xvals)
1792 {
1793 const size_t xsize = gridptr->x.size;
1794 if (xsize > 1) gridptr->x.inc = grid_calc_increment_in_meter(xsize, xvals);
1795 }
1796
1797 return gridptr->x.inc;
1798 }
1799
1800 /*
1801 @Function
1802 @Title
1803
1804 @Prototype
1805 @Parameter
1806 @Item Grid identifier
1807
1808 @EndFunction
1809 */
gridInqYinc(int gridID)1810 double gridInqYinc(int gridID)
1811 {
1812 grid_t *gridptr = grid_to_pointer(gridID);
1813 const double *yvals = gridptr->vtable->inqYValsPtr(gridptr);
1814
1815 if (fabs(gridptr->y.inc) <= 0 && yvals)
1816 {
1817 const size_t ysize = gridptr->y.size;
1818 if (ysize > 1) gridptr->y.inc = grid_calc_increment(ysize, yvals);
1819 }
1820
1821 return gridptr->y.inc;
1822 }
1823
gridInqYincInMeter(int gridID)1824 double gridInqYincInMeter(int gridID)
1825 {
1826 grid_t *gridptr = grid_to_pointer(gridID);
1827 const double *yvals = gridptr->vtable->inqYValsPtr(gridptr);
1828
1829 if (fabs(gridptr->y.inc) <= 0 && yvals)
1830 {
1831 const size_t ysize = gridptr->y.size;
1832 if (ysize > 1) gridptr->y.inc = grid_calc_increment_in_meter(ysize, yvals);
1833 }
1834
1835 return gridptr->y.inc;
1836 }
1837
1838 /*
1839 @Function
1840 @Title
1841
1842 @Prototype
1843 @Parameter
1844 @Item Grid identifier
1845
1846 @EndFunction
1847 */
gridInqParamRLL(int gridID,double * xpole,double * ypole,double * angle)1848 void gridInqParamRLL(int gridID, double *xpole, double *ypole, double *angle)
1849 {
1850 *xpole = 0; *ypole = 0; *angle = 0;
1851
1852 const char *projection = "rotated_latitude_longitude";
1853 char gmapname[CDI_MAX_NAME];
1854 int length = CDI_MAX_NAME;
1855 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname, &length);
1856 if (gmapname[0] && strIsEqual(gmapname, projection))
1857 {
1858 int atttype, attlen;
1859 char attname[CDI_MAX_NAME+1];
1860
1861 int natts;
1862 cdiInqNatts(gridID, CDI_GLOBAL, &natts);
1863
1864 for (int iatt = 0; iatt < natts; ++iatt)
1865 {
1866 cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
1867 if (attlen != 1) continue;
1868
1869 double attflt;
1870 if (cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, &attflt))
1871 {
1872 if ( strIsEqual(attname, "grid_north_pole_longitude") ) *xpole = attflt;
1873 else if ( strIsEqual(attname, "grid_north_pole_latitude") ) *ypole = attflt;
1874 else if ( strIsEqual(attname, "north_pole_grid_longitude") ) *angle = attflt;
1875 }
1876 }
1877 }
1878 else
1879 Warning("%s mapping parameter missing!", projection);
1880 }
1881
1882 /*
1883 @Function
1884 @Title
1885
1886 @Prototype
1887 @Parameter
1888 @Item Grid identifier
1889
1890 @EndFunction
1891 */
gridDefParamRLL(int gridID,double xpole,double ypole,double angle)1892 void gridDefParamRLL(int gridID, double xpole, double ypole, double angle)
1893 {
1894 cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, "degrees");
1895 cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, "degrees");
1896
1897 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARNAME, "rotated_pole");
1898
1899 const char *gmapname = "rotated_latitude_longitude";
1900 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname);
1901 cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", (int)(strlen(gmapname)), gmapname);
1902 cdiDefAttFlt(gridID, CDI_GLOBAL, "grid_north_pole_longitude", CDI_DATATYPE_FLT64, 1, &xpole);
1903 cdiDefAttFlt(gridID, CDI_GLOBAL, "grid_north_pole_latitude", CDI_DATATYPE_FLT64, 1, &ypole);
1904 if ( IS_NOT_EQUAL(angle, 0) ) cdiDefAttFlt(gridID, CDI_GLOBAL, "north_pole_grid_longitude", CDI_DATATYPE_FLT64, 1, &angle);
1905
1906 grid_t *gridptr = grid_to_pointer(gridID);
1907 gridptr->projtype = CDI_PROJ_RLL;
1908
1909 gridVerifyProj(gridID);
1910 }
1911
1912 /*
1913 @Function
1914 @Title
1915
1916 @Prototype
1917 @Parameter
1918 @Item Grid identifier
1919
1920 @EndFunction
1921 */
gridInqParamGME(int gridID,int * nd,int * ni,int * ni2,int * ni3)1922 void gridInqParamGME(int gridID, int *nd, int *ni, int *ni2, int *ni3)
1923 {
1924 grid_t *gridptr = grid_to_pointer(gridID);
1925
1926 *nd = gridptr->gme.nd;
1927 *ni = gridptr->gme.ni;
1928 *ni2 = gridptr->gme.ni2;
1929 *ni3 = gridptr->gme.ni3;
1930 }
1931
1932 /*
1933 @Function
1934 @Title
1935
1936 @Prototype
1937 @Parameter
1938 @Item Grid identifier
1939
1940 @EndFunction
1941 */
gridDefParamGME(int gridID,int nd,int ni,int ni2,int ni3)1942 void gridDefParamGME(int gridID, int nd, int ni, int ni2, int ni3)
1943 {
1944 grid_t *gridptr = grid_to_pointer(gridID);
1945
1946 if ( gridptr->gme.nd != nd )
1947 {
1948 gridptr->gme.nd = nd;
1949 gridptr->gme.ni = ni;
1950 gridptr->gme.ni2 = ni2;
1951 gridptr->gme.ni3 = ni3;
1952 gridMark4Update(gridID);
1953 }
1954 }
1955
1956 /*
1957 @Function
1958 @Title
1959
1960 @Prototype
1961 @Parameter
1962 @Item Grid identifier
1963
1964 @EndFunction
1965 */
gridChangeType(int gridID,int gridtype)1966 void gridChangeType(int gridID, int gridtype)
1967 {
1968 grid_t *gridptr = grid_to_pointer(gridID);
1969
1970 if ( CDI_Debug )
1971 Message("Changed grid type from %s to %s", gridNamePtr(gridptr->type), gridNamePtr(gridtype));
1972
1973 if (gridptr->type != gridtype)
1974 {
1975 gridptr->type = gridtype;
1976 gridMark4Update(gridID);
1977 }
1978 }
1979
1980 static
grid_check_cyclic(grid_t * gridptr)1981 void grid_check_cyclic(grid_t *gridptr)
1982 {
1983 gridptr->isCyclic = 0;
1984 enum { numVertices = 4 };
1985 size_t xsize = gridptr->x.size,
1986 ysize = gridptr->y.size;
1987 const double *xvals = gridptr->vtable->inqXValsPtr(gridptr),
1988 *yvals = gridptr->vtable->inqYValsPtr(gridptr),
1989 (*xbounds)[numVertices]
1990 = (const double (*)[numVertices])gridptr->vtable->inqXBoundsPtr(gridptr);
1991
1992 if ( gridptr->type == GRID_GAUSSIAN || gridptr->type == GRID_LONLAT )
1993 {
1994 if ( xvals && xsize > 1 )
1995 {
1996 double xval1 = xvals[0];
1997 double xval2 = xvals[1];
1998 double xvaln = xvals[xsize-1];
1999 if ( xval2 < xval1 ) xval2 += 360;
2000 if ( xvaln < xval1 ) xvaln += 360;
2001
2002 if ( IS_NOT_EQUAL(xval1, xvaln) )
2003 {
2004 double xinc = xval2 - xval1;
2005 if ( IS_EQUAL(xinc, 0) ) xinc = (xvaln - xval1)/(xsize-1);
2006
2007 const double x0 = xvaln + xinc - 360;
2008
2009 if ( fabs(x0 - xval1) < 0.01*xinc ) gridptr->isCyclic = 1;
2010 }
2011 }
2012 }
2013 else if ( gridptr->type == GRID_CURVILINEAR )
2014 {
2015 bool lcheck = true;
2016 if ( yvals && xvals )
2017 {
2018 if ( (fabs(yvals[0] - yvals[xsize-1]) > fabs(yvals[0] - yvals[xsize*ysize-xsize])) &&
2019 (fabs(yvals[xsize*ysize-xsize] - yvals[xsize*ysize-1]) > fabs(yvals[xsize-1] - yvals[xsize*ysize-1])) )
2020 lcheck = false;
2021 }
2022 else lcheck = false;
2023
2024 if ( lcheck && xvals && xsize > 1 )
2025 {
2026 size_t nc = 0;
2027 for ( size_t j = 0; j < ysize; ++j )
2028 {
2029 size_t i1 = j*xsize,
2030 i2 = j*xsize+1,
2031 in = j*xsize+(xsize-1);
2032 double val1 = xvals[i1],
2033 val2 = xvals[i2],
2034 valn = xvals[in];
2035 double xinc = fabs(val2-val1);
2036
2037 if ( val1 < 1 && valn > 300 ) val1 += 360;
2038 if ( valn < 1 && val1 > 300 ) valn += 360;
2039 if ( val1 < -179 && valn > 120 ) val1 += 360;
2040 if ( valn < -179 && val1 > 120 ) valn += 360;
2041 if ( fabs(valn-val1) > 180 ) val1 += 360;
2042
2043 double x0 = valn + copysign(xinc, val1 - valn);
2044
2045 nc += fabs(x0-val1) < 0.5*xinc;
2046 }
2047 gridptr->isCyclic = nc > ysize/2;
2048 }
2049
2050 if ( lcheck && xbounds && xsize > 1 )
2051 {
2052 bool isCyclic = true;
2053 for ( size_t j = 0; j < ysize; ++j )
2054 {
2055 size_t i1 = j*xsize,
2056 i2 = j*xsize+(xsize-1);
2057 for (size_t k1 = 0; k1 < numVertices; ++k1 )
2058 {
2059 double val1 = xbounds[i1][k1];
2060 for (size_t k2 = 0; k2 < numVertices; ++k2 )
2061 {
2062 double val2 = xbounds[i2][k2];
2063
2064 if ( val1 < 1 && val2 > 300 ) val1 += 360;
2065 if ( val2 < 1 && val1 > 300 ) val2 += 360;
2066 if ( val1 < -179 && val2 > 120 ) val1 += 360;
2067 if ( val2 < -179 && val1 > 120 ) val2 += 360;
2068 if ( fabs(val2-val1) > 180 ) val1 += 360;
2069
2070 if ( fabs(val1-val2) < 0.001 )
2071 goto foundCloseVertices;
2072 }
2073 }
2074 // all vertices more than 0.001 degrees apart
2075 isCyclic = false;
2076 break;
2077 foundCloseVertices:
2078 ;
2079 }
2080 gridptr->isCyclic = isCyclic;
2081 }
2082 }
2083 }
2084
2085
gridIsCircular(int gridID)2086 int gridIsCircular(int gridID)
2087 {
2088 grid_t *gridptr = grid_to_pointer(gridID);
2089
2090 if ( gridptr->isCyclic == CDI_UNDEFID ) grid_check_cyclic(gridptr);
2091
2092 return gridptr->isCyclic;
2093 }
2094
2095 static
compareXYvals(grid_t * gridRef,grid_t * gridTest)2096 bool compareXYvals(grid_t *gridRef, grid_t *gridTest)
2097 {
2098 bool differ = false;
2099
2100 size_t xsizeTest = gridTest->x.size, ysizeTest = gridTest->y.size;
2101 if (!differ && xsizeTest > 0 && xsizeTest == (size_t)gridRef->vtable->inqXVals(gridRef, NULL))
2102 {
2103 const double *restrict xvalsRef = gridRef->vtable->inqXValsPtr(gridRef),
2104 *restrict xvalsTest = gridTest->vtable->inqXValsPtr(gridTest);
2105
2106 for (size_t i = 0; i < xsizeTest; ++i)
2107 if (fabs(xvalsTest[i] - xvalsRef[i]) > 1.e-10)
2108 {
2109 differ = true;
2110 break;
2111 }
2112 }
2113
2114 if (!differ && ysizeTest > 0 && ysizeTest == (size_t)gridRef->vtable->inqYVals(gridRef, NULL))
2115 {
2116 const double *restrict yvalsRef = gridRef->vtable->inqYValsPtr(gridRef),
2117 *restrict yvalsTest = gridTest->vtable->inqYValsPtr(gridTest);
2118 for (size_t i = 0; i < ysizeTest; ++i)
2119 if (fabs(yvalsTest[i] - yvalsRef[i]) > 1.e-10)
2120 {
2121 differ = true;
2122 break;
2123 }
2124 }
2125
2126 return differ;
2127 }
2128
2129 static
compareXYvals2(grid_t * gridRef,grid_t * gridTest)2130 bool compareXYvals2(grid_t *gridRef, grid_t *gridTest)
2131 {
2132 size_t gridsize = gridTest->size;
2133 bool differ = ((gridTest->x.vals == NULL) ^ (gridRef->x.vals == NULL))
2134 || ((gridTest->y.vals == NULL) ^ (gridRef->y.vals == NULL))
2135 || ((gridTest->x.bounds == NULL) ^ (gridRef->x.bounds == NULL))
2136 || ((gridTest->y.bounds == NULL) ^ (gridRef->y.bounds == NULL));
2137
2138 typedef double (*inqVal)(grid_t *grid, SizeType index);
2139 inqVal inqXValRef = gridRef->vtable->inqXVal,
2140 inqYValRef = gridRef->vtable->inqYVal,
2141 inqXValTest = gridTest->vtable->inqXVal,
2142 inqYValTest = gridTest->vtable->inqYVal;
2143
2144 if ( !differ && gridTest->x.vals )
2145 differ = fabs(inqXValTest(gridTest, 0) - inqXValRef(gridRef, 0)) > 1.e-9
2146 || fabs(inqXValTest(gridTest, gridsize-1) - inqXValRef(gridRef, gridsize-1)) > 1.e-9;
2147
2148 if ( !differ && gridTest->y.vals )
2149 differ = fabs(inqYValTest(gridTest, 0) - inqYValRef(gridRef, 0)) > 1.e-9
2150 || fabs(inqYValTest(gridTest, gridsize-1) - inqYValRef(gridRef, gridsize-1)) > 1.e-9;
2151
2152 return differ;
2153 }
2154
2155 static
compare_bounds(const grid_t * grid,const grid_t * gridRef)2156 bool compare_bounds(const grid_t *grid, const grid_t *gridRef)
2157 {
2158 bool differ = false;
2159
2160 if ((grid->x.bounds && !gridRef->x.bounds) || (!grid->x.bounds && gridRef->x.bounds) ||
2161 (grid->y.bounds && !gridRef->y.bounds) || (!grid->y.bounds && gridRef->y.bounds))
2162 differ = true;
2163
2164 return differ;
2165 }
2166
2167 static
compare_lonlat(int gridID,const grid_t * grid,const grid_t * gridRef)2168 bool compare_lonlat(int gridID, const grid_t *grid, const grid_t *gridRef)
2169 {
2170 bool differ = false;
2171 /*
2172 printf("gridID %d\n", gridID);
2173 printf("grid.xdef %d\n", grid->x.flag);
2174 printf("grid.ydef %d\n", grid->y.flag);
2175 printf("grid.xsize %zu\n", grid->x.size);
2176 printf("grid.ysize %zu\n", grid->y.size);
2177 printf("grid.xfirst %f\n", grid->x.first);
2178 printf("grid.yfirst %f\n", grid->y.first);
2179 printf("grid.xfirst %f\n", gridInqXval(gridID, 0));
2180 printf("grid.yfirst %f\n", gridInqYval(gridID, 0));
2181 printf("grid.xinc %f\n", grid->x.inc);
2182 printf("grid.yinc %f\n", grid->y.inc);
2183 printf("grid.xinc %f\n", gridInqXinc(gridID));
2184 printf("grid.yinc %f\n", gridInqYinc(gridID));
2185 */
2186 if ( grid->x.size == gridRef->x.size && grid->y.size == gridRef->y.size )
2187 {
2188 if ( grid->x.flag == 2 && grid->y.flag == 2 )
2189 {
2190 if ( ! (IS_EQUAL(grid->x.first, 0) && IS_EQUAL(grid->x.last, 0) && IS_EQUAL(grid->x.inc, 0)) &&
2191 ! (IS_EQUAL(grid->y.first, 0) && IS_EQUAL(grid->y.last, 0) && IS_EQUAL(grid->y.inc, 0)) &&
2192 IS_NOT_EQUAL(grid->x.first, grid->x.last) && IS_NOT_EQUAL(grid->y.first, grid->y.last) )
2193 {
2194 if ( IS_NOT_EQUAL(grid->x.first, gridInqXval(gridID, 0)) ||
2195 IS_NOT_EQUAL(grid->y.first, gridInqYval(gridID, 0)))
2196 {
2197 differ = true;
2198 }
2199 if ( !differ && fabs(grid->x.inc) > 0 &&
2200 fabs(fabs(grid->x.inc) - fabs(gridRef->x.inc)) > fabs(grid->x.inc/1000))
2201 {
2202 differ = true;
2203 }
2204 if ( !differ && fabs(grid->y.inc) > 0 &&
2205 fabs(fabs(grid->y.inc) - fabs(gridRef->y.inc)) > fabs(grid->y.inc/1000))
2206 {
2207 differ = true;
2208 }
2209 }
2210 }
2211 else if ( grid->x.vals && grid->y.vals )
2212 differ = gridRef->vtable->compareXYFull((grid_t *)gridRef, (grid_t *)grid);
2213
2214 if (!differ) differ = compare_bounds(grid, gridRef);
2215 }
2216 else
2217 differ = true;
2218
2219 return differ;
2220 }
2221
2222
2223 static
compare_projection(int gridID,const grid_t * grid,const grid_t * gridRef)2224 bool compare_projection(int gridID, const grid_t *grid, const grid_t *gridRef)
2225 {
2226 bool differ = compare_lonlat(gridID, grid, gridRef);
2227
2228 if (!differ)
2229 {
2230 //printf(">%s< >%s<\n", cdiInqVarKeyString(&grid->keys, CDI_KEY_GRIDMAP_VARNAME), cdiInqVarKeyString(&gridRef->keys, CDI_KEY_GRIDMAP_VARNAME));
2231 //printf(">%s< >%s<\n", cdiInqVarKeyString(&grid->keys, CDI_KEY_GRIDMAP_NAME), cdiInqVarKeyString(&gridRef->keys, CDI_KEY_GRIDMAP_NAME));
2232 // if (strcmp(cdiInqVarKeyString(&grid->keys, CDI_KEY_GRIDMAP_VARNAME), cdiInqVarKeyString(&gridRef->keys, CDI_KEY_GRIDMAP_VARNAME))) return true;
2233 // if (strcmp(cdiInqVarKeyString(&grid->keys, CDI_KEY_GRIDMAP_NAME), cdiInqVarKeyString(&gridRef->keys, CDI_KEY_GRIDMAP_NAME))) return true;
2234 }
2235
2236 return differ;
2237 }
2238
2239 static
compare_generic(const grid_t * grid,const grid_t * gridRef)2240 bool compare_generic(const grid_t *grid, const grid_t *gridRef)
2241 {
2242 bool differ = false;
2243
2244 if ( grid->x.size == gridRef->x.size && grid->y.size == gridRef->y.size )
2245 {
2246 if ( grid->x.flag == 1 && grid->y.flag == 1
2247 && grid->x.vals && grid->y.vals )
2248 differ = gridRef->vtable->compareXYFull((grid_t *)gridRef, (grid_t *)grid);
2249 }
2250 else if ( (grid->y.size == 0 || grid->y.size == 1) &&
2251 grid->x.size == gridRef->x.size*gridRef->y.size )
2252 {
2253 }
2254 else
2255 differ = true;
2256
2257 return differ;
2258 }
2259
2260 static
compare_gaussian(int gridID,const grid_t * grid,const grid_t * gridRef)2261 bool compare_gaussian(int gridID, const grid_t *grid, const grid_t *gridRef)
2262 {
2263 const double cmp_eps = 0.0015;
2264 bool differ = false;
2265
2266 if ( grid->x.size == gridRef->x.size && grid->y.size == gridRef->y.size )
2267 {
2268 if ( grid->x.flag == 2 && grid->y.flag == 2 )
2269 {
2270 if ( ! (IS_EQUAL(grid->x.first, 0) && IS_EQUAL(grid->x.last, 0) && IS_EQUAL(grid->x.inc, 0)) &&
2271 ! (IS_EQUAL(grid->y.first, 0) && IS_EQUAL(grid->y.last, 0)) )
2272 if ( fabs(grid->x.first - gridInqXval(gridID, 0)) > cmp_eps ||
2273 fabs(grid->y.first - gridInqYval(gridID, 0)) > cmp_eps ||
2274 (fabs(grid->x.inc)>0 && fabs(fabs(grid->x.inc) - fabs(gridRef->x.inc)) > fabs(grid->x.inc/1000)) )
2275 {
2276 differ = true;
2277 }
2278 }
2279 else if ( grid->x.vals && grid->y.vals )
2280 differ = gridRef->vtable->compareXYFull((grid_t *)gridRef, (grid_t *)grid);
2281
2282 if (!differ) differ = compare_bounds(grid, gridRef);
2283 }
2284 else
2285 differ = true;
2286
2287 return differ;
2288 }
2289
2290 static
compare_curvilinear(const grid_t * grid,const grid_t * gridRef)2291 bool compare_curvilinear(const grid_t *grid, const grid_t *gridRef)
2292 {
2293 bool differ = false;
2294
2295 /*
2296 printf("gridID %d\n", gridID);
2297 printf("grid.xsize %d\n", grid->x.size);
2298 printf("grid.ysize %d\n", grid->y.size);
2299 printf("grid.xfirst %f\n", grid->x.vals[0]);
2300 printf("grid.yfirst %f\n", grid->y.vals[0]);
2301 printf("grid xfirst %f\n", gridInqXval(gridID, 0));
2302 printf("grid yfirst %f\n", gridInqYval(gridID, 0));
2303 printf("grid.xlast %f\n", grid->x.vals[grid->size-1]);
2304 printf("grid.ylast %f\n", grid->y.vals[grid->size-1]);
2305 printf("grid xlast %f\n", gridInqXval(gridID, grid->size-1));
2306 printf("grid ylast %f\n", gridInqYval(gridID, grid->size-1));
2307 printf("grid.nv %d\n", grid->nvertex);
2308 printf("grid nv %d\n", gridInqNvertex(gridID));
2309 */
2310 if ( grid->x.size == gridRef->x.size && grid->y.size == gridRef->y.size )
2311 differ = gridRef->vtable->compareXYAO((grid_t *)gridRef, (grid_t *)grid);
2312
2313 return differ;
2314 }
2315
2316 static
compare_unstructured(const grid_t * grid,const grid_t * gridRef,bool coord_compare)2317 bool compare_unstructured(const grid_t *grid, const grid_t *gridRef, bool coord_compare)
2318 {
2319 bool differ = false;
2320
2321 unsigned char uuid1[CDI_UUID_SIZE] = { 0 };
2322 unsigned char uuid2[CDI_UUID_SIZE] = { 0 };
2323 int length = CDI_UUID_SIZE;
2324 cdiInqVarKeyBytes(&gridRef->keys, CDI_KEY_UUID, uuid1, &length);
2325 length = CDI_UUID_SIZE;
2326 cdiInqVarKeyBytes(&grid->keys, CDI_KEY_UUID, uuid2, &length);
2327 differ = ((!cdiUUIDIsNull(uuid1) || !cdiUUIDIsNull(uuid2)) && memcmp(uuid1, uuid2, CDI_UUID_SIZE));
2328 if (!differ)
2329 {
2330 int numberA = cdiInqVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDUSED);
2331 int numberB = cdiInqVarKeyInt(&gridRef->keys, CDI_KEY_NUMBEROFGRIDUSED);
2332 int positionA = cdiInqVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE);
2333 int positionB = cdiInqVarKeyInt(&gridRef->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE);
2334 if ( coord_compare )
2335 {
2336 differ = grid->nvertex != gridRef->nvertex
2337 || (numberA > 0 && positionA != positionB)
2338 || gridRef->vtable->compareXYFull((grid_t *)gridRef, (grid_t *)grid);
2339 }
2340 else
2341 {
2342 if ( ((grid->x.vals == NULL) ^ (gridRef->x.vals == NULL)) &&
2343 ((grid->y.vals == NULL) ^ (gridRef->y.vals == NULL)) )
2344 {
2345 int nvertexA = grid->nvertex, nvertexB = gridRef->nvertex;
2346 differ = ( nvertexA && nvertexB && (nvertexA != nvertexB) )
2347 || (( numberA && numberB && (numberA != numberB) )
2348 || ( numberA && numberB && positionA != positionB ) );
2349 }
2350 else
2351 {
2352 differ = grid->nvertex != gridRef->nvertex
2353 || numberA != numberB
2354 || (numberA > 0 && positionA != positionB)
2355 || gridRef->vtable->compareXYAO((grid_t *)gridRef, (grid_t *)grid);
2356 }
2357 }
2358 }
2359
2360 return differ;
2361 }
2362
2363 static
gridCompare(int gridID,const grid_t * grid,bool coord_compare)2364 bool gridCompare(int gridID, const grid_t *grid, bool coord_compare)
2365 {
2366 bool differ = true;
2367 const grid_t *gridRef = grid_to_pointer(gridID);
2368
2369 if ( grid->type == gridRef->type || grid->type == GRID_GENERIC )
2370 {
2371 if ( grid->size == gridRef->size )
2372 {
2373 differ = false;
2374 if ( grid->type == GRID_LONLAT )
2375 {
2376 differ = compare_lonlat(gridID, grid, gridRef);
2377 }
2378 else if ( grid->type == GRID_PROJECTION )
2379 {
2380 differ = compare_projection(gridID, grid, gridRef);
2381 }
2382 else if ( grid->type == GRID_GENERIC )
2383 {
2384 differ = compare_generic(grid, gridRef);
2385 }
2386 else if ( grid->type == GRID_GAUSSIAN )
2387 {
2388 differ = compare_gaussian(gridID, grid, gridRef);
2389 }
2390 else if ( grid->type == GRID_CURVILINEAR )
2391 {
2392 differ = compare_curvilinear(grid, gridRef);
2393 }
2394 else if ( grid->type == GRID_UNSTRUCTURED )
2395 {
2396 differ = compare_unstructured(grid, gridRef, coord_compare);
2397 }
2398 }
2399 }
2400
2401 int scanningModeA = cdiInqVarKeyInt(&grid->keys, CDI_KEY_SCANNINGMODE);
2402 int scanningModeB = cdiInqVarKeyInt(&gridRef->keys, CDI_KEY_SCANNINGMODE);
2403 if ( scanningModeA != scanningModeB )
2404 {
2405 // often grid definition may differ in UV-relativeToGrid
2406 differ = 1;
2407 #ifdef HIRLAM_EXTENSIONS
2408 if ( cdiDebugExt>=200 )
2409 printf("gridCompare(gridID=%d): Differs: scanningModeA [%d] != scanningModeB(gridID) [%d]\n",
2410 gridID, scanningModeA, scanningModeB);
2411 #endif // HIRLAM_EXTENSIONS
2412 }
2413
2414 return differ;
2415 }
2416
2417
cmp_key_int(const cdi_keys_t * keysp1,const cdi_keys_t * keysp2,int key)2418 int cmp_key_int(const cdi_keys_t *keysp1, const cdi_keys_t *keysp2, int key)
2419 {
2420 int v1 = cdiInqVarKeyInt(keysp1, key);
2421 int v2 = cdiInqVarKeyInt(keysp2, key);
2422 if (v1 != v2) return 1;
2423 return 0;
2424 }
2425
2426
gridCompareP(void * gridptr1,void * gridptr2)2427 int gridCompareP(void *gridptr1, void *gridptr2)
2428 {
2429 grid_t *g1 = ( grid_t * ) gridptr1;
2430 grid_t *g2 = ( grid_t * ) gridptr2;
2431 enum { equal = 0,
2432 differ = -1 };
2433 size_t i, size;
2434
2435 xassert ( g1 );
2436 xassert ( g2 );
2437
2438 if (cdiInqVarKeyInt(&g1->keys, CDI_KEY_DATATYPE) != cdiInqVarKeyInt(&g2->keys, CDI_KEY_DATATYPE)) return differ;
2439 if ( g1->type != g2->type ) return differ;
2440 if ( g1->isCyclic != g2->isCyclic ) return differ;
2441 if ( g1->x.flag != g2->x.flag ) return differ;
2442 if ( g1->y.flag != g2->y.flag ) return differ;
2443 if ( g1->gme.nd != g2->gme.nd ) return differ;
2444 if ( g1->gme.ni != g2->gme.ni ) return differ;
2445 if ( g1->gme.ni2 != g2->gme.ni2 ) return differ;
2446 if ( g1->gme.ni3 != g2->gme.ni3 ) return differ;
2447 if ( cmp_key_int(&g1->keys, &g2->keys, CDI_KEY_NUMBEROFGRIDUSED) ) return differ;
2448 if ( cmp_key_int(&g1->keys, &g2->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE) ) return differ;
2449 if ( g1->trunc != g2->trunc ) return differ;
2450 if ( g1->nvertex != g2->nvertex ) return differ;
2451 if ( g1->reducedPointsSize != g2->reducedPointsSize ) return differ;
2452 if ( g1->size != g2->size ) return differ;
2453 if ( g1->x.size != g2->x.size ) return differ;
2454 if ( g1->y.size != g2->y.size ) return differ;
2455 if ( g1->lcomplex != g2->lcomplex ) return differ;
2456
2457 if ( IS_NOT_EQUAL(g1->x.first , g2->x.first) ) return differ;
2458 if ( IS_NOT_EQUAL(g1->y.first , g2->y.first) ) return differ;
2459 if ( IS_NOT_EQUAL(g1->x.last , g2->x.last) ) return differ;
2460 if ( IS_NOT_EQUAL(g1->y.last , g2->y.last) ) return differ;
2461 if ( IS_NOT_EQUAL(g1->x.inc , g2->x.inc) ) return differ;
2462 if ( IS_NOT_EQUAL(g1->y.inc , g2->y.inc) ) return differ;
2463 if ( cmp_key_int(&g1->keys, &g2->keys, CDI_KEY_SCANNINGMODE) ) return differ;
2464
2465 const double *restrict g1_xvals = g1->vtable->inqXValsPtr(g1),
2466 *restrict g2_xvals = g2->vtable->inqXValsPtr(g2);
2467 if ( g1_xvals )
2468 {
2469 if ( g1->type == GRID_UNSTRUCTURED || g1->type == GRID_CURVILINEAR )
2470 size = g1->size;
2471 else
2472 size = g1->x.size;
2473 xassert ( size );
2474
2475 if ( !g2_xvals ) return differ;
2476
2477 for ( i = 0; i < size; i++ )
2478 if ( IS_NOT_EQUAL(g1_xvals[i], g2_xvals[i]) ) return differ;
2479 }
2480 else if ( g2_xvals )
2481 return differ;
2482
2483 const double *restrict g1_yvals = g1->vtable->inqYValsPtr(g1),
2484 *restrict g2_yvals = g2->vtable->inqYValsPtr(g2);
2485 if ( g1_yvals )
2486 {
2487 if ( g1->type == GRID_UNSTRUCTURED || g1->type == GRID_CURVILINEAR )
2488 size = g1->size;
2489 else
2490 size = g1->y.size;
2491 xassert ( size );
2492
2493 if ( !g2_yvals ) return differ;
2494
2495 for ( i = 0; i < size; i++ )
2496 if ( IS_NOT_EQUAL(g1_yvals[i], g2_yvals[i]) ) return differ;
2497 }
2498 else if ( g2_yvals )
2499 return differ;
2500
2501 const double *restrict g1_area = g1->vtable->inqAreaPtr(g1),
2502 *restrict g2_area = g2->vtable->inqAreaPtr(g2);
2503 if ( g1_area )
2504 {
2505 xassert ( g1->size );
2506
2507 if ( !g2_area ) return differ;
2508
2509 for ( i = 0; i < g1->size; i++ )
2510 if ( IS_NOT_EQUAL(g1_area[i], g2_area[i]) ) return differ;
2511 }
2512 else if ( g2_area )
2513 return differ;
2514
2515 {
2516 const double *restrict g1_xbounds, *restrict g2_xbounds;
2517 if ( (g1_xbounds = g1->vtable->inqXBoundsPtr(g1)) )
2518 {
2519 xassert ( g1->nvertex );
2520 if ( g1->type == GRID_CURVILINEAR || g1->type == GRID_UNSTRUCTURED )
2521 size = g1->nvertex * g1->size;
2522 else
2523 size = g1->nvertex * g1->x.size;
2524 xassert ( size );
2525
2526 if ( !(g2_xbounds = g2->vtable->inqXBoundsPtr(g2)) ) return differ;
2527
2528 for ( i = 0; i < size; i++ )
2529 if ( IS_NOT_EQUAL(g1_xbounds[i], g2_xbounds[i]) ) return differ;
2530 }
2531 else if ( g2->vtable->inqXBoundsPtr(g2) )
2532 return differ;
2533 }
2534
2535 {
2536 const double *restrict g1_ybounds, *restrict g2_ybounds;
2537 if ( (g1_ybounds = g1->vtable->inqYBoundsPtr(g1)) )
2538 {
2539 xassert ( g1->nvertex );
2540 if ( g1->type == GRID_CURVILINEAR || g1->type == GRID_UNSTRUCTURED )
2541 size = g1->nvertex * g1->size;
2542 else
2543 size = g1->nvertex * g1->y.size;
2544 xassert ( size );
2545
2546 if ( ! (g2_ybounds = g2->vtable->inqYBoundsPtr(g2)) ) return differ;
2547
2548 for ( i = 0; i < size; i++ )
2549 if ( IS_NOT_EQUAL(g1->y.bounds[i], g2->y.bounds[i]) ) return differ;
2550 }
2551 else if ( g2->vtable->inqYBoundsPtr(g2) )
2552 return differ;
2553 }
2554
2555 if (strcmp(cdiInqVarKeyString(&g1->x.keys, CDI_KEY_NAME), cdiInqVarKeyString(&g2->x.keys, CDI_KEY_NAME))) return differ;
2556 if (strcmp(cdiInqVarKeyString(&g1->y.keys, CDI_KEY_NAME), cdiInqVarKeyString(&g2->y.keys, CDI_KEY_NAME))) return differ;
2557 if (strcmp(cdiInqVarKeyString(&g1->x.keys, CDI_KEY_LONGNAME), cdiInqVarKeyString(&g2->x.keys, CDI_KEY_LONGNAME))) return differ;
2558 if (strcmp(cdiInqVarKeyString(&g1->y.keys, CDI_KEY_LONGNAME), cdiInqVarKeyString(&g2->y.keys, CDI_KEY_LONGNAME))) return differ;
2559 if (strcmp(cdiInqVarKeyString(&g1->x.keys, CDI_KEY_UNITS), cdiInqVarKeyString(&g2->x.keys, CDI_KEY_UNITS))) return differ;
2560 if (strcmp(cdiInqVarKeyString(&g1->y.keys, CDI_KEY_UNITS), cdiInqVarKeyString(&g2->y.keys, CDI_KEY_UNITS))) return differ;
2561 if (strcmp(cdiInqVarKeyString(&g1->x.keys, CDI_KEY_STDNAME), cdiInqVarKeyString(&g2->x.keys, CDI_KEY_STDNAME))) return differ;
2562 if (strcmp(cdiInqVarKeyString(&g1->y.keys, CDI_KEY_STDNAME), cdiInqVarKeyString(&g2->y.keys, CDI_KEY_STDNAME))) return differ;
2563
2564 if (strcmp(cdiInqVarKeyString(&g1->y.keys, CDI_KEY_REFERENCEURI), cdiInqVarKeyString(&g2->y.keys, CDI_KEY_REFERENCEURI))) return differ;
2565
2566 if ( g1->mask )
2567 {
2568 xassert ( g1->size );
2569 if ( !g2->mask ) return differ;
2570 if ( memcmp ( g1->mask, g2->mask, g1->size * sizeof(mask_t)) ) return differ;
2571 }
2572 else if ( g2->mask )
2573 return differ;
2574
2575 if ( g1->mask_gme )
2576 {
2577 xassert ( g1->size );
2578 if ( !g2->mask_gme ) return differ;
2579 if ( memcmp ( g1->mask_gme, g2->mask_gme, g1->size * sizeof(mask_t)) ) return differ;
2580 }
2581 else if ( g2->mask_gme )
2582 return differ;
2583
2584 unsigned char uuid1[CDI_UUID_SIZE] = { 0 };
2585 unsigned char uuid2[CDI_UUID_SIZE] = { 0 };
2586 int length = CDI_UUID_SIZE;
2587 cdiInqVarKeyBytes(&g1->keys, CDI_KEY_UUID, uuid1, &length);
2588 length = CDI_UUID_SIZE;
2589 cdiInqVarKeyBytes(&g2->keys, CDI_KEY_UUID, uuid2, &length);
2590 if ( memcmp(uuid1, uuid2, CDI_UUID_SIZE) ) return differ;
2591
2592 return equal;
2593 }
2594
2595 static
gridComplete(grid_t * grid)2596 void gridComplete(grid_t *grid)
2597 {
2598 const int gridID = grid->self;
2599
2600 if (grid->datatype != CDI_UNDEFID) cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_DATATYPE, grid->datatype);
2601
2602 const int gridtype = grid->type;
2603 switch (gridtype)
2604 {
2605 case GRID_LONLAT:
2606 case GRID_GAUSSIAN:
2607 case GRID_UNSTRUCTURED:
2608 case GRID_CURVILINEAR:
2609 case GRID_GENERIC:
2610 case GRID_PROJECTION:
2611 case GRID_CHARXY:
2612 {
2613 if ( grid->x.size > 0 ) gridDefXsize(gridID, grid->x.size);
2614 if ( grid->y.size > 0 ) gridDefYsize(gridID, grid->y.size);
2615
2616 if ( gridtype == GRID_GAUSSIAN ) gridDefNP(gridID, grid->np);
2617
2618 if ( grid->nvertex > 0 )
2619 gridDefNvertex(gridID, grid->nvertex);
2620
2621 if ( grid->x.flag == 2 )
2622 {
2623 assert(gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR);
2624 double *xvals = (double *) Malloc(grid->x.size * sizeof (double));
2625 gridGenXvals(grid->x.size, grid->x.first, grid->x.last, grid->x.inc, xvals);
2626 grid->x.vals = xvals;
2627 // gridDefXinc(gridID, grid->x.inc);
2628 }
2629
2630 if ( grid->y.flag == 2 )
2631 {
2632 assert(gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR);
2633 double *yvals = (double *) Malloc(grid->y.size * sizeof (double));
2634 gridGenYvals(gridtype, grid->y.size, grid->y.first, grid->y.last, grid->y.inc, yvals);
2635 grid->y.vals = yvals;
2636 // gridDefYinc(gridID, grid->y.inc);
2637 }
2638
2639 if ( grid->projtype == CDI_PROJ_RLL )
2640 {
2641 const char *name = cdiInqVarKeyString(&grid->x.keys, CDI_KEY_NAME);
2642 if ( name[0] == 0 || name[0] == 'x' ) cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, "rlon");
2643 name = cdiInqVarKeyString(&grid->y.keys, CDI_KEY_NAME);
2644 if ( name[0] == 0 || name[0] == 'y' ) cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, "rlat");
2645 name = cdiInqVarKeyString(&grid->x.keys, CDI_KEY_LONGNAME);
2646 if ( name[0] == 0 ) cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_LONGNAME, "longitude in rotated pole grid");
2647 name = cdiInqVarKeyString(&grid->y.keys, CDI_KEY_LONGNAME);
2648 if ( name[0] == 0 ) cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_LONGNAME, "latitude in rotated pole grid");
2649 name = cdiInqVarKeyString(&grid->x.keys, CDI_KEY_UNITS);
2650 if ( name[0] == 0 ) cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, "degrees");
2651 name = cdiInqVarKeyString(&grid->y.keys, CDI_KEY_UNITS);
2652 if ( name[0] == 0 ) cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, "degrees");
2653 cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_grid_latlon][0]);
2654 cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_STDNAME, xystdname_tab[grid_xystdname_grid_latlon][1]);
2655 }
2656
2657 if (gridtype == GRID_UNSTRUCTURED)
2658 {
2659 int number = cdiInqVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDUSED);
2660 if (number > 0)
2661 {
2662 cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, number);
2663 int position = cdiInqVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE);
2664 if (position > 0) cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, position);
2665 }
2666 }
2667
2668 break;
2669 }
2670 case GRID_GAUSSIAN_REDUCED:
2671 {
2672 gridDefNP(gridID, grid->np);
2673 gridDefYsize(gridID, grid->y.size);
2674 if ( grid->x.flag == 2 )
2675 {
2676 double xvals[2] = { grid->x.first, grid->x.last };
2677 gridDefXsize(gridID, 2);
2678 gridDefXvals(gridID, xvals);
2679 }
2680
2681 if ( grid->y.flag == 2 )
2682 {
2683 double *yvals = (double *) Malloc(grid->y.size * sizeof (double));
2684 gridGenYvals(gridtype, grid->y.size, grid->y.first, grid->y.last, grid->y.inc, yvals);
2685 grid->y.vals = yvals;
2686 // gridDefYinc(gridID, grid->y.inc);
2687 }
2688 break;
2689 }
2690 case GRID_SPECTRAL:
2691 {
2692 gridDefTrunc(gridID, grid->trunc);
2693 if ( grid->lcomplex ) gridDefComplexPacking(gridID, 1);
2694 break;
2695 }
2696 case GRID_FOURIER:
2697 {
2698 gridDefTrunc(gridID, grid->trunc);
2699 break;
2700 }
2701 case GRID_GME:
2702 {
2703 gridDefParamGME(gridID, grid->gme.nd, grid->gme.ni, grid->gme.ni2, grid->gme.ni3);
2704 break;
2705 }
2706 /*
2707 case GRID_GENERIC:
2708 {
2709 if ( grid->x.size > 0 && grid->y.size > 0 )
2710 {
2711 gridDefXsize(gridID, grid->x.size);
2712 gridDefYsize(gridID, grid->y.size);
2713 if ( grid->x.vals ) gridDefXvals(gridID, grid->x.vals);
2714 if ( grid->y.vals ) gridDefYvals(gridID, grid->y.vals);
2715 }
2716 break;
2717 }
2718 */
2719 case GRID_TRAJECTORY:
2720 {
2721 gridDefXsize(gridID, 1);
2722 gridDefYsize(gridID, 1);
2723 break;
2724 }
2725 default:
2726 {
2727 Error("Gridtype %s unsupported!", gridNamePtr(gridtype));
2728 break;
2729 }
2730 }
2731 }
2732
2733 // Used only in iterator_grib.c
gridGenerate(const grid_t * grid)2734 int gridGenerate(const grid_t *grid)
2735 {
2736 int gridtype = grid->type;
2737 int gridID = gridCreate(gridtype, grid->size);
2738 grid_t *restrict gridptr = grid_to_pointer(gridID);
2739 cdiCopyVarKey(&grid->keys, CDI_KEY_DATATYPE, &gridptr->keys);
2740 gridptr->x.size = grid->x.size;
2741 gridptr->y.size = grid->y.size;
2742 gridptr->np = grid->np;
2743 gridptr->nvertex = grid->nvertex;
2744 gridptr->x.flag = grid->x.flag;
2745 int valdef_group1 = 0;
2746 static const int valdef_group1_tab[] = {
2747 GRID_LONLAT, GRID_GAUSSIAN, GRID_UNSTRUCTURED, GRID_CURVILINEAR, GRID_GENERIC, GRID_PROJECTION
2748 };
2749 for ( size_t i = 0; i < sizeof (valdef_group1_tab) / sizeof (valdef_group1_tab[0]); ++i)
2750 valdef_group1 |= (gridtype == valdef_group1_tab[i]);
2751 if ( valdef_group1 && grid->x.flag == 1 )
2752 {
2753 gridDefXvals(gridID, grid->x.vals);
2754 if ( grid->x.bounds )
2755 gridDefXbounds(gridID, grid->x.bounds);
2756 }
2757 gridptr->x.first = grid->x.first;
2758 gridptr->x.last = grid->x.last;
2759 gridptr->x.inc = grid->x.inc;
2760 gridptr->y.flag = grid->y.flag;
2761 if ( (valdef_group1 || gridtype == GRID_GAUSSIAN_REDUCED) && grid->y.flag == 1)
2762 {
2763 gridDefYvals(gridID, grid->y.vals);
2764 if ( grid->y.bounds )
2765 gridDefYbounds(gridID, grid->y.bounds);
2766 }
2767 gridptr->y.first = grid->y.first;
2768 gridptr->y.last = grid->y.last;
2769 gridptr->y.inc = grid->y.inc;
2770 if ( valdef_group1 && grid->area) gridDefArea(gridID, grid->area);
2771
2772 cdiCopyVarKey(&grid->keys, CDI_KEY_NUMBEROFGRIDUSED, &gridptr->keys);
2773 cdiCopyVarKey(&grid->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE, &gridptr->keys);
2774 cdiCopyVarKey(&grid->keys, CDI_KEY_REFERENCEURI, &gridptr->keys);
2775
2776 cdiCopyVarKey(&grid->keys, CDI_KEY_SCANNINGMODE, &gridptr->keys);
2777
2778 if (gridtype == GRID_PROJECTION) gridptr->name = strdup(grid->name);
2779 if (gridtype == GRID_GAUSSIAN_REDUCED) gridDefReducedPoints(gridID, grid->y.size, grid->reducedPoints);
2780 gridptr->trunc = grid->trunc;
2781 gridptr->lcomplex = grid->lcomplex;
2782 gridptr->gme.nd = grid->gme.nd;
2783 gridptr->gme.ni = grid->gme.ni;
2784 gridptr->gme.ni2 = grid->gme.ni2;
2785 gridptr->gme.ni3 = grid->gme.ni3;
2786
2787 gridComplete(gridptr);
2788
2789 return gridID;
2790 }
2791
2792 static void
grid_copy_base_array_fields(grid_t * gridptrOrig,grid_t * gridptrDup)2793 grid_copy_base_array_fields(grid_t *gridptrOrig, grid_t *gridptrDup)
2794 {
2795 size_t reducedPointsSize = (SizeType)gridptrOrig->reducedPointsSize;
2796 size_t gridsize = gridptrOrig->size;
2797 int gridtype = gridptrOrig->type;
2798 int irregular = gridtype == GRID_CURVILINEAR || gridtype == GRID_UNSTRUCTURED;
2799 if ( reducedPointsSize )
2800 {
2801 gridptrDup->reducedPoints = (int*) Malloc(reducedPointsSize * sizeof(int));
2802 memcpy(gridptrDup->reducedPoints, gridptrOrig->reducedPoints, reducedPointsSize * sizeof(int));
2803 }
2804
2805 if ( gridptrOrig->x.vals != NULL )
2806 {
2807 size_t size = irregular ? gridsize : gridptrOrig->x.size;
2808 gridptrDup->x.vals = (double *)Malloc(size * sizeof (double));
2809 memcpy(gridptrDup->x.vals, gridptrOrig->x.vals, size * sizeof (double));
2810 }
2811
2812 if ( gridptrOrig->y.vals != NULL )
2813 {
2814 size_t size = irregular ? gridsize : gridptrOrig->y.size;
2815 gridptrDup->y.vals = (double *)Malloc(size * sizeof (double));
2816 memcpy(gridptrDup->y.vals, gridptrOrig->y.vals, size * sizeof (double));
2817 }
2818
2819 if ( gridptrOrig->x.bounds != NULL )
2820 {
2821 size_t size = (irregular ? gridsize : gridptrOrig->x.size) * gridptrOrig->nvertex;
2822 gridptrDup->x.bounds = (double *)Malloc(size * sizeof (double));
2823 memcpy(gridptrDup->x.bounds, gridptrOrig->x.bounds, size * sizeof (double));
2824 }
2825
2826 if ( gridptrOrig->y.bounds != NULL )
2827 {
2828 size_t size = (irregular ? gridsize : gridptrOrig->y.size) * gridptrOrig->nvertex;
2829 gridptrDup->y.bounds = (double *)Malloc(size * sizeof (double));
2830 memcpy(gridptrDup->y.bounds, gridptrOrig->y.bounds, size * sizeof (double));
2831 }
2832
2833 {
2834 const double *gridptrOrig_area = gridptrOrig->vtable->inqAreaPtr(gridptrOrig);
2835 if ( gridptrOrig_area != NULL )
2836 {
2837 size_t size = gridsize;
2838 gridptrDup->area = (double *)Malloc(size * sizeof (double));
2839 memcpy(gridptrDup->area, gridptrOrig_area, size * sizeof (double));
2840 }
2841 }
2842
2843 if ( gridptrOrig->mask != NULL )
2844 {
2845 size_t size = gridsize;
2846 gridptrDup->mask = (mask_t *)Malloc(size * sizeof(mask_t));
2847 memcpy(gridptrDup->mask, gridptrOrig->mask, size * sizeof (mask_t));
2848 }
2849
2850 if ( gridptrOrig->mask_gme != NULL )
2851 {
2852 size_t size = gridsize;
2853 gridptrDup->mask_gme = (mask_t *)Malloc(size * sizeof (mask_t));
2854 memcpy(gridptrDup->mask_gme, gridptrOrig->mask_gme, size * sizeof(mask_t));
2855 }
2856 }
2857
2858
2859 /*
2860 @Function gridDuplicate
2861 @Title Duplicate a horizontal Grid
2862
2863 @Prototype int gridDuplicate(int gridID)
2864 @Parameter
2865 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
2866
2867 @Description
2868 The function @func{gridDuplicate} duplicates a horizontal Grid.
2869
2870 @Result
2871 @func{gridDuplicate} returns an identifier to the duplicated Grid.
2872
2873 @EndFunction
2874 */
gridDuplicate(int gridID)2875 int gridDuplicate(int gridID)
2876 {
2877 grid_t *gridptr = grid_to_pointer(gridID);
2878 grid_t *gridptrnew = gridptr->vtable->copy(gridptr);
2879 int gridIDnew = reshPut(gridptrnew, &gridOps);
2880 gridptrnew->self = gridIDnew;
2881 return gridIDnew;
2882 }
2883
2884
gridCompress(int gridID)2885 void gridCompress(int gridID)
2886 {
2887 grid_t *gridptr = grid_to_pointer(gridID);
2888
2889 int gridtype = gridInqType(gridID);
2890 if ( gridtype == GRID_UNSTRUCTURED )
2891 {
2892 if ( gridptr->mask_gme != NULL )
2893 {
2894 size_t gridsize = gridInqSize(gridID);
2895 size_t nv = (size_t)gridptr->nvertex;
2896 double *restrict area
2897 = (double *)gridptr->vtable->inqAreaPtr(gridptr),
2898 *restrict xvals = (double *)gridptr->vtable->inqXValsPtr(gridptr),
2899 *restrict yvals = (double *)gridptr->vtable->inqYValsPtr(gridptr),
2900 *restrict xbounds = (double *)gridptr->vtable->inqXBoundsPtr(gridptr),
2901 *restrict ybounds = (double *)gridptr->vtable->inqYBoundsPtr(gridptr);
2902 mask_t *restrict mask_gme = gridptr->mask_gme;
2903 size_t *restrict selection = (size_t *)Malloc(gridsize * sizeof (selection[0]));
2904 size_t nselect;
2905 {
2906 size_t j = 0;
2907 for (size_t i = 0; i < gridsize; i++ )
2908 selection[j] = i, j += (mask_gme[i] != 0);
2909 nselect = j;
2910 }
2911 selection = (size_t *)Realloc(selection, nselect * sizeof (selection[0]));
2912 if (xvals)
2913 for (size_t i = 0; i < nselect; i++ )
2914 xvals[i] = xvals[selection[i]];
2915 if (yvals)
2916 for (size_t i = 0; i < nselect; i++ )
2917 yvals[i] = yvals[selection[i]];
2918 if (area)
2919 for (size_t i = 0; i < nselect; i++ )
2920 area[i] = area[selection[i]];
2921 if (xbounds)
2922 for (size_t i = 0; i < nselect; i++ )
2923 for (size_t iv = 0; iv < nv; iv++)
2924 xbounds[i * nv + iv] = xbounds[selection[i] * nv + iv];
2925 if (ybounds)
2926 for (size_t i = 0; i < nselect; i++ )
2927 for (size_t iv = 0; iv < nv; iv++)
2928 ybounds[i * nv + iv] = ybounds[selection[i] * nv + iv];
2929 Free(selection);
2930
2931 /* fprintf(stderr, "grid compress %d %d %d\n", i, j, gridsize); */
2932 gridsize = nselect;
2933 gridptr->size = (int)gridsize;
2934 gridptr->x.size = (int)gridsize;
2935 gridptr->y.size = (int)gridsize;
2936
2937 double **resizeP[] = { &gridptr->x.vals, &gridptr->y.vals,
2938 &gridptr->area,
2939 &gridptr->x.bounds, &gridptr->y.bounds };
2940 size_t newSize[] = { gridsize, gridsize, gridsize, nv*gridsize,
2941 nv*gridsize };
2942 for ( size_t i = 0; i < sizeof (resizeP) / sizeof (resizeP[0]); ++i)
2943 if ( *(resizeP[i]) )
2944 *(resizeP[i]) = (double *)Realloc(*(resizeP[i]), newSize[i]*sizeof(double));
2945
2946 Free(gridptr->mask_gme);
2947 gridptr->mask_gme = NULL;
2948 gridMark4Update(gridID);
2949 }
2950 }
2951 else
2952 Warning("Unsupported grid type: %s", gridNamePtr(gridtype));
2953 }
2954
2955 static void
gridDefAreaSerial(grid_t * gridptr,const double * area)2956 gridDefAreaSerial(grid_t *gridptr, const double *area)
2957 {
2958 size_t size = gridptr->size;
2959
2960 if (size == 0)
2961 Error("size undefined for gridID = %d", gridptr->self);
2962
2963 if (gridptr->area == NULL)
2964 gridptr->area = (double *) Malloc(size*sizeof(double));
2965 else if (CDI_Debug)
2966 Warning("values already defined!");
2967
2968 memcpy(gridptr->area, area, size * sizeof(double));
2969 }
2970
2971
gridDefArea(int gridID,const double * area)2972 void gridDefArea(int gridID, const double *area)
2973 {
2974 grid_t *gridptr = grid_to_pointer(gridID);
2975 gridptr->vtable->defArea(gridptr, area);
2976 gridMark4Update(gridID);
2977 }
2978
2979 static void
gridInqAreaSerial(grid_t * gridptr,double * area)2980 gridInqAreaSerial(grid_t *gridptr, double *area)
2981 {
2982 if (gridptr->area)
2983 memcpy(area, gridptr->area, gridptr->size * sizeof (double));
2984 }
2985
2986
gridInqArea(int gridID,double * area)2987 void gridInqArea(int gridID, double *area)
2988 {
2989 grid_t *gridptr = grid_to_pointer(gridID);
2990 gridptr->vtable->inqArea(gridptr, area);
2991 }
2992
2993 static int
gridHasAreaBase(grid_t * gridptr)2994 gridHasAreaBase(grid_t *gridptr)
2995 {
2996 return gridptr->area != NULL;
2997 }
2998
gridHasArea(int gridID)2999 int gridHasArea(int gridID)
3000 {
3001 grid_t *gridptr = grid_to_pointer(gridID);
3002 return gridptr->vtable->hasArea(gridptr);
3003 }
3004
3005
gridInqAreaPtrBase(grid_t * gridptr)3006 static const double *gridInqAreaPtrBase(grid_t *gridptr)
3007 {
3008 return gridptr->area;
3009 }
3010
gridInqAreaPtr(int gridID)3011 const double *gridInqAreaPtr(int gridID)
3012 {
3013 grid_t *gridptr = grid_to_pointer(gridID);
3014 return gridptr->vtable->inqAreaPtr(gridptr);
3015 }
3016
3017
gridDefNvertex(int gridID,int nvertex)3018 void gridDefNvertex(int gridID, int nvertex)
3019 {
3020 grid_t *gridptr = grid_to_pointer(gridID);
3021
3022 if (gridptr->nvertex != nvertex)
3023 {
3024 gridptr->nvertex = nvertex;
3025 gridMark4Update(gridID);
3026 }
3027 }
3028
3029
gridInqNvertex(int gridID)3030 int gridInqNvertex(int gridID)
3031 {
3032 grid_t *gridptr = grid_to_pointer(gridID);
3033 return gridptr->nvertex;
3034 }
3035
3036 static void
gridDefBoundsGeneric(grid_t * gridptr,const double * bounds,size_t regularSize,double ** field)3037 gridDefBoundsGeneric(grid_t *gridptr, const double *bounds, size_t regularSize, double **field)
3038 {
3039 const int irregular = gridptr->type == GRID_CURVILINEAR || gridptr->type == GRID_UNSTRUCTURED;
3040 const size_t nvertex = (size_t)gridptr->nvertex;
3041 if (nvertex == 0)
3042 {
3043 Warning("nvertex undefined for gridID = %d. Cannot define bounds!", gridptr->self);
3044 return;
3045 }
3046
3047 const size_t size = nvertex * (irregular ? gridptr->size : regularSize);
3048 if (size == 0)
3049 Error("size undefined for gridID = %d", gridptr->self);
3050
3051 if (*field == NULL && size)
3052 *field = (double *)Malloc(size * sizeof (double));
3053 else if (CDI_Debug)
3054 Warning("values already defined!");
3055
3056 memcpy(*field, bounds, size * sizeof (double));
3057 }
3058
3059
3060 static void
gridDefXBoundsSerial(grid_t * gridptr,const double * xbounds)3061 gridDefXBoundsSerial(grid_t *gridptr, const double *xbounds)
3062 {
3063 gridDefBoundsGeneric(gridptr, xbounds, gridptr->x.size, &gridptr->x.bounds);
3064 }
3065
3066 /*
3067 @Function gridDefXbounds
3068 @Title Define the bounds of a X-axis
3069
3070 @Prototype void gridDefXbounds(int gridID, const double *xbounds)
3071 @Parameter
3072 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3073 @Item xbounds X-bounds of the grid.
3074
3075 @Description
3076 The function @func{gridDefXbounds} defines all bounds of the X-axis.
3077
3078 @EndFunction
3079 */
gridDefXbounds(int gridID,const double * xbounds)3080 void gridDefXbounds(int gridID, const double *xbounds)
3081 {
3082 grid_t *gridptr = grid_to_pointer(gridID);
3083 gridptr->vtable->defXBounds(gridptr, xbounds);
3084 gridMark4Update(gridID);
3085 }
3086
3087 static SizeType
gridInqXBoundsSerial(grid_t * gridptr,double * xbounds)3088 gridInqXBoundsSerial(grid_t *gridptr, double *xbounds)
3089 {
3090 const size_t nvertex = (size_t)gridptr->nvertex;
3091
3092 const int irregular = (gridptr->type == GRID_CURVILINEAR || gridptr->type == GRID_UNSTRUCTURED);
3093 size_t size = nvertex * (irregular ? gridptr->size : gridptr->x.size);
3094
3095 const double *gridptr_xbounds = gridptr->vtable->inqXBoundsPtr(gridptr);
3096 if (gridptr_xbounds)
3097 {
3098 if (size && xbounds) memcpy(xbounds, gridptr_xbounds, size * sizeof (double));
3099 }
3100 else
3101 size = 0;
3102
3103 return (SizeType)size;
3104 }
3105
3106 /*
3107 @Function gridInqXbounds
3108 @Title Get the bounds of a X-axis
3109
3110 @Prototype SizeType gridInqXbounds(int gridID, double *xbounds)
3111 @Parameter
3112 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
3113 @Item xbounds Pointer to the location into which the X-bounds are read.
3114 The caller must allocate space for the returned values.
3115
3116 @Description
3117 The function @func{gridInqXbounds} returns the bounds of the X-axis.
3118
3119 @Result
3120 Upon successful completion @func{gridInqXbounds} returns the number of bounds and
3121 the bounds are stored in @func{xbounds}.
3122 Otherwise, 0 is returned and @func{xbounds} is empty.
3123
3124 @EndFunction
3125 */
gridInqXbounds(int gridID,double * xbounds)3126 SizeType gridInqXbounds(int gridID, double *xbounds)
3127 {
3128 grid_t *gridptr = grid_to_pointer(gridID);
3129 return gridptr->vtable->inqXBounds(gridptr, xbounds);
3130 }
3131
3132 static const double *
gridInqXBoundsPtrSerial(grid_t * gridptr)3133 gridInqXBoundsPtrSerial(grid_t *gridptr)
3134 {
3135 return gridptr->x.bounds;
3136 }
3137
3138
gridInqXboundsPtr(int gridID)3139 const double *gridInqXboundsPtr(int gridID)
3140 {
3141 grid_t *gridptr = grid_to_pointer(gridID);
3142 return gridptr->vtable->inqXBoundsPtr(gridptr);
3143 }
3144
3145 static void
gridDefYBoundsSerial(grid_t * gridptr,const double * ybounds)3146 gridDefYBoundsSerial(grid_t *gridptr, const double *ybounds)
3147 {
3148 gridDefBoundsGeneric(gridptr, ybounds, gridptr->y.size, &gridptr->y.bounds);
3149 }
3150
3151 //----------------------------------------------------------------------------
3152 // Parallel Version
3153 //----------------------------------------------------------------------------
3154
gridInqXboundsPart(int gridID,int start,SizeType size,double * xbounds)3155 SizeType gridInqXboundsPart(int gridID, int start, SizeType size, double *xbounds)
3156 {
3157 grid_t *gridptr = grid_to_pointer(gridID);
3158 const double *gridptr_xbounds = gridptr->vtable->inqXBoundsPtr(gridptr);
3159 if (gridptr_xbounds && size && xbounds)
3160 memcpy(xbounds, gridptr_xbounds+start, size * sizeof (double));
3161
3162 return size;
3163 }
3164
gridInqYboundsPart(int gridID,int start,SizeType size,double * ybounds)3165 SizeType gridInqYboundsPart(int gridID, int start, SizeType size, double *ybounds)
3166 {
3167 grid_t *gridptr = grid_to_pointer(gridID);
3168 const double *gridptr_ybounds = gridptr->vtable->inqYBoundsPtr(gridptr);
3169 if (gridptr_ybounds && size && ybounds)
3170 memcpy(ybounds, gridptr_ybounds+start, size * sizeof (double));
3171
3172 return size;
3173 }
3174
3175 /*
3176 @Function gridDefYbounds
3177 @Title Define the bounds of a Y-axis
3178
3179 @Prototype void gridDefYbounds(int gridID, const double *ybounds)
3180 @Parameter
3181 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3182 @Item ybounds Y-bounds of the grid.
3183
3184 @Description
3185 The function @func{gridDefYbounds} defines all bounds of the Y-axis.
3186
3187 @EndFunction
3188 */
gridDefYbounds(int gridID,const double * ybounds)3189 void gridDefYbounds(int gridID, const double *ybounds)
3190 {
3191 grid_t *gridptr = grid_to_pointer(gridID);
3192 gridptr->vtable->defYBounds(gridptr, ybounds);
3193 gridMark4Update(gridID);
3194 }
3195
3196 static SizeType
gridInqYBoundsSerial(grid_t * gridptr,double * ybounds)3197 gridInqYBoundsSerial(grid_t *gridptr, double *ybounds)
3198 {
3199 const size_t nvertex = (size_t)gridptr->nvertex;
3200
3201 const int irregular = (gridptr->type == GRID_CURVILINEAR || gridptr->type == GRID_UNSTRUCTURED);
3202 size_t size = nvertex * (irregular ? gridptr->size : gridptr->y.size);
3203
3204 const double *gridptr_ybounds = gridptr->vtable->inqYBoundsPtr(gridptr);
3205 if (gridptr_ybounds)
3206 {
3207 if (size && ybounds) memcpy(ybounds, gridptr_ybounds, size * sizeof (double));
3208 }
3209 else
3210 size = 0;
3211
3212 return (SizeType)size;
3213 }
3214
3215
3216 /*
3217 @Function gridInqYbounds
3218 @Title Get the bounds of a Y-axis
3219
3220 @Prototype SizeType gridInqYbounds(int gridID, double *ybounds)
3221 @Parameter
3222 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
3223 @Item ybounds Pointer to the location into which the Y-bounds are read.
3224 The caller must allocate space for the returned values.
3225
3226 @Description
3227 The function @func{gridInqYbounds} returns the bounds of the Y-axis.
3228
3229 @Result
3230 Upon successful completion @func{gridInqYbounds} returns the number of bounds and
3231 the bounds are stored in @func{ybounds}.
3232 Otherwise, 0 is returned and @func{ybounds} is empty.
3233
3234 @EndFunction
3235 */
gridInqYbounds(int gridID,double * ybounds)3236 SizeType gridInqYbounds(int gridID, double *ybounds)
3237 {
3238 grid_t *gridptr = grid_to_pointer(gridID);
3239 return gridptr->vtable->inqYBounds(gridptr, ybounds);
3240 }
3241
3242 static const double *
gridInqYBoundsPtrSerial(grid_t * gridptr)3243 gridInqYBoundsPtrSerial(grid_t *gridptr)
3244 {
3245 return gridptr->y.bounds;
3246 }
3247
3248
gridInqYboundsPtr(int gridID)3249 const double *gridInqYboundsPtr(int gridID)
3250 {
3251 grid_t *gridptr = grid_to_pointer(gridID);
3252 return gridptr->vtable->inqYBoundsPtr(gridptr);
3253 }
3254
3255 static void
printDblsPrefixAutoBrk(FILE * fp,int dig,const char prefix[],size_t nbyte0,size_t n,const double vals[])3256 printDblsPrefixAutoBrk(FILE *fp, int dig, const char prefix[], size_t nbyte0,
3257 size_t n, const double vals[])
3258 {
3259 fputs(prefix, fp);
3260 size_t nbyte = nbyte0;
3261 for ( size_t i = 0; i < n; i++ )
3262 {
3263 if ( nbyte > 80 )
3264 {
3265 fprintf(fp, "\n%*s", (int)nbyte0, "");
3266 nbyte = nbyte0;
3267 }
3268 nbyte += (size_t)fprintf(fp, "%.*g ", dig, vals[i]);
3269 }
3270 fputs("\n", fp);
3271 }
3272
3273 static inline
resizeBuffer(void ** buf,size_t * bufSize,size_t reqSize)3274 void *resizeBuffer(void **buf, size_t *bufSize, size_t reqSize)
3275 {
3276 if (reqSize > *bufSize)
3277 {
3278 *buf = Realloc(*buf, reqSize);
3279 *bufSize = reqSize;
3280 }
3281 return *buf;
3282 }
3283
3284 static
gridPrintAttributes(FILE * fp,int gridID)3285 void gridPrintAttributes(FILE *fp, int gridID)
3286 {
3287 int cdiID = gridID;
3288 int varID = CDI_GLOBAL;
3289 int atttype, attlen;
3290 char attname[CDI_MAX_NAME+1];
3291 void *attBuf = NULL;
3292 size_t attBufSize = 0;
3293
3294 int natts;
3295 cdiInqNatts(cdiID, varID, &natts);
3296
3297 for ( int iatt = 0; iatt < natts; ++iatt )
3298 {
3299 cdiInqAtt(cdiID, varID, iatt, attname, &atttype, &attlen);
3300
3301 if ( attlen == 0 ) continue;
3302
3303 if ( atttype == CDI_DATATYPE_TXT )
3304 {
3305 size_t attSize = (size_t)(attlen+1)*sizeof(char);
3306 char *atttxt = (char *)resizeBuffer(&attBuf, &attBufSize, attSize);
3307 cdiInqAttTxt(cdiID, varID, attname, attlen, atttxt);
3308 atttxt[attlen] = 0;
3309 fprintf(fp, "ATTR_TXT: %s = \"%s\"\n", attname, atttxt);
3310 }
3311 else if ( atttype == CDI_DATATYPE_INT8 || atttype == CDI_DATATYPE_UINT8 ||
3312 atttype == CDI_DATATYPE_INT16 || atttype == CDI_DATATYPE_UINT16 ||
3313 atttype == CDI_DATATYPE_INT32 || atttype == CDI_DATATYPE_UINT32 )
3314 {
3315 size_t attSize = (size_t)attlen*sizeof(int);
3316 int *attint = (int *)resizeBuffer(&attBuf, &attBufSize, attSize);
3317 cdiInqAttInt(cdiID, varID, attname, attlen, &attint[0]);
3318 if ( attlen == 1 )
3319 fprintf(fp, "ATTR_INT: %s =", attname);
3320 else
3321 fprintf(fp, "ATTR_INT_%d: %s =", attlen, attname);
3322 for ( int i = 0; i < attlen; ++i ) fprintf(fp, " %d", attint[i]);
3323 fprintf(fp, "\n");
3324 }
3325 else if ( atttype == CDI_DATATYPE_FLT32 || atttype == CDI_DATATYPE_FLT64 )
3326 {
3327 size_t attSize = (size_t)attlen * sizeof(double);
3328 double *attflt = (double *)resizeBuffer(&attBuf, &attBufSize, attSize);
3329 int dig = (atttype == CDI_DATATYPE_FLT64) ? 15 : 7;
3330 cdiInqAttFlt(cdiID, varID, attname, attlen, attflt);
3331 if ( attlen == 1 )
3332 fprintf(fp, "ATTR_FLT: %s =", attname);
3333 else
3334 fprintf(fp, "ATTR_FLT_%d: %s =", attlen, attname);
3335 for ( int i = 0; i < attlen; ++i ) fprintf(fp, " %.*g", dig, attflt[i]);
3336 fprintf(fp, "\n");
3337 }
3338 }
3339
3340 Free(attBuf);
3341 }
3342
3343 static
gridPrintKernel(int gridID,int opt,FILE * fp)3344 void gridPrintKernel(int gridID, int opt, FILE *fp)
3345 {
3346 char attstr[CDI_MAX_NAME];
3347 char attstr2[CDI_MAX_NAME];
3348 size_t nxvals = gridInqXvals(gridID, NULL);
3349 size_t nyvals = gridInqYvals(gridID, NULL);
3350
3351 int type = gridInqType(gridID);
3352 size_t gridsize = gridInqSize(gridID);
3353 size_t xsize = gridInqXsize(gridID);
3354 size_t ysize = gridInqYsize(gridID);
3355 int nvertex = gridInqNvertex(gridID);
3356 int datatype;
3357 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
3358
3359 int dig = (datatype == CDI_DATATYPE_FLT64) ? 15 : 7;
3360
3361 fprintf(fp, "gridtype = %s\n" "gridsize = %zu\n", gridNamePtr(type), gridsize);
3362
3363 if ( type != GRID_GME )
3364 {
3365 if ( type != GRID_UNSTRUCTURED && type != GRID_SPECTRAL && type != GRID_FOURIER )
3366 {
3367 if ( xsize > 0 ) fprintf(fp, "xsize = %zu\n", xsize);
3368 if ( ysize > 0 ) fprintf(fp, "ysize = %zu\n", ysize);
3369 }
3370
3371 if ( nxvals > 0 )
3372 {
3373 int length = CDI_MAX_NAME;
3374 cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_NAME, attstr, &length);
3375 if ( attstr[0] ) fprintf(fp, "xname = %s\n", attstr);
3376 length = CDI_MAX_NAME;
3377 cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_DIMNAME, attstr2, &length);
3378 if ( attstr2[0] && strcmp(attstr, attstr2) ) fprintf(fp, "xdimname = %s\n", attstr2);
3379 length = CDI_MAX_NAME;
3380 cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_LONGNAME, attstr, &length);
3381 if ( attstr[0] ) fprintf(fp, "xlongname = %s\n", attstr);
3382 length = CDI_MAX_NAME;
3383 cdiInqKeyString(gridID, CDI_XAXIS, CDI_KEY_UNITS, attstr, &length);
3384 if ( attstr[0] ) fprintf(fp, "xunits = %s\n", attstr);
3385 }
3386
3387 if ( nyvals > 0 )
3388 {
3389 int length = CDI_MAX_NAME;
3390 cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_NAME, attstr, &length);
3391 if ( attstr[0] ) fprintf(fp, "yname = %s\n", attstr);
3392 length = CDI_MAX_NAME;
3393 cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_DIMNAME, attstr2, &length);
3394 if ( attstr2[0] && strcmp(attstr, attstr2) ) fprintf(fp, "ydimname = %s\n", attstr2);
3395 length = CDI_MAX_NAME;
3396 cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_LONGNAME, attstr, &length);
3397 if ( attstr[0] ) fprintf(fp, "ylongname = %s\n", attstr);
3398 length = CDI_MAX_NAME;
3399 cdiInqKeyString(gridID, CDI_YAXIS, CDI_KEY_UNITS, attstr, &length);
3400 if ( attstr[0] ) fprintf(fp, "yunits = %s\n", attstr);
3401 }
3402
3403 if ( type == GRID_UNSTRUCTURED && nvertex > 0 ) fprintf(fp, "nvertex = %d\n", nvertex);
3404 }
3405
3406 switch (type)
3407 {
3408 case GRID_LONLAT:
3409 case GRID_GAUSSIAN:
3410 case GRID_GAUSSIAN_REDUCED:
3411 case GRID_GENERIC:
3412 case GRID_PROJECTION:
3413 case GRID_CURVILINEAR:
3414 case GRID_UNSTRUCTURED:
3415 case GRID_CHARXY:
3416 {
3417 if ( type == GRID_GAUSSIAN || type == GRID_GAUSSIAN_REDUCED ) fprintf(fp, "np = %d\n", gridInqNP(gridID));
3418
3419 if ( type == GRID_UNSTRUCTURED )
3420 {
3421 int number = 0;
3422 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
3423 if ( number > 0 )
3424 {
3425 fprintf(fp, "number = %d\n", number);
3426 int position = 0;
3427 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, &position);
3428 if ( position >= 0 ) fprintf(fp, "position = %d\n", position);
3429 }
3430
3431 int length;
3432 if (CDI_NOERR == cdiInqKeyLen(gridID, CDI_GLOBAL, CDI_KEY_REFERENCEURI, &length))
3433 {
3434 char reference_link[8192];
3435 length = sizeof(reference_link);
3436 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_REFERENCEURI, reference_link, &length);
3437 fprintf(fp, "uri = %s\n", reference_link);
3438 }
3439 }
3440
3441 if ( nxvals > 0 )
3442 {
3443 double xfirst = 0.0, xinc = 0.0;
3444
3445 if ( type == GRID_LONLAT || type == GRID_GAUSSIAN ||
3446 type == GRID_PROJECTION || type == GRID_GENERIC )
3447 {
3448 xfirst = gridInqXval(gridID, 0);
3449 xinc = gridInqXinc(gridID);
3450 }
3451
3452 if ( IS_NOT_EQUAL(xinc, 0) && opt )
3453 {
3454 fprintf(fp, "xfirst = %.*g\n"
3455 "xinc = %.*g\n", dig, xfirst, dig, xinc);
3456 }
3457 else
3458 {
3459 double *xvals = (double*) Malloc(nxvals*sizeof(double));
3460 gridInqXvals(gridID, xvals);
3461 static const char prefix[] = "xvals = ";
3462 printDblsPrefixAutoBrk(fp, dig, prefix, sizeof(prefix)-1, nxvals, xvals);
3463 Free(xvals);
3464 }
3465 }
3466
3467 if ( nyvals > 0 )
3468 {
3469 double yfirst = 0.0, yinc = 0.0;
3470
3471 if ( type == GRID_LONLAT || type == GRID_GENERIC ||
3472 type == GRID_PROJECTION || type == GRID_GENERIC )
3473 {
3474 yfirst = gridInqYval(gridID, 0);
3475 yinc = gridInqYinc(gridID);
3476 }
3477
3478 if ( IS_NOT_EQUAL(yinc, 0) && opt )
3479 {
3480 fprintf(fp, "yfirst = %.*g\n"
3481 "yinc = %.*g\n", dig, yfirst, dig, yinc);
3482 }
3483 else
3484 {
3485 double *yvals = (double*) Malloc(nyvals*sizeof(double));
3486 gridInqYvals(gridID, yvals);
3487 static const char prefix[] = "yvals = ";
3488 printDblsPrefixAutoBrk(fp, dig, prefix, sizeof(prefix)-1, nyvals, yvals);
3489 Free(yvals);
3490 }
3491 }
3492
3493 if ( type == GRID_PROJECTION ) gridPrintAttributes(fp, gridID);
3494
3495 break;
3496 }
3497 case GRID_SPECTRAL:
3498 {
3499 fprintf(fp, "truncation = %d\n"
3500 "complexpacking = %d\n", gridInqTrunc(gridID), gridInqComplexPacking(gridID) );
3501 break;
3502 }
3503 case GRID_FOURIER:
3504 {
3505 fprintf(fp, "truncation = %d\n", gridInqTrunc(gridID));
3506 break;
3507 }
3508 case GRID_GME:
3509 {
3510 int nd, ni, ni2, ni3;
3511 gridInqParamGME(gridID, &nd, &ni, &ni2, &ni3);
3512 fprintf(fp, "ni = %d\n", ni );
3513 break;
3514 }
3515 default:
3516 {
3517 fprintf(stderr, "Unsupported grid type: %s\n", gridNamePtr(type));
3518 break;
3519 }
3520 }
3521 }
3522
3523
gridPrintP(void * voidptr,FILE * fp)3524 void gridPrintP(void *voidptr, FILE *fp)
3525 {
3526 grid_t *gridptr = (grid_t *) voidptr;
3527 int gridID = gridptr->self;
3528
3529 xassert( gridptr );
3530
3531 gridPrintKernel(gridID, 0, fp);
3532
3533 int datatype = CDI_UNDEFID;
3534 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
3535
3536 fprintf(fp,
3537 "datatype = %d\n"
3538 "nd = %d\n"
3539 "ni = %d\n"
3540 "ni2 = %d\n"
3541 "ni3 = %d\n"
3542 "trunc = %d\n"
3543 "lcomplex = %d\n"
3544 "reducedPointsSize = %d\n",
3545 datatype, gridptr->gme.nd, gridptr->gme.ni, gridptr->gme.ni2,
3546 gridptr->gme.ni3, gridptr->trunc,
3547 gridptr->lcomplex, gridptr->reducedPointsSize );
3548 }
3549
gridInqXValsPtrSerial(grid_t * gridptr)3550 static const double *gridInqXValsPtrSerial(grid_t *gridptr)
3551 {
3552 return gridptr->x.vals;
3553 }
3554
3555 #ifndef USE_MPI
gridInqXCvalsPtrSerial(grid_t * gridptr)3556 static const char **gridInqXCvalsPtrSerial(grid_t *gridptr)
3557 {
3558 return (const char **) gridptr->x.cvals;
3559 }
3560 #endif
3561
gridInqXvalsPtr(int gridID)3562 const double *gridInqXvalsPtr(int gridID)
3563 {
3564 grid_t *gridptr = grid_to_pointer(gridID);
3565 return gridptr->vtable->inqXValsPtr(gridptr);
3566 }
3567
3568 #ifndef USE_MPI
gridInqXCvalsPtr(int gridID)3569 const char **gridInqXCvalsPtr(int gridID)
3570 {
3571 grid_t *gridptr = grid_to_pointer(gridID);
3572 return gridptr->vtable->inqXCvalsPtr(gridptr);
3573 }
3574 #endif
3575
gridInqYValsPtrSerial(grid_t * gridptr)3576 static const double *gridInqYValsPtrSerial(grid_t *gridptr)
3577 {
3578 return gridptr->y.vals;
3579 }
3580
3581 #ifndef USE_MPI
gridInqYCvalsPtrSerial(grid_t * gridptr)3582 static const char **gridInqYCvalsPtrSerial(grid_t *gridptr)
3583 {
3584 return (const char **) gridptr->y.cvals;
3585 }
3586 #endif
3587
gridInqYvalsPtr(int gridID)3588 const double *gridInqYvalsPtr(int gridID)
3589 {
3590 grid_t *gridptr = grid_to_pointer(gridID);
3591 return gridptr->vtable->inqYValsPtr(gridptr);
3592 }
3593
3594 #ifndef USE_MPI
gridInqYCvalsPtr(int gridID)3595 const char **gridInqYCvalsPtr(int gridID)
3596 {
3597 grid_t *gridptr = grid_to_pointer(gridID);
3598 return gridptr->vtable->inqYCvalsPtr(gridptr);
3599 }
3600 #endif
3601
3602
gridProjParamsInit(struct CDI_GridProjParams * gpp)3603 void gridProjParamsInit(struct CDI_GridProjParams *gpp)
3604 {
3605 // clang-format off
3606 gpp->mv = CDI_Grid_Missval; // Missing value
3607 gpp->lon_0 = CDI_Grid_Missval; // The East longitude of the meridian which is parallel to the Y-axis
3608 gpp->lat_0 = CDI_Grid_Missval; // Latitude of the projection origin
3609 gpp->lat_1 = CDI_Grid_Missval; // First latitude from the pole at which the secant cone cuts the sphere
3610 gpp->lat_2 = CDI_Grid_Missval; // Second latitude at which the secant cone cuts the sphere
3611 gpp->a = CDI_Grid_Missval; // Semi-major axis or earth radius in metres (optional)
3612 gpp->b = CDI_Grid_Missval; // Semi-minor axis in metres (optional)
3613 gpp->rf = CDI_Grid_Missval; // Inverse flattening (1/f) (optional)
3614 gpp->xval_0 = CDI_Grid_Missval; // Longitude of the first grid point in degree (optional)
3615 gpp->yval_0 = CDI_Grid_Missval; // Latitude of the first grid point in degree (optional)
3616 gpp->x_0 = CDI_Grid_Missval; // False easting (optional)
3617 gpp->y_0 = CDI_Grid_Missval; // False northing (optional)
3618 // clang-format on
3619 }
3620
3621 static
gridDefParamsCommon(int gridID,struct CDI_GridProjParams gpp)3622 void gridDefParamsCommon(int gridID, struct CDI_GridProjParams gpp)
3623 {
3624 if (IS_NOT_EQUAL(gpp.a, gpp.mv))
3625 {
3626 if (IS_NOT_EQUAL(gpp.b, gpp.mv))
3627 {
3628 cdiDefAttFlt(gridID, CDI_GLOBAL, "semi_major_axis", CDI_DATATYPE_FLT64, 1, &gpp.a);
3629 cdiDefAttFlt(gridID, CDI_GLOBAL, "semi_minor_axis", CDI_DATATYPE_FLT64, 1, &gpp.b);
3630 }
3631 else
3632 {
3633 cdiDefAttFlt(gridID, CDI_GLOBAL, "earth_radius", CDI_DATATYPE_FLT64, 1, &gpp.a);
3634 }
3635 }
3636 if (IS_NOT_EQUAL(gpp.rf, gpp.mv)) cdiDefAttFlt(gridID, CDI_GLOBAL, "inverse_flattening", CDI_DATATYPE_FLT64, 1, &gpp.rf);
3637 if (IS_NOT_EQUAL(gpp.x_0, gpp.mv)) cdiDefAttFlt(gridID, CDI_GLOBAL, "false_easting", CDI_DATATYPE_FLT64, 1, &gpp.x_0);
3638 if (IS_NOT_EQUAL(gpp.y_0, gpp.mv)) cdiDefAttFlt(gridID, CDI_GLOBAL, "false_northing", CDI_DATATYPE_FLT64, 1, &gpp.y_0);
3639 if (IS_NOT_EQUAL(gpp.xval_0, gpp.mv)) cdiDefAttFlt(gridID, CDI_GLOBAL, "longitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.xval_0);
3640 if (IS_NOT_EQUAL(gpp.yval_0, gpp.mv)) cdiDefAttFlt(gridID, CDI_GLOBAL, "latitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.yval_0);
3641 }
3642
3643 /*
3644 @Function gridDefParamsLCC
3645 @Title Define the parameters of a Lambert Conformal Conic grid
3646
3647 @Prototype void gridDefParamsLCC(int gridID, struct CDI_GridProjParams gridProjParams)
3648 @Parameter
3649 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3650 @Item gridProjParams Grid projection parameters.
3651
3652 @Description
3653 The function @func{gridDefParamsLCC} defines the parameters of a Lambert Conformal Conic grid.
3654
3655 @EndFunction
3656 */
gridDefParamsLCC(int gridID,struct CDI_GridProjParams gpp)3657 void gridDefParamsLCC(int gridID, struct CDI_GridProjParams gpp)
3658 {
3659 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARNAME, "Lambert_Conformal");
3660
3661 const char *gmapname = "lambert_conformal_conic";
3662 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname);
3663 cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", (int)(strlen(gmapname)), gmapname);
3664 int nlats = 0;
3665 double lats[2];
3666 lats[nlats++] = gpp.lat_1;
3667 if (IS_NOT_EQUAL(gpp.lat_1, gpp.lat_2)) lats[nlats++] = gpp.lat_2;
3668 cdiDefAttFlt(gridID, CDI_GLOBAL, "standard_parallel", CDI_DATATYPE_FLT64, nlats, lats);
3669 cdiDefAttFlt(gridID, CDI_GLOBAL, "longitude_of_central_meridian", CDI_DATATYPE_FLT64, 1, &gpp.lon_0);
3670 cdiDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", CDI_DATATYPE_FLT64, 1, &gpp.lat_0);
3671
3672 gridDefParamsCommon(gridID, gpp);
3673
3674 grid_t *gridptr = grid_to_pointer(gridID);
3675 gridptr->projtype = CDI_PROJ_LCC;
3676
3677 if (gridptr->type != GRID_PROJECTION) gridptr->type = GRID_PROJECTION;
3678
3679 gridVerifyProj(gridID);
3680 }
3681
3682 /*
3683 @Function gridInqParamsLCC
3684 @Title Get the parameter of a Lambert Conformal Conic grid
3685
3686 @Prototype void gridInqParamsLCC(int gridID, struct CDI_GridProjParams *gpp)
3687 @Parameter
3688 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
3689 @Item gridProjParams Grid projection parameters.
3690
3691 @Description
3692 The function @func{gridInqParamsLCC} returns the parameter of a Lambert Conformal Conic grid.
3693
3694 @EndFunction
3695 */
gridInqParamsLCC(int gridID,struct CDI_GridProjParams * gpp)3696 int gridInqParamsLCC(int gridID, struct CDI_GridProjParams *gpp)
3697 {
3698 int status = -1;
3699 if (gridInqType(gridID) != GRID_PROJECTION) return status;
3700
3701 gridProjParamsInit(gpp);
3702
3703 status = -2;
3704 const char *projection = "lambert_conformal_conic";
3705 char gmapname[CDI_MAX_NAME];
3706 int length = CDI_MAX_NAME;
3707 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname, &length);
3708 if (gmapname[0] && strIsEqual(gmapname, projection))
3709 {
3710 char attname[CDI_MAX_NAME+1];
3711
3712 int natts;
3713 cdiInqNatts(gridID, CDI_GLOBAL, &natts);
3714
3715 if (natts) status = 0;
3716
3717 for (int iatt = 0; iatt < natts; ++iatt)
3718 {
3719 int atttype, attlen;
3720 cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
3721 if (attlen > 2) continue;
3722
3723 double attflt[2];
3724 if (cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, attflt))
3725 {
3726 // clang-format off
3727 if (strIsEqual(attname, "earth_radius")) gpp->a = attflt[0];
3728 else if (strIsEqual(attname, "semi_major_axis")) gpp->a = attflt[0];
3729 else if (strIsEqual(attname, "semi_minor_axis")) gpp->b = attflt[0];
3730 else if (strIsEqual(attname, "inverse_flattening")) gpp->rf = attflt[0];
3731 else if (strIsEqual(attname, "longitude_of_central_meridian")) gpp->lon_0 = attflt[0];
3732 else if (strIsEqual(attname, "latitude_of_projection_origin")) gpp->lat_0 = attflt[0];
3733 else if (strIsEqual(attname, "false_easting")) gpp->x_0 = attflt[0];
3734 else if (strIsEqual(attname, "false_northing")) gpp->y_0 = attflt[0];
3735 else if (strIsEqual(attname, "longitudeOfFirstGridPointInDegrees")) gpp->xval_0 = attflt[0];
3736 else if (strIsEqual(attname, "latitudeOfFirstGridPointInDegrees")) gpp->yval_0 = attflt[0];
3737 else if (strIsEqual(attname, "standard_parallel"))
3738 {
3739 gpp->lat_1 = attflt[0];
3740 gpp->lat_2 = (attlen == 2) ? attflt[1] : attflt[0];
3741 }
3742 // clang-format on
3743 }
3744 }
3745 }
3746
3747 return status;
3748 }
3749
3750
gridVerifyProjParamsLCC(struct CDI_GridProjParams * gpp)3751 int gridVerifyProjParamsLCC(struct CDI_GridProjParams *gpp)
3752 {
3753 static bool lwarn = true;
3754
3755 if (lwarn)
3756 {
3757 // lwarn = false;
3758 const char *projection = "lambert_conformal_conic";
3759 if (IS_EQUAL(gpp->lon_0, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "longitude_of_central_meridian"); }
3760 if (IS_EQUAL(gpp->lat_0, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "latitude_of_central_meridian"); }
3761 if (IS_EQUAL(gpp->lat_1, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "standard_parallel"); }
3762 if (IS_NOT_EQUAL(gpp->x_0, gpp->mv) && IS_NOT_EQUAL(gpp->y_0, gpp->mv) && (IS_EQUAL(gpp->xval_0, gpp->mv) || IS_EQUAL(gpp->yval_0, gpp->mv)))
3763 {
3764 if (proj_lcc_to_lonlat_func)
3765 {
3766 gpp->xval_0 = -gpp->x_0; gpp->yval_0 = -gpp->y_0;
3767 proj_lcc_to_lonlat_func(gpp, 0.0, 0.0, (SizeType)1, &gpp->xval_0, &gpp->yval_0);
3768 }
3769 if (IS_EQUAL(gpp->xval_0, gpp->mv) || IS_EQUAL(gpp->yval_0, gpp->mv))
3770 Warning("%s mapping parameter %s missing!", projection, "longitudeOfFirstGridPointInDegrees and latitudeOfFirstGridPointInDegrees");
3771 }
3772 }
3773
3774 return 0;
3775 }
3776
3777
gridVerifyProjParamsSTERE(struct CDI_GridProjParams * gpp)3778 int gridVerifyProjParamsSTERE(struct CDI_GridProjParams *gpp)
3779 {
3780 static bool lwarn = true;
3781
3782 if (lwarn)
3783 {
3784 // lwarn = false;
3785 const char *projection = "lambert_conformal_conic";
3786 if (IS_EQUAL(gpp->lon_0, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "straight_vertical_longitude_from_pole"); }
3787 if (IS_EQUAL(gpp->lat_0, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "latitude_of_projection_origin"); }
3788 if (IS_EQUAL(gpp->lat_1, gpp->mv)) { Warning("%s mapping parameter %s missing!", projection, "standard_parallel"); }
3789 if (IS_NOT_EQUAL(gpp->x_0, gpp->mv) && IS_NOT_EQUAL(gpp->y_0, gpp->mv) && (IS_EQUAL(gpp->xval_0, gpp->mv) || IS_EQUAL(gpp->yval_0, gpp->mv)))
3790 {
3791 if (proj_stere_to_lonlat_func)
3792 {
3793 gpp->xval_0 = -gpp->x_0; gpp->xval_0 = -gpp->y_0;
3794 proj_stere_to_lonlat_func(gpp, 0.0, 0.0, (SizeType)1, &gpp->xval_0, &gpp->yval_0);
3795 }
3796 if (IS_EQUAL(gpp->xval_0, gpp->mv) || IS_EQUAL(gpp->yval_0, gpp->mv))
3797 Warning("%s mapping parameter %s missing!", projection, "longitudeOfFirstGridPointInDegrees and latitudeOfFirstGridPointInDegrees");
3798 }
3799 }
3800
3801 return 0;
3802 }
3803
3804 /*
3805 @Function gridDefParamsSTERE
3806 @Title Define the parameter of a Polar stereographic grid
3807
3808 @Prototype void gridDefParamsSTERE(int gridID, struct CDI_GridProjParams gridProjParams)
3809 @Parameter
3810 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3811 @Item gridProjParams Grid projection parameters.
3812
3813 @Description
3814 The function @func{gridDefParamsSTERE} defines the parameter of a Polar stereographic grid.
3815
3816 @EndFunction
3817 */
gridDefParamsSTERE(int gridID,struct CDI_GridProjParams gpp)3818 void gridDefParamsSTERE(int gridID, struct CDI_GridProjParams gpp)
3819 {
3820 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARNAME, "Polar_Stereographic");
3821
3822 const char *gmapname = "polar_stereographic";
3823 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname);
3824 cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", (int)(strlen(gmapname)), gmapname);
3825 cdiDefAttFlt(gridID, CDI_GLOBAL, "standard_parallel", CDI_DATATYPE_FLT64, 1, &gpp.lat_1);
3826 cdiDefAttFlt(gridID, CDI_GLOBAL, "straight_vertical_longitude_from_pole", CDI_DATATYPE_FLT64, 1, &gpp.lon_0);
3827 cdiDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", CDI_DATATYPE_FLT64, 1, &gpp.lat_0);
3828
3829 gridDefParamsCommon(gridID, gpp);
3830
3831 grid_t *gridptr = grid_to_pointer(gridID);
3832 gridptr->projtype = CDI_PROJ_STERE;
3833
3834 gridVerifyProj(gridID);
3835 }
3836
3837 /*
3838 @Function gridInqParamsSTERE
3839 @Title Get the parameter of a Polar stereographic grid
3840
3841 @Prototype void gridInqParamsSTERE(int gridID, struct CDI_GridProjParams *gpp)
3842 @Parameter
3843 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
3844 @Item gridProjParams Grid projection parameters.
3845
3846 @Description
3847 The function @func{gridInqParamsSTERE} returns the parameter of a Polar stereographic grid.
3848
3849 @EndFunction
3850 */
gridInqParamsSTERE(int gridID,struct CDI_GridProjParams * gpp)3851 int gridInqParamsSTERE(int gridID, struct CDI_GridProjParams *gpp)
3852 {
3853 int status = -1;
3854 if ( gridInqType(gridID) != GRID_PROJECTION ) return status;
3855
3856 gridProjParamsInit(gpp);
3857
3858 status = -2;
3859 const char *projection = "polar_stereographic";
3860 char gmapname[CDI_MAX_NAME];
3861 int length = CDI_MAX_NAME;
3862 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname, &length);
3863 if (gmapname[0] && strIsEqual(gmapname, projection))
3864 {
3865 int atttype, attlen;
3866 char attname[CDI_MAX_NAME+1];
3867
3868 int natts;
3869 cdiInqNatts(gridID, CDI_GLOBAL, &natts);
3870
3871 if (natts) status = 0;
3872
3873 for (int iatt = 0; iatt < natts; ++iatt)
3874 {
3875 cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
3876 if (attlen > 2) continue;
3877
3878 double attflt[2];
3879 if (cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, attflt))
3880 {
3881 if (strIsEqual(attname, "earth_radius")) gpp->a = attflt[0];
3882 else if (strIsEqual(attname, "semi_major_axis")) gpp->a = attflt[0];
3883 else if (strIsEqual(attname, "semi_minor_axis")) gpp->b = attflt[0];
3884 else if (strIsEqual(attname, "inverse_flattening")) gpp->rf = attflt[0];
3885 else if (strIsEqual(attname, "standard_parallel")) gpp->lat_1 = attflt[0];
3886 else if (strIsEqual(attname, "straight_vertical_longitude_from_pole")) gpp->lon_0 = attflt[0];
3887 else if (strIsEqual(attname, "latitude_of_projection_origin")) gpp->lat_0 = attflt[0];
3888 else if (strIsEqual(attname, "false_easting")) gpp->x_0 = attflt[0];
3889 else if (strIsEqual(attname, "false_northing")) gpp->y_0 = attflt[0];
3890 else if (strIsEqual(attname, "longitudeOfFirstGridPointInDegrees")) gpp->xval_0 = attflt[0];
3891 else if (strIsEqual(attname, "latitudeOfFirstGridPointInDegrees")) gpp->yval_0 = attflt[0];
3892 }
3893 }
3894 }
3895
3896 return status;
3897 }
3898
3899
gridDefComplexPacking(int gridID,int lcomplex)3900 void gridDefComplexPacking(int gridID, int lcomplex)
3901 {
3902 grid_t *gridptr = grid_to_pointer(gridID);
3903
3904 if (gridptr->lcomplex != lcomplex)
3905 {
3906 gridptr->lcomplex = lcomplex != 0;
3907 gridMark4Update(gridID);
3908 }
3909 }
3910
3911
gridInqComplexPacking(int gridID)3912 int gridInqComplexPacking(int gridID)
3913 {
3914 grid_t* gridptr = grid_to_pointer(gridID);
3915
3916 return (int)gridptr->lcomplex;
3917 }
3918
3919
gridDefHasDims(int gridID,int hasdims)3920 void gridDefHasDims(int gridID, int hasdims)
3921 {
3922 grid_t* gridptr = grid_to_pointer(gridID);
3923
3924 if ( gridptr->hasdims != (hasdims != 0) )
3925 {
3926 gridptr->hasdims = hasdims != 0;
3927 gridMark4Update(gridID);
3928 }
3929 }
3930
3931
gridInqHasDims(int gridID)3932 int gridInqHasDims(int gridID)
3933 {
3934 grid_t* gridptr = grid_to_pointer(gridID);
3935
3936 return (int)gridptr->hasdims;
3937 }
3938
3939 /*
3940 @Function gridDefNumber
3941 @Title Define the reference number for an unstructured grid
3942
3943 @Prototype void gridDefNumber(int gridID, const int number)
3944 @Parameter
3945 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3946 @Item number Reference number for an unstructured grid.
3947
3948 @Description
3949 The function @func{gridDefNumber} defines the reference number for an unstructured grid.
3950
3951 @EndFunction
3952 */
gridDefNumber(int gridID,int number)3953 void gridDefNumber(int gridID, int number)
3954 {
3955 cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, number);
3956 }
3957
3958 /*
3959 @Function gridInqNumber
3960 @Title Get the reference number to an unstructured grid
3961
3962 @Prototype int gridInqNumber(int gridID)
3963 @Parameter
3964 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
3965
3966 @Description
3967 The function @func{gridInqNumber} returns the reference number to an unstructured grid.
3968
3969 @Result
3970 @func{gridInqNumber} returns the reference number to an unstructured grid.
3971 @EndFunction
3972 */
gridInqNumber(int gridID)3973 int gridInqNumber(int gridID)
3974 {
3975 int number = 0;
3976 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
3977 return number;
3978 }
3979
3980 /*
3981 @Function gridDefPosition
3982 @Title Define the position of grid in the reference file
3983
3984 @Prototype void gridDefPosition(int gridID, const int position)
3985 @Parameter
3986 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
3987 @Item position Position of grid in the reference file.
3988
3989 @Description
3990 The function @func{gridDefPosition} defines the position of grid in the reference file.
3991
3992 @EndFunction
3993 */
gridDefPosition(int gridID,int position)3994 void gridDefPosition(int gridID, int position)
3995 {
3996 cdiDefKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, position);
3997 }
3998
3999 /*
4000 @Function gridInqPosition
4001 @Title Get the position of grid in the reference file
4002
4003 @Prototype int gridInqPosition(int gridID)
4004 @Parameter
4005 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
4006
4007 @Description
4008 The function @func{gridInqPosition} returns the position of grid in the reference file.
4009
4010 @Result
4011 @func{gridInqPosition} returns the position of grid in the reference file.
4012 @EndFunction
4013 */
gridInqPosition(int gridID)4014 int gridInqPosition(int gridID)
4015 {
4016 int position = 0;
4017 cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDINREFERENCE, &position);
4018 return position;
4019 }
4020
4021 /*
4022 @Function gridDefReference
4023 @Title Define the reference URI for an unstructured grid
4024
4025 @Prototype void gridDefReference(int gridID, const char *reference)
4026 @Parameter
4027 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
4028 @Item reference Reference URI for an unstructured grid.
4029
4030 @Description
4031 The function @func{gridDefReference} defines the reference URI for an unstructured grid.
4032
4033 @EndFunction
4034 */
gridDefReference(int gridID,const char * reference)4035 void gridDefReference(int gridID, const char *reference)
4036 {
4037 if (reference)
4038 {
4039 cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_REFERENCEURI, reference);
4040 gridMark4Update(gridID);
4041 }
4042 }
4043
4044 /*
4045 @Function gridInqReference
4046 @Title Get the reference URI to an unstructured grid
4047
4048 @Prototype char *gridInqReference(int gridID, char *reference)
4049 @Parameter
4050 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
4051
4052 @Description
4053 The function @func{gridInqReference} returns the reference URI to an unstructured grid.
4054
4055 @Result
4056 @func{gridInqReference} returns the reference URI to an unstructured grid.
4057 @EndFunction
4058 */
gridInqReference(int gridID,char * reference)4059 int gridInqReference(int gridID, char *reference)
4060 {
4061 int length = 0;
4062 if (CDI_NOERR == cdiInqKeyLen(gridID, CDI_GLOBAL, CDI_KEY_REFERENCEURI, &length))
4063 {
4064 if (reference)
4065 cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_REFERENCEURI, reference, &length);
4066 }
4067
4068 return length;
4069 }
4070
4071 /*
4072 @Function gridDefUUID
4073 @Title Define the UUID for an unstructured grid
4074
4075 @Prototype void gridDefUUID(int gridID, const char *uuid)
4076 @Parameter
4077 @Item gridID Grid ID, from a previous call to @fref{gridCreate}.
4078 @Item uuid UUID for an unstructured grid.
4079
4080 @Description
4081 The function @func{gridDefUUID} defines the UUID for an unstructured grid.
4082
4083 @EndFunction
4084 */
gridDefUUID(int gridID,const unsigned char uuid[CDI_UUID_SIZE])4085 void gridDefUUID(int gridID, const unsigned char uuid[CDI_UUID_SIZE])
4086 {
4087 cdiDefKeyBytes(gridID, CDI_GLOBAL, CDI_KEY_UUID, uuid, CDI_UUID_SIZE);
4088
4089 gridMark4Update(gridID);
4090 }
4091
4092 /*
4093 @Function gridInqUUID
4094 @Title Get the UUID to an unstructured grid
4095
4096 @Prototype void gridInqUUID(int gridID, char *uuid)
4097 @Parameter
4098 @Item gridID Grid ID, from a previous call to @fref{gridCreate} or @fref{vlistInqVarGrid}.
4099
4100 @Description
4101 The function @func{gridInqUUID} returns the UUID to an unstructured grid.
4102
4103 @Result
4104 @func{gridInqUUID} returns the UUID to an unstructured grid to the parameter uuid.
4105 @EndFunction
4106 */
gridInqUUID(int gridID,unsigned char uuid[CDI_UUID_SIZE])4107 void gridInqUUID(int gridID, unsigned char uuid[CDI_UUID_SIZE])
4108 {
4109 memset(uuid, 0, CDI_UUID_SIZE);
4110 int length = CDI_UUID_SIZE;
4111 cdiInqKeyBytes(gridID, CDI_GLOBAL, CDI_KEY_UUID, uuid, &length);
4112 }
4113
4114
cdiGridGetIndexList(unsigned ngrids,int * gridIndexList)4115 void cdiGridGetIndexList(unsigned ngrids, int * gridIndexList)
4116 {
4117 reshGetResHListOfType(ngrids, gridIndexList, &gridOps);
4118 }
4119
4120 static int
gridTxCode()4121 gridTxCode ()
4122 {
4123 return GRID;
4124 }
4125
4126 enum {
4127 GRID_PACK_INT_IDX_SELF,
4128 GRID_PACK_INT_IDX_TYPE,
4129 GRID_PACK_INT_IDX_IS_CYCLIC,
4130 GRID_PACK_INT_IDX_X_FLAG,
4131 GRID_PACK_INT_IDX_Y_FLAG,
4132 GRID_PACK_INT_IDX_GME_ND,
4133 GRID_PACK_INT_IDX_GME_NI,
4134 GRID_PACK_INT_IDX_GME_NI2,
4135 GRID_PACK_INT_IDX_GME_NI3,
4136 GRID_PACK_INT_IDX_TRUNC,
4137 GRID_PACK_INT_IDX_NVERTEX,
4138 REDUCEDPOINTSSIZE,
4139 GRID_PACK_INT_IDX_SIZE,
4140 GRID_PACK_INT_IDX_X_SIZE,
4141 GRID_PACK_INT_IDX_Y_SIZE,
4142 GRID_PACK_INT_IDX_LCOMPLEX,
4143 GRID_PACK_INT_IDX_MEMBERMASK,
4144 GRID_PACK_INT_IDX_XTSTDNNAME,
4145 GRID_PACK_INT_IDX_YTSTDNNAME,
4146 /*
4147 GRID_PACK_INT_IDX_ISCANSNEGATIVELY,
4148 GRID_PACK_INT_IDX_JSCANSPOSITIVELY,
4149 GRID_PACK_INT_IDX_JPOINTSARECONSECUTIVE,
4150 */
4151 gridNint
4152 };
4153
4154 enum {
4155 GRID_PACK_DBL_IDX_X_FIRST,
4156 GRID_PACK_DBL_IDX_Y_FIRST,
4157 GRID_PACK_DBL_IDX_X_LAST,
4158 GRID_PACK_DBL_IDX_Y_LAST,
4159 GRID_PACK_DBL_IDX_X_INC,
4160 GRID_PACK_DBL_IDX_Y_INC,
4161 gridNdouble
4162 };
4163
4164 enum {
4165 gridHasMaskFlag = 1 << 0,
4166 gridHasGMEMaskFlag = 1 << 1,
4167 gridHasXValsFlag = 1 << 2,
4168 gridHasYValsFlag = 1 << 3,
4169 gridHasAreaFlag = 1 << 4,
4170 gridHasXBoundsFlag = 1 << 5,
4171 gridHasYBoundsFlag = 1 << 6,
4172 gridHasReducedPointsFlag = 1 << 7,
4173 };
4174
4175
gridGetComponentFlags(const grid_t * gridP)4176 static int gridGetComponentFlags(const grid_t * gridP)
4177 {
4178 int flags = (gridHasMaskFlag & (int)((unsigned)(gridP->mask == NULL) - 1U))
4179 | (gridHasGMEMaskFlag & (int)((unsigned)(gridP->mask_gme == NULL) - 1U))
4180 | (gridHasXValsFlag & (int)((unsigned)(gridP->vtable->inqXValsPtr((grid_t *)gridP) == NULL) - 1U))
4181 | (gridHasYValsFlag & (int)((unsigned)(gridP->vtable->inqYValsPtr((grid_t *)gridP) == NULL) - 1U))
4182 | (gridHasAreaFlag & (int)((unsigned)(gridP->vtable->inqAreaPtr((grid_t *)gridP) == NULL) - 1U))
4183 | (gridHasXBoundsFlag & (int)((unsigned)(gridP->x.bounds == NULL) - 1U))
4184 | (gridHasYBoundsFlag & (int)((unsigned)(gridP->y.bounds == NULL) - 1U))
4185 | (gridHasReducedPointsFlag & (int)((unsigned)(gridP->reducedPoints == NULL) - 1U));
4186 return flags;
4187 }
4188
4189 static int
gridGetPackSize(void * voidP,void * context)4190 gridGetPackSize(void * voidP, void *context)
4191 {
4192 grid_t * gridP = ( grid_t * ) voidP;
4193 int packBuffSize = 0, count;
4194
4195 packBuffSize += serializeGetSize(gridNint, CDI_DATATYPE_INT, context)
4196 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4197
4198 if (gridP->reducedPoints)
4199 {
4200 xassert(gridP->reducedPointsSize);
4201 packBuffSize += serializeGetSize(gridP->reducedPointsSize, CDI_DATATYPE_INT, context)
4202 + serializeGetSize( 1, CDI_DATATYPE_UINT32, context);
4203 }
4204
4205 packBuffSize += serializeGetSize(gridNdouble, CDI_DATATYPE_FLT64, context);
4206
4207 if (gridP->vtable->inqXValsPtr(gridP))
4208 {
4209 if (gridP->type == GRID_UNSTRUCTURED || gridP->type == GRID_CURVILINEAR)
4210 count = gridP->size;
4211 else
4212 count = gridP->x.size;
4213 xassert(count);
4214 packBuffSize += serializeGetSize(count, CDI_DATATYPE_FLT64, context)
4215 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4216 }
4217
4218 if (gridP->vtable->inqYValsPtr(gridP))
4219 {
4220 if (gridP->type == GRID_UNSTRUCTURED || gridP->type == GRID_CURVILINEAR)
4221 count = gridP->size;
4222 else
4223 count = gridP->y.size;
4224 xassert(count);
4225 packBuffSize += serializeGetSize(count, CDI_DATATYPE_FLT64, context)
4226 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4227 }
4228
4229 if (gridP->vtable->inqAreaPtr(gridP))
4230 {
4231 xassert(gridP->size);
4232 packBuffSize +=
4233 serializeGetSize(gridP->size, CDI_DATATYPE_FLT64, context)
4234 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4235 }
4236
4237 if (gridP->x.bounds)
4238 {
4239 xassert(gridP->nvertex);
4240 if (gridP->type == GRID_CURVILINEAR || gridP->type == GRID_UNSTRUCTURED)
4241 count = gridP->size;
4242 else
4243 count = gridP->x.size;
4244 xassert(count);
4245 packBuffSize
4246 += (serializeGetSize(gridP->nvertex * count, CDI_DATATYPE_FLT64, context)
4247 + serializeGetSize(1, CDI_DATATYPE_UINT32, context));
4248 }
4249
4250 if (gridP->y.bounds)
4251 {
4252 xassert(gridP->nvertex);
4253 if (gridP->type == GRID_CURVILINEAR || gridP->type == GRID_UNSTRUCTURED)
4254 count = gridP->size;
4255 else
4256 count = gridP->y.size;
4257 xassert(count);
4258 packBuffSize
4259 += (serializeGetSize(gridP->nvertex * count, CDI_DATATYPE_FLT64, context)
4260 + serializeGetSize(1, CDI_DATATYPE_UINT32, context));
4261 }
4262
4263 packBuffSize += serializeKeysGetPackSize(&gridP->keys, context);
4264 packBuffSize += serializeKeysGetPackSize(&gridP->x.keys, context);
4265 packBuffSize += serializeKeysGetPackSize(&gridP->y.keys, context);
4266
4267 if (gridP->mask)
4268 {
4269 xassert(gridP->size);
4270 packBuffSize
4271 += serializeGetSize(gridP->size, CDI_DATATYPE_UCHAR, context)
4272 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4273 }
4274
4275 if (gridP->mask_gme)
4276 {
4277 xassert(gridP->size);
4278 packBuffSize += serializeGetSize(gridP->size, CDI_DATATYPE_UCHAR, context)
4279 + serializeGetSize(1, CDI_DATATYPE_UINT32, context);
4280 }
4281
4282 return packBuffSize;
4283 }
4284
4285 void
gridUnpack(char * unpackBuffer,int unpackBufferSize,int * unpackBufferPos,int originNamespace,void * context,int force_id)4286 gridUnpack(char * unpackBuffer, int unpackBufferSize,
4287 int * unpackBufferPos, int originNamespace, void *context,
4288 int force_id)
4289 {
4290 grid_t * gridP;
4291 uint32_t d;
4292 int memberMask, size;
4293
4294 gridInit();
4295
4296 {
4297 int intBuffer[gridNint];
4298 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4299 intBuffer, gridNint, CDI_DATATYPE_INT, context);
4300 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4301 &d, 1, CDI_DATATYPE_UINT32, context);
4302
4303 xassert(cdiCheckSum(CDI_DATATYPE_INT, gridNint, intBuffer) == d);
4304 int targetID = namespaceAdaptKey(intBuffer[0], originNamespace);
4305 gridP = gridNewEntry(force_id?targetID:CDI_UNDEFID);
4306
4307 xassert(!force_id || targetID == gridP->self);
4308
4309 gridP->type = intBuffer[GRID_PACK_INT_IDX_TYPE];
4310 gridP->isCyclic = (signed char)intBuffer[GRID_PACK_INT_IDX_IS_CYCLIC];
4311 gridP->x.flag = (short)intBuffer[GRID_PACK_INT_IDX_X_FLAG];
4312 gridP->y.flag = (short)intBuffer[GRID_PACK_INT_IDX_Y_FLAG];
4313 gridP->gme.nd = intBuffer[GRID_PACK_INT_IDX_GME_ND];
4314 gridP->gme.ni = intBuffer[GRID_PACK_INT_IDX_GME_NI];
4315 gridP->gme.ni2 = intBuffer[GRID_PACK_INT_IDX_GME_NI2];
4316 gridP->gme.ni3 = intBuffer[GRID_PACK_INT_IDX_GME_NI3];
4317 gridP->trunc = intBuffer[GRID_PACK_INT_IDX_TRUNC];
4318 gridP->nvertex = intBuffer[GRID_PACK_INT_IDX_NVERTEX];
4319 gridP->reducedPointsSize = intBuffer[REDUCEDPOINTSSIZE];
4320 gridP->size = intBuffer[GRID_PACK_INT_IDX_SIZE];
4321 gridP->x.size = intBuffer[GRID_PACK_INT_IDX_X_SIZE];
4322 gridP->y.size = intBuffer[GRID_PACK_INT_IDX_Y_SIZE];
4323 gridP->lcomplex = (bool)intBuffer[GRID_PACK_INT_IDX_LCOMPLEX];
4324 memberMask = intBuffer[GRID_PACK_INT_IDX_MEMBERMASK];
4325 }
4326
4327 if (memberMask & gridHasReducedPointsFlag)
4328 {
4329 xassert(gridP->reducedPointsSize);
4330 gridP->reducedPoints = (int *) Malloc((size_t)gridP->reducedPointsSize * sizeof (int));
4331 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4332 gridP->reducedPoints, gridP->reducedPointsSize , CDI_DATATYPE_INT, context);
4333 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4334 &d, 1, CDI_DATATYPE_UINT32, context);
4335 xassert(cdiCheckSum(CDI_DATATYPE_INT, gridP->reducedPointsSize, gridP->reducedPoints) == d);
4336 }
4337
4338 {
4339 double doubleBuffer[gridNdouble];
4340 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4341 doubleBuffer, gridNdouble, CDI_DATATYPE_FLT64, context);
4342 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4343 &d, 1, CDI_DATATYPE_UINT32, context);
4344 xassert(d == cdiCheckSum(CDI_DATATYPE_FLT, gridNdouble, doubleBuffer));
4345
4346 gridP->x.first = doubleBuffer[GRID_PACK_DBL_IDX_X_FIRST];
4347 gridP->y.first = doubleBuffer[GRID_PACK_DBL_IDX_Y_FIRST];
4348 gridP->x.last = doubleBuffer[GRID_PACK_DBL_IDX_X_LAST];
4349 gridP->y.last = doubleBuffer[GRID_PACK_DBL_IDX_Y_LAST];
4350 gridP->x.inc = doubleBuffer[GRID_PACK_DBL_IDX_X_INC];
4351 gridP->y.inc = doubleBuffer[GRID_PACK_DBL_IDX_Y_INC];
4352 }
4353
4354 bool irregular = gridP->type == GRID_UNSTRUCTURED
4355 || gridP->type == GRID_CURVILINEAR;
4356 if (memberMask & gridHasXValsFlag)
4357 {
4358 size = irregular ? gridP->size : gridP->x.size;
4359
4360 gridP->x.vals = (double *) Malloc(size * sizeof (double));
4361 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4362 gridP->x.vals, size, CDI_DATATYPE_FLT64, context);
4363 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4364 &d, 1, CDI_DATATYPE_UINT32, context);
4365 xassert(cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->x.vals) == d );
4366 }
4367
4368 if (memberMask & gridHasYValsFlag)
4369 {
4370 size = irregular ? gridP->size : gridP->y.size;
4371
4372 gridP->y.vals = (double *) Malloc(size * sizeof (double));
4373 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4374 gridP->y.vals, size, CDI_DATATYPE_FLT64, context);
4375 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4376 &d, 1, CDI_DATATYPE_UINT32, context);
4377 xassert(cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->y.vals) == d);
4378 }
4379
4380 if (memberMask & gridHasAreaFlag)
4381 {
4382 size = gridP->size;
4383 xassert(size);
4384 gridP->area = (double *) Malloc(size * sizeof (double));
4385 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4386 gridP->area, size, CDI_DATATYPE_FLT64, context);
4387 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4388 &d, 1, CDI_DATATYPE_UINT32, context);
4389 xassert(cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->area) == d);
4390 }
4391
4392 if (memberMask & gridHasXBoundsFlag)
4393 {
4394 size = gridP->nvertex * (irregular ? gridP->size : gridP->x.size);
4395 xassert(size);
4396
4397 gridP->x.bounds = (double *) Malloc(size * sizeof (double));
4398 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4399 gridP->x.bounds, size, CDI_DATATYPE_FLT64, context);
4400 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4401 &d, 1, CDI_DATATYPE_UINT32, context);
4402 xassert(cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->x.bounds) == d);
4403 }
4404
4405 if (memberMask & gridHasYBoundsFlag)
4406 {
4407 size = gridP->nvertex * (irregular ? gridP->size : gridP->y.size);
4408 xassert(size);
4409
4410 gridP->y.bounds = (double *) Malloc(size * sizeof (double));
4411 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4412 gridP->y.bounds, size, CDI_DATATYPE_FLT64, context);
4413 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4414 &d, 1, CDI_DATATYPE_UINT32, context);
4415 xassert(cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->y.bounds) == d);
4416 }
4417
4418 serializeKeysUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos, &gridP->keys, context);
4419 serializeKeysUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos, &gridP->x.keys, context);
4420 serializeKeysUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos, &gridP->y.keys, context);
4421
4422 if (memberMask & gridHasMaskFlag)
4423 {
4424 xassert((size = gridP->size));
4425 gridP->mask = (mask_t *) Malloc(size * sizeof (mask_t));
4426 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4427 gridP->mask, gridP->size, CDI_DATATYPE_UCHAR, context);
4428 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4429 &d, 1, CDI_DATATYPE_UINT32, context);
4430 xassert(cdiCheckSum(CDI_DATATYPE_UCHAR, gridP->size, gridP->mask) == d);
4431 }
4432
4433 if (memberMask & gridHasGMEMaskFlag)
4434 {
4435 xassert((size = gridP->size));
4436 gridP->mask_gme = (mask_t *) Malloc(size * sizeof (mask_t));
4437 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4438 gridP->mask_gme, gridP->size, CDI_DATATYPE_UCHAR, context);
4439 serializeUnpack(unpackBuffer, unpackBufferSize, unpackBufferPos,
4440 &d, 1, CDI_DATATYPE_UINT32, context);
4441 xassert(cdiCheckSum(CDI_DATATYPE_UCHAR, gridP->size, gridP->mask_gme) == d);
4442 }
4443
4444 reshSetStatus(gridP->self, &gridOps,
4445 reshGetStatus(gridP->self, &gridOps) & ~RESH_SYNC_BIT);
4446 }
4447
4448
4449 static void
gridPack(void * voidP,void * packBuffer,int packBufferSize,int * packBufferPos,void * context)4450 gridPack(void * voidP, void * packBuffer, int packBufferSize,
4451 int * packBufferPos, void *context)
4452 {
4453 grid_t * gridP = ( grid_t * ) voidP;
4454 int size;
4455 uint32_t d;
4456 int memberMask;
4457
4458 {
4459 int intBuffer[gridNint];
4460
4461 intBuffer[GRID_PACK_INT_IDX_SELF] = gridP->self;
4462 intBuffer[GRID_PACK_INT_IDX_TYPE] = gridP->type;
4463 intBuffer[GRID_PACK_INT_IDX_IS_CYCLIC] = gridP->isCyclic;
4464 intBuffer[GRID_PACK_INT_IDX_X_FLAG] = gridP->x.flag;
4465 intBuffer[GRID_PACK_INT_IDX_Y_FLAG] = gridP->y.flag;
4466 intBuffer[GRID_PACK_INT_IDX_GME_ND] = gridP->gme.nd;
4467 intBuffer[GRID_PACK_INT_IDX_GME_NI] = gridP->gme.ni;
4468 intBuffer[GRID_PACK_INT_IDX_GME_NI2] = gridP->gme.ni2;
4469 intBuffer[GRID_PACK_INT_IDX_GME_NI3] = gridP->gme.ni3;
4470 intBuffer[GRID_PACK_INT_IDX_TRUNC] = gridP->trunc;
4471 intBuffer[GRID_PACK_INT_IDX_NVERTEX] = gridP->nvertex;
4472 intBuffer[REDUCEDPOINTSSIZE] = gridP->reducedPointsSize;
4473 intBuffer[GRID_PACK_INT_IDX_SIZE] = gridP->size;
4474 intBuffer[GRID_PACK_INT_IDX_X_SIZE] = gridP->x.size;
4475 intBuffer[GRID_PACK_INT_IDX_Y_SIZE] = gridP->y.size;
4476 intBuffer[GRID_PACK_INT_IDX_LCOMPLEX] = gridP->lcomplex;
4477 intBuffer[GRID_PACK_INT_IDX_MEMBERMASK] = memberMask
4478 = gridGetComponentFlags(gridP);
4479
4480 serializePack(intBuffer, gridNint, CDI_DATATYPE_INT,
4481 packBuffer, packBufferSize, packBufferPos, context);
4482 d = cdiCheckSum(CDI_DATATYPE_INT, gridNint, intBuffer);
4483 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4484 packBuffer, packBufferSize, packBufferPos, context);
4485 }
4486
4487 if (memberMask & gridHasReducedPointsFlag)
4488 {
4489 size = gridP->reducedPointsSize;
4490 xassert(size > 0);
4491 serializePack(gridP->reducedPoints, size, CDI_DATATYPE_INT,
4492 packBuffer, packBufferSize, packBufferPos, context);
4493 d = cdiCheckSum(CDI_DATATYPE_INT , size, gridP->reducedPoints);
4494 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4495 packBuffer, packBufferSize, packBufferPos, context);
4496 }
4497
4498 {
4499 double doubleBuffer[gridNdouble];
4500
4501 doubleBuffer[GRID_PACK_DBL_IDX_X_FIRST] = gridP->x.first;
4502 doubleBuffer[GRID_PACK_DBL_IDX_Y_FIRST] = gridP->y.first;
4503 doubleBuffer[GRID_PACK_DBL_IDX_X_LAST] = gridP->x.last;
4504 doubleBuffer[GRID_PACK_DBL_IDX_Y_LAST] = gridP->y.last;
4505 doubleBuffer[GRID_PACK_DBL_IDX_X_INC] = gridP->x.inc;
4506 doubleBuffer[GRID_PACK_DBL_IDX_Y_INC] = gridP->y.inc;
4507
4508 serializePack(doubleBuffer, gridNdouble, CDI_DATATYPE_FLT64,
4509 packBuffer, packBufferSize, packBufferPos, context);
4510 d = cdiCheckSum(CDI_DATATYPE_FLT, gridNdouble, doubleBuffer);
4511 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4512 packBuffer, packBufferSize, packBufferPos, context);
4513 }
4514
4515 if (memberMask & gridHasXValsFlag)
4516 {
4517 if (gridP->type == GRID_UNSTRUCTURED || gridP->type == GRID_CURVILINEAR)
4518 size = gridP->size;
4519 else
4520 size = gridP->x.size;
4521 xassert(size);
4522
4523 const double *gridP_xvals = gridP->vtable->inqXValsPtr(gridP);
4524 serializePack(gridP_xvals, size, CDI_DATATYPE_FLT64,
4525 packBuffer, packBufferSize, packBufferPos, context);
4526 d = cdiCheckSum(CDI_DATATYPE_FLT, size, gridP_xvals);
4527 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4528 packBuffer, packBufferSize, packBufferPos, context);
4529 }
4530
4531 if (memberMask & gridHasYValsFlag)
4532 {
4533 if (gridP->type == GRID_UNSTRUCTURED || gridP->type == GRID_CURVILINEAR )
4534 size = gridP->size;
4535 else
4536 size = gridP->y.size;
4537 xassert(size);
4538 const double *gridP_yvals = gridP->vtable->inqYValsPtr(gridP);
4539 serializePack(gridP_yvals, size, CDI_DATATYPE_FLT64,
4540 packBuffer, packBufferSize, packBufferPos, context);
4541 d = cdiCheckSum(CDI_DATATYPE_FLT, size, gridP_yvals);
4542 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4543 packBuffer, packBufferSize, packBufferPos, context);
4544 }
4545
4546 if (memberMask & gridHasAreaFlag)
4547 {
4548 xassert(gridP->size);
4549
4550 serializePack(gridP->area, gridP->size, CDI_DATATYPE_FLT64,
4551 packBuffer, packBufferSize, packBufferPos, context);
4552 d = cdiCheckSum(CDI_DATATYPE_FLT, gridP->size, gridP->area);
4553 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4554 packBuffer, packBufferSize, packBufferPos, context);
4555 }
4556
4557 if (memberMask & gridHasXBoundsFlag)
4558 {
4559 xassert ( gridP->nvertex );
4560 if (gridP->type == GRID_CURVILINEAR || gridP->type == GRID_UNSTRUCTURED)
4561 size = gridP->nvertex * gridP->size;
4562 else
4563 size = gridP->nvertex * gridP->x.size;
4564 xassert ( size );
4565
4566 serializePack(gridP->x.bounds, size, CDI_DATATYPE_FLT64,
4567 packBuffer, packBufferSize, packBufferPos, context);
4568 d = cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->x.bounds);
4569 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4570 packBuffer, packBufferSize, packBufferPos, context);
4571 }
4572
4573 if (memberMask & gridHasYBoundsFlag)
4574 {
4575 xassert(gridP->nvertex);
4576 if (gridP->type == GRID_CURVILINEAR || gridP->type == GRID_UNSTRUCTURED)
4577 size = gridP->nvertex * gridP->size;
4578 else
4579 size = gridP->nvertex * gridP->y.size;
4580 xassert ( size );
4581
4582 serializePack(gridP->y.bounds, size, CDI_DATATYPE_FLT64,
4583 packBuffer, packBufferSize, packBufferPos, context);
4584 d = cdiCheckSum(CDI_DATATYPE_FLT, size, gridP->y.bounds);
4585 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4586 packBuffer, packBufferSize, packBufferPos, context);
4587 }
4588
4589 serializeKeysPack(&gridP->keys, packBuffer, packBufferSize, packBufferPos, context);
4590 serializeKeysPack(&gridP->x.keys, packBuffer, packBufferSize, packBufferPos, context);
4591 serializeKeysPack(&gridP->y.keys, packBuffer, packBufferSize, packBufferPos, context);
4592
4593 if (memberMask & gridHasMaskFlag)
4594 {
4595 xassert((size = gridP->size));
4596 serializePack(gridP->mask, size, CDI_DATATYPE_UCHAR,
4597 packBuffer, packBufferSize, packBufferPos, context);
4598 d = cdiCheckSum(CDI_DATATYPE_UCHAR, size, gridP->mask);
4599 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4600 packBuffer, packBufferSize, packBufferPos, context);
4601 }
4602
4603 if (memberMask & gridHasGMEMaskFlag)
4604 {
4605 xassert((size = gridP->size));
4606
4607 serializePack(gridP->mask_gme, size, CDI_DATATYPE_UCHAR,
4608 packBuffer, packBufferSize, packBufferPos, context);
4609 d = cdiCheckSum(CDI_DATATYPE_UCHAR, size, gridP->mask_gme);
4610 serializePack(&d, 1, CDI_DATATYPE_UINT32,
4611 packBuffer, packBufferSize, packBufferPos, context);
4612 }
4613 }
4614
4615
4616 struct gridCompareSearchState
4617 {
4618 int resIDValue;
4619 const grid_t *queryKey;
4620 };
4621
4622 static enum cdiApplyRet
gridCompareSearch(int id,void * res,void * data)4623 gridCompareSearch(int id, void *res, void *data)
4624 {
4625 struct gridCompareSearchState *state = (struct gridCompareSearchState*)data;
4626 (void)res;
4627 if ( gridCompare(id, state->queryKey, true) == false )
4628 {
4629 state->resIDValue = id;
4630 return CDI_APPLY_STOP;
4631 }
4632 else
4633 return CDI_APPLY_GO_ON;
4634 }
4635
4636 // Add grid (which must be Malloc'ed to vlist if not already found)
cdiVlistAddGridIfNew(int vlistID,grid_t * grid,int mode)4637 struct addIfNewRes cdiVlistAddGridIfNew(int vlistID, grid_t *grid, int mode)
4638 {
4639 /*
4640 mode: 0 search in vlist and grid table
4641 1 search in grid table only
4642 2 search in grid table only and don't store the grid in vlist
4643 */
4644 bool gridglobdefined = false;
4645 bool griddefined = false;
4646 int gridID = CDI_UNDEFID;
4647 vlist_t *vlistptr = vlist_to_pointer(vlistID);
4648
4649 unsigned ngrids = (unsigned)vlistptr->ngrids;
4650
4651 if ( mode == 0 )
4652 for ( unsigned index = 0; index < ngrids; index++ )
4653 {
4654 if ( (gridID = vlistptr->gridIDs[index]) != CDI_UNDEFID )
4655 {
4656 if ( gridCompare(gridID, grid, false) == false )
4657 {
4658 griddefined = true;
4659 break;
4660 }
4661 }
4662 else
4663 Error("Internal problem: undefined gridID in vlist %d, position %u!", vlistID, index);
4664 }
4665
4666 if ( ! griddefined )
4667 {
4668 struct gridCompareSearchState query;
4669 query.queryKey = grid;// = { .queryKey = grid };
4670 if ( (gridglobdefined = (cdiResHFilterApply(&gridOps, gridCompareSearch, &query) == CDI_APPLY_STOP)) )
4671 gridID = query.resIDValue;
4672
4673 if ( mode == 1 && gridglobdefined )
4674 for ( unsigned index = 0; index < ngrids; index++ )
4675 if ( vlistptr->gridIDs[index] == gridID )
4676 {
4677 gridglobdefined = false;
4678 break;
4679 }
4680 }
4681
4682 if ( ! griddefined )
4683 {
4684 if ( ! gridglobdefined )
4685 {
4686 grid->self = gridID = reshPut(grid, &gridOps);
4687 gridComplete(grid);
4688 }
4689 if ( mode < 2 )
4690 {
4691 if (ngrids >= MAX_GRIDS_PS)
4692 Error("Internal limit exceeded, MAX_GRIDS_PS=%d needs to be increased!", MAX_GRIDS_PS);
4693 vlistptr->gridIDs[ngrids] = gridID;
4694 vlistptr->ngrids++;
4695 }
4696 }
4697
4698 return (struct addIfNewRes){ .Id = gridID, .isNew = !griddefined && !gridglobdefined };
4699 }
4700
4701
4702 const struct gridVirtTable cdiGridVtable
4703 = {
4704 .destroy = gridDestroyKernel,
4705 .copy = grid_copy_base,
4706 .copyScalarFields = grid_copy_base_scalar_fields,
4707 .copyArrayFields = grid_copy_base_array_fields,
4708 .defXVals = gridDefXValsSerial,
4709 .defYVals = gridDefYValsSerial,
4710 .defMask = gridDefMaskSerial,
4711 .defMaskGME = gridDefMaskGMESerial,
4712 .defXBounds = gridDefXBoundsSerial,
4713 .defYBounds = gridDefYBoundsSerial,
4714 .defArea = gridDefAreaSerial,
4715 .inqXVal = gridInqXValSerial,
4716 .inqYVal = gridInqYValSerial,
4717 .inqXVals = gridInqXValsSerial,
4718 .inqXValsPart = gridInqXValsPartSerial,
4719 .inqYVals = gridInqYValsSerial,
4720 .inqYValsPart = gridInqYValsPartSerial,
4721 .inqXValsPtr = gridInqXValsPtrSerial,
4722 .inqYValsPtr = gridInqYValsPtrSerial,
4723 #ifndef USE_MPI
4724 .inqXIsc = gridInqXIscSerial,
4725 .inqYIsc = gridInqYIscSerial,
4726 .inqXCvals = gridInqXCvalsSerial,
4727 .inqYCvals = gridInqYCvalsSerial,
4728 .inqXCvalsPtr = gridInqXCvalsPtrSerial,
4729 .inqYCvalsPtr = gridInqYCvalsPtrSerial,
4730 #endif
4731 .compareXYFull = compareXYvals,
4732 .compareXYAO = compareXYvals2,
4733 .inqArea = gridInqAreaSerial,
4734 .inqAreaPtr = gridInqAreaPtrBase,
4735 .hasArea = gridHasAreaBase,
4736 .inqMask = gridInqMaskSerial,
4737 .inqMaskGME = gridInqMaskGMESerial,
4738 .inqXBounds = gridInqXBoundsSerial,
4739 .inqYBounds = gridInqYBoundsSerial,
4740 .inqXBoundsPtr = gridInqXBoundsPtrSerial,
4741 .inqYBoundsPtr = gridInqYBoundsPtrSerial,
4742 };
4743
4744 /*
4745 * Local Variables:
4746 * c-file-style: "Java"
4747 * c-basic-offset: 2
4748 * indent-tabs-mode: nil
4749 * show-trailing-whitespace: t
4750 * require-trailing-newline: t
4751 * End:
4752 */
4753