1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #ifdef HAVE_LIBNETCDF
6
7 #include <ctype.h>
8 #include <limits.h>
9
10 #include "dmemory.h"
11 #include "cdi_int.h"
12 #include "cdi_uuid.h"
13 #include "stream_cdf.h"
14 #include "cdf_int.h"
15 #include "varscan.h"
16 #include "vlist.h"
17 #include "cdf_util.h"
18 #include "cdf_lazy_grid.h"
19
20
21 enum VarStatus {
22 UndefVar = -1,
23 CoordVar = 0,
24 DataVar = 1,
25 };
26
27 enum AxisType {
28 X_AXIS = 1,
29 Y_AXIS = 2,
30 Z_AXIS = 3,
31 E_AXIS = 4,
32 T_AXIS = 5,
33 };
34
35 typedef struct {
36 int dimid;
37 int ncvarid;
38 int dimtype;
39 size_t len;
40 char name[CDI_MAX_NAME];
41 }
42 ncdim_t;
43
44 #define MAX_COORDVARS 5
45 #define MAX_AUXVARS 4
46 #define MAX_DIMS_CDF 8
47
48 typedef struct {
49 int ncid;
50 int varStatus;
51 bool ignoreVar;
52 bool isCubeSphere;
53 bool isCharAxis;
54 bool isXaxis;
55 bool isYaxis;
56 bool isZaxis;
57 bool isTaxis;
58 bool isLon;
59 bool isLat;
60 bool isClimatology;
61 bool hasCalendar;
62 bool hasFormulaterms;
63 bool printWarning;
64 int timetype;
65 int param;
66 int code;
67 int tabnum;
68 int bounds;
69 int gridID;
70 int zaxisID;
71 int gridtype;
72 int zaxistype;
73 int xdim;
74 int ydim;
75 int zdim;
76 int xvarid;
77 int yvarid;
78 int rpvarid;
79 int zvarid;
80 int tvarid;
81 int psvarid;
82 int p0varid;
83 int ncoordvars;
84 int cvarids[MAX_COORDVARS];
85 int coordvarids[MAX_COORDVARS];
86 int auxvarids[MAX_AUXVARS];
87 int nauxvars;
88 int cellarea;
89 int tableID;
90 int truncation;
91 int position;
92 int numLPE;
93 bool defmissval;
94 bool deffillval;
95 int xtype;
96 int gmapid;
97 int positive;
98 int ndims;
99 int dimids[MAX_DIMS_CDF];
100 int dimtype[MAX_DIMS_CDF];
101 size_t chunks[MAX_DIMS_CDF];
102 bool chunked;
103 int chunktype;
104 int natts;
105 int *atts;
106 size_t vctsize;
107 double *vct;
108 double missval;
109 double fillval;
110 double addoffset;
111 double scalefactor;
112 bool ldeflate;
113 bool lszip;
114 bool lunsigned;
115 bool lvalidrange;
116 double validrange[2];
117 int typeOfEnsembleForecast;
118 int numberOfForecastsInEnsemble;
119 int perturbationNumber;
120 char name[CDI_MAX_NAME];
121 char longname[CDI_MAX_NAME];
122 char stdname[CDI_MAX_NAME];
123 char units[CDI_MAX_NAME];
124 char extra[CDI_MAX_NAME];
125 }
126 ncvar_t;
127
128 typedef struct {
129 char gridfile[8912];
130 unsigned char uuid[CDI_UUID_SIZE];
131 int number_of_grid_used;
132 int timedimid;
133 }
134 GridInfo;
135
136 static
scanTimeString(const char * ptu,int64_t * rdate,int * rtime)137 void scanTimeString(const char *ptu, int64_t *rdate, int *rtime)
138 {
139 int year = 1, month = 1, day = 1;
140 int hour = 0, minute = 0, second = 0;
141 int v1 = 1, v2 = 1, v3 = 1;
142 char ch = ' ';
143
144 if ( *ptu ) sscanf(ptu, "%d-%d-%d%c%d:%d:%d", &v1, &v2, &v3, &ch, &hour, &minute, &second);
145
146 if ( v3 > 999 && v1 < 32 )
147 { year = v3; month = v2; day = v1; }
148 else
149 { year = v1; month = v2; day = v3; }
150
151 *rdate = cdiEncodeDate(year, month, day);
152 *rtime = cdiEncodeTime(hour, minute, second);
153 }
154
155 static
scanTimeUnit(const char * unitstr)156 int scanTimeUnit(const char *unitstr)
157 {
158 const size_t len = strlen(unitstr);
159 const int timeunit = get_timeunit(len, unitstr);
160 if ( timeunit == -1 ) Message("Unsupported TIMEUNIT: %s!", unitstr);
161 return timeunit;
162 }
163
164 static
setForecastTime(const char * timestr,taxis_t * taxis)165 void setForecastTime(const char *timestr, taxis_t *taxis)
166 {
167 const size_t len = strlen(timestr);
168 if ( len != 0 )
169 scanTimeString(timestr, &taxis->fdate, &taxis->ftime);
170 else
171 taxis->fdate = taxis->ftime = 0;
172 }
173
174 static
setBaseTime(const char * timeunits,taxis_t * taxis)175 int setBaseTime(const char *timeunits, taxis_t *taxis)
176 {
177 int taxistype = TAXIS_ABSOLUTE;
178
179 size_t len = strlen(timeunits);
180 while ( isspace(*timeunits) && len ) { timeunits++; len--; }
181
182 char *tu = (char *)Malloc((len+1) * sizeof(char));
183
184 for ( size_t i = 0; i < len; i++ ) tu[i] = (char)tolower((int)timeunits[i]);
185 tu[len] = 0;
186
187 int timeunit = get_timeunit(len, tu);
188 if ( timeunit == -1 )
189 {
190 Message("Unsupported TIMEUNIT: %s!", timeunits);
191 return 1;
192 }
193
194 size_t pos = 0;
195 while ( pos < len && !isspace(tu[pos]) ) ++pos;
196 if ( tu[pos] )
197 {
198 while ( isspace(tu[pos]) ) ++pos;
199
200 if ( strStartsWith(tu+pos, "since") ) taxistype = TAXIS_RELATIVE;
201
202 while ( pos < len && !isspace(tu[pos]) ) ++pos;
203 if ( tu[pos] )
204 {
205 while ( isspace(tu[pos]) ) ++pos;
206
207 if ( taxistype == TAXIS_ABSOLUTE )
208 {
209 if ( timeunit == TUNIT_DAY )
210 {
211 if ( !strStartsWith(tu+pos, "%y%m%d.%f") )
212 {
213 Warning("Unsupported format %s for TIMEUNIT day!", tu+pos);
214 timeunit = -1;
215 }
216 }
217 else if ( timeunit == TUNIT_MONTH )
218 {
219 if ( !strStartsWith(tu+pos, "%y%m.%f") )
220 {
221 Warning("Unsupported format %s for TIMEUNIT month!", tu+pos);
222 timeunit = -1;
223 }
224 }
225 else if ( timeunit == TUNIT_YEAR )
226 {
227 if ( !strStartsWith(tu+pos, "%y.%f") )
228 {
229 Warning("Unsupported format %s for TIMEUNIT year!", tu+pos);
230 timeunit = -1;
231 }
232 }
233 else
234 {
235 Warning("Unsupported format for time units: %s!", tu);
236 }
237 }
238 else if ( taxistype == TAXIS_RELATIVE )
239 {
240 int64_t rdate = -1;
241 int rtime = -1;
242 scanTimeString(tu+pos, &rdate, &rtime);
243 taxis->rdate = rdate;
244 taxis->rtime = rtime;
245
246 if ( CDI_Debug ) Message("rdate = %lld rtime = %d", rdate, rtime);
247 }
248 }
249 }
250
251 taxis->type = taxistype;
252 taxis->unit = timeunit;
253
254 Free(tu);
255
256 if ( CDI_Debug ) Message("taxistype = %d unit = %d", taxistype, timeunit);
257
258 return 0;
259 }
260
261 static
xtypeIsText(int xtype)262 bool xtypeIsText(int xtype)
263 {
264 const bool isText = (xtype == NC_CHAR)
265 #ifdef HAVE_NETCDF4
266 || (xtype == NC_STRING)
267 #endif
268 ;
269 return isText;
270 }
271
272 static
xtypeIsFloat(nc_type xtype)273 bool xtypeIsFloat(nc_type xtype)
274 {
275 return (xtype == NC_FLOAT || xtype == NC_DOUBLE);
276 }
277
278 static
xtypeIsInt(nc_type xtype)279 bool xtypeIsInt(nc_type xtype)
280 {
281 const bool isInt = xtype == NC_SHORT || xtype == NC_INT || xtype == NC_BYTE
282 #ifdef HAVE_NETCDF4
283 || xtype == NC_USHORT || xtype == NC_UINT || xtype == NC_UBYTE
284 #endif
285 ;
286
287 return isInt;
288 }
289
290 static
cdfInqDatatype(stream_t * streamptr,int xtype,bool lunsigned)291 int cdfInqDatatype(stream_t *streamptr, int xtype, bool lunsigned)
292 {
293 int datatype = -1;
294
295 #ifdef HAVE_NETCDF4
296 if ( xtype == NC_BYTE && lunsigned ) xtype = NC_UBYTE;
297 #endif
298
299 if ( xtype == NC_BYTE ) datatype = CDI_DATATYPE_INT8;
300 else if ( xtype == NC_CHAR ) datatype = CDI_DATATYPE_UINT8;
301 else if ( xtype == NC_SHORT ) datatype = CDI_DATATYPE_INT16;
302 else if ( xtype == NC_INT ) datatype = CDI_DATATYPE_INT32;
303 else if ( xtype == NC_FLOAT ) datatype = CDI_DATATYPE_FLT32;
304 else if ( xtype == NC_DOUBLE ) datatype = CDI_DATATYPE_FLT64;
305 #ifdef HAVE_NETCDF4
306 else if ( xtype == NC_UBYTE ) datatype = CDI_DATATYPE_UINT8;
307 else if ( xtype == NC_LONG ) datatype = CDI_DATATYPE_INT32;
308 else if ( xtype == NC_USHORT ) datatype = CDI_DATATYPE_UINT16;
309 else if ( xtype == NC_UINT ) datatype = CDI_DATATYPE_UINT32;
310 else if ( xtype == NC_INT64 ) datatype = CDI_DATATYPE_FLT64;
311 else if ( xtype == NC_UINT64 ) datatype = CDI_DATATYPE_FLT64;
312 else
313 {
314 if (xtype != streamptr->nc_complex_float_id && xtype != streamptr->nc_complex_double_id)
315 {
316 bool isUserDefinedType = false;
317 #ifdef NC_FIRSTUSERTYPEID
318 isUserDefinedType = (xtype >= NC_FIRSTUSERTYPEID);
319 #endif
320 if (isUserDefinedType)
321 {
322 int fileID = streamptr->fileID;
323 size_t nfields = 0, compoundsize = 0;
324 int status = nc_inq_compound(fileID, xtype, NULL, &compoundsize, &nfields);
325 if (status == NC_NOERR && nfields == 2 && (compoundsize == 8 || compoundsize == 16))
326 {
327 nc_type field_type1 = -1, field_type2 = -1;
328 int field_dims1 = 0, field_dims2 = 0;
329 nc_inq_compound_field(fileID, xtype, 0, NULL, NULL, &field_type1, &field_dims1, NULL);
330 nc_inq_compound_field(fileID, xtype, 1, NULL, NULL, &field_type2, &field_dims2, NULL);
331 if (field_type1 == field_type2 && field_dims1 == 0 && field_dims2 == 0)
332 {
333 if (field_type1 == NC_FLOAT) streamptr->nc_complex_float_id = xtype;
334 else if (field_type1 == NC_DOUBLE) streamptr->nc_complex_double_id = xtype;
335 }
336 }
337 }
338 }
339 if ( xtype == streamptr->nc_complex_float_id ) datatype = CDI_DATATYPE_CPX32;
340 else if ( xtype == streamptr->nc_complex_double_id ) datatype = CDI_DATATYPE_CPX64;
341 }
342 #endif
343
344 return datatype;
345 }
346
347 static
cdfGetAttInt(int fileID,int ncvarid,const char * attname,size_t attlen,int * attint)348 void cdfGetAttInt(int fileID, int ncvarid, const char *attname, size_t attlen, int *attint)
349 {
350 *attint = 0;
351
352 nc_type atttype;
353 size_t nc_attlen;
354 cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
355 cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
356
357 if (xtypeIsFloat(atttype) || xtypeIsInt(atttype))
358 {
359 const bool lalloc = (nc_attlen > attlen);
360 int *pintatt = lalloc ? (int *)(Malloc(nc_attlen*sizeof(int))) : attint;
361 cdf_get_att_int(fileID, ncvarid, attname, pintatt);
362 if (lalloc)
363 {
364 memcpy(attint, pintatt, attlen*sizeof(int));
365 Free(pintatt);
366 }
367 }
368 }
369
370 static
cdfGetAttDouble(int fileID,int ncvarid,char * attname,size_t attlen,double * attdouble)371 void cdfGetAttDouble(int fileID, int ncvarid, char *attname, size_t attlen, double *attdouble)
372 {
373 *attdouble = 0.0;
374
375 nc_type atttype;
376 size_t nc_attlen;
377 cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
378 cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
379
380 if (xtypeIsFloat(atttype) || xtypeIsInt(atttype))
381 {
382 const bool lalloc = (nc_attlen > attlen);
383 double *pdoubleatt = lalloc ? (double*)Malloc(nc_attlen*sizeof(double)) : attdouble;
384 cdf_get_att_double(fileID, ncvarid, attname, pdoubleatt);
385 if (lalloc)
386 {
387 memcpy(attdouble, pdoubleatt, attlen*sizeof(double));
388 Free(pdoubleatt);
389 }
390 }
391 }
392
393 static
cdfCheckAttText(int fileID,int ncvarid,const char * attname)394 bool cdfCheckAttText(int fileID, int ncvarid, const char *attname)
395 {
396 bool status = false;
397 nc_type atttype;
398
399 const int status_nc = nc_inq_atttype(fileID, ncvarid, attname, &atttype);
400 if ( status_nc == NC_NOERR && (atttype == NC_CHAR
401 #ifdef HAVE_NETCDF4
402 || atttype == NC_STRING
403 #endif
404 ) )
405 {
406 status = true;
407 }
408
409 return status;
410 }
411
412 static
cdfGetAttText(int fileID,int ncvarid,const char * attname,size_t attlen,char * atttext)413 void cdfGetAttText(int fileID, int ncvarid, const char *attname, size_t attlen, char *atttext)
414 {
415 atttext[0] = 0;
416
417 nc_type atttype;
418 size_t nc_attlen;
419 cdf_inq_atttype(fileID, ncvarid, attname, &atttype);
420 cdf_inq_attlen(fileID, ncvarid, attname, &nc_attlen);
421
422 if ( atttype == NC_CHAR )
423 {
424 char attbuf[65636];
425 if ( nc_attlen < sizeof(attbuf) )
426 {
427 cdf_get_att_text(fileID, ncvarid, attname, attbuf);
428
429 if ( nc_attlen > (attlen-1) ) nc_attlen = (attlen-1);
430
431 attbuf[nc_attlen++] = 0;
432 memcpy(atttext, attbuf, nc_attlen);
433 }
434 }
435 #ifdef HAVE_NETCDF4
436 else if ( atttype == NC_STRING )
437 {
438 if ( nc_attlen == 1 )
439 {
440 char *attbuf = NULL;
441 cdf_get_att_string(fileID, ncvarid, attname, &attbuf);
442
443 size_t ssize = strlen(attbuf) + 1;
444 if ( ssize > attlen ) ssize = attlen;
445 memcpy(atttext, attbuf, ssize);
446 atttext[ssize - 1] = 0;
447 Free(attbuf);
448 }
449 }
450 #endif
451 }
452
453
cdf_scale_add(size_t size,double * data,double addoffset,double scalefactor)454 void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor)
455 {
456 const bool laddoffset = IS_NOT_EQUAL(addoffset, 0);
457 const bool lscalefactor = IS_NOT_EQUAL(scalefactor, 1);
458
459 if ( laddoffset && lscalefactor )
460 {
461 for (size_t i = 0; i < size; ++i) data[i] = data[i] * scalefactor + addoffset;
462 }
463 else if (lscalefactor)
464 {
465 for (size_t i = 0; i < size; ++i) data[i] *= scalefactor;
466 }
467 else if (laddoffset)
468 {
469 for (size_t i = 0; i < size; ++i) data[i] += addoffset;
470 }
471 }
472
473 static
cdfCreateRecords(stream_t * streamptr,int tsID)474 void cdfCreateRecords(stream_t *streamptr, int tsID)
475 {
476 if ( tsID < 0 || (tsID >= streamptr->ntsteps && tsID > 0) ) return;
477
478 if ( streamptr->tsteps[tsID].nallrecs > 0 ) return;
479
480 int vlistID = streamptr->vlistID;
481
482 tsteps_t *sourceTstep = streamptr->tsteps;
483 tsteps_t *destTstep = sourceTstep + tsID;
484
485 const int nvars = vlistNvars(vlistID);
486 const int nrecs = vlistNrecs(vlistID);
487 if ( nrecs <= 0 ) return;
488
489 if ( tsID == 0 )
490 {
491 const int nvrecs = nrecs; /* use all records at first timestep */
492
493 streamptr->nrecs += nrecs;
494
495 destTstep->records = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
496 destTstep->nrecs = nrecs;
497 destTstep->nallrecs = nrecs;
498 destTstep->recordSize = nrecs;
499 destTstep->curRecID = CDI_UNDEFID;
500 destTstep->recIDs = (int *) Malloc((size_t)nvrecs*sizeof (int));;
501 for ( int recID = 0; recID < nvrecs; recID++ ) destTstep->recIDs[recID] = recID;
502
503 record_t *records = destTstep->records;
504
505 for ( int varID = 0, recID = 0; varID < nvars; varID++ )
506 {
507 const int zaxisID = vlistInqVarZaxis(vlistID, varID);
508 const int nlev = zaxisInqSize(zaxisID);
509 for ( int levelID = 0; levelID < nlev; levelID++ )
510 {
511 recordInitEntry(&records[recID]);
512 records[recID].varID = (short)varID;
513 records[recID].levelID = levelID;
514 recID++;
515 }
516 }
517 }
518 else if ( tsID == 1 )
519 {
520 int nvrecs = 0;
521 for ( int varID = 0; varID < nvars; varID++ )
522 {
523 if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
524 {
525 int zaxisID = vlistInqVarZaxis(vlistID, varID);
526 nvrecs += zaxisInqSize(zaxisID);
527 }
528 }
529
530 streamptr->nrecs += nvrecs;
531
532 destTstep->records = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
533 destTstep->nrecs = nvrecs;
534 destTstep->nallrecs = nrecs;
535 destTstep->recordSize = nrecs;
536 destTstep->curRecID = CDI_UNDEFID;
537
538 memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));
539
540 if ( nvrecs )
541 {
542 destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof (int));
543 for ( int recID = 0, vrecID = 0; recID < nrecs; recID++ )
544 {
545 int varID = destTstep->records[recID].varID;
546 if ( vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT )
547 {
548 destTstep->recIDs[vrecID++] = recID;
549 }
550 }
551 }
552 }
553 else
554 {
555 if ( streamptr->tsteps[1].records == 0 ) cdfCreateRecords(streamptr, 1);
556
557 const int nvrecs = streamptr->tsteps[1].nrecs;
558
559 streamptr->nrecs += nvrecs;
560
561 destTstep->records = (record_t *) Malloc((size_t)nrecs*sizeof(record_t));
562 destTstep->nrecs = nvrecs;
563 destTstep->nallrecs = nrecs;
564 destTstep->recordSize = nrecs;
565 destTstep->curRecID = CDI_UNDEFID;
566
567 memcpy(destTstep->records, sourceTstep->records, (size_t)nrecs*sizeof(record_t));
568
569 destTstep->recIDs = (int *) Malloc((size_t)nvrecs * sizeof(int));
570
571 memcpy(destTstep->recIDs, streamptr->tsteps[1].recIDs, (size_t)nvrecs*sizeof(int));
572 }
573 }
574
575 static
cdf_time_dimid(int fileID,int ndims,int nvars,ncdim_t * ncdims)576 int cdf_time_dimid(int fileID, int ndims, int nvars, ncdim_t *ncdims)
577 {
578 char dimname[CDI_MAX_NAME];
579
580 for ( int dimid = 0; dimid < ndims; ++dimid )
581 {
582 dimname[0] = 0;
583 cdf_inq_dimname(fileID, ncdims[dimid].dimid, dimname);
584 if ( strIsEqual(dimname, "time") || strIsEqual(dimname, "Time") ) return dimid;
585 }
586
587 for ( int varid = 0; varid < nvars; ++varid )
588 {
589 nc_type xtype;
590 int nvdims, nvatts, dimids[9];
591 cdf_inq_var(fileID, varid, NULL, &xtype, &nvdims, dimids, &nvatts);
592 for ( int i = 0; i < nvdims; ++i )
593 for ( int dimid = 0; dimid < ndims; ++dimid )
594 if ( ncdims[dimid].dimid == dimids[i] )
595 {
596 dimids[i] = dimid;
597 break;
598 }
599
600 if ( nvdims == 1 )
601 {
602 char sbuf[CDI_MAX_NAME];
603 for ( int iatt = 0; iatt < nvatts; ++iatt )
604 {
605 sbuf[0] = 0;
606 cdf_inq_attname(fileID, varid, iatt, sbuf);
607 if (strStartsWith(sbuf, "units"))
608 {
609 cdfGetAttText(fileID, varid, "units", sizeof(sbuf), sbuf);
610 strToLower(sbuf);
611
612 if ( is_time_units(sbuf) ) return dimids[0];
613 }
614 }
615 }
616 }
617
618 return CDI_UNDEFID;
619 }
620
621 static
init_ncdims(int ndims,ncdim_t * ncdims)622 void init_ncdims(int ndims, ncdim_t *ncdims)
623 {
624 for (int ncdimid = 0; ncdimid < ndims; ncdimid++)
625 {
626 ncdim_t *ncdim = &ncdims[ncdimid];
627 ncdim->dimid = CDI_UNDEFID;
628 ncdim->ncvarid = CDI_UNDEFID;
629 ncdim->dimtype = CDI_UNDEFID;
630 ncdim->len = 0;
631 ncdim->name[0] = 0;
632 }
633 }
634
635 static
init_ncvars(int nvars,ncvar_t * ncvars)636 void init_ncvars(int nvars, ncvar_t *ncvars)
637 {
638 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
639 {
640 ncvar_t *ncvar = &ncvars[ncvarid];
641 ncvar->ncid = CDI_UNDEFID;
642 ncvar->varStatus = UndefVar;
643 ncvar->ignoreVar = false;
644 ncvar->isCubeSphere = false;
645 ncvar->isCharAxis = false;
646 ncvar->isXaxis = false;
647 ncvar->isYaxis = false;
648 ncvar->isZaxis = false;
649 ncvar->isTaxis = false;
650 ncvar->isLon = false;
651 ncvar->isLat = false;
652 ncvar->isClimatology = false;
653 ncvar->hasCalendar = false;
654 ncvar->hasFormulaterms = false;
655 ncvar->printWarning = true;
656 ncvar->timetype = TIME_CONSTANT;
657 ncvar->param = CDI_UNDEFID;
658 ncvar->code = CDI_UNDEFID;
659 ncvar->tabnum = 0;
660 ncvar->bounds = CDI_UNDEFID;
661 ncvar->gridID = CDI_UNDEFID;
662 ncvar->zaxisID = CDI_UNDEFID;
663 ncvar->gridtype = CDI_UNDEFID;
664 ncvar->zaxistype = CDI_UNDEFID;
665 ncvar->xdim = CDI_UNDEFID;
666 ncvar->ydim = CDI_UNDEFID;
667 ncvar->zdim = CDI_UNDEFID;
668 ncvar->xvarid = CDI_UNDEFID;
669 ncvar->yvarid = CDI_UNDEFID;
670 ncvar->rpvarid = CDI_UNDEFID;
671 ncvar->zvarid = CDI_UNDEFID;
672 ncvar->tvarid = CDI_UNDEFID;
673 ncvar->psvarid = CDI_UNDEFID;
674 ncvar->p0varid = CDI_UNDEFID;
675 ncvar->ncoordvars = 0;
676 for (int i = 0; i < MAX_COORDVARS; ++i) ncvar->cvarids[i] = CDI_UNDEFID;
677 for (int i = 0; i < MAX_COORDVARS; ++i) ncvar->coordvarids[i] = CDI_UNDEFID;
678 for (int i = 0; i < MAX_AUXVARS; ++i) ncvar->auxvarids[i] = CDI_UNDEFID;
679 ncvar->nauxvars = 0;
680 ncvar->cellarea = CDI_UNDEFID;
681 ncvar->tableID = CDI_UNDEFID;
682 ncvar->truncation = 0;
683 ncvar->position = 0;
684 ncvar->numLPE = 0;
685 ncvar->defmissval = false;
686 ncvar->deffillval = false;
687 ncvar->xtype = 0;
688 ncvar->gmapid = CDI_UNDEFID;
689 ncvar->positive = 0;
690 ncvar->ndims = 0;
691 for (int i = 0; i < MAX_DIMS_CDF; ++i) ncvar->dimids[i] = CDI_UNDEFID;
692 for (int i = 0; i < MAX_DIMS_CDF; ++i) ncvar->dimtype[i] = CDI_UNDEFID;
693 for (int i = 0; i < MAX_DIMS_CDF; ++i) ncvar->chunks[i] = 0;
694 ncvar->chunked = false;
695 ncvar->chunktype = CDI_UNDEFID;
696 ncvar->natts = 0;
697 ncvar->atts = NULL;
698 ncvar->vctsize = 0;
699 ncvar->vct = NULL;
700 ncvar->missval = 0;
701 ncvar->fillval = 0;
702 ncvar->addoffset = 0;
703 ncvar->scalefactor = 1;
704 ncvar->ldeflate = false;
705 ncvar->lszip = false;
706 ncvar->lunsigned = false;
707 ncvar->lvalidrange = false;
708 ncvar->validrange[0] = VALIDMISS;
709 ncvar->validrange[1] = VALIDMISS;
710 ncvar->typeOfEnsembleForecast = -1;
711 ncvar->numberOfForecastsInEnsemble = -1;
712 ncvar->perturbationNumber = -1;
713 memset(ncvar->name, 0, CDI_MAX_NAME);
714 memset(ncvar->longname, 0, CDI_MAX_NAME);
715 memset(ncvar->stdname, 0, CDI_MAX_NAME);
716 memset(ncvar->units, 0, CDI_MAX_NAME);
717 memset(ncvar->extra, 0, CDI_MAX_NAME);
718 }
719 }
720
721 static
cdf_set_var(ncvar_t * ncvar,int varStatus)722 void cdf_set_var(ncvar_t *ncvar, int varStatus)
723 {
724 if (ncvar->varStatus != UndefVar && ncvar->varStatus != varStatus && ncvar->printWarning)
725 {
726 if (!ncvar->ignoreVar) Warning("Inconsistent variable definition for %s!", ncvar->name);
727
728 ncvar->printWarning = false;
729 varStatus = CoordVar;
730 }
731
732 ncvar->varStatus = varStatus;
733 }
734
735 static
cdf_set_dim(ncvar_t * ncvar,int dimid,int dimtype)736 void cdf_set_dim(ncvar_t *ncvar, int dimid, int dimtype)
737 {
738 if (ncvar->dimtype[dimid] != CDI_UNDEFID && ncvar->dimtype[dimid] != dimtype)
739 {
740 Warning("Inconsistent dimension definition for %s! dimid = %d; type = %d; newtype = %d",
741 ncvar->name, dimid, ncvar->dimtype[dimid], dimtype);
742 }
743
744 ncvar->dimtype[dimid] = dimtype;
745 }
746
747 static
scan_hybrid_formulaterms(int ncid,int ncfvarid,int * avarid,int * bvarid,int * psvarid,int * p0varid)748 void scan_hybrid_formulaterms(int ncid, int ncfvarid, int *avarid, int *bvarid, int *psvarid, int *p0varid)
749 {
750 *avarid = -1;
751 *bvarid = -1;
752 *psvarid = -1;
753 *p0varid = -1;
754
755 char attstring[1024];
756 cdfGetAttText(ncid, ncfvarid, "formula_terms", sizeof(attstring), attstring);
757 char *pstring = attstring;
758
759 bool lstop = false;
760 for ( int i = 0; i < 4; i++ )
761 {
762 while ( isspace((int) *pstring) ) pstring++;
763 if ( *pstring == 0 ) break;
764 char *tagname = pstring;
765 while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
766 if ( *pstring == 0 ) lstop = true;
767 *pstring++ = 0;
768
769 while ( isspace((int) *pstring) ) pstring++;
770 if ( *pstring == 0 ) break;
771 char *varname = pstring;
772 while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
773 if ( *pstring == 0 ) lstop = true;
774 *pstring++ = 0;
775
776 int dimvarid;
777 int status_nc = nc_inq_varid(ncid, varname, &dimvarid);
778 if ( status_nc == NC_NOERR )
779 {
780 if ( strIsEqual(tagname, "ap:") ) *avarid = dimvarid;
781 else if ( strIsEqual(tagname, "a:") ) *avarid = dimvarid;
782 else if ( strIsEqual(tagname, "b:") ) *bvarid = dimvarid;
783 else if ( strIsEqual(tagname, "ps:") ) *psvarid = dimvarid;
784 else if ( strIsEqual(tagname, "p0:") ) *p0varid = dimvarid;
785 }
786 else if ( !strIsEqual(tagname, "ps:") )
787 {
788 Warning("%s - %s", nc_strerror(status_nc), varname);
789 }
790
791 if ( lstop ) break;
792 }
793 }
794
795 static
readVCT(int ncid,int ndims2,size_t dimlen,size_t dimlen2,int avarid2,int bvarid2,double * vct)796 void readVCT(int ncid, int ndims2, size_t dimlen, size_t dimlen2, int avarid2, int bvarid2, double *vct)
797 {
798 double *abuf = (double*) Malloc(dimlen*2*sizeof(double));
799 double *bbuf = (double*) Malloc(dimlen*2*sizeof(double));
800 cdf_get_var_double(ncid, avarid2, abuf);
801 cdf_get_var_double(ncid, bvarid2, bbuf);
802
803 if ( ndims2 == 2 )
804 {
805 for ( size_t i = 0; i < dimlen; ++i )
806 {
807 vct[i] = abuf[i*2];
808 vct[i+dimlen+1] = bbuf[i*2];
809 }
810 vct[dimlen] = abuf[dimlen*2-1];
811 vct[dimlen*2+1] = bbuf[dimlen*2-1];
812 }
813 else
814 {
815 for ( size_t i = 0; i < dimlen2; ++i )
816 {
817 vct[i] = abuf[i];
818 vct[i+dimlen+1] = bbuf[i];
819 }
820 }
821
822 Free(abuf);
823 Free(bbuf);
824 }
825
826 static
isHybridSigmaPressureCoordinate(int ncid,int ncvarid,ncvar_t * ncvars,const ncdim_t * ncdims)827 bool isHybridSigmaPressureCoordinate(int ncid, int ncvarid, ncvar_t *ncvars, const ncdim_t *ncdims)
828 {
829 bool status = false;
830 ncvar_t *ncvar = &ncvars[ncvarid];
831
832 if ( strIsEqual(ncvar->stdname, "atmosphere_hybrid_sigma_pressure_coordinate") )
833 {
834 CDI_Convention = CDI_CONVENTION_CF;
835
836 status = true;
837 ncvar->zaxistype = ZAXIS_HYBRID;
838 //int ndims = ncvar->ndims;
839 int dimid = ncvar->dimids[0];
840 size_t dimlen = ncdims[dimid].len;
841 int avarid1 = -1, bvarid1 = -1, psvarid1 = -1, p0varid1 = -1;
842 int ncfvarid = ncvarid;
843 if ( ncvars[ncfvarid].hasFormulaterms )
844 scan_hybrid_formulaterms(ncid, ncfvarid, &avarid1, &bvarid1, &psvarid1, &p0varid1);
845 // printf("avarid1, bvarid1, psvarid1, p0varid1 %d %d %d %d\n", avarid1, bvarid1, psvarid1, p0varid1);
846 if ( avarid1 != -1 ) ncvars[avarid1].varStatus = CoordVar;
847 if ( bvarid1 != -1 ) ncvars[bvarid1].varStatus = CoordVar;
848 if ( psvarid1 != -1 ) ncvar->psvarid = psvarid1;
849 if ( p0varid1 != -1 ) ncvar->p0varid = p0varid1;
850
851 if ( ncvar->bounds != CDI_UNDEFID && ncvars[ncvar->bounds].hasFormulaterms )
852 {
853 ncfvarid = ncvar->bounds;
854 int avarid2 = -1, bvarid2 = -1, psvarid2 = -1, p0varid2 = -1;
855 if ( ncvars[ncfvarid].hasFormulaterms )
856 scan_hybrid_formulaterms(ncid, ncfvarid, &avarid2, &bvarid2, &psvarid2, &p0varid2);
857 // printf("avarid2, bvarid2, psvarid2, p0varid2 %d %d %d %d\n", avarid2, bvarid2, psvarid2, p0varid2);
858 if ( avarid2 != -1 && bvarid2 != -1 )
859 {
860 ncvars[avarid2].varStatus = CoordVar;
861 ncvars[bvarid2].varStatus = CoordVar;
862
863 const int ndims2 = ncvars[avarid2].ndims;
864 const int dimid2 = ncvars[avarid2].dimids[0];
865 const size_t dimlen2 = ncdims[dimid2].len;
866
867 if ( (ndims2 == 2 && dimid == ncvars[avarid2].dimids[0] ) ||
868 (ndims2 == 1 && dimlen == dimlen2-1 ) )
869 {
870 double px = 1;
871 if ( p0varid1 != -1 && p0varid1 == p0varid2 )
872 cdf_get_var_double(ncid, p0varid2, &px);
873
874 const size_t vctsize = (dimlen+1)*2;
875 double *vct = (double *) Malloc(vctsize*sizeof(double));
876
877 readVCT(ncid, ndims2, dimlen, dimlen2, avarid2, bvarid2, vct);
878
879 if ( p0varid1 != -1 && IS_NOT_EQUAL(px, 1) )
880 for ( size_t i = 0; i < dimlen+1; ++i ) vct[i] *= px;
881
882 ncvar->vct = vct;
883 ncvar->vctsize = vctsize;
884 }
885 }
886 }
887 }
888
889 return status;
890 }
891
892 static
cdf_set_cdi_attr(int ncid,int ncvarid,int attnum,int cdiID,int varID)893 void cdf_set_cdi_attr(int ncid, int ncvarid, int attnum, int cdiID, int varID)
894 {
895 nc_type atttype;
896 size_t attlen;
897 char attname[CDI_MAX_NAME];
898
899 cdf_inq_attname(ncid, ncvarid, attnum, attname);
900 cdf_inq_attlen(ncid, ncvarid, attname, &attlen);
901 cdf_inq_atttype(ncid, ncvarid, attname, &atttype);
902 if (xtypeIsInt(atttype))
903 {
904 int attint;
905 int *pattint = (attlen > 1) ? (int*) malloc(attlen*sizeof(int)) : &attint;
906 cdfGetAttInt(ncid, ncvarid, attname, attlen, pattint);
907 int datatype = (atttype == NC_SHORT) ? CDI_DATATYPE_INT16 :
908 (atttype == NC_BYTE) ? CDI_DATATYPE_INT8 :
909 #ifdef HAVE_NETCDF4
910 (atttype == NC_UBYTE) ? CDI_DATATYPE_UINT8 :
911 (atttype == NC_USHORT) ? CDI_DATATYPE_UINT16 :
912 (atttype == NC_UINT) ? CDI_DATATYPE_UINT32 :
913 #endif
914 CDI_DATATYPE_INT32;
915 cdiDefAttInt(cdiID, varID, attname, datatype, (int)attlen, pattint);
916 if (attlen > 1) free(pattint);
917 }
918 else if (xtypeIsFloat(atttype))
919 {
920 double attflt;
921 double *pattflt = (attlen > 1) ? (double*) malloc(attlen*sizeof(double)) : &attflt;
922 cdfGetAttDouble(ncid, ncvarid, attname, attlen, pattflt);
923 const int datatype = (atttype == NC_FLOAT) ? CDI_DATATYPE_FLT32 : CDI_DATATYPE_FLT64;
924 cdiDefAttFlt(cdiID, varID, attname, datatype, (int)attlen, pattflt);
925 if (attlen > 1) free(pattflt);
926 }
927 else if (xtypeIsText(atttype))
928 {
929 char attstring[8192];
930 cdfGetAttText(ncid, ncvarid, attname, sizeof(attstring), attstring);
931 cdiDefAttTxt(cdiID, varID, attname, strlen(attstring), attstring);
932 }
933 }
934
935 static
cdf_print_vars(const ncvar_t * ncvars,int nvars,const char * oname)936 void cdf_print_vars(const ncvar_t *ncvars, int nvars, const char *oname)
937 {
938 char axis[7];
939 enum { TAXIS = 't', ZAXIS = 'z', EAXIS = 'e', YAXIS = 'y', XAXIS = 'x' };
940
941 fprintf(stderr, "%s:\n", oname);
942
943 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
944 {
945 const ncvar_t *ncvar = &ncvars[ncvarid];
946 int ndim = 0;
947 if (ncvar->varStatus == DataVar)
948 {
949 axis[ndim++] = 'v';
950 axis[ndim++] = ':';
951 for (int i = 0; i < ncvar->ndims; i++)
952 {
953 if (ncvar->dimtype[i] == T_AXIS) axis[ndim++] = TAXIS;
954 else if (ncvar->dimtype[i] == Z_AXIS) axis[ndim++] = ZAXIS;
955 else if (ncvar->dimtype[i] == E_AXIS) axis[ndim++] = EAXIS;
956 else if (ncvar->dimtype[i] == Y_AXIS) axis[ndim++] = YAXIS;
957 else if (ncvar->dimtype[i] == X_AXIS) axis[ndim++] = XAXIS;
958 else axis[ndim++] = '?';
959 }
960 }
961 else
962 {
963 axis[ndim++] = 'c';
964 axis[ndim++] = ':';
965 if (ncvar->isTaxis) axis[ndim++] = TAXIS;
966 else if (ncvar->isZaxis) axis[ndim++] = ZAXIS;
967 else if (ncvar->isLat ) axis[ndim++] = YAXIS;
968 else if (ncvar->isYaxis) axis[ndim++] = YAXIS;
969 else if (ncvar->isLon ) axis[ndim++] = XAXIS;
970 else if (ncvar->isXaxis) axis[ndim++] = XAXIS;
971 else axis[ndim++] = '?';
972 }
973
974 axis[ndim++] = 0;
975
976 fprintf(stderr, "%3d %3d %-6s %s\n", ncvarid, ndim-3, axis, ncvar->name);
977 }
978 }
979
980 static
cdfScanAttrAxis(ncvar_t * ncvars,ncdim_t * ncdims,int ncvarid,const char * attstring,int nvdims,const int * dimidsp)981 void cdfScanAttrAxis(ncvar_t *ncvars, ncdim_t *ncdims, int ncvarid, const char *attstring, int nvdims, const int *dimidsp)
982 {
983 ncvar_t *ncvar = &ncvars[ncvarid];
984 int attlen = (int) strlen(attstring);
985
986 if (nvdims == 0 && attlen == 1)
987 {
988 if (attstring[0] == 'z')
989 {
990 cdf_set_var(ncvar, CoordVar);
991 ncvar->isZaxis = true;
992 return;
993 }
994 }
995
996 if (attlen != nvdims) return;
997
998 for (int i = 0; i < attlen; ++i)
999 if (attstring[i] != '-' && attstring[i] != 't' && attstring[i] != 'z' && attstring[i] != 'y' && attstring[i] != 'x')
1000 return;
1001
1002 while (attlen--)
1003 {
1004 const int attchr = attstring[attlen];
1005 if (attchr == 't')
1006 {
1007 if (attlen != 0) Warning("axis attribute 't' not on first position!");
1008 cdf_set_dim(ncvar, attlen, T_AXIS);
1009 }
1010 else if (attchr == 'z')
1011 {
1012 ncvar->zdim = dimidsp[attlen];
1013 cdf_set_dim(ncvar, attlen, Z_AXIS);
1014
1015 if (ncvar->ndims == 1)
1016 {
1017 cdf_set_var(ncvar, CoordVar);
1018 ncdims[ncvar->dimids[0]].dimtype = Z_AXIS;
1019 }
1020 }
1021 else if (attchr == 'y')
1022 {
1023 ncvar->ydim = dimidsp[attlen];
1024 cdf_set_dim(ncvar, attlen, Y_AXIS);
1025
1026 if (ncvar->ndims == 1)
1027 {
1028 cdf_set_var(ncvar, CoordVar);
1029 ncdims[ncvar->dimids[0]].dimtype = Y_AXIS;
1030 }
1031 }
1032 else if (attchr == 'x')
1033 {
1034 ncvar->xdim = dimidsp[attlen];
1035 cdf_set_dim(ncvar, attlen, X_AXIS);
1036
1037 if (ncvar->ndims == 1)
1038 {
1039 cdf_set_var(ncvar, CoordVar);
1040 ncdims[ncvar->dimids[0]].dimtype = X_AXIS;
1041 }
1042 }
1043 }
1044 }
1045
1046 static
cdf_get_cell_varid(char * attstring,int ncid)1047 int cdf_get_cell_varid(char *attstring, int ncid)
1048 {
1049 int nc_cell_id = CDI_UNDEFID;
1050
1051 char *pstring = attstring;
1052 while ( isspace((int) *pstring) ) pstring++;
1053 char *cell_measures = pstring;
1054 while ( isalnum((int) *pstring) ) pstring++;
1055 *pstring++ = 0;
1056 while ( isspace((int) *pstring) ) pstring++;
1057 char *cell_var = pstring;
1058 while ( ! isspace((int) *pstring) && *pstring != 0 ) pstring++;
1059 *pstring++ = 0;
1060 /*
1061 printf("cell_measures >%s<\n", cell_measures);
1062 printf("cell_var >%s<\n", cell_var);
1063 */
1064 if ( strStartsWith(cell_measures, "area") )
1065 {
1066 int nc_var_id;
1067 int status = nc_inq_varid(ncid, cell_var, &nc_var_id);
1068 if ( status == NC_NOERR )
1069 nc_cell_id = nc_var_id;
1070 /*
1071 else
1072 Warning("%s - %s", nc_strerror(status), cell_var);
1073 */
1074 }
1075
1076 return nc_cell_id;
1077 }
1078
1079 static
cdfScanVarAttr(int nvars,ncvar_t * ncvars,int ndims,ncdim_t * ncdims,int timedimid,int modelID,int format)1080 void cdfScanVarAttr(int nvars, ncvar_t *ncvars, int ndims, ncdim_t *ncdims, int timedimid, int modelID, int format)
1081 {
1082 int ncdimid;
1083 int nvdims, nvatts;
1084 int iatt;
1085 nc_type xtype, atttype;
1086 size_t attlen;
1087 char name[CDI_MAX_NAME];
1088 char attname[CDI_MAX_NAME];
1089 char attstring[8192];
1090
1091 int nchecked_vars = 0;
1092 enum { max_check_vars = 9 };
1093 char *checked_vars[max_check_vars];
1094 for ( int i = 0; i < max_check_vars; ++i ) checked_vars[i] = NULL;
1095
1096 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1097 {
1098 ncvar_t *ncvar = &ncvars[ncvarid];
1099 const int ncid = ncvar->ncid;
1100 int *dimidsp = ncvar->dimids;
1101
1102 cdf_inq_var(ncid, ncvarid, name, &xtype, &nvdims, dimidsp, &nvatts);
1103 for ( int i = 0; i < nvdims; ++i )
1104 for ( int dimid = 0; dimid < ndims; ++dimid )
1105 if ( ncdims[dimid].dimid == dimidsp[i] )
1106 {
1107 dimidsp[i] = dimid;
1108 break;
1109 }
1110 strcpy(ncvar->name, name);
1111
1112 for ( ncdimid = 0; ncdimid < nvdims; ncdimid++ )
1113 ncvar->dimtype[ncdimid] = -1;
1114
1115 ncvar->xtype = xtype;
1116 ncvar->ndims = nvdims;
1117
1118 #ifdef HAVE_NETCDF4
1119 if ( format == NC_FORMAT_NETCDF4_CLASSIC || format == NC_FORMAT_NETCDF4 )
1120 {
1121 int shuffle = 0, deflate = 0, deflate_level = 0;
1122 nc_inq_var_deflate(ncid, ncvarid, &shuffle, &deflate, &deflate_level);
1123 if (deflate > 0) ncvar->ldeflate = true;
1124
1125 #ifdef HAVE_NC_DEF_VAR_SZIP
1126 int options_mask = 0, pixels_per_block = 0;
1127 nc_inq_var_szip(ncid, ncvarid, &options_mask, &pixels_per_block);
1128 if (options_mask && pixels_per_block) ncvar->lszip = true;
1129 #endif
1130 /*
1131 size_t cache_size, nelems;
1132 float preemption;
1133 nc_get_chunk_cache(&cache_size, &nelems, &preemption);
1134 printf("cache_size %lu nelems %lu preemption %g\n", cache_size, nelems, preemption);
1135 nc_get_var_chunk_cache(ncid, ncvarid, &cache_size, &nelems, &preemption);
1136 printf("varid %d cache_size %lu nelems %lu preemption %g\n", ncvarid, cache_size, nelems, preemption);
1137 */
1138 size_t chunks[nvdims];
1139 int storage_in;
1140 if (nc_inq_var_chunking(ncid, ncvarid, &storage_in, chunks) == NC_NOERR)
1141 {
1142 if (storage_in == NC_CHUNKED)
1143 {
1144 ncvar->chunked = true;
1145 for (int i = 0; i < nvdims; ++i) ncvar->chunks[i] = chunks[i];
1146 if (CDI_Debug)
1147 {
1148 fprintf(stderr, "%s: chunking %d %d %d chunks ", name, storage_in, NC_CONTIGUOUS, NC_CHUNKED);
1149 for (int i = 0; i < nvdims; ++i) fprintf(stderr, "%zu ", chunks[i]);
1150 fprintf(stderr, "\n");
1151 }
1152 {
1153 char *buf = ncvar->extra;
1154 size_t pos = strlen(buf);
1155 static const char prefix[] = "chunks=";
1156 memcpy(buf + pos, prefix, sizeof (prefix));
1157 pos += sizeof (prefix) - 1;
1158 for ( int i = nvdims-1; i >= 0; --i )
1159 {
1160 pos += (size_t)(sprintf(buf + pos, "%zu%s", chunks[i], i > 0 ? "x" : ""));
1161 }
1162 buf[pos] = ' '; buf[pos + 1] = 0;
1163 }
1164 }
1165 }
1166 }
1167 #endif
1168
1169 if ( nvdims > 0 )
1170 {
1171 if ( timedimid == dimidsp[0] )
1172 {
1173 ncvar->timetype = TIME_VARYING;
1174 cdf_set_dim(ncvar, 0, T_AXIS);
1175 }
1176 else
1177 {
1178 for ( ncdimid = 1; ncdimid < nvdims; ncdimid++ )
1179 {
1180 if ( timedimid == dimidsp[ncdimid] )
1181 {
1182 Warning("Time must be the first dimension! Unsupported array structure, skipped variable %s!", ncvar->name);
1183 ncvar->varStatus = CoordVar;
1184 }
1185 }
1186 }
1187 }
1188
1189 if (ncvar->natts == 0 && nvatts > 0)
1190 ncvar->atts = (int*) Malloc((size_t)nvatts*sizeof(int));
1191
1192 for ( iatt = 0; iatt < nvatts; iatt++ )
1193 {
1194 int nc_cell_id = CDI_UNDEFID;
1195
1196 cdf_inq_attname(ncid, ncvarid, iatt, attname);
1197 cdf_inq_atttype(ncid, ncvarid, attname, &atttype);
1198 cdf_inq_attlen(ncid, ncvarid, attname, &attlen);
1199
1200 size_t attstringsize = sizeof(attstring);
1201 const bool isText = xtypeIsText(atttype);
1202 const bool isNumber = (xtypeIsFloat(atttype) || xtypeIsInt(atttype));
1203 if (isText)
1204 {
1205 cdfGetAttText(ncid, ncvarid, attname, sizeof(attstring), attstring);
1206 attstringsize = strlen(attstring) + 1;
1207 if ( attstringsize > CDI_MAX_NAME ) attstringsize = CDI_MAX_NAME;
1208 }
1209
1210 if ( isText && strIsEqual(attname, "long_name") )
1211 {
1212 memcpy(ncvar->longname, attstring, attstringsize);
1213 }
1214 else if ( isText && strIsEqual(attname, "standard_name") )
1215 {
1216 memcpy(ncvar->stdname, attstring, attstringsize);
1217 }
1218 else if ( isText && strIsEqual(attname, "units") )
1219 {
1220 memcpy(ncvar->units, attstring, attstringsize);
1221 }
1222 else if ( isText && strIsEqual(attname, "calendar") )
1223 {
1224 ncvar->hasCalendar = true;
1225 }
1226 else if ( isText && strIsEqual(attname, "param") )
1227 {
1228 int pnum = 0, pcat = 255, pdis = 255;
1229 sscanf(attstring, "%d.%d.%d", &pnum, &pcat, &pdis);
1230 ncvar->param = cdiEncodeParam(pnum, pcat, pdis);
1231 cdf_set_var(ncvar, DataVar);
1232 }
1233 else if ( isNumber && strIsEqual(attname, "code") )
1234 {
1235 cdfGetAttInt(ncid, ncvarid, attname, 1, &ncvar->code);
1236 cdf_set_var(ncvar, DataVar);
1237 }
1238 else if ( isNumber && strIsEqual(attname, "table") )
1239 {
1240 int tablenum;
1241 cdfGetAttInt(ncid, ncvarid, attname, 1, &tablenum);
1242 if ( tablenum > 0 )
1243 {
1244 ncvar->tabnum = tablenum;
1245 ncvar->tableID = tableInq(modelID, tablenum, NULL);
1246 if ( ncvar->tableID == CDI_UNDEFID )
1247 ncvar->tableID = tableDef(modelID, tablenum, NULL);
1248 }
1249 cdf_set_var(ncvar, DataVar);
1250 }
1251 else if ( isText && strIsEqual(attname, "trunc_type") )
1252 {
1253 if ( strIsEqual(attstring, "Triangular") )
1254 ncvar->gridtype = GRID_SPECTRAL;
1255 }
1256 else if ( isText && (strIsEqual(attname, "grid_type") || strIsEqual(attname, "CDI_grid_type")) )
1257 {
1258 strToLower(attstring);
1259 cdf_set_gridtype(attstring, &ncvar->gridtype);
1260 cdf_set_var(ncvar, DataVar);
1261 }
1262 else if ( isText && strIsEqual(attname, "CDI_grid_latitudes") )
1263 {
1264 int nc_yvar_id;
1265 int status = nc_inq_varid(ncid, attstring, &nc_yvar_id);
1266 if ( status == NC_NOERR )
1267 {
1268 ncvar->yvarid = nc_yvar_id;
1269 cdf_set_var(&ncvars[ncvar->yvarid], CoordVar);
1270 }
1271 else
1272 Warning("%s - %s", nc_strerror(status), attstring);
1273
1274 cdf_set_var(ncvar, DataVar);
1275 }
1276 else if ( isText && strIsEqual(attname, "CDI_grid_reduced_points") )
1277 {
1278 int nc_rpvar_id;
1279 int status = nc_inq_varid(ncid, attstring, &nc_rpvar_id);
1280 if ( status == NC_NOERR )
1281 {
1282 ncvar->rpvarid = nc_rpvar_id;
1283 cdf_set_var(&ncvars[ncvar->rpvarid], CoordVar);
1284 }
1285 else
1286 Warning("%s - %s", nc_strerror(status), attstring);
1287
1288 cdf_set_var(ncvar, DataVar);
1289 }
1290 else if ( isNumber && strIsEqual(attname, "CDI_grid_num_LPE") )
1291 {
1292 cdfGetAttInt(ncid, ncvarid, attname, 1, &ncvar->numLPE);
1293 }
1294 else if ( isText && strIsEqual(attname, "level_type") )
1295 {
1296 strToLower(attstring);
1297 cdf_set_zaxistype(attstring, &ncvar->zaxistype);
1298 cdf_set_var(ncvar, DataVar);
1299 }
1300 else if ( isNumber && strIsEqual(attname, "trunc_count") )
1301 {
1302 cdfGetAttInt(ncid, ncvarid, attname, 1, &ncvar->truncation);
1303 }
1304 else if ( isNumber && strIsEqual(attname, "truncation") )
1305 {
1306 cdfGetAttInt(ncid, ncvarid, attname, 1, &ncvar->truncation);
1307 }
1308 else if ( isNumber && strIsEqual(attname, "number_of_grid_in_reference") )
1309 {
1310 cdfGetAttInt(ncid, ncvarid, attname, 1, &ncvar->position);
1311 }
1312 else if ( isNumber && strIsEqual(attname, "add_offset") )
1313 {
1314 cdfGetAttDouble(ncid, ncvarid, attname, 1, &ncvar->addoffset);
1315 /*
1316 if ( atttype != NC_BYTE && atttype != NC_SHORT && atttype != NC_INT )
1317 if ( ncvar->addoffset != 0 )
1318 Warning("attribute add_offset not supported for atttype %d", atttype);
1319 */
1320 // (also used for lon/lat) cdf_set_var(ncvar, DataVar);
1321 }
1322 else if ( isNumber && strIsEqual(attname, "scale_factor") )
1323 {
1324 cdfGetAttDouble(ncid, ncvarid, attname, 1, &ncvar->scalefactor);
1325 /*
1326 if ( atttype != NC_BYTE && atttype != NC_SHORT && atttype != NC_INT )
1327 if ( ncvar->scalefactor != 1 )
1328 Warning("attribute scale_factor not supported for atttype %d", atttype);
1329 */
1330 // (also used for lon/lat) cdf_set_var(ncvar, DataVar);
1331 }
1332 else if ( isText && strIsEqual(attname, "climatology") )
1333 {
1334 int ncboundsid;
1335 int status = nc_inq_varid(ncid, attstring, &ncboundsid);
1336 if ( status == NC_NOERR )
1337 {
1338 ncvar->isClimatology = true;
1339 ncvar->bounds = ncboundsid;
1340 cdf_set_var(&ncvars[ncvar->bounds], CoordVar);
1341 cdf_set_var(ncvar, CoordVar);
1342 }
1343 else
1344 Warning("%s - %s", nc_strerror(status), attstring);
1345 }
1346 else if ( isText && strIsEqual(attname, "bounds") )
1347 {
1348 int ncboundsid;
1349 int status = nc_inq_varid(ncid, attstring, &ncboundsid);
1350 if ( status == NC_NOERR )
1351 {
1352 ncvar->bounds = ncboundsid;
1353 cdf_set_var(&ncvars[ncvar->bounds], CoordVar);
1354 cdf_set_var(ncvar, CoordVar);
1355 }
1356 else
1357 Warning("%s - %s", nc_strerror(status), attstring);
1358 }
1359 else if ( isText && strIsEqual(attname, "formula_terms") )
1360 {
1361 ncvar->hasFormulaterms = true;
1362 }
1363 else if ( isText && strIsEqual(attname, "cell_measures") && (nc_cell_id=cdf_get_cell_varid(attstring, ncid)) != CDI_UNDEFID )
1364 {
1365 ncvar->cellarea = nc_cell_id;
1366 ncvars[nc_cell_id].varStatus = CoordVar;
1367 cdf_set_var(ncvar, DataVar);
1368 }
1369 else if ( isText && (strIsEqual(attname, "associate") || strIsEqual(attname, "coordinates")) )
1370 {
1371 bool lstop = false;
1372 char *pstring = attstring;
1373
1374 for ( int i = 0; i < MAX_COORDVARS; i++ )
1375 {
1376 while ( isspace((int) *pstring) ) pstring++;
1377 if ( *pstring == 0 ) break;
1378 char *varname = pstring;
1379 while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
1380 if ( *pstring == 0 ) lstop = true;
1381 if ( *(pstring-1) == ',' ) *(pstring-1) = 0;
1382 *pstring++ = 0;
1383
1384 int dimvarid;
1385 int status = nc_inq_varid(ncid, varname, &dimvarid);
1386 if ( status == NC_NOERR )
1387 {
1388 cdf_set_var(&ncvars[dimvarid], CoordVar);
1389 if ( !CDI_Ignore_Att_Coordinates )
1390 {
1391 ncvar->coordvarids[i] = dimvarid;
1392 ncvar->ncoordvars++;
1393 }
1394 }
1395 else
1396 {
1397 if ( !CDI_Ignore_Att_Coordinates ) ncvar->ncoordvars++;
1398
1399 int k;
1400 for ( k = 0; k < nchecked_vars; ++k )
1401 if ( strIsEqual(checked_vars[k], varname) ) break;
1402
1403 if ( k == nchecked_vars )
1404 {
1405 if ( nchecked_vars < max_check_vars ) checked_vars[nchecked_vars++] = strdup(varname);
1406 Warning("%s - >%s<", nc_strerror(status), varname);
1407 }
1408 }
1409
1410 if ( lstop ) break;
1411 }
1412
1413 cdf_set_var(ncvar, DataVar);
1414 }
1415 else if ( isText && strIsEqual(attname, "auxiliary_variable") )
1416 {
1417 bool lstop = false;
1418 char *pstring = attstring;
1419
1420 for ( int i = 0; i < MAX_AUXVARS; i++ )
1421 {
1422 while ( isspace((int) *pstring) ) pstring++;
1423 if ( *pstring == 0 ) break;
1424 char *varname = pstring;
1425 while ( !isspace((int) *pstring) && *pstring != 0 ) pstring++;
1426 if ( *pstring == 0 ) lstop = true;
1427 *pstring++ = 0;
1428
1429 int dimvarid;
1430 int status = nc_inq_varid(ncid, varname, &dimvarid);
1431 if ( status == NC_NOERR )
1432 {
1433 cdf_set_var(&ncvars[dimvarid], CoordVar);
1434 // if ( !CDI_Ignore_Att_Coordinates )
1435 {
1436 ncvar->auxvarids[i] = dimvarid;
1437 ncvar->nauxvars++;
1438 }
1439 }
1440 else
1441 Warning("%s - %s", nc_strerror(status), varname);
1442
1443 if ( lstop ) break;
1444 }
1445
1446 cdf_set_var(ncvar, DataVar);
1447 }
1448 else if ( isText && strIsEqual(attname, "grid_mapping") )
1449 {
1450 int nc_gmap_id;
1451 int status = nc_inq_varid(ncid, attstring, &nc_gmap_id);
1452 if ( status == NC_NOERR )
1453 {
1454 ncvar->gmapid = nc_gmap_id;
1455 cdf_set_var(&ncvars[ncvar->gmapid], CoordVar);
1456 }
1457 else
1458 Warning("%s - %s", nc_strerror(status), attstring);
1459
1460 cdf_set_var(ncvar, DataVar);
1461 }
1462 else if ( isText && strIsEqual(attname, "positive") )
1463 {
1464 strToLower(attstring);
1465
1466 if ( strIsEqual(attstring, "down") ) ncvar->positive = POSITIVE_DOWN;
1467 else if ( strIsEqual(attstring, "up") ) ncvar->positive = POSITIVE_UP;
1468
1469 if ( ncvar->ndims == 1 )
1470 {
1471 cdf_set_var(ncvar, CoordVar);
1472 cdf_set_dim(ncvar, 0, Z_AXIS);
1473 ncdims[ncvar->dimids[0]].dimtype = Z_AXIS;
1474 }
1475 else if ( ncvar->ndims == 0 )
1476 {
1477 cdf_set_var(ncvar, CoordVar);
1478 ncvar->isZaxis = true;
1479 }
1480 else
1481 {
1482 ncvar->atts[ncvar->natts++] = iatt;
1483 }
1484 }
1485 else if ( isNumber && strIsEqual(attname, "_FillValue") )
1486 {
1487 cdfGetAttDouble(ncid, ncvarid, attname, 1, &ncvar->fillval);
1488 ncvar->deffillval = true;
1489 // cdf_set_var(ncvar, DataVar);
1490 }
1491 else if ( isNumber && strIsEqual(attname, "missing_value") )
1492 {
1493 cdfGetAttDouble(ncid, ncvarid, attname, 1, &ncvar->missval);
1494 ncvar->defmissval = true;
1495 // cdf_set_var(ncvar, DataVar);
1496 }
1497 else if ( isNumber && strIsEqual(attname, "valid_range") && attlen == 2 )
1498 {
1499 if ( ncvar->lvalidrange == false )
1500 {
1501 bool lignore = (xtypeIsFloat(atttype) != xtypeIsFloat(xtype));
1502 if ( !CDI_Ignore_Valid_Range && lignore == false )
1503 {
1504 cdfGetAttDouble(ncid, ncvarid, attname, 2, ncvar->validrange);
1505 ncvar->lvalidrange = (ncvar->validrange[0] <= ncvar->validrange[1]);
1506 if ( ((int)ncvar->validrange[0]) == 0 && ((int)ncvar->validrange[1]) == 255 )
1507 ncvar->lunsigned = true;
1508 // cdf_set_var(ncvar, DataVar);
1509 }
1510 else if ( lignore )
1511 {
1512 Warning("Inconsistent data type for attribute %s:valid_range, ignored!", name);
1513 }
1514 }
1515 }
1516 else if ( isNumber && strIsEqual(attname, "valid_min") && attlen == 1 )
1517 {
1518 bool lignore = (xtypeIsFloat(atttype) != xtypeIsFloat(xtype));
1519 if ( !CDI_Ignore_Valid_Range && lignore == false )
1520 {
1521 cdfGetAttDouble(ncid, ncvarid, attname, 1, &(ncvar->validrange)[0]);
1522 ncvar->lvalidrange = true;
1523 }
1524 else if ( lignore )
1525 {
1526 Warning("Inconsistent data type for attribute %s:valid_min, ignored!", name);
1527 }
1528 }
1529 else if ( isNumber && strIsEqual(attname, "valid_max") && attlen == 1 )
1530 {
1531 bool lignore = (xtypeIsFloat(atttype) != xtypeIsFloat(xtype));
1532 if ( !CDI_Ignore_Valid_Range && lignore == false )
1533 {
1534 cdfGetAttDouble(ncid, ncvarid, attname, 1, &(ncvar->validrange)[1]);
1535 ncvar->lvalidrange = true;
1536 }
1537 else if ( lignore )
1538 {
1539 Warning("Inconsistent data type for attribute %s:valid_max, ignored!", name);
1540 }
1541 }
1542 else if ( isText && strIsEqual(attname, "_Unsigned") )
1543 {
1544 strToLower(attstring);
1545 if ( strIsEqual(attstring, "true") )
1546 {
1547 ncvar->lunsigned = true;
1548 /*
1549 ncvar->lvalidrange = true;
1550 ncvar->validrange[0] = 0;
1551 ncvar->validrange[1] = 255;
1552 */
1553 }
1554 // cdf_set_var(ncvar, DataVar);
1555 }
1556 else if ( isText && strIsEqual(attname, "cdi") )
1557 {
1558 strToLower(attstring);
1559
1560 if ( strIsEqual(attstring, "ignore") )
1561 {
1562 ncvar->ignoreVar = true;
1563 cdf_set_var(ncvar, CoordVar);
1564 }
1565 }
1566 else if ( isNumber &&
1567 (strIsEqual(attname, "realization") ||
1568 strIsEqual(attname, "ensemble_members") ||
1569 strIsEqual(attname, "forecast_init_type")) )
1570 {
1571 int temp;
1572 cdfGetAttInt(ncid, ncvarid, attname, 1, &temp);
1573
1574 if ( strIsEqual(attname, "realization") ) ncvar->perturbationNumber = temp;
1575 else if ( strIsEqual(attname, "ensemble_members") ) ncvar->numberOfForecastsInEnsemble = temp;
1576 else if ( strIsEqual(attname, "forecast_init_type") ) ncvar->typeOfEnsembleForecast = temp;
1577
1578 cdf_set_var(ncvar, DataVar);
1579 }
1580 else
1581 {
1582 ncvar->atts[ncvar->natts++] = iatt;
1583 }
1584 }
1585 }
1586
1587 for ( int i = 0; i < max_check_vars; ++i ) if ( checked_vars[i] ) Free(checked_vars[i]);
1588 }
1589
1590 static
cdfVerifyVarAttr(int nvars,ncvar_t * ncvars,ncdim_t * ncdims)1591 void cdfVerifyVarAttr(int nvars, ncvar_t *ncvars, ncdim_t *ncdims)
1592 {
1593 nc_type atttype;
1594 size_t attlen;
1595 char attname[CDI_MAX_NAME];
1596 char attstring[8192];
1597
1598 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1599 {
1600 ncvar_t *ncvar = &ncvars[ncvarid];
1601 const int ncid = ncvar->ncid;
1602 const int *dimidsp = ncvar->dimids;
1603 const int nvdims = ncvar->ndims;
1604 const int nvatts = ncvar->natts;
1605
1606 for (int i = 0; i < nvatts; i++)
1607 {
1608 const int attnum = ncvar->atts[i];
1609 cdf_inq_attname(ncid, ncvarid, attnum, attname);
1610 cdf_inq_attlen(ncid, ncvarid, attname, &attlen);
1611 cdf_inq_atttype(ncid, ncvarid, attname, &atttype);
1612
1613 size_t attstringsize = sizeof(attstring);
1614 const bool isText = xtypeIsText(atttype);
1615 // const bool isNumber = (xtypeIsFloat(atttype) || xtypeIsInt(atttype));
1616
1617 if (isText)
1618 {
1619 cdfGetAttText(ncid, ncvarid, attname, sizeof(attstring), attstring);
1620 attstringsize = strlen(attstring) + 1;
1621 if (attstringsize > CDI_MAX_NAME) attstringsize = CDI_MAX_NAME;
1622
1623 if (strIsEqual(attname, "axis"))
1624 {
1625 cdfScanAttrAxis(ncvars, ncdims, ncvarid, strToLower(attstring), nvdims, dimidsp);
1626 }
1627 /*
1628 else if (strIsEqual(attname, "standard_name"))
1629 {
1630 memcpy(ncvar->stdname, attstring, attstringsize);
1631 }
1632 */
1633 }
1634 }
1635 }
1636 }
1637
1638
1639 static
cdf_set_dimtype(int nvars,ncvar_t * ncvars,ncdim_t * ncdims)1640 void cdf_set_dimtype(int nvars, ncvar_t *ncvars, ncdim_t *ncdims)
1641 {
1642 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1643 {
1644 ncvar_t *ncvar = &ncvars[ncvarid];
1645 if (ncvar->varStatus == DataVar)
1646 {
1647 int ndims = ncvar->ndims;
1648 for (int i = 0; i < ndims; i++)
1649 {
1650 int ncdimid = ncvar->dimids[i];
1651 int dimtype = ncdims[ncdimid].dimtype;
1652 if (dimtype >= X_AXIS && dimtype <= T_AXIS)
1653 cdf_set_dim(ncvar, i, dimtype);
1654 }
1655
1656 if ( CDI_Debug )
1657 {
1658 Message("var %d %s", ncvarid, ncvar->name);
1659 for ( int i = 0; i < ndims; i++ )
1660 printf(" dim%d type=%d ", i, ncvar->dimtype[i]);
1661 printf("\n");
1662 }
1663 }
1664 }
1665
1666 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1667 {
1668 ncvar_t *ncvar = &ncvars[ncvarid];
1669 if (ncvar->varStatus == DataVar)
1670 {
1671 bool lxdim = false, lydim = false, lzdim = false/* , ltdim = false */;
1672 int lcdim = 0;
1673 int ndims = ncvar->ndims;
1674 for (int i = 0; i < ndims; i++)
1675 {
1676 int dimtype = ncvar->dimtype[i];
1677 lxdim |= (dimtype == X_AXIS);
1678 lydim |= (dimtype == Y_AXIS);
1679 lzdim |= (dimtype == Z_AXIS);
1680 if (ncvar->cvarids[i] != CDI_UNDEFID) lcdim++;
1681 // ltdim |= (dimtype == T_AXIS);
1682 }
1683
1684 int allcdims = lcdim;
1685
1686 if ( !lxdim && ncvar->xvarid != CDI_UNDEFID
1687 && ncvars[ncvar->xvarid].ndims == 0 ) lxdim = true;
1688
1689 if ( !lydim && ncvar->yvarid != CDI_UNDEFID
1690 && ncvars[ncvar->yvarid].ndims == 0 ) lydim = true;
1691
1692 if ( lxdim && (lydim || ncvar->gridtype == GRID_UNSTRUCTURED) )
1693 for ( int i = ndims-1; i >= 0; i-- )
1694 {
1695 if ( ncvar->dimtype[i] == -1 && !lzdim )
1696 {
1697 if ( lcdim )
1698 {
1699 int cdimvar = ncvar->cvarids[allcdims-lcdim];
1700 ncvar->zvarid = cdimvar;
1701 lcdim--;
1702 ncvars[cdimvar].zaxistype = ZAXIS_CHAR;
1703 }
1704 cdf_set_dim(ncvar, i, Z_AXIS);
1705 lzdim = true;
1706 int ncdimid = ncvar->dimids[i];
1707 if ( ncdims[ncdimid].dimtype == CDI_UNDEFID )
1708 ncdims[ncdimid].dimtype = Z_AXIS;
1709 }
1710 }
1711 }
1712 }
1713
1714 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1715 {
1716 ncvar_t *ncvar = &ncvars[ncvarid];
1717 int ndims = ncvar->ndims;
1718 for ( int i = 0; i < ndims; i++ )
1719 {
1720 if ( ncvar->dimtype[i] == CDI_UNDEFID )
1721 {
1722 int ncdimid = ncvar->dimids[i];
1723 if ( ncdims[ncdimid].dimtype == Z_AXIS )
1724 {
1725 ncvar->isZaxis = true;
1726 cdf_set_dim(ncvar, i, Z_AXIS);
1727 }
1728 }
1729 }
1730 }
1731
1732 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
1733 {
1734 ncvar_t *ncvar = &ncvars[ncvarid];
1735 if (ncvar->varStatus == DataVar)
1736 {
1737 bool lxdim = false, lydim = false, lzdim = false /*, ltdim = false */;
1738 int lcdim = 0;
1739 const int ndims = ncvar->ndims;
1740 for (int i = 0; i < ndims; i++)
1741 {
1742 int dimtype = ncvar->dimtype[i];
1743 lxdim |= (dimtype == X_AXIS);
1744 lydim |= (dimtype == Y_AXIS);
1745 lzdim |= (dimtype == Z_AXIS);
1746 if ( ncvar->cvarids[i] != CDI_UNDEFID ) lcdim++;
1747 // ltdim |= (dimtype == T_AXIS);
1748 }
1749
1750 const int allcdims = lcdim;
1751
1752 if (!lxdim && ncvar->xvarid != CDI_UNDEFID && ncvars[ncvar->xvarid].ndims == 0) lxdim = true;
1753 if (!lydim && ncvar->yvarid != CDI_UNDEFID && ncvars[ncvar->yvarid].ndims == 0) lydim = true;
1754
1755 // if ( ndims > 1 )
1756 for ( int i = ndims-1; i >= 0; i-- )
1757 {
1758 if ( ncvar->dimtype[i] == -1 )
1759 {
1760 if ( !lxdim )
1761 {
1762 if ( lcdim && ncvar->xvarid == CDI_UNDEFID )
1763 {
1764 int cdimvar = ncvar->cvarids[allcdims-lcdim];
1765 ncvar->xvarid = cdimvar;
1766 lcdim--;
1767 }
1768 cdf_set_dim(ncvar, i, X_AXIS);
1769 lxdim = true;
1770 }
1771 else if ( !lydim && ncvar->gridtype != GRID_UNSTRUCTURED )
1772 // else if ( !lydim && ! (ncvars[ncvar->xvarid].dimids[0] == ncvars[ncvar->yvarid].dimids[0] &&
1773 // ncvars[ncvar->xvarid].ndims == 1 && ncvars[ncvar->yvarid].ndims == 1))
1774 {
1775 if ( lcdim && ncvar->yvarid == CDI_UNDEFID )
1776 {
1777 int cdimvar = ncvar->cvarids[allcdims-lcdim];
1778 ncvar->yvarid = cdimvar;
1779 lcdim--;
1780 }
1781 cdf_set_dim(ncvar, i, Y_AXIS);
1782 lydim = true;
1783 }
1784 else if ( !lzdim )
1785 {
1786 if ( lcdim > 0 )
1787 {
1788 int cdimvar = ncvar->cvarids[allcdims-lcdim];
1789 ncvar->zvarid = cdimvar;
1790 lcdim--;
1791 ncvars[cdimvar].zaxistype = ZAXIS_CHAR;
1792 }
1793 cdf_set_dim(ncvar, i, Z_AXIS);
1794 lzdim = true;
1795 }
1796 }
1797 }
1798 if ( lcdim > 0 )
1799 Warning("Could not assign all character coordinates to data variable!");
1800 }
1801 }
1802 }
1803
1804 // verify coordinates vars - first scan (dimname == varname)
1805 static
verify_coordinates_vars_1(int ncid,int ndims,ncdim_t * ncdims,ncvar_t * ncvars,int timedimid,bool * isHybridCF)1806 void verify_coordinates_vars_1(int ncid, int ndims, ncdim_t *ncdims, ncvar_t *ncvars, int timedimid, bool *isHybridCF)
1807 {
1808 for ( int ncdimid = 0; ncdimid < ndims; ncdimid++ )
1809 {
1810 int ncvarid = ncdims[ncdimid].ncvarid;
1811 if ( ncvarid != -1 )
1812 {
1813 ncvar_t *ncvar = &ncvars[ncvarid];
1814 if ( ncvar->dimids[0] == timedimid )
1815 {
1816 ncvar->isTaxis = true;
1817 ncdims[ncdimid].dimtype = T_AXIS;
1818 continue;
1819 }
1820
1821 if ( isHybridSigmaPressureCoordinate(ncid, ncvarid, ncvars, ncdims) )
1822 {
1823 *isHybridCF = true;
1824 continue;
1825 }
1826
1827 if ( ncvar->units[0] != 0 )
1828 {
1829 if ( is_lon_axis(ncvar->units, ncvar->stdname) )
1830 {
1831 ncvar->isLon = true;
1832 cdf_set_var(ncvar, CoordVar);
1833 cdf_set_dim(ncvar, 0, X_AXIS);
1834 ncdims[ncdimid].dimtype = X_AXIS;
1835 }
1836 else if ( is_lat_axis(ncvar->units, ncvar->stdname) )
1837 {
1838 ncvar->isLat = true;
1839 cdf_set_var(ncvar, CoordVar);
1840 cdf_set_dim(ncvar, 0, Y_AXIS);
1841 ncdims[ncdimid].dimtype = Y_AXIS;
1842 }
1843 else if ( is_x_axis(ncvar->units, ncvar->stdname) )
1844 {
1845 ncvar->isXaxis = true;
1846 cdf_set_var(ncvar, CoordVar);
1847 cdf_set_dim(ncvar, 0, X_AXIS);
1848 ncdims[ncdimid].dimtype = X_AXIS;
1849 }
1850 else if ( is_y_axis(ncvar->units, ncvar->stdname) )
1851 {
1852 ncvar->isYaxis = true;
1853 cdf_set_var(ncvar, CoordVar);
1854 cdf_set_dim(ncvar, 0, Y_AXIS);
1855 ncdims[ncdimid].dimtype = Y_AXIS;
1856 }
1857 else if ( is_pressure_units(ncvar->units) )
1858 {
1859 ncvar->zaxistype = ZAXIS_PRESSURE;
1860 }
1861 else if (strIsEqual(ncvar->units, "level") || strIsEqual(ncvar->units, "1"))
1862 {
1863 if (strIsEqual(ncvar->longname, "hybrid level at layer midpoints"))
1864 ncvar->zaxistype = ZAXIS_HYBRID;
1865 else if (strStartsWith(ncvar->longname, "hybrid level at midpoints"))
1866 ncvar->zaxistype = ZAXIS_HYBRID;
1867 else if (strIsEqual(ncvar->longname, "hybrid level at layer interfaces"))
1868 ncvar->zaxistype = ZAXIS_HYBRID_HALF;
1869 else if (strStartsWith(ncvar->longname, "hybrid level at interfaces"))
1870 ncvar->zaxistype = ZAXIS_HYBRID_HALF;
1871 else if (strIsEqual(ncvar->units, "level"))
1872 ncvar->zaxistype = ZAXIS_GENERIC;
1873 }
1874 else if ( is_DBL_axis(ncvar->longname) )
1875 {
1876 ncvar->zaxistype = ZAXIS_DEPTH_BELOW_LAND;
1877 }
1878 else if ( is_height_units(ncvar->units) )
1879 {
1880 if ( is_depth_axis(ncvar->stdname, ncvar->longname) )
1881 ncvar->zaxistype = ZAXIS_DEPTH_BELOW_SEA;
1882 else if ( is_height_axis(ncvar->stdname, ncvar->longname) )
1883 ncvar->zaxistype = ZAXIS_HEIGHT;
1884 else if ( is_altitude_axis(ncvar->stdname, ncvar->longname) )
1885 ncvar->zaxistype = ZAXIS_ALTITUDE;
1886 }
1887 }
1888 else
1889 {
1890 if ( (strIsEqual(ncvar->longname, "generalized_height") ||
1891 strIsEqual(ncvar->longname, "generalized height")) &&
1892 strIsEqual(ncvar->stdname, "height"))
1893 ncvar->zaxistype = ZAXIS_REFERENCE;
1894 else if (strIsEqual(ncvar->stdname, "air_pressure"))
1895 ncvar->zaxistype = ZAXIS_PRESSURE;
1896 }
1897
1898 if ( !ncvar->isLon && (ncvar->longname[0] != 0) &&
1899 !ncvar->isLat && (ncvar->longname[1] != 0) )
1900 {
1901 if (strStartsWith(ncvar->longname+1, "ongitude"))
1902 {
1903 ncvar->isLon = true;
1904 cdf_set_var(ncvar, CoordVar);
1905 cdf_set_dim(ncvar, 0, X_AXIS);
1906 ncdims[ncdimid].dimtype = X_AXIS;
1907 continue;
1908 }
1909 else if (strStartsWith(ncvar->longname+1, "atitude"))
1910 {
1911 ncvar->isLat = true;
1912 cdf_set_var(ncvar, CoordVar);
1913 cdf_set_dim(ncvar, 0, Y_AXIS);
1914 ncdims[ncdimid].dimtype = Y_AXIS;
1915 continue;
1916 }
1917 }
1918
1919 if ( ncvar->zaxistype != CDI_UNDEFID )
1920 {
1921 ncvar->isZaxis = true;
1922 cdf_set_var(ncvar, CoordVar);
1923 cdf_set_dim(ncvar, 0, Z_AXIS);
1924 ncdims[ncdimid].dimtype = Z_AXIS;
1925 }
1926 }
1927 }
1928 }
1929
1930 // verify coordinates vars - second scan (all other variables)
1931 static
verify_coordinates_vars_2(stream_t * streamptr,int nvars,ncvar_t * ncvars)1932 void verify_coordinates_vars_2(stream_t *streamptr, int nvars, ncvar_t *ncvars)
1933 {
1934 for ( int ncvarid = 0; ncvarid < nvars; ncvarid++ )
1935 {
1936 ncvar_t *ncvar = &ncvars[ncvarid];
1937 if (ncvar->varStatus == CoordVar)
1938 {
1939 if ( ncvar->units[0] != 0 )
1940 {
1941 if ( is_lon_axis(ncvar->units, ncvar->stdname) )
1942 {
1943 ncvar->isLon = true;
1944 continue;
1945 }
1946 else if ( is_lat_axis(ncvar->units, ncvar->stdname) )
1947 {
1948 ncvar->isLat = true;
1949 continue;
1950 }
1951 else if ( is_x_axis(ncvar->units, ncvar->stdname) )
1952 {
1953 ncvar->isXaxis = true;
1954 continue;
1955 }
1956 else if ( is_y_axis(ncvar->units, ncvar->stdname) )
1957 {
1958 ncvar->isYaxis = true;
1959 continue;
1960 }
1961 else if ( ncvar->zaxistype == CDI_UNDEFID &&
1962 (strIsEqual(ncvar->units, "level") || strIsEqual(ncvar->units, "1")) )
1963 {
1964 if (strIsEqual(ncvar->longname, "hybrid level at layer midpoints"))
1965 ncvar->zaxistype = ZAXIS_HYBRID;
1966 else if (strStartsWith(ncvar->longname, "hybrid level at midpoints"))
1967 ncvar->zaxistype = ZAXIS_HYBRID;
1968 else if (strIsEqual(ncvar->longname, "hybrid level at layer interfaces"))
1969 ncvar->zaxistype = ZAXIS_HYBRID_HALF;
1970 else if (strStartsWith(ncvar->longname, "hybrid level at interfaces"))
1971 ncvar->zaxistype = ZAXIS_HYBRID_HALF;
1972 else if (strIsEqual(ncvar->units, "level"))
1973 ncvar->zaxistype = ZAXIS_GENERIC;
1974 continue;
1975 }
1976 else if ( ncvar->zaxistype == CDI_UNDEFID && is_pressure_units(ncvar->units) )
1977 {
1978 ncvar->zaxistype = ZAXIS_PRESSURE;
1979 continue;
1980 }
1981 else if ( is_DBL_axis(ncvar->longname) )
1982 {
1983 ncvar->zaxistype = ZAXIS_DEPTH_BELOW_LAND;
1984 continue;
1985 }
1986 else if ( is_height_units(ncvar->units) )
1987 {
1988 if ( is_depth_axis(ncvar->stdname, ncvar->longname) )
1989 ncvar->zaxistype = ZAXIS_DEPTH_BELOW_SEA;
1990 else if ( is_height_axis(ncvar->stdname, ncvar->longname) )
1991 ncvar->zaxistype = ZAXIS_HEIGHT;
1992 continue;
1993 }
1994 }
1995 else if (strIsEqual(ncvar->stdname, "region") ||
1996 strIsEqual(ncvar->stdname, "area_type") ||
1997 cdfInqDatatype(streamptr, ncvar->xtype, ncvar->lunsigned) == CDI_DATATYPE_UINT8 )
1998 {
1999 ncvar->isCharAxis = true;
2000 }
2001 else if (strIsEqual(ncvar->stdname, "air_pressure"))
2002 ncvar->zaxistype = ZAXIS_PRESSURE;
2003
2004 // not needed anymore for rotated grids
2005 if ( !ncvar->isLon && (ncvar->longname[0] != 0) &&
2006 !ncvar->isLat && (ncvar->longname[1] != 0) )
2007 {
2008 if (strStartsWith(ncvar->longname+1, "ongitude"))
2009 {
2010 ncvar->isLon = true;
2011 continue;
2012 }
2013 else if (strStartsWith(ncvar->longname+1, "atitude"))
2014 {
2015 ncvar->isLat = true;
2016 continue;
2017 }
2018 }
2019 }
2020 }
2021 }
2022
2023 static
grid_set_chunktype(grid_t * grid,ncvar_t * ncvar)2024 void grid_set_chunktype(grid_t *grid, ncvar_t *ncvar)
2025 {
2026 if (ncvar->chunked)
2027 {
2028 int ndims = ncvar->ndims;
2029
2030 if ( grid->type == GRID_UNSTRUCTURED )
2031 {
2032 ncvar->chunktype = ncvar->chunks[ndims-1] == grid->size
2033 ? CDI_CHUNK_GRID : CDI_CHUNK_AUTO;
2034 }
2035 else
2036 {
2037 if ( grid->x.size > 1 && grid->y.size > 1 && ndims > 1 &&
2038 grid->x.size == ncvar->chunks[ndims-1] &&
2039 grid->y.size == ncvar->chunks[ndims-2] )
2040 ncvar->chunktype = CDI_CHUNK_GRID;
2041 else if ( grid->x.size > 1 && grid->x.size == ncvar->chunks[ndims-1] )
2042 ncvar->chunktype = CDI_CHUNK_LINES;
2043 else
2044 ncvar->chunktype = CDI_CHUNK_AUTO;
2045 }
2046 }
2047 }
2048
2049 // define all input grids
2050 static
cdf_load_vals(size_t size,int ndims,int varid,ncvar_t * ncvar,double ** gridvals,struct xyValGet * valsGet,int ntdims,size_t * start,size_t * count)2051 void cdf_load_vals(size_t size, int ndims, int varid, ncvar_t *ncvar, double **gridvals, struct xyValGet *valsGet,
2052 int ntdims, size_t *start, size_t *count)
2053 {
2054 if (CDI_Netcdf_Lazy_Grid_Load)
2055 {
2056 *valsGet = (struct xyValGet){
2057 .scalefactor = ncvar->scalefactor,
2058 .addoffset = ncvar->addoffset,
2059 .start = { start[0], start[1], start[2] },
2060 .count = { count[0], count[1], count[2] },
2061 .size = size,
2062 .datasetNCId = ncvar->ncid,
2063 .varNCId = varid,
2064 .ndims = (short)ndims,
2065 };
2066 *gridvals = cdfPendingLoad;
2067 }
2068 else
2069 {
2070 *gridvals = (double*) Malloc(size*sizeof(double));
2071 if (ntdims == 1)
2072 cdf_get_vara_double(ncvar->ncid, varid, start, count, *gridvals);
2073 else
2074 cdf_get_var_double(ncvar->ncid, varid, *gridvals);
2075 cdf_scale_add(size, *gridvals, ncvar->addoffset, ncvar->scalefactor);
2076 }
2077 }
2078
2079 #ifndef USE_MPI
2080 static
cdf_load_cvals(size_t size,int varid,ncvar_t * ncvar,char *** gridvals,size_t dimlength)2081 void cdf_load_cvals(size_t size, int varid, ncvar_t *ncvar, char ***gridvals, size_t dimlength)
2082 {
2083 size_t startc[] = {0, 0};
2084 size_t countc[] = {1, size/dimlength};
2085 *gridvals = (char **) Malloc(dimlength * sizeof(char *));
2086 for ( size_t i = 0; i < dimlength; i++ )
2087 {
2088 (*gridvals)[i] = (char*) Malloc((size/dimlength) * sizeof(char));
2089 cdf_get_vara_text(ncvar->ncid, varid, startc, countc, (*gridvals)[i]);
2090 startc[0] = i+1;
2091 }
2092 }
2093 #endif
2094
2095 static
cdf_load_bounds(size_t size,ncvar_t * ncvar,double ** gridbounds,struct cdfLazyGridIds * cellBoundsGet)2096 void cdf_load_bounds(size_t size, ncvar_t *ncvar, double **gridbounds, struct cdfLazyGridIds *cellBoundsGet)
2097 {
2098 if (CDI_Netcdf_Lazy_Grid_Load)
2099 {
2100 cellBoundsGet->datasetNCId = ncvar->ncid;
2101 cellBoundsGet->varNCId = ncvar->bounds;
2102 *gridbounds = cdfPendingLoad;
2103 }
2104 else
2105 {
2106 *gridbounds = (double*) Malloc(size*sizeof(double));
2107 cdf_get_var_double(ncvar->ncid, ncvar->bounds, *gridbounds);
2108 }
2109 }
2110
2111 static
cdf_load_bounds_cube_sphere(size_t bxsize,size_t bysize,size_t size,ncvar_t * ncvar,double ** gridbounds,struct cdfLazyGridIds * cellBoundsGet)2112 void cdf_load_bounds_cube_sphere(size_t bxsize, size_t bysize, size_t size, ncvar_t *ncvar, double **gridbounds, struct cdfLazyGridIds *cellBoundsGet)
2113 {
2114 if (CDI_Netcdf_Lazy_Grid_Load)
2115 {
2116 cellBoundsGet->datasetNCId = ncvar->ncid;
2117 cellBoundsGet->varNCId = ncvar->bounds;
2118 *gridbounds = cdfPendingLoad;
2119 }
2120 else
2121 {
2122 float *bounds = (float*) Malloc(6*bxsize*bysize*sizeof(float));
2123 cdf_get_var_float(ncvar->ncid, ncvar->bounds, bounds);
2124
2125 *gridbounds = (double*) Malloc(size*sizeof(double));
2126 double *pbounds = *gridbounds;
2127
2128 size_t m = 0;
2129 for (size_t k = 0; k < 6; ++k)
2130 for (size_t j = 0; j < (bysize - 1); ++j)
2131 for (size_t i = 0; i < (bxsize - 1); ++i)
2132 {
2133 const size_t offset = k * bysize * bxsize;
2134 pbounds[m + 0] = bounds[offset + (j + 1) * bxsize + i];
2135 pbounds[m + 1] = bounds[offset + j * bxsize + i];
2136 pbounds[m + 2] = bounds[offset + j * bxsize + (i + 1)];
2137 pbounds[m + 3] = bounds[offset + (j + 1) * bxsize + (i + 1)];
2138 m += 4;
2139 }
2140
2141 Free(bounds);
2142 }
2143 }
2144
2145 static
cdf_load_cellarea(size_t size,ncvar_t * ncvar,double ** gridarea,struct cdfLazyGridIds * cellAreaGet)2146 void cdf_load_cellarea(size_t size, ncvar_t *ncvar, double **gridarea, struct cdfLazyGridIds *cellAreaGet)
2147 {
2148 if ( CDI_Netcdf_Lazy_Grid_Load )
2149 {
2150 cellAreaGet->datasetNCId = ncvar->ncid;
2151 cellAreaGet->varNCId = ncvar->cellarea;
2152 *gridarea = cdfPendingLoad;
2153 }
2154 else
2155 {
2156 *gridarea = (double*) Malloc(size*sizeof(double));
2157 cdf_get_var_double(ncvar->ncid, ncvar->cellarea, *gridarea);
2158 }
2159 }
2160
2161 static
cdf_copy_grid_axis_attr(ncvar_t * ncvar,struct gridaxis_t * gridaxis)2162 void cdf_copy_grid_axis_attr(ncvar_t *ncvar, struct gridaxis_t *gridaxis)
2163 {
2164 cdiDefVarKeyBytes(&gridaxis->keys, CDI_KEY_NAME, (const unsigned char*)ncvar->name, (int)strlen(ncvar->name)+1);
2165 if (ncvar->longname[0])
2166 cdiDefVarKeyBytes(&gridaxis->keys, CDI_KEY_LONGNAME, (const unsigned char*)ncvar->longname, (int)strlen(ncvar->longname)+1);
2167 if (ncvar->units[0])
2168 cdiDefVarKeyBytes(&gridaxis->keys, CDI_KEY_UNITS, (const unsigned char*)ncvar->units, (int)strlen(ncvar->units)+1);
2169 #ifndef USE_MPI
2170 if ( gridaxis->cvals )
2171 if (ncvar->stdname[0])
2172 cdiDefVarKeyBytes(&gridaxis->keys, CDI_KEY_STDNAME, (const unsigned char*)ncvar->stdname, (int)strlen(ncvar->stdname)+1);
2173 #endif
2174 }
2175
2176 static
cdf_get_xydimid(int ndims,int * dimids,int * dimtype,int * xdimid,int * ydimid)2177 int cdf_get_xydimid(int ndims, int *dimids, int *dimtype, int *xdimid, int *ydimid)
2178 {
2179 int nxdims = 0, nydims = 0;
2180 int xdimids[2] = {-1, -1}, ydimids[2] = {-1, -1};
2181
2182 for ( int i = 0; i < ndims; i++ )
2183 {
2184 if ( dimtype[i] == X_AXIS && nxdims < 2 )
2185 {
2186 xdimids[nxdims] = dimids[i];
2187 nxdims++;
2188 }
2189 else if ( dimtype[i] == Y_AXIS && nydims < 2 )
2190 {
2191 ydimids[nydims] = dimids[i];
2192 nydims++;
2193 }
2194 }
2195
2196 if ( nxdims == 2 )
2197 {
2198 *xdimid = xdimids[1];
2199 *ydimid = xdimids[0];
2200 }
2201 else if ( nydims == 2 )
2202 {
2203 *xdimid = ydimids[1];
2204 *ydimid = ydimids[0];
2205 }
2206 else
2207 {
2208 *xdimid = xdimids[0];
2209 *ydimid = ydimids[0];
2210 }
2211
2212 return nydims;
2213 }
2214
2215 static
cdf_check_gridtype(int * gridtype,bool isLon,bool isLat,size_t xsize,size_t ysize,grid_t * grid)2216 void cdf_check_gridtype(int *gridtype, bool isLon, bool isLat, size_t xsize, size_t ysize, grid_t *grid)
2217 {
2218 if ( isLat && (isLon || xsize == 0) )
2219 {
2220 double yinc = 0;
2221 if ( isLon && ysize > 1 )
2222 {
2223 yinc = fabs(grid->y.vals[0] - grid->y.vals[1]);
2224 for ( size_t i = 2; i < ysize; i++ )
2225 if ( (fabs(grid->y.vals[i-1] - grid->y.vals[i]) - yinc) > (yinc/1000) )
2226 {
2227 yinc = 0;
2228 break;
2229 }
2230 }
2231 if ( ysize < 10000 && IS_EQUAL(yinc, 0.0) && isGaussianLatitudes(ysize, grid->y.vals) )
2232 {
2233 *gridtype = GRID_GAUSSIAN;
2234 grid->np = (int)(ysize/2);
2235 }
2236 else
2237 *gridtype = GRID_LONLAT;
2238 }
2239 else if ( isLon && !isLat && ysize == 0 )
2240 {
2241 *gridtype = GRID_LONLAT;
2242 }
2243 else
2244 *gridtype = GRID_GENERIC;
2245 }
2246
2247 static
cdf_read_xcoord(stream_t * streamptr,struct cdfLazyGrid * lazyGrid,ncdim_t * ncdims,ncvar_t * ncvar,int xvarid,ncvar_t * axisvar,size_t * xsize,size_t ysize,int ntdims,size_t * start,size_t * count,bool * isLon)2248 bool cdf_read_xcoord(stream_t *streamptr, struct cdfLazyGrid *lazyGrid, ncdim_t *ncdims, ncvar_t *ncvar, int xvarid, ncvar_t *axisvar,
2249 size_t *xsize, size_t ysize, int ntdims, size_t *start, size_t *count, bool *isLon)
2250 {
2251 grid_t *grid = &lazyGrid->base;
2252 bool skipvar = true;
2253 *isLon = axisvar->isLon;
2254 int ndims = axisvar->ndims;
2255 size_t size = 0;
2256
2257 if (xtypeIsText(axisvar->xtype))
2258 {
2259 ncvar->varStatus = CoordVar;
2260 Warning("Unsupported coordinate type (char/string), skipped variable %s!", ncvar->name);
2261 return true;
2262 }
2263
2264 const int datatype = cdfInqDatatype(streamptr, axisvar->xtype, axisvar->lunsigned);
2265
2266 if ((ndims - ntdims) == 2)
2267 {
2268 // Check size of 2 dimensional coordinates variables
2269 int dimid = axisvar->dimids[ndims-2];
2270 const size_t dimsize1 = ncdims[dimid].len;
2271 dimid = axisvar->dimids[ndims-1];
2272 const size_t dimsize2 = ncdims[dimid].len;
2273
2274 if (datatype == CDI_DATATYPE_UINT8)
2275 {
2276 ncvar->gridtype = GRID_CHARXY;
2277 size = dimsize1 * dimsize2;
2278 skipvar = (dimsize1 != *xsize);
2279 }
2280 else
2281 {
2282 ncvar->gridtype = GRID_CURVILINEAR;
2283 size = (*xsize) * ysize;
2284 skipvar = (dimsize1 * dimsize2 != size);
2285 }
2286 }
2287 else if ((ndims - ntdims) == 1)
2288 {
2289 size = *xsize;
2290 // Check size of 1 dimensional coordinates variables
2291 const int dimid = axisvar->dimids[ndims-1];
2292 const size_t dimsize = ncdims[dimid].len;
2293 skipvar = (dimsize != size);
2294 }
2295 else if (ndims == 0 && *xsize == 0)
2296 {
2297 size = *xsize = 1;
2298 skipvar = false;
2299 }
2300 else if (ncvar->isCubeSphere)
2301 {
2302 size = *xsize;
2303 skipvar = false;
2304 }
2305
2306 if (skipvar)
2307 {
2308 Warning("Unsupported array structure, skipped variable %s!", ncvar->name);
2309 ncvar->varStatus = UndefVar;
2310 return true;
2311 }
2312
2313 if (datatype != -1) grid->datatype = datatype;
2314
2315 if (datatype == CDI_DATATYPE_UINT8 && !CDI_Netcdf_Lazy_Grid_Load)
2316 {
2317 #ifndef USE_MPI
2318 cdf_load_cvals(size, xvarid, axisvar, &grid->x.cvals, *xsize);
2319 grid->x.clength = size / (*xsize) ;
2320 #endif
2321 }
2322 else
2323 cdf_load_vals(size, ndims, xvarid, axisvar, &grid->x.vals, &lazyGrid->xValsGet, ntdims, start, count);
2324
2325 cdf_copy_grid_axis_attr(axisvar, &grid->x);
2326
2327 return false;
2328 }
2329
2330 static
cdf_read_ycoord(stream_t * streamptr,struct cdfLazyGrid * lazyGrid,ncdim_t * ncdims,ncvar_t * ncvar,int yvarid,ncvar_t * axisvar,size_t xsize,size_t * ysize,int ntdims,size_t * start,size_t * count,bool * isLat)2331 bool cdf_read_ycoord(stream_t *streamptr, struct cdfLazyGrid *lazyGrid, ncdim_t *ncdims, ncvar_t *ncvar, int yvarid, ncvar_t *axisvar,
2332 size_t xsize, size_t *ysize, int ntdims, size_t *start, size_t *count, bool *isLat)
2333 {
2334 grid_t *grid = &lazyGrid->base;
2335 bool skipvar = true;
2336 *isLat = axisvar->isLat;
2337 int ndims = axisvar->ndims;
2338 size_t size = 0;
2339
2340 if (xtypeIsText(axisvar->xtype))
2341 {
2342 ncvar->varStatus = CoordVar;
2343 Warning("Unsupported coordinate type (char/string), skipped variable %s!", ncvar->name);
2344 return true;
2345 }
2346
2347 const int datatype = cdfInqDatatype(streamptr, axisvar->xtype, axisvar->lunsigned);
2348
2349 if ((ndims - ntdims) == 2)
2350 {
2351 // Check size of 2 dimensional coordinates variables
2352 int dimid = axisvar->dimids[ndims-2];
2353 const size_t dimsize1 = ncdims[dimid].len;
2354 dimid = axisvar->dimids[ndims-1];
2355 const size_t dimsize2 = ncdims[dimid].len;
2356
2357 if (datatype == CDI_DATATYPE_UINT8)
2358 {
2359 ncvar->gridtype = GRID_CHARXY;
2360 size = dimsize1 * dimsize2;
2361 skipvar = (dimsize1 != *ysize);
2362 }
2363 else
2364 {
2365 ncvar->gridtype = GRID_CURVILINEAR;
2366 size = xsize * (*ysize);
2367 skipvar = (dimsize1 * dimsize2 != size);
2368 }
2369 }
2370 else if ( (ndims - ntdims) == 1 )
2371 {
2372 size = (*ysize) ? *ysize : xsize;
2373 const int dimid = axisvar->dimids[ndims-1];
2374 const size_t dimsize = ncdims[dimid].len;
2375 skipvar = (dimsize != size);
2376 }
2377 else if ( ndims == 0 && *ysize == 0 )
2378 {
2379 size = *ysize = 1;
2380 skipvar = false;
2381 }
2382 else if (ncvar->isCubeSphere)
2383 {
2384 size = *ysize;
2385 skipvar = false;
2386 }
2387
2388 if (skipvar)
2389 {
2390 Warning("Unsupported array structure, skipped variable %s!", ncvar->name);
2391 ncvar->varStatus = UndefVar;
2392 return true;
2393 }
2394
2395 if (datatype != -1) grid->datatype = datatype;
2396
2397 if (datatype == CDI_DATATYPE_UINT8 && !CDI_Netcdf_Lazy_Grid_Load)
2398 {
2399 #ifndef USE_MPI
2400 cdf_load_cvals(size, yvarid, axisvar, &grid->y.cvals, *ysize);
2401 grid->y.clength = size / (*ysize) ;
2402 #endif
2403 }
2404 else
2405 cdf_load_vals(size, ndims, yvarid, axisvar, &grid->y.vals, &lazyGrid->yValsGet, ntdims, start, count);
2406
2407 cdf_copy_grid_axis_attr(axisvar, &grid->y);
2408
2409 return false;
2410 }
2411
2412 static
cdf_read_coordinates(stream_t * streamptr,struct cdfLazyGrid * lazyGrid,ncvar_t * ncvar,ncvar_t * ncvars,ncdim_t * ncdims,int timedimid,int xvarid,int yvarid,size_t xsize,size_t ysize,int * vdimid)2413 bool cdf_read_coordinates(stream_t *streamptr, struct cdfLazyGrid *lazyGrid, ncvar_t *ncvar, ncvar_t *ncvars, ncdim_t *ncdims,
2414 int timedimid, int xvarid, int yvarid, size_t xsize, size_t ysize, int *vdimid)
2415 {
2416 grid_t *grid = &lazyGrid->base;
2417 size_t size = 0;
2418
2419 grid->datatype = CDI_DATATYPE_FLT64;
2420
2421 if (ncvar->gridtype == GRID_TRAJECTORY)
2422 {
2423 if (ncvar->xvarid == CDI_UNDEFID) Error("Longitude coordinates undefined for %s!", ncvar->name);
2424 if (ncvar->yvarid == CDI_UNDEFID) Error("Latitude coordinates undefined for %s!", ncvar->name);
2425 }
2426 else
2427 {
2428 size_t start[3] = {0, 0, 0}, count[3] = {1, 1, 1};
2429 int ntdims = 0;
2430
2431 if (xvarid != CDI_UNDEFID && yvarid != CDI_UNDEFID)
2432 {
2433 const int ndims = ncvars[xvarid].ndims;
2434 if ( ndims != ncvars[yvarid].ndims && !ncvars[xvarid].isCharAxis && !ncvars[yvarid].isCharAxis )
2435 {
2436 Warning("Inconsistent grid structure for variable %s!", ncvar->name);
2437 ncvar->xvarid = xvarid = CDI_UNDEFID;
2438 ncvar->yvarid = yvarid = CDI_UNDEFID;
2439 }
2440 if (ndims > 1)
2441 {
2442 if (ndims <= 3)
2443 {
2444 if (ncvars[xvarid].dimids[0] == timedimid && ncvars[yvarid].dimids[0] == timedimid)
2445 {
2446 static bool ltwarn = true;
2447 size_t ntsteps = 0;
2448 cdf_inq_dimlen(ncvar->ncid, ncdims[timedimid].dimid, &ntsteps);
2449 if (ltwarn && ntsteps > 1) Warning("Time varying grids unsupported, using grid at time step 1!");
2450 ltwarn = false;
2451 ntdims = 1;
2452 count[1] = ysize; count[2] = xsize;
2453 }
2454 }
2455 else
2456 {
2457 Warning("Unsupported grid structure for variable %s (grid dims > 2)!", ncvar->name);
2458 ncvar->xvarid = xvarid = CDI_UNDEFID;
2459 ncvar->yvarid = yvarid = CDI_UNDEFID;
2460 }
2461 }
2462 }
2463
2464 if (xvarid != CDI_UNDEFID)
2465 {
2466 if (!ncvar->isCubeSphere && (ncvars[xvarid].ndims - ntdims) > 2)
2467 {
2468 Warning("Coordinates variable %s has too many dimensions (%d), skipped!", ncvars[xvarid].name, ncvars[xvarid].ndims);
2469 //ncvar->xvarid = CDI_UNDEFID;
2470 xvarid = CDI_UNDEFID;
2471 }
2472 }
2473
2474 if (yvarid != CDI_UNDEFID)
2475 {
2476 if (!ncvar->isCubeSphere && (ncvars[yvarid].ndims - ntdims) > 2)
2477 {
2478 Warning("Coordinates variable %s has too many dimensions (%d), skipped!", ncvars[yvarid].name, ncvars[yvarid].ndims);
2479 //ncvar->yvarid = CDI_UNDEFID;
2480 yvarid = CDI_UNDEFID;
2481 }
2482 }
2483
2484 bool isLon = false, isLat = false;
2485
2486 if (xvarid != CDI_UNDEFID)
2487 if (cdf_read_xcoord(streamptr, lazyGrid, ncdims, ncvar, xvarid, &ncvars[xvarid],
2488 &xsize, ysize, ntdims, start, count, &isLon))
2489 return true;
2490
2491 if (yvarid != CDI_UNDEFID)
2492 if (cdf_read_ycoord(streamptr, lazyGrid, ncdims, ncvar, yvarid, &ncvars[yvarid],
2493 xsize, &ysize, ntdims, start, count, &isLat))
2494 return true;
2495
2496 if (ncvar->gridtype == GRID_UNSTRUCTURED) size = xsize;
2497 else if (ncvar->gridtype == GRID_GAUSSIAN_REDUCED) size = xsize;
2498 else if (ysize == 0) size = xsize;
2499 else if (xsize == 0) size = ysize;
2500 else size = xsize*ysize;
2501
2502 if (ncvar->gridtype == CDI_UNDEFID || ncvar->gridtype == GRID_GENERIC)
2503 cdf_check_gridtype(&ncvar->gridtype, isLon, isLat, xsize, ysize, grid);
2504 }
2505
2506 int gridtype = grid->type;
2507 if (gridtype != GRID_PROJECTION) gridtype = ncvar->gridtype;
2508 else if (gridtype == GRID_PROJECTION && ncvar->gridtype == GRID_LONLAT)
2509 {
2510 const int gmapvarid = ncvar->gmapid;
2511 if (gmapvarid != CDI_UNDEFID && cdfCheckAttText(ncvar->ncid, gmapvarid, "grid_mapping_name"))
2512 {
2513 char attstring[CDI_MAX_NAME];
2514 cdfGetAttText(ncvar->ncid, gmapvarid, "grid_mapping_name", CDI_MAX_NAME, attstring);
2515 if (strIsEqual(attstring, "latitude_longitude")) gridtype = ncvar->gridtype;
2516 }
2517 }
2518
2519 switch (gridtype)
2520 {
2521 case GRID_GENERIC:
2522 case GRID_LONLAT:
2523 case GRID_GAUSSIAN:
2524 case GRID_UNSTRUCTURED:
2525 case GRID_CURVILINEAR:
2526 case GRID_PROJECTION:
2527 {
2528 grid->size = size;
2529 grid->x.size = xsize;
2530 grid->y.size = ysize;
2531 if (xvarid != CDI_UNDEFID && CDI_Read_Cell_Corners)
2532 {
2533 grid->x.flag = 1;
2534 const int bvarid = ncvars[xvarid].bounds;
2535 if ( bvarid != CDI_UNDEFID )
2536 {
2537 const int nbdims = ncvars[bvarid].ndims;
2538 if (nbdims == 2 || nbdims == 3)
2539 {
2540 if (ncvar->isCubeSphere)
2541 {
2542 grid->nvertex = 4;
2543 size_t bxsize = ncdims[ncvars[bvarid].dimids[2]].len;
2544 size_t bysize = ncdims[ncvars[bvarid].dimids[1]].len;
2545 cdf_load_bounds_cube_sphere(bxsize, bysize, size * (size_t)grid->nvertex, &ncvars[xvarid], &grid->x.bounds, &lazyGrid->xBoundsGet);
2546 }
2547 else
2548 {
2549 *vdimid = ncvars[bvarid].dimids[nbdims-1];
2550 grid->nvertex = (int)ncdims[*vdimid].len;
2551 cdf_load_bounds(size * (size_t)grid->nvertex, &ncvars[xvarid], &grid->x.bounds, &lazyGrid->xBoundsGet);
2552 }
2553 }
2554 }
2555 }
2556 if (yvarid != CDI_UNDEFID && CDI_Read_Cell_Corners)
2557 {
2558 grid->y.flag = 1;
2559 const int bvarid = ncvars[yvarid].bounds;
2560 if (bvarid != CDI_UNDEFID)
2561 {
2562 const int nbdims = ncvars[bvarid].ndims;
2563 if (nbdims == 2 || nbdims == 3)
2564 {
2565 if (ncvar->isCubeSphere)
2566 {
2567 grid->nvertex = 4;
2568 size_t bxsize = ncdims[ncvars[bvarid].dimids[2]].len;
2569 size_t bysize = ncdims[ncvars[bvarid].dimids[1]].len;
2570 cdf_load_bounds_cube_sphere(bxsize, bysize, size * (size_t)grid->nvertex, &ncvars[yvarid], &grid->y.bounds, &lazyGrid->yBoundsGet);
2571 }
2572 else
2573 {
2574 if (*vdimid == CDI_UNDEFID)
2575 {
2576 *vdimid = ncvars[bvarid].dimids[nbdims-1];
2577 grid->nvertex = (int)ncdims[*vdimid].len;
2578 }
2579 cdf_load_bounds(size * (size_t)grid->nvertex, &ncvars[yvarid], &grid->y.bounds, &lazyGrid->yBoundsGet);
2580 }
2581 }
2582 }
2583 }
2584
2585 if (ncvar->cellarea != CDI_UNDEFID) cdf_load_cellarea(size, ncvar, &grid->area, &lazyGrid->cellAreaGet);
2586
2587 if (gridtype == GRID_GAUSSIAN)
2588 {
2589 if (ncvar->numLPE > 0) grid->np = ncvar->numLPE;
2590 }
2591
2592 break;
2593 }
2594 case GRID_GAUSSIAN_REDUCED:
2595 {
2596 if (ncvar->numLPE > 0 && ncvar->rpvarid != CDI_UNDEFID)
2597 {
2598 if (ncvars[ncvar->rpvarid].ndims == 1)
2599 {
2600 grid->size = size;
2601 const int dimid = ncvars[ncvar->rpvarid].dimids[0];
2602 const size_t len = ncdims[dimid].len;
2603 grid->y.size = len;
2604 grid->reducedPointsSize = len;
2605 grid->reducedPoints = (int*) Malloc(len*sizeof(int));
2606 cdf_get_var_int(ncvar->ncid, ncvar->rpvarid, grid->reducedPoints);
2607 grid->np = ncvar->numLPE;
2608
2609 const int bvarid = (yvarid == CDI_UNDEFID) ? CDI_UNDEFID : ncvars[yvarid].bounds;
2610 if (bvarid != CDI_UNDEFID)
2611 {
2612 const int nbdims = ncvars[bvarid].ndims;
2613 if (nbdims == 2 || nbdims == 3)
2614 {
2615 if (*vdimid == CDI_UNDEFID)
2616 {
2617 *vdimid = ncvars[bvarid].dimids[nbdims-1];
2618 grid->nvertex = (int)ncdims[*vdimid].len;
2619 }
2620 cdf_load_bounds(size*(size_t)grid->nvertex, &ncvars[yvarid], &grid->y.bounds, &lazyGrid->yBoundsGet);
2621 }
2622 }
2623 }
2624 }
2625 break;
2626 }
2627 case GRID_SPECTRAL:
2628 {
2629 grid->size = size;
2630 grid->lcomplex = 1;
2631 grid->trunc = ncvar->truncation;
2632 break;
2633 }
2634 case GRID_FOURIER:
2635 {
2636 grid->size = size;
2637 grid->trunc = ncvar->truncation;
2638 break;
2639 }
2640 case GRID_TRAJECTORY:
2641 {
2642 grid->size = 1;
2643 break;
2644 }
2645 case GRID_CHARXY:
2646 {
2647 grid->size = size;
2648 grid->x.size = xsize;
2649 grid->y.size = ysize;
2650 break;
2651 }
2652 }
2653
2654 // if ( grid->type != GRID_PROJECTION && grid->type != ncvar->gridtype )
2655 if ( grid->type != gridtype )
2656 {
2657 // int gridtype = ncvar->gridtype;
2658 grid->type = gridtype;
2659 cdiGridTypeInit(grid, gridtype, grid->size);
2660 }
2661
2662 if ( grid->size == 0 )
2663 {
2664 const int ndims = ncvar->ndims;
2665 int *dimtype = ncvar->dimtype;
2666 if ( ndims == 0 ||
2667 (ndims == 1 && dimtype[0] == T_AXIS) ||
2668 (ndims == 1 && dimtype[0] == Z_AXIS) ||
2669 (ndims == 2 && dimtype[0] == T_AXIS && dimtype[1] == Z_AXIS) )
2670 {
2671 grid->type = GRID_GENERIC;
2672 grid->size = 1;
2673 grid->x.size = 0;
2674 grid->y.size = 0;
2675 }
2676 else
2677 {
2678 Warning("Unsupported grid, skipped variable %s!", ncvar->name);
2679 ncvar->varStatus = UndefVar;
2680 return true;
2681 }
2682 }
2683
2684 return false;
2685 }
2686
2687 static
cdf_set_unstructured_par(ncvar_t * ncvar,grid_t * grid,int * xdimid,int * ydimid,GridInfo * gridInfo)2688 bool cdf_set_unstructured_par(ncvar_t *ncvar, grid_t *grid, int *xdimid, int *ydimid, GridInfo *gridInfo)
2689 {
2690 int ndims = ncvar->ndims;
2691 int *dimtype = ncvar->dimtype;
2692
2693 int zdimid = CDI_UNDEFID;
2694 int xdimidx = CDI_UNDEFID, ydimidx = CDI_UNDEFID;
2695
2696 for (int i = 0; i < ndims; i++)
2697 {
2698 if (dimtype[i] == X_AXIS) xdimidx = i;
2699 else if (dimtype[i] == Y_AXIS) ydimidx = i;
2700 else if (dimtype[i] == Z_AXIS) zdimid = ncvar->dimids[i];
2701 }
2702
2703 if (*xdimid != CDI_UNDEFID && *ydimid != CDI_UNDEFID && zdimid == CDI_UNDEFID)
2704 {
2705 if (grid->x.size > grid->y.size && grid->y.size < 1000)
2706 {
2707 dimtype[ydimidx] = Z_AXIS;
2708 *ydimid = CDI_UNDEFID;
2709 grid->size = grid->x.size;
2710 grid->y.size = 0;
2711 }
2712 else if (grid->y.size > grid->x.size && grid->x.size < 1000)
2713 {
2714 dimtype[xdimidx] = Z_AXIS;
2715 *xdimid = *ydimid;
2716 *ydimid = CDI_UNDEFID;
2717 grid->size = grid->y.size;
2718 grid->x.size = grid->y.size;
2719 grid->y.size = 0;
2720 }
2721 }
2722
2723 if (grid->size != grid->x.size)
2724 {
2725 Warning("Unsupported array structure, skipped variable %s!", ncvar->name);
2726 ncvar->varStatus = UndefVar;
2727 return true;
2728 }
2729
2730 if (gridInfo->number_of_grid_used != CDI_UNDEFID) cdiDefVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDUSED, gridInfo->number_of_grid_used);
2731 if (ncvar->position > 0) cdiDefVarKeyInt(&grid->keys, CDI_KEY_NUMBEROFGRIDINREFERENCE, ncvar->position);
2732 if (!cdiUUIDIsNull(gridInfo->uuid)) cdiDefVarKeyBytes(&grid->keys, CDI_KEY_UUID, gridInfo->uuid, CDI_UUID_SIZE);
2733
2734 return false;
2735 }
2736
2737 static
cdf_read_mapping_atts(int ncid,int gmapvarid,int projID)2738 void cdf_read_mapping_atts(int ncid, int gmapvarid, int projID)
2739 {
2740 if (cdfCheckAttText(ncid, gmapvarid, "grid_mapping_name"))
2741 {
2742 char attstring[CDI_MAX_NAME];
2743 cdfGetAttText(ncid, gmapvarid, "grid_mapping_name", CDI_MAX_NAME, attstring);
2744 cdiDefKeyString(projID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, attstring);
2745 }
2746
2747 int nvatts;
2748 cdf_inq_varnatts(ncid, gmapvarid, &nvatts);
2749 for (int attnum = 0; attnum < nvatts; ++attnum)
2750 cdf_set_cdi_attr(ncid, gmapvarid, attnum, projID, CDI_GLOBAL);
2751 }
2752
2753 static
cdf_set_grid_to_similar_vars(ncvar_t * ncvar1,ncvar_t * ncvar2,int gridtype,int xdimid,int ydimid)2754 void cdf_set_grid_to_similar_vars(ncvar_t *ncvar1, ncvar_t *ncvar2, int gridtype, int xdimid, int ydimid)
2755 {
2756 if (ncvar2->varStatus == DataVar && ncvar2->gridID == CDI_UNDEFID)
2757 {
2758 int xdimid2 = CDI_UNDEFID, ydimid2 = CDI_UNDEFID, zdimid2 = CDI_UNDEFID;
2759 int xdimidx = CDI_UNDEFID, ydimidx = CDI_UNDEFID;
2760 int ndims2 = ncvar2->ndims;
2761
2762 int *dimtype2 = ncvar2->dimtype;
2763 int *dimids2 = ncvar2->dimids;
2764 for (int i = 0; i < ndims2; i++)
2765 {
2766 if (dimtype2[i] == X_AXIS) { xdimid2 = dimids2[i]; xdimidx = i; }
2767 else if (dimtype2[i] == Y_AXIS) { ydimid2 = dimids2[i]; ydimidx = i; }
2768 else if (dimtype2[i] == Z_AXIS) { zdimid2 = dimids2[i]; }
2769 }
2770
2771 if (!ncvar2->isCubeSphere && ncvar2->gridtype == CDI_UNDEFID && gridtype == GRID_UNSTRUCTURED)
2772 {
2773 if (xdimid == xdimid2 && ydimid2 != CDI_UNDEFID && zdimid2 == CDI_UNDEFID)
2774 {
2775 ncvar2->dimtype[ydimidx] = Z_AXIS;
2776 ydimid2 = CDI_UNDEFID;
2777 }
2778
2779 if (xdimid == ydimid2 && xdimid2 != CDI_UNDEFID && zdimid2 == CDI_UNDEFID)
2780 {
2781 ncvar2->dimtype[xdimidx] = Z_AXIS;
2782 xdimid2 = ydimid2;
2783 ydimid2 = CDI_UNDEFID;
2784 }
2785 }
2786
2787 if (xdimid == xdimid2 && (ydimid == ydimid2 || (xdimid == ydimid && ydimid2 == CDI_UNDEFID)))
2788 {
2789 bool is_same_grid = (ncvar1->xvarid == ncvar2->xvarid
2790 && ncvar1->yvarid == ncvar2->yvarid
2791 && ncvar1->position == ncvar2->position);
2792
2793 // if ( xvarid != -1 && ncvar2->xvarid != CDI_UNDEFID && xvarid != ncvar2->xvarid ) is_same_grid = false;
2794 // if ( yvarid != -1 && ncvar2->yvarid != CDI_UNDEFID && yvarid != ncvar2->yvarid ) is_same_grid = false;
2795
2796 if (is_same_grid)
2797 {
2798 if (CDI_Debug) Message("Same gridID %d %s", ncvar1->gridID, ncvar2->name);
2799 ncvar2->gridID = ncvar1->gridID;
2800 ncvar2->chunktype = ncvar1->chunktype;
2801 }
2802 }
2803 }
2804 }
2805
2806 static
cdf_define_all_grids(stream_t * streamptr,ncgrid_t * ncgrid,int vlistID,ncdim_t * ncdims,int nvars,ncvar_t * ncvars,GridInfo * gridInfo)2807 int cdf_define_all_grids(stream_t *streamptr, ncgrid_t *ncgrid, int vlistID, ncdim_t *ncdims, int nvars, ncvar_t *ncvars, GridInfo *gridInfo)
2808 {
2809 for (int ncvarid = 0; ncvarid < nvars; ++ncvarid)
2810 {
2811 ncvar_t *ncvar = &ncvars[ncvarid];
2812 if (ncvar->varStatus == DataVar && ncvar->gridID == CDI_UNDEFID)
2813 {
2814 int ndims = ncvar->ndims;
2815 int *dimtype = ncvar->dimtype;
2816 int vdimid = CDI_UNDEFID;
2817 struct addIfNewRes projAdded = { .Id = CDI_UNDEFID, .isNew = 0 },
2818 gridAdded = { .Id = CDI_UNDEFID, .isNew = 0 };
2819 int xdimid = CDI_UNDEFID, ydimid = CDI_UNDEFID;
2820 int nydims = cdf_get_xydimid(ndims, ncvar->dimids, dimtype, &xdimid, &ydimid);
2821
2822 int xaxisid = (xdimid != CDI_UNDEFID) ? ncdims[xdimid].ncvarid : CDI_UNDEFID;
2823 int yaxisid = (ydimid != CDI_UNDEFID) ? ncdims[ydimid].ncvarid : CDI_UNDEFID;
2824 int xvarid = (ncvar->xvarid != CDI_UNDEFID) ? ncvar->xvarid : xaxisid;
2825 int yvarid = (ncvar->yvarid != CDI_UNDEFID) ? ncvar->yvarid : yaxisid;
2826
2827 size_t xsize = (xdimid != CDI_UNDEFID) ? ncdims[xdimid].len : 0;
2828 size_t ysize = (ydimid != CDI_UNDEFID) ? ncdims[ydimid].len : 0;
2829
2830 if (ydimid == CDI_UNDEFID && yvarid != CDI_UNDEFID)
2831 {
2832 if (ncvars[yvarid].ndims == 1)
2833 {
2834 ydimid = ncvars[yvarid].dimids[0];
2835 ysize = ncdims[ydimid].len;
2836 }
2837 }
2838
2839 const int gmapvarid = ncvar->gmapid;
2840 bool lproj = (gmapvarid != CDI_UNDEFID);
2841
2842 if (!lproj && xaxisid != CDI_UNDEFID && xaxisid != xvarid && yaxisid != CDI_UNDEFID && yaxisid != yvarid) lproj = true;
2843
2844 if (ncvar->isCubeSphere && lproj && xvarid != CDI_UNDEFID && yvarid != CDI_UNDEFID && ncvars[xvarid].ndims == 3 && ncvars[yvarid].ndims == 3)
2845 {
2846 lproj = false;
2847 ncvar->gridtype = GRID_UNSTRUCTURED;
2848 xsize = xsize * ysize * 6;
2849 ysize = xsize;
2850 }
2851
2852 const bool lgrid = !(lproj && ncvar->xvarid == CDI_UNDEFID);
2853
2854 const bool lunstructured = (xdimid != CDI_UNDEFID && xdimid == ydimid && nydims == 0);
2855 if ((ncvar->gridtype == CDI_UNDEFID || ncvar->gridtype == GRID_GENERIC) && lunstructured)
2856 ncvar->gridtype = GRID_UNSTRUCTURED;
2857
2858 struct cdfLazyGrid *lazyGrid = NULL, *lazyProj = NULL;
2859
2860 {
2861 const int gridtype = (!lgrid) ? GRID_PROJECTION : ncvar->gridtype;
2862 if (CDI_Netcdf_Lazy_Grid_Load)
2863 {
2864 cdfLazyGridRenew(&lazyGrid, gridtype);
2865 if (lgrid && lproj) cdfLazyGridRenew(&lazyProj, GRID_PROJECTION);
2866 }
2867 else
2868 {
2869 cdfBaseGridRenew(&lazyGrid, gridtype);
2870 if (lgrid && lproj) cdfBaseGridRenew(&lazyProj, GRID_PROJECTION);
2871 }
2872 }
2873 grid_t *grid = &lazyGrid->base;
2874 grid_t *proj = (lgrid && lproj) ? &lazyProj->base : NULL;
2875
2876 xaxisid = (xdimid != CDI_UNDEFID) ? ncdims[xdimid].ncvarid : CDI_UNDEFID;
2877 yaxisid = (ydimid != CDI_UNDEFID) ? ncdims[ydimid].ncvarid : CDI_UNDEFID;
2878
2879 if (cdf_read_coordinates(streamptr, lazyGrid, ncvar, ncvars, ncdims, gridInfo->timedimid,
2880 xvarid, yvarid, xsize, ysize, &vdimid))
2881 continue;
2882
2883 if (gridInfo->number_of_grid_used != CDI_UNDEFID &&
2884 (grid->type == CDI_UNDEFID || grid->type == GRID_GENERIC) &&
2885 xdimid != CDI_UNDEFID && xsize > 999)
2886 grid->type = GRID_UNSTRUCTURED;
2887
2888 if (!ncvar->isCubeSphere && grid->type == GRID_UNSTRUCTURED)
2889 if (cdf_set_unstructured_par(ncvar, grid, &xdimid, &ydimid, gridInfo))
2890 continue;
2891
2892 if (lgrid && lproj)
2893 {
2894 int dimid;
2895 cdf_read_coordinates(streamptr, lazyProj, ncvar, ncvars, ncdims, gridInfo->timedimid,
2896 xaxisid, yaxisid, xsize, ysize, &dimid);
2897 }
2898
2899 if (CDI_Debug)
2900 {
2901 Message("grid: type = %d, size = %zu, nx = %zu, ny = %zu",
2902 grid->type, grid->size, grid->x.size, grid->y.size);
2903 if (proj)
2904 Message("proj: type = %d, size = %zu, nx = %zu, ny = %zu",
2905 proj->type, proj->size, proj->x.size, proj->y.size);
2906 }
2907
2908
2909 if (lgrid && lproj)
2910 {
2911 projAdded = cdiVlistAddGridIfNew(vlistID, proj, 2);
2912 grid->proj = projAdded.Id;
2913 }
2914
2915 gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 1);
2916 ncvar->gridID = gridAdded.Id;
2917
2918 const int gridID = ncvar->gridID;
2919
2920 if (lproj && gmapvarid != CDI_UNDEFID)
2921 {
2922 const bool gridIsNew = lgrid ? projAdded.isNew : gridAdded.isNew;
2923 if (gridIsNew)
2924 {
2925 const int projID = lgrid ? grid->proj : gridID;
2926 const int ncid = ncvars[gmapvarid].ncid;
2927 const int gmapvartype = ncvars[gmapvarid].xtype;
2928 cdiDefKeyInt(projID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARTYPE, gmapvartype);
2929 const char *gmapvarname = ncvars[gmapvarid].name;
2930 cdf_read_mapping_atts(ncid, gmapvarid, projID);
2931 cdiDefKeyString(projID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARNAME, gmapvarname);
2932 gridVerifyProj(projID);
2933 }
2934 }
2935
2936 if (grid->type == GRID_UNSTRUCTURED && gridInfo->gridfile[0] != 0)
2937 gridDefReference(gridID, gridInfo->gridfile);
2938
2939 if (ncvar->chunked) grid_set_chunktype(grid, ncvar);
2940
2941 const int gridindex = vlistGridIndex(vlistID, gridID);
2942 ncgrid[gridindex].gridID = gridID;
2943 if (grid->type == GRID_TRAJECTORY)
2944 {
2945 ncgrid[gridindex].ncIDs[CDF_VARID_X] = xvarid;
2946 ncgrid[gridindex].ncIDs[CDF_VARID_Y] = yvarid;
2947 }
2948 else
2949 {
2950 if (xdimid != CDI_UNDEFID) ncgrid[gridindex].ncIDs[CDF_DIMID_X] = ncdims[xdimid].dimid;
2951 if (ydimid != CDI_UNDEFID) ncgrid[gridindex].ncIDs[CDF_DIMID_Y] = ncdims[ydimid].dimid;
2952 if (ncvar->isCubeSphere) ncgrid[gridindex].ncIDs[CDF_DIMID_E] = ncdims[ncvar->dimids[ndims-3]].dimid;
2953 }
2954
2955 if (xdimid == CDI_UNDEFID && ydimid == CDI_UNDEFID && grid->size == 1) gridDefHasDims(gridID, CoordVar);
2956
2957 if (xdimid != CDI_UNDEFID) cdiDefKeyString(gridID, CDI_XAXIS, CDI_KEY_DIMNAME, ncdims[xdimid].name);
2958 if (ydimid != CDI_UNDEFID) cdiDefKeyString(gridID, CDI_YAXIS, CDI_KEY_DIMNAME, ncdims[ydimid].name);
2959 if (vdimid != CDI_UNDEFID) cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_VDIMNAME, ncdims[vdimid].name);
2960
2961 if (CDI_Debug) Message("gridID %d %d %s", gridID, ncvarid, ncvar->name);
2962
2963 for (int ncvarid2 = ncvarid+1; ncvarid2 < nvars; ncvarid2++)
2964 cdf_set_grid_to_similar_vars(ncvar, &ncvars[ncvarid2], grid->type, xdimid, ydimid);
2965
2966 if (gridAdded.isNew) lazyGrid = NULL;
2967 if (projAdded.isNew) lazyProj = NULL;
2968
2969 if (lazyGrid)
2970 {
2971 if (CDI_Netcdf_Lazy_Grid_Load) cdfLazyGridDestroy(lazyGrid);
2972 if (grid) { grid_free(grid); Free(grid); }
2973 }
2974
2975 if (lazyProj)
2976 {
2977 if (CDI_Netcdf_Lazy_Grid_Load) cdfLazyGridDestroy(lazyProj);
2978 if (proj) { grid_free(proj); Free(proj); }
2979 }
2980 }
2981 }
2982
2983 return 0;
2984 }
2985
2986 // define all input zaxes
2987 static
cdf_define_all_zaxes(stream_t * streamptr,int vlistID,ncdim_t * ncdims,int nvars,ncvar_t * ncvars,size_t vctsize_echam,double * vct_echam,unsigned char * uuidOfVGrid)2988 int cdf_define_all_zaxes(stream_t *streamptr, int vlistID, ncdim_t *ncdims, int nvars, ncvar_t *ncvars,
2989 size_t vctsize_echam, double *vct_echam, unsigned char *uuidOfVGrid)
2990 {
2991 char *pname, *plongname, *punits, *pstdname;
2992 size_t vctsize = vctsize_echam;
2993 double *vct = vct_echam;
2994
2995 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
2996 {
2997 ncvar_t *ncvar = &ncvars[ncvarid];
2998 if ( ncvar->varStatus == DataVar && ncvar->zaxisID == CDI_UNDEFID )
2999 {
3000 bool is_scalar = false;
3001 bool with_bounds = false;
3002 int zdimid = CDI_UNDEFID;
3003 int zvarid = CDI_UNDEFID;
3004 size_t zsize = 1;
3005 int psvarid = -1;
3006 int p0varid = -1;
3007
3008 int positive = 0;
3009 int ndims = ncvar->ndims;
3010
3011 if ( ncvar->zvarid != -1 && ncvars[ncvar->zvarid].ndims == 0 )
3012 {
3013 zvarid = ncvar->zvarid;
3014 is_scalar = true;
3015 }
3016 else
3017 {
3018 for ( int i = 0; i < ndims; i++ )
3019 {
3020 if (ncvar->dimtype[i] == Z_AXIS) zdimid = ncvar->dimids[i];
3021 }
3022
3023 if (zdimid != CDI_UNDEFID)
3024 {
3025 // zvarid = ncdims[zdimid].ncvarid;
3026 zvarid = (ncvar->zvarid != CDI_UNDEFID) ? ncvar->zvarid : ncdims[zdimid].ncvarid;
3027 zsize = ncdims[zdimid].len;
3028 }
3029 }
3030
3031 if (CDI_Debug) Message("nlevs = %zu", zsize);
3032
3033 double *zvar = NULL;
3034 char **zcvals = NULL;
3035 size_t zclength = 0;
3036
3037 int zaxisType = CDI_UNDEFID;
3038 if (zvarid != CDI_UNDEFID) zaxisType = ncvars[zvarid].zaxistype;
3039 if (zaxisType == CDI_UNDEFID) zaxisType = ZAXIS_GENERIC;
3040
3041 int zdatatype = CDI_DATATYPE_FLT64;
3042 double *lbounds = NULL;
3043 double *ubounds = NULL;
3044
3045 if (zvarid != CDI_UNDEFID)
3046 {
3047 positive = ncvars[zvarid].positive;
3048 pname = ncvars[zvarid].name;
3049 plongname = ncvars[zvarid].longname;
3050 punits = ncvars[zvarid].units;
3051 pstdname = ncvars[zvarid].stdname;
3052 if (ncvars[zvarid].xtype == NC_FLOAT) zdatatype = CDI_DATATYPE_FLT32;
3053 else if (ncvars[zvarid].xtype == NC_INT) zdatatype = CDI_DATATYPE_INT32;
3054 else if (ncvars[zvarid].xtype == NC_SHORT) zdatatype = CDI_DATATYPE_INT16;
3055 // don't change the name !!!
3056 /*
3057 if ((len = strlen(pname)) > 2)
3058 if (pname[len-2] == '_' && isdigit((int) pname[len-1])) pname[len-2] = 0;
3059 */
3060 #ifndef USE_MPI
3061 if (zaxisType == ZAXIS_CHAR && ncvars[zvarid].ndims == 2)
3062 {
3063 zdatatype = CDI_DATATYPE_UINT8;
3064 zclength = ncdims[ncvars[zvarid].dimids[1]].len;
3065 cdf_load_cvals(zsize*zclength, zvarid, ncvar, &zcvals, zsize);
3066 }
3067 #endif
3068
3069 if (zaxisType == ZAXIS_HYBRID && ncvars[zvarid].vct)
3070 {
3071 vct = ncvars[zvarid].vct;
3072 vctsize = ncvars[zvarid].vctsize;
3073
3074 if (ncvars[zvarid].psvarid != -1) psvarid = ncvars[zvarid].psvarid;
3075 if (ncvars[zvarid].p0varid != -1) p0varid = ncvars[zvarid].p0varid;
3076 }
3077
3078 if (zaxisType != ZAXIS_CHAR)
3079 {
3080 zvar = (double*) Malloc(zsize*sizeof(double));
3081 cdf_get_var_double(ncvars[zvarid].ncid, zvarid, zvar);
3082 }
3083
3084 if ( ncvars[zvarid].bounds != CDI_UNDEFID )
3085 {
3086 const int nbdims = ncvars[ncvars[zvarid].bounds].ndims;
3087 if ( nbdims == 2 || is_scalar )
3088 {
3089 const size_t nlevel = is_scalar ? 1 : ncdims[ncvars[ncvars[zvarid].bounds].dimids[0]].len;
3090 const int nvertex = (int)ncdims[ncvars[ncvars[zvarid].bounds].dimids[1-is_scalar]].len;
3091 if ( nlevel == zsize && nvertex == 2 )
3092 {
3093 with_bounds = true;
3094 lbounds = (double *) Malloc(4 * nlevel*sizeof(double));
3095 ubounds = lbounds + nlevel;
3096 double *zbounds = lbounds + 2 * nlevel;
3097 cdf_get_var_double(ncvars[zvarid].ncid, ncvars[zvarid].bounds, zbounds);
3098 for ( size_t i = 0; i < nlevel; ++i )
3099 {
3100 lbounds[i] = zbounds[i*2];
3101 ubounds[i] = zbounds[i*2+1];
3102 }
3103 }
3104 }
3105 }
3106 }
3107 else
3108 {
3109 pname = (zdimid != CDI_UNDEFID) ? ncdims[zdimid].name : NULL;
3110 plongname = NULL;
3111 punits = NULL;
3112 pstdname = NULL;
3113
3114 if ( zsize == 1 && zdimid == CDI_UNDEFID )
3115 {
3116 zaxisType = (ncvar->zaxistype != CDI_UNDEFID) ? ncvar->zaxistype : ZAXIS_SURFACE;
3117 // if ( pname )
3118 {
3119 zvar = (double*) Malloc(sizeof(double));
3120 zvar[0] = 0;
3121 }
3122 }
3123 }
3124
3125 if ( zsize > INT_MAX )
3126 {
3127 Warning("Size limit exceeded for z-axis dimension (limit=%d)!", INT_MAX);
3128 return CDI_EDIMSIZE;
3129 }
3130
3131 ncvar->zaxisID = varDefZaxis(vlistID, zaxisType, (int) zsize, zvar, (const char **)zcvals, zclength, with_bounds, lbounds, ubounds,
3132 (int)vctsize, vct, pname, plongname, punits, zdatatype, 1, 0, -1);
3133
3134 const int zaxisID = ncvar->zaxisID;
3135
3136 if (CDI_CMOR_Mode && zsize == 1 && zaxisType != ZAXIS_HYBRID) zaxisDefScalar(zaxisID);
3137
3138 if (pstdname && *pstdname)
3139 cdiDefKeyBytes(zaxisID, CDI_GLOBAL, CDI_KEY_STDNAME, (const unsigned char*)pstdname, (int)strlen(pstdname)+1);
3140
3141 if (!cdiUUIDIsNull(uuidOfVGrid))
3142 cdiDefKeyBytes(zaxisID, CDI_GLOBAL, CDI_KEY_UUID, uuidOfVGrid, CDI_UUID_SIZE);
3143
3144 if (zaxisType == ZAXIS_HYBRID)
3145 {
3146 if (psvarid != -1) cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_PSNAME, ncvars[psvarid].name);
3147 if (p0varid != -1)
3148 {
3149 double px = 1;
3150 cdf_get_var_double(ncvars[p0varid].ncid, p0varid, &px);
3151 cdiDefKeyFloat(zaxisID, CDI_GLOBAL, CDI_KEY_P0VALUE, px);
3152 cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_P0NAME, ncvars[p0varid].name);
3153 }
3154 }
3155
3156 if (positive > 0) zaxisDefPositive(zaxisID, positive);
3157 if (is_scalar) zaxisDefScalar(zaxisID);
3158
3159 if (zdimid != CDI_UNDEFID)
3160 cdiDefKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_DIMNAME, ncdims[zdimid].name);
3161
3162 if (zvar) Free(zvar);
3163 if (zcvals)
3164 {
3165 for (size_t i = 0; i < zsize; i++) Free(zcvals[i]);
3166 Free(zcvals);
3167 }
3168 if (lbounds) Free(lbounds);
3169
3170 if (zvarid != CDI_UNDEFID)
3171 {
3172 const int ncid = ncvars[zvarid].ncid;
3173 const int nvatts = ncvars[zvarid].natts;
3174 for (int iatt = 0; iatt < nvatts; ++iatt)
3175 {
3176 const int attnum = ncvars[zvarid].atts[iatt];
3177 cdf_set_cdi_attr(ncid, zvarid, attnum, zaxisID, CDI_GLOBAL);
3178 }
3179 }
3180
3181 const int zaxisindex = vlistZaxisIndex(vlistID, zaxisID);
3182 streamptr->zaxisID[zaxisindex] = zdimid >= 0 ? ncdims[zdimid].dimid : zdimid;
3183
3184 if (CDI_Debug) Message("zaxisID %d %d %s", zaxisID, ncvarid, ncvar->name);
3185
3186 for (int ncvarid2 = ncvarid+1; ncvarid2 < nvars; ncvarid2++)
3187 if (ncvars[ncvarid2].varStatus == DataVar && ncvars[ncvarid2].zaxisID == CDI_UNDEFID /*&& ncvars[ncvarid2].zaxistype == CDI_UNDEFID*/)
3188 {
3189 int zvarid2 = CDI_UNDEFID;
3190 if (ncvars[ncvarid2].zvarid != CDI_UNDEFID && ncvars[ncvars[ncvarid2].zvarid].ndims == 0)
3191 zvarid2 = ncvars[ncvarid2].zvarid;
3192
3193 int zdimid2 = CDI_UNDEFID;
3194 ndims = ncvars[ncvarid2].ndims;
3195 for (int i = 0; i < ndims; i++)
3196 {
3197 if (ncvars[ncvarid2].dimtype[i] == Z_AXIS)
3198 zdimid2 = ncvars[ncvarid2].dimids[i];
3199 }
3200
3201 if (zdimid == zdimid2 /* && zvarid == zvarid2 */)
3202 {
3203 if ((zdimid != CDI_UNDEFID && ncvars[ncvarid2].zaxistype == CDI_UNDEFID) ||
3204 (zdimid == CDI_UNDEFID && zvarid != CDI_UNDEFID && zvarid == zvarid2) ||
3205 (zdimid == CDI_UNDEFID && zaxisType == ncvars[ncvarid2].zaxistype) ||
3206 (zdimid == CDI_UNDEFID && zvarid2 == CDI_UNDEFID && ncvars[ncvarid2].zaxistype == CDI_UNDEFID))
3207 {
3208 if (CDI_Debug) Message("zaxisID %d %d %s", zaxisID, ncvarid2, ncvars[ncvarid2].name);
3209 ncvars[ncvarid2].zaxisID = zaxisID;
3210 }
3211 }
3212 }
3213 }
3214 }
3215
3216 return 0;
3217 }
3218
3219
3220 struct cdf_varinfo
3221 {
3222 int varid;
3223 const char *name;
3224 };
3225
3226 static
cdf_cmp_varname(const void * s1,const void * s2)3227 int cdf_cmp_varname(const void *s1, const void *s2)
3228 {
3229 const struct cdf_varinfo *x = (const struct cdf_varinfo *)s1,
3230 *y = (const struct cdf_varinfo *)s2;
3231 return strcmp(x->name, y->name);
3232 }
3233
3234 static
cdf_sort_varnames(int * varids,const int nvars,const ncvar_t * ncvars)3235 void cdf_sort_varnames(int *varids, const int nvars, const ncvar_t *ncvars)
3236 {
3237 struct cdf_varinfo *varInfo
3238 = (struct cdf_varinfo *) Malloc((size_t)nvars * sizeof(struct cdf_varinfo));
3239
3240 for ( int varID = 0; varID < nvars; varID++ )
3241 {
3242 const int ncvarid = varids[varID];
3243 varInfo[varID].varid = ncvarid;
3244 varInfo[varID].name = ncvars[ncvarid].name;
3245 }
3246 qsort(varInfo, (size_t)nvars, sizeof(varInfo[0]), cdf_cmp_varname);
3247 for ( int varID = 0; varID < nvars; varID++ )
3248 {
3249 varids[varID] = varInfo[varID].varid;
3250 }
3251 Free(varInfo);
3252 if ( CDI_Debug )
3253 for ( int i = 0; i < nvars; i++ ) Message("sorted varids[%d] = %d", i, varids[i]);
3254 }
3255
3256 static
cdf_define_code_and_param(int vlistID,int varID)3257 void cdf_define_code_and_param(int vlistID, int varID)
3258 {
3259 if ( vlistInqVarCode(vlistID, varID) == -varID-1 )
3260 {
3261 char name[CDI_MAX_NAME]; name[0] = 0;
3262 vlistInqVarName(vlistID, varID, name);
3263 const size_t len = strlen(name);
3264 if ( len > 3 && isdigit((int) name[3]) )
3265 {
3266 if ( strStartsWith(name, "var") )
3267 {
3268 vlistDefVarCode(vlistID, varID, atoi(name+3));
3269 }
3270 }
3271 else if ( len > 4 && isdigit((int) name[4]) )
3272 {
3273 if ( strStartsWith(name, "code") )
3274 {
3275 vlistDefVarCode(vlistID, varID, atoi(name+4));
3276 }
3277 }
3278 else if ( len > 5 && isdigit((int) name[5]) )
3279 {
3280 if ( strStartsWith(name, "param") )
3281 {
3282 int pnum = -1, pcat = 255, pdis = 255;
3283 sscanf(name+5, "%d.%d.%d", &pnum, &pcat, &pdis);
3284 vlistDefVarParam(vlistID, varID, cdiEncodeParam(pnum, pcat, pdis));
3285 }
3286 }
3287 }
3288 }
3289
3290 static
cdf_define_institut_and_model_id(int vlistID,int varID)3291 void cdf_define_institut_and_model_id(int vlistID, int varID)
3292 {
3293 int varInstID = vlistInqVarInstitut(vlistID, varID);
3294 int varModelID = vlistInqVarModel(vlistID, varID);
3295 int varTableID = vlistInqVarTable(vlistID, varID);
3296 const int code = vlistInqVarCode(vlistID, varID);
3297 if ( CDI_Default_TableID != CDI_UNDEFID )
3298 {
3299 char name[CDI_MAX_NAME]; name[0] = 0;
3300 char longname[CDI_MAX_NAME]; longname[0] = 0;
3301 char units[CDI_MAX_NAME]; units[0] = 0;
3302 tableInqEntry(CDI_Default_TableID, code, -1, name, longname, units);
3303 if ( name[0] )
3304 {
3305 cdiDeleteKey(vlistID, varID, CDI_KEY_NAME);
3306 cdiDeleteKey(vlistID, varID, CDI_KEY_LONGNAME);
3307 cdiDeleteKey(vlistID, varID, CDI_KEY_UNITS);
3308
3309 if ( varTableID != CDI_UNDEFID )
3310 {
3311 cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, name);
3312 if ( longname[0] ) cdiDefKeyString(vlistID, varID, CDI_KEY_LONGNAME, longname);
3313 if ( units[0] ) cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, units);
3314 }
3315 else
3316 {
3317 varTableID = CDI_Default_TableID;
3318 }
3319 }
3320
3321 if ( CDI_Default_ModelID != CDI_UNDEFID ) varModelID = CDI_Default_ModelID;
3322 if ( CDI_Default_InstID != CDI_UNDEFID ) varInstID = CDI_Default_InstID;
3323 }
3324 if ( varInstID != CDI_UNDEFID ) vlistDefVarInstitut(vlistID, varID, varInstID);
3325 if ( varModelID != CDI_UNDEFID ) vlistDefVarModel(vlistID, varID, varModelID);
3326 if ( varTableID != CDI_UNDEFID ) vlistDefVarTable(vlistID, varID, varTableID);
3327 }
3328
3329 // define all input data variables
3330 static
cdf_define_all_vars(stream_t * streamptr,int vlistID,int instID,int modelID,int nvars,int num_ncvars,ncvar_t * ncvars,ncdim_t * ncdims)3331 void cdf_define_all_vars(stream_t *streamptr, int vlistID, int instID, int modelID, int nvars, int num_ncvars, ncvar_t *ncvars, ncdim_t *ncdims)
3332 {
3333 int *varids = (int *) Malloc((size_t)nvars * sizeof (int));
3334 int n = 0;
3335 for (int ncvarid = 0; ncvarid < num_ncvars; ncvarid++)
3336 if (ncvars[ncvarid].varStatus == DataVar) varids[n++] = ncvarid;
3337
3338 if (CDI_Debug) for (int i = 0; i < nvars; i++) Message("varids[%d] = %d", i, varids[i]);
3339
3340 if (streamptr->sortname) cdf_sort_varnames(varids, nvars, ncvars);
3341
3342 for (int varID1 = 0; varID1 < nvars; varID1++)
3343 {
3344 const int ncvarid = varids[varID1];
3345 ncvar_t *ncvar = &ncvars[ncvarid];
3346 const int gridID = ncvar->gridID;
3347 const int zaxisID = ncvar->zaxisID;
3348
3349 stream_new_var(streamptr, gridID, zaxisID, CDI_UNDEFID);
3350 const int varID = vlistDefVar(vlistID, gridID, zaxisID, ncvar->timetype);
3351
3352 #ifdef HAVE_NETCDF4
3353 if (ncvar->ldeflate) vlistDefVarCompType(vlistID, varID, CDI_COMPRESS_ZIP);
3354 if (ncvar->lszip) vlistDefVarCompType(vlistID, varID, CDI_COMPRESS_SZIP);
3355 if (ncvar->chunked && ncvar->chunktype != CDI_UNDEFID) vlistDefVarChunkType(vlistID, varID, ncvar->chunktype);
3356 #endif
3357
3358 streamptr->vars[varID1].defmiss = false;
3359 streamptr->vars[varID1].ncvarid = ncvarid;
3360
3361 cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, ncvar->name);
3362 if (ncvar->param != CDI_UNDEFID) vlistDefVarParam(vlistID, varID, ncvar->param);
3363 if (ncvar->code != CDI_UNDEFID) vlistDefVarCode(vlistID, varID, ncvar->code);
3364 if (ncvar->code != CDI_UNDEFID)
3365 vlistDefVarParam(vlistID, varID, cdiEncodeParam(ncvar->code, ncvar->tabnum, 255));
3366 if (ncvar->longname[0]) cdiDefKeyString(vlistID, varID, CDI_KEY_LONGNAME, ncvar->longname);
3367 if (ncvar->stdname[0]) cdiDefKeyString(vlistID, varID, CDI_KEY_STDNAME, ncvar->stdname);
3368 if (ncvar->units[0]) cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, ncvar->units);
3369
3370 if (ncvar->lvalidrange) vlistDefVarValidrange(vlistID, varID, ncvar->validrange);
3371
3372 if (IS_NOT_EQUAL(ncvar->addoffset, 0)) vlistDefVarAddoffset(vlistID, varID, ncvar->addoffset);
3373 if (IS_NOT_EQUAL(ncvar->scalefactor, 1)) vlistDefVarScalefactor(vlistID, varID, ncvar->scalefactor);
3374
3375 vlistDefVarDatatype(vlistID, varID, cdfInqDatatype(streamptr, ncvar->xtype, ncvar->lunsigned));
3376
3377 vlistDefVarInstitut(vlistID, varID, instID);
3378 vlistDefVarModel(vlistID, varID, modelID);
3379 if (ncvar->tableID != CDI_UNDEFID) vlistDefVarTable(vlistID, varID, ncvar->tableID);
3380
3381 if (ncvar->deffillval == false && ncvar->defmissval)
3382 {
3383 ncvar->deffillval = true;
3384 ncvar->fillval = ncvar->missval;
3385 }
3386
3387 if (ncvar->deffillval) vlistDefVarMissval(vlistID, varID, ncvar->fillval);
3388
3389 if ( CDI_Debug )
3390 Message("varID = %d gridID = %d zaxisID = %d", varID,
3391 vlistInqVarGrid(vlistID, varID), vlistInqVarZaxis(vlistID, varID));
3392
3393 const int gridindex = vlistGridIndex(vlistID, gridID);
3394 const int xdimid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
3395 const int ydimid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_Y];
3396
3397 const int zaxisindex = vlistZaxisIndex(vlistID, zaxisID);
3398 const int zdimid = streamptr->zaxisID[zaxisindex];
3399
3400 const int ndims = ncvar->ndims;
3401 static const int ipow10[4] = {1, 10, 100, 1000};
3402
3403 const int iodim = (ncvar->timetype != TIME_CONSTANT);
3404
3405 const int *dimids = ncvar->dimids;
3406
3407 int ixyz = 0;
3408 if (gridInqType(gridID) == GRID_UNSTRUCTURED && ndims-iodim <= 2 && (ydimid == xdimid || ydimid == CDI_UNDEFID))
3409 {
3410 ixyz = (xdimid == ncdims[dimids[ndims-1]].dimid) ? 321 : 213;
3411 }
3412 else
3413 {
3414 for (int idim = iodim; idim < ndims; idim++)
3415 {
3416 const int dimid = ncdims[dimids[idim]].dimid;
3417 if (xdimid == dimid) ixyz += 1*ipow10[ndims-idim-1];
3418 else if (ydimid == dimid) ixyz += 2*ipow10[ndims-idim-1];
3419 else if (zdimid == dimid) ixyz += 3*ipow10[ndims-idim-1];
3420 }
3421 }
3422
3423 if (ncvar->isCubeSphere) ixyz = 0;
3424 vlistDefVarXYZ(vlistID, varID, ixyz);
3425 /*
3426 printf("ixyz %d\n", ixyz);
3427 printf("ndims %d\n", ncvar->ndims);
3428 for ( int i = 0; i < ncvar->ndims; ++i )
3429 printf("dimids: %d %d\n", i, dimids[i]);
3430 printf("xdimid, ydimid %d %d\n", xdimid, ydimid);
3431 */
3432 if (ncvar->numberOfForecastsInEnsemble != -1)
3433 {
3434 cdiDefKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, ncvar->typeOfEnsembleForecast);
3435 cdiDefKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, ncvar->numberOfForecastsInEnsemble);
3436 cdiDefKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, ncvar->perturbationNumber);
3437 }
3438
3439 if (ncvar->extra[0] != 0) vlistDefVarExtra(vlistID, varID, ncvar->extra);
3440 }
3441
3442 for (int varID = 0; varID < nvars; varID++)
3443 {
3444 const int ncvarid = varids[varID];
3445 ncvar_t *ncvar = &ncvars[ncvarid];
3446 const int ncid = ncvar->ncid;
3447 const int nvatts = ncvar->natts;
3448 for (int iatt = 0; iatt < nvatts; ++iatt) cdf_set_cdi_attr(ncid, ncvarid, ncvar->atts[iatt], vlistID, varID);
3449
3450 if (ncvar->atts) Free(ncvar->atts);
3451 if (ncvar->vct) Free(ncvar->vct);
3452
3453 ncvar->atts = NULL;
3454 ncvar->vct = NULL;
3455 }
3456
3457 // release mem of not freed attributes
3458 for (int ncvarid = 0; ncvarid < num_ncvars; ncvarid++)
3459 if (ncvars[ncvarid].atts) Free(ncvars[ncvarid].atts);
3460
3461 if (varids) Free(varids);
3462
3463 for (int varID = 0; varID < nvars; varID++) cdf_define_code_and_param(vlistID, varID);
3464 for (int varID = 0; varID < nvars; varID++) cdf_define_institut_and_model_id(vlistID, varID);
3465 }
3466
3467 static
cdf_copy_attint(int fileID,int vlistID,nc_type xtype,size_t attlen,char * attname)3468 void cdf_copy_attint(int fileID, int vlistID, nc_type xtype, size_t attlen, char *attname)
3469 {
3470 int attint[8];
3471 int *pattint = (attlen > 8) ? (int*) malloc(attlen*sizeof(int)) : attint;
3472 cdfGetAttInt(fileID, NC_GLOBAL, attname, attlen, pattint);
3473 const int datatype = (xtype == NC_SHORT) ? CDI_DATATYPE_INT16 : CDI_DATATYPE_INT32;
3474 cdiDefAttInt(vlistID, CDI_GLOBAL, attname, datatype, (int)attlen, pattint);
3475 if (attlen > 8) free(pattint);
3476 }
3477
3478 static
cdf_copy_attflt(int fileID,int vlistID,nc_type xtype,size_t attlen,char * attname)3479 void cdf_copy_attflt(int fileID, int vlistID, nc_type xtype, size_t attlen, char *attname)
3480 {
3481 double attflt[8];
3482 double *pattflt = (attlen > 8) ? (double*) malloc(attlen*sizeof(double)) : attflt;
3483 cdfGetAttDouble(fileID, NC_GLOBAL, attname, attlen, pattflt);
3484 const int datatype = (xtype == NC_FLOAT) ? CDI_DATATYPE_FLT32 : CDI_DATATYPE_FLT64;
3485 cdiDefAttFlt(vlistID, CDI_GLOBAL, attname, datatype, (int)attlen, pattflt);
3486 if (attlen > 8) free(pattflt);
3487 }
3488
3489 static
check_cube_sphere(int vlistID,int nvars,ncvar_t * ncvars,ncdim_t * ncdims)3490 bool check_cube_sphere(int vlistID, int nvars, ncvar_t *ncvars, ncdim_t *ncdims)
3491 {
3492 bool isGeosData = false;
3493 const int numNames = 4;
3494 const char *attnames[] = { "additional_vars", "file_format_version", "gridspec_file", "grid_mapping_name" };
3495 const char *grid_mapping = "gnomonic cubed-sphere";
3496 char attstring[256];
3497 int nf_dimid = -1, ncontact_dimid = -1;
3498
3499 int numFound = 0;
3500 for (int i = 0; i < numNames; ++i)
3501 if (0 == cdiInqAttTxt(vlistID, CDI_GLOBAL, attnames[i], (int)sizeof(attstring), attstring)) numFound++;
3502
3503 if (numFound == numNames && strStartsWith(attstring, grid_mapping))
3504 {
3505 for (int i = 0; i < numNames; ++i) cdiDelAtt(vlistID, CDI_GLOBAL, attnames[i]);
3506
3507 const char *nf_name = "nf";
3508 const char *ncontact_name = "ncontact";
3509 for (int varid = 0; varid < nvars; ++varid)
3510 {
3511 ncvar_t *ncvar = &ncvars[varid];
3512 if (ncvar->ndims == 1)
3513 {
3514 int dimid = ncvar->dimids[0];
3515 if (ncdims[dimid].len == 6 && strIsEqual(nf_name, ncvar->name))
3516 {
3517 isGeosData = true;
3518 nf_dimid = ncvar->dimids[0];
3519 }
3520 if (ncdims[dimid].len == 4 && strIsEqual(ncontact_name, ncvar->name)) ncontact_dimid = ncvar->dimids[0];
3521 }
3522 if (isGeosData && ncontact_dimid != -1) break;
3523 }
3524 }
3525
3526 if (isGeosData)
3527 {
3528 ncdims[nf_dimid].dimtype = E_AXIS;
3529 for (int varid = 0; varid < nvars; ++varid)
3530 {
3531 ncvar_t *ncvar = &ncvars[varid];
3532 if (strIsEqual("orientation", ncvar->name) ||
3533 strIsEqual("anchor", ncvar->name) ||
3534 strIsEqual("contacts", ncvar->name))
3535 cdf_set_var(ncvar, CoordVar);
3536 }
3537
3538 for (int varid = 0; varid < nvars; ++varid)
3539 {
3540 ncvar_t *ncvar = &ncvars[varid];
3541 const int ndims = ncvar->ndims;
3542 if (ndims >= 3 && ncvar->dimids[ndims - 3] == nf_dimid && ncvar->ncoordvars == 2 && ncvar->gmapid != -1)
3543 ncvar->isCubeSphere = true;
3544 }
3545
3546 int xvarid = -1, yvarid = -1;
3547 int xboundsid = -1, yboundsid = -1;
3548 for (int varid = 0; varid < nvars; ++varid)
3549 {
3550 ncvar_t *ncvar = &ncvars[varid];
3551 const int ndims = ncvar->ndims;
3552 if (ndims == 3)
3553 {
3554 if (strIsEqual("lons", ncvar->name)) xvarid = varid;
3555 else if (strIsEqual("lats", ncvar->name)) yvarid = varid;
3556 else if (strIsEqual("corner_lons", ncvar->name)) xboundsid = varid;
3557 else if (strIsEqual("corner_lats", ncvar->name)) yboundsid = varid;
3558 }
3559 if (xvarid != -1 && xboundsid != -1 && yvarid != -1 && yboundsid != -1)
3560 {
3561 cdf_set_var(&ncvars[xboundsid], CoordVar);
3562 cdf_set_var(&ncvars[yboundsid], CoordVar);
3563 ncvars[xvarid].bounds = xboundsid;
3564 ncvars[yvarid].bounds = yboundsid;
3565 break;
3566 }
3567 }
3568 }
3569
3570 return isGeosData;
3571 }
3572
3573 static
cdf_scan_global_attr(int fileID,int vlistID,int ngatts,int * instID,int * modelID,bool * ucla_les,unsigned char * uuidOfVGrid,GridInfo * gridInfo)3574 void cdf_scan_global_attr(int fileID, int vlistID, int ngatts, int *instID, int *modelID, bool *ucla_les, unsigned char *uuidOfVGrid, GridInfo *gridInfo)
3575 {
3576 nc_type xtype;
3577 size_t attlen;
3578 char attname[CDI_MAX_NAME];
3579 char attstring[65636];
3580
3581 for ( int iatt = 0; iatt < ngatts; iatt++ )
3582 {
3583 cdf_inq_attname(fileID, NC_GLOBAL, iatt, attname);
3584 cdf_inq_atttype(fileID, NC_GLOBAL, attname, &xtype);
3585 cdf_inq_attlen(fileID, NC_GLOBAL, attname, &attlen);
3586
3587 if ( xtypeIsText(xtype) )
3588 {
3589 cdfGetAttText(fileID, NC_GLOBAL, attname, sizeof(attstring), attstring);
3590
3591 size_t attstrlen = strlen(attstring);
3592
3593 if ( attlen > 0 && attstring[0] != 0 )
3594 {
3595 if ( strIsEqual(attname, "institution") )
3596 {
3597 *instID = institutInq(0, 0, NULL, attstring);
3598 if ( *instID == CDI_UNDEFID )
3599 *instID = institutDef(0, 0, NULL, attstring);
3600 }
3601 else if ( strIsEqual(attname, "source") )
3602 {
3603 *modelID = modelInq(-1, 0, attstring);
3604 if ( *modelID == CDI_UNDEFID )
3605 *modelID = modelDef(-1, 0, attstring);
3606 }
3607 else if ( strIsEqual(attname, "Source") )
3608 {
3609 if (strStartsWith(attstring, "UCLA-LES")) *ucla_les = true;
3610 }
3611 /*
3612 else if ( strIsEqual(attname, "Conventions") )
3613 {
3614 }
3615 */
3616 else if ( strIsEqual(attname, "_NCProperties") )
3617 {
3618 }
3619 else if ( strIsEqual(attname, "CDI") )
3620 {
3621 }
3622 else if ( strIsEqual(attname, "CDO") )
3623 {
3624 }
3625 /*
3626 else if ( strIsEqual(attname, "forecast_reference_time") )
3627 {
3628 memcpy(fcreftime, attstring, attstrlen+1);
3629 }
3630 */
3631 else if ( strIsEqual(attname, "grid_file_uri") )
3632 {
3633 memcpy(gridInfo->gridfile, attstring, attstrlen+1);
3634 }
3635 else if ( strIsEqual(attname, "uuidOfHGrid") && attstrlen == 36 )
3636 {
3637 attstring[36] = 0;
3638 cdiStr2UUID(attstring, gridInfo->uuid);
3639 }
3640 else if ( strIsEqual(attname, "uuidOfVGrid") && attstrlen == 36 )
3641 {
3642 attstring[36] = 0;
3643 cdiStr2UUID(attstring, uuidOfVGrid);
3644 }
3645 else
3646 {
3647 if ( strIsEqual(attname, "ICON_grid_file_uri") && gridInfo->gridfile[0] == 0 )
3648 memcpy(gridInfo->gridfile, attstring, attstrlen+1);
3649
3650 cdiDefAttTxt(vlistID, CDI_GLOBAL, attname, (int)attstrlen, attstring);
3651 }
3652 }
3653 }
3654 else if ( xtype == NC_SHORT || xtype == NC_INT )
3655 {
3656 if ( strIsEqual(attname, "number_of_grid_used") )
3657 {
3658 gridInfo->number_of_grid_used = CDI_UNDEFID;
3659 cdfGetAttInt(fileID, NC_GLOBAL, attname, 1, &gridInfo->number_of_grid_used);
3660 }
3661 else
3662 {
3663 cdf_copy_attint(fileID, vlistID, xtype, attlen, attname);
3664 }
3665 }
3666 else if ( xtype == NC_FLOAT || xtype == NC_DOUBLE )
3667 {
3668 cdf_copy_attflt(fileID, vlistID, xtype, attlen, attname);
3669 }
3670 }
3671 }
3672
3673 static
find_leadtime(int nvars,ncvar_t * ncvars,int timedimid)3674 int find_leadtime(int nvars, ncvar_t *ncvars, int timedimid)
3675 {
3676 int leadtime_id = CDI_UNDEFID;
3677
3678 for ( int ncvarid = 0; ncvarid < nvars; ncvarid++ )
3679 {
3680 ncvar_t *ncvar = &ncvars[ncvarid];
3681 if ( ncvar->ndims == 1 && timedimid == ncvar->dimids[0])
3682 if ( ncvar->stdname[0] && strIsEqual(ncvar->stdname, "forecast_period") )
3683 {
3684 leadtime_id = ncvarid;
3685 break;
3686 }
3687 }
3688
3689 return leadtime_id;
3690 }
3691
3692 static
find_time_vars(int nvars,ncvar_t * ncvars,ncdim_t * ncdims,int timedimid,stream_t * streamptr,bool * time_has_units,bool * time_has_bounds,bool * time_climatology)3693 void find_time_vars(int nvars, ncvar_t *ncvars, ncdim_t *ncdims, int timedimid, stream_t *streamptr,
3694 bool *time_has_units, bool *time_has_bounds, bool *time_climatology)
3695 {
3696 int ncvarid;
3697
3698 if ( timedimid == CDI_UNDEFID )
3699 {
3700 char timeunits[CDI_MAX_NAME];
3701
3702 for ( ncvarid = 0; ncvarid < nvars; ncvarid++ )
3703 {
3704 ncvar_t *ncvar = &ncvars[ncvarid];
3705 if ( ncvar->ndims == 0 && strIsEqual(ncvar->name, "time") )
3706 {
3707 if ( ncvar->units[0] )
3708 {
3709 strcpy(timeunits, ncvar->units);
3710 strToLower(timeunits);
3711
3712 if ( is_time_units(timeunits) )
3713 {
3714 streamptr->basetime.ncvarid = ncvarid;
3715 break;
3716 }
3717 }
3718 }
3719 }
3720 }
3721 else
3722 {
3723 bool ltimevar = false;
3724
3725 if ( ncdims[timedimid].ncvarid != CDI_UNDEFID )
3726 {
3727 streamptr->basetime.ncvarid = ncdims[timedimid].ncvarid;
3728 ltimevar = true;
3729 }
3730
3731 for ( ncvarid = 0; ncvarid < nvars; ncvarid++ )
3732 {
3733 ncvar_t *ncvar = &ncvars[ncvarid];
3734 if ( ncvarid != streamptr->basetime.ncvarid &&
3735 ncvar->ndims == 1 &&
3736 timedimid == ncvar->dimids[0] &&
3737 !xtypeIsText(ncvar->xtype) &&
3738 is_timeaxis_units(ncvar->units) )
3739 {
3740 ncvar->varStatus = CoordVar;
3741
3742 if ( !ltimevar )
3743 {
3744 streamptr->basetime.ncvarid = ncvarid;
3745 ltimevar = true;
3746 if ( CDI_Debug ) fprintf(stderr, "timevar %s\n", ncvar->name);
3747 }
3748 else
3749 {
3750 Warning("Found more than one time variable, skipped variable %s!", ncvar->name);
3751 }
3752 }
3753 }
3754
3755 if ( ltimevar == false ) // search for WRF time description
3756 {
3757 for ( ncvarid = 0; ncvarid < nvars; ncvarid++ )
3758 {
3759 ncvar_t *ncvar = &ncvars[ncvarid];
3760 if ( ncvarid != streamptr->basetime.ncvarid &&
3761 ncvar->ndims == 2 &&
3762 timedimid == ncvar->dimids[0] &&
3763 xtypeIsText(ncvar->xtype) &&
3764 (ncdims[ncvar->dimids[1]].len == 19 || ncdims[ncvar->dimids[1]].len == 64))
3765 {
3766 ncvar->isTaxis = true;
3767 streamptr->basetime.ncvarid = ncvarid;
3768 streamptr->basetime.lwrf = true;
3769 break;
3770 }
3771 }
3772 }
3773
3774 // time varID
3775 ncvarid = streamptr->basetime.ncvarid;
3776
3777 if ( ncvarid == CDI_UNDEFID ) Warning("Time variable >%s< not found!", ncdims[timedimid].name);
3778 }
3779
3780 // time varID
3781 ncvarid = streamptr->basetime.ncvarid;
3782
3783 if ( ncvarid != CDI_UNDEFID && streamptr->basetime.lwrf == false )
3784 {
3785 ncvar_t *ncvar = &ncvars[ncvarid];
3786 if ( ncvar->units[0] != 0 ) *time_has_units = true;
3787
3788 if ( ncvar->bounds != CDI_UNDEFID )
3789 {
3790 int nbdims = ncvars[ncvar->bounds].ndims;
3791 if ( nbdims == 2 )
3792 {
3793 int len = (int) ncdims[ncvars[ncvar->bounds].dimids[nbdims-1]].len;
3794 if ( len == 2 && timedimid == ncvars[ncvar->bounds].dimids[0] )
3795 {
3796 *time_has_bounds = true;
3797 streamptr->basetime.ncvarboundsid = ncvar->bounds;
3798 if ( ncvar->isClimatology ) *time_climatology = true;
3799 }
3800 }
3801 }
3802 }
3803 }
3804
3805 static
read_vct_echam(int fileID,int nvars,ncvar_t * ncvars,ncdim_t * ncdims,double ** vct,size_t * pvctsize)3806 void read_vct_echam(int fileID, int nvars, ncvar_t *ncvars, ncdim_t *ncdims, double **vct, size_t *pvctsize)
3807 {
3808 // find ECHAM VCT
3809 int nvcth_id = CDI_UNDEFID, vcta_id = CDI_UNDEFID, vctb_id = CDI_UNDEFID;
3810 // int p0_id = CDI_UNDEFID;
3811
3812 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
3813 {
3814 ncvar_t *ncvar = &ncvars[ncvarid];
3815 const char *name = ncvar->name;
3816 if (ncvar->ndims == 1)
3817 {
3818 const size_t len = strlen(name);
3819 if (len == 4 && name[0] == 'h' && name[1] == 'y')
3820 {
3821 if (name[2] == 'a' && name[3] == 'i') // hyai
3822 {
3823 vcta_id = ncvarid;
3824 nvcth_id = ncvar->dimids[0];
3825 ncvar->varStatus = CoordVar;
3826 }
3827 else if (name[2] == 'b' && name[3] == 'i') // hybi
3828 {
3829 vctb_id = ncvarid;
3830 nvcth_id = ncvar->dimids[0];
3831 ncvar->varStatus = CoordVar;
3832 }
3833 else if ((name[2] == 'a' || name[2] == 'b') && name[3] == 'm' )
3834 {
3835 ncvar->varStatus = CoordVar; // hyam or hybm
3836 }
3837 }
3838 }
3839 /*
3840 else if (ncvar->ndims == 0)
3841 {
3842 const size_t len = strlen(name);
3843 if (len == 2 && name[0] == 'P' && name[1] == '0') p0_id = ncvarid;
3844 }
3845 */
3846 }
3847
3848 // read VCT
3849 if (nvcth_id != CDI_UNDEFID && vcta_id != CDI_UNDEFID && vctb_id != CDI_UNDEFID)
3850 {
3851 const size_t vctsize = 2 * ncdims[nvcth_id].len;
3852 *vct = (double *) Malloc(vctsize * sizeof(double));
3853 cdf_get_var_double(fileID, vcta_id, *vct);
3854 cdf_get_var_double(fileID, vctb_id, *vct + vctsize/2);
3855 *pvctsize = vctsize;
3856 /*
3857 if (p0_id != CDI_UNDEFID)
3858 {
3859 double p0;
3860 cdf_get_var_double(fileID, p0_id, &p0);
3861 }
3862 */
3863 }
3864 }
3865
3866 static
cdf_set_ucla_dimtype(int ndims,ncdim_t * ncdims,ncvar_t * ncvars)3867 void cdf_set_ucla_dimtype(int ndims, ncdim_t *ncdims, ncvar_t *ncvars)
3868 {
3869 for ( int ncdimid = 0; ncdimid < ndims; ncdimid++ )
3870 {
3871 int ncvarid = ncdims[ncdimid].ncvarid;
3872 if (ncvarid != -1)
3873 {
3874 ncvar_t *ncvar = &ncvars[ncvarid];
3875 if ( ncdims[ncdimid].dimtype == CDI_UNDEFID && ncvar->units[0] == 'm' )
3876 {
3877 if ( ncvar->name[0] == 'x' ) ncdims[ncdimid].dimtype = X_AXIS;
3878 else if ( ncvar->name[0] == 'y' ) ncdims[ncdimid].dimtype = Y_AXIS;
3879 else if ( ncvar->name[0] == 'z' ) ncdims[ncdimid].dimtype = Z_AXIS;
3880 }
3881 }
3882 }
3883 }
3884
3885 static
cdf_check_variables(stream_t * streamptr,int nvars,ncvar_t * ncvars,size_t ntsteps,int timedimid)3886 int cdf_check_variables(stream_t *streamptr, int nvars, ncvar_t *ncvars, size_t ntsteps, int timedimid)
3887 {
3888 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
3889 {
3890 ncvar_t *ncvar = &ncvars[ncvarid];
3891 if (ncvar->isTaxis && ncvar->ndims == 2)
3892 {
3893 ncvar->varStatus = CoordVar;
3894 continue;
3895 }
3896
3897 if (ncvar->varStatus == UndefVar && ncvar->ndims > 1 && timedimid != CDI_UNDEFID && timedimid == ncvar->dimids[0])
3898 cdf_set_var(ncvar, DataVar);
3899
3900 if (ncvar->varStatus == UndefVar)
3901 {
3902 if (ncvar->ndims == 0)
3903 cdf_set_var(ncvar, nvars == 1 ? DataVar : CoordVar);
3904 else if (ncvar->ndims > 0)
3905 cdf_set_var(ncvar, DataVar);
3906 else
3907 {
3908 ncvar->varStatus = CoordVar;
3909 Warning("Variable %s has an unknown type, skipped!", ncvar->name);
3910 }
3911 }
3912
3913 if (ncvar->varStatus == CoordVar) continue;
3914
3915 if ((ncvar->ndims > 4 && !ncvar->isCubeSphere) || ncvar->ndims > 5)
3916 {
3917 ncvar->varStatus = CoordVar;
3918 Warning("%d dimensional variables are not supported, skipped variable %s!",
3919 ncvar->ndims, ncvar->name);
3920 continue;
3921 }
3922
3923 if (((ncvar->ndims == 4 && !ncvar->isCubeSphere) || ncvar->ndims == 5) && timedimid == CDI_UNDEFID)
3924 {
3925 ncvar->varStatus = CoordVar;
3926 Warning("%d dimensional variables without time dimension are not supported, skipped variable %s!",
3927 ncvar->ndims, ncvar->name);
3928 continue;
3929 }
3930
3931 if (xtypeIsText(ncvar->xtype))
3932 {
3933 ncvar->varStatus = CoordVar;
3934 Warning("Unsupported data type (char/string), skipped variable %s!", ncvar->name);
3935 continue;
3936 }
3937
3938 if (cdfInqDatatype(streamptr, ncvar->xtype, ncvar->lunsigned) == -1)
3939 {
3940 ncvar->varStatus = CoordVar;
3941 Warning("Unsupported data type, skipped variable %s!", ncvar->name);
3942 continue;
3943 }
3944
3945 if (timedimid != CDI_UNDEFID && ntsteps == 0 && ncvar->ndims > 0)
3946 {
3947 if (timedimid == ncvar->dimids[0])
3948 {
3949 ncvar->varStatus = CoordVar;
3950 Warning("Number of time steps undefined, skipped variable %s!", ncvar->name);
3951 continue;
3952 }
3953 }
3954 }
3955
3956 return timedimid;
3957 }
3958
3959 static
cdfVerifyVars(int nvars,ncvar_t * ncvars,ncdim_t * ncdims)3960 void cdfVerifyVars(int nvars, ncvar_t *ncvars, ncdim_t *ncdims)
3961 {
3962 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
3963 {
3964 ncvar_t *ncvar = &ncvars[ncvarid];
3965 if (ncvar->varStatus == DataVar && ncvar->ndims > 0)
3966 {
3967 int ndims = 0;
3968 for (int i = 0; i < ncvar->ndims; ++i)
3969 {
3970 if (ncvar->dimtype[i] == T_AXIS) ndims++;
3971 else if (ncvar->dimtype[i] == E_AXIS) ndims++;
3972 else if (ncvar->dimtype[i] == Z_AXIS) ndims++;
3973 else if (ncvar->dimtype[i] == Y_AXIS) ndims++;
3974 else if (ncvar->dimtype[i] == X_AXIS) ndims++;
3975 }
3976
3977 if ( ncvar->ndims != ndims )
3978 {
3979 ncvar->varStatus = CoordVar;
3980 Warning("Inconsistent number of dimensions, skipped variable %s!", ncvar->name);
3981 }
3982
3983 int zdimid = -1;
3984 for (int i = 0; i < ncvar->ndims; ++i)
3985 {
3986 if (ncvar->dimtype[i] == Z_AXIS) zdimid = ncvar->dimids[i];
3987 }
3988
3989 int zaxisID = ncvar->zaxisID;
3990 if (zaxisInqScalar(zaxisID) && zdimid != -1)
3991 {
3992 ncvar->varStatus = CoordVar;
3993 Warning("Unsupported dimension >%s<, skipped variable %s!", ncdims[zdimid].name, ncvar->name);
3994 }
3995 }
3996 }
3997 }
3998
3999
cdfInqContents(stream_t * streamptr)4000 int cdfInqContents(stream_t *streamptr)
4001 {
4002 int ndims, nvars, ngatts, unlimdimid;
4003 bool time_has_units = false;
4004 bool time_has_bounds = false;
4005 bool time_climatology = false;
4006 int leadtime_id = CDI_UNDEFID;
4007 int format = 0;
4008 char fcreftime[CDI_MAX_NAME];
4009 GridInfo gridInfo;
4010 gridInfo.gridfile[0] = 0;
4011 memset(gridInfo.uuid, 0, CDI_UUID_SIZE);
4012 gridInfo.number_of_grid_used = CDI_UNDEFID;
4013 gridInfo.timedimid = CDI_UNDEFID;
4014
4015 fcreftime[0] = 0;
4016
4017 const int vlistID = streamptr->vlistID;
4018 const int fileID = streamptr->fileID;
4019
4020 if (CDI_Debug) Message("streamID = %d, fileID = %d", streamptr->self, fileID);
4021
4022 #ifdef HAVE_NETCDF4
4023 nc_inq_format(fileID, &format);
4024 #endif
4025
4026 cdf_inq(fileID, &ndims , &nvars, &ngatts, &unlimdimid);
4027
4028 if (CDI_Debug) Message("root: ndims %d, nvars %d, ngatts %d", ndims, nvars, ngatts);
4029 /*
4030 if ( ndims == 0 )
4031 {
4032 Warning("No dimensions found!");
4033 return CDI_EUFSTRUCT;
4034 }
4035 */
4036 // alloc ncdims
4037 ncdim_t *ncdims = ndims ? (ncdim_t *) Malloc((size_t)ndims * sizeof(ncdim_t)) : NULL;
4038 init_ncdims(ndims, ncdims);
4039
4040 if (ndims)
4041 {
4042 int ncdimid = 0;
4043 for (int dimid = 0; dimid < NC_MAX_DIMS; ++dimid)
4044 {
4045 size_t dimlen;
4046 int status = nc_inq_dimlen(fileID, dimid, &dimlen);
4047 if (status == NC_NOERR)
4048 {
4049 ncdims[ncdimid++].dimid = dimid;
4050 if (ncdimid == ndims) break;
4051 }
4052 }
4053 }
4054
4055 #ifdef HAVE_NETCDF4
4056 if (format == NC_FORMAT_NETCDF4)
4057 {
4058 int numgrps = 0;
4059 int ncids[NC_MAX_VARS];
4060 char gname[CDI_MAX_NAME];
4061 nc_inq_grps(fileID, &numgrps, ncids);
4062 for (int i = 0; i < numgrps; ++i)
4063 {
4064 int ncid = ncids[i];
4065 nc_inq_grpname(ncid, gname);
4066 int gndims, gnvars, gngatts, gunlimdimid;
4067 cdf_inq(ncid, &gndims , &gnvars, &gngatts, &gunlimdimid);
4068
4069 if (CDI_Debug) Message("%s: ndims %d, nvars %d, ngatts %d", gname, gndims, gnvars, gngatts);
4070
4071 if (gndims == 0)
4072 {
4073 }
4074 }
4075 if (numgrps) Warning("NetCDF4 groups not supported! Found %d root group%s.", numgrps, (numgrps > 1) ? "s" : "");
4076 }
4077 #endif
4078
4079 if (nvars == 0)
4080 {
4081 Warning("No arrays found!");
4082 return CDI_EUFSTRUCT;
4083 }
4084
4085 // alloc ncvars
4086 ncvar_t *ncvars = (ncvar_t *) Malloc((size_t)nvars * sizeof (ncvar_t));
4087 init_ncvars(nvars, ncvars);
4088
4089 for (int ncvarid = 0; ncvarid < nvars; ++ncvarid) ncvars[ncvarid].ncid = fileID;
4090
4091
4092 // scan global attributes
4093 int instID = CDI_UNDEFID;
4094 int modelID = CDI_UNDEFID;
4095 bool ucla_les = false;
4096 unsigned char uuidOfVGrid[CDI_UUID_SIZE] = { 0 };
4097 cdf_scan_global_attr(fileID, vlistID, ngatts, &instID, &modelID, &ucla_les, uuidOfVGrid, &gridInfo);
4098
4099 // find time dim
4100 int timedimid = (unlimdimid >= 0) ? unlimdimid : cdf_time_dimid(fileID, ndims, nvars, ncdims);
4101
4102 streamptr->basetime.ncdimid = timedimid;
4103
4104 size_t ntsteps = 0;
4105 if ( timedimid != CDI_UNDEFID ) cdf_inq_dimlen(fileID, ncdims[timedimid].dimid, &ntsteps);
4106 if ( ntsteps > INT_MAX )
4107 {
4108 Warning("Size limit exceeded for time dimension (limit=%d)!", INT_MAX);
4109 return CDI_EDIMSIZE;
4110 }
4111
4112 if (CDI_Debug) Message("Number of timesteps = %zu", ntsteps);
4113 if (CDI_Debug) Message("Time dimid = %d", streamptr->basetime.ncdimid);
4114
4115 // read ncdims
4116 for (int ncdimid = 0; ncdimid < ndims; ncdimid++)
4117 {
4118 cdf_inq_dimlen(fileID, ncdims[ncdimid].dimid, &ncdims[ncdimid].len);
4119 cdf_inq_dimname(fileID, ncdims[ncdimid].dimid, ncdims[ncdimid].name);
4120 if ( timedimid == ncdimid )
4121 ncdims[ncdimid].dimtype = T_AXIS;
4122 }
4123
4124 if (CDI_Debug) cdf_print_vars(ncvars, nvars, "cdfScanVarAttr");
4125
4126 // scan attributes of all variables
4127 cdfScanVarAttr(nvars, ncvars, ndims, ncdims, timedimid, modelID, format);
4128
4129 cdfVerifyVarAttr(nvars, ncvars, ncdims);
4130
4131
4132 if (CDI_Convert_Cubesphere)
4133 {
4134 bool isGeosData = check_cube_sphere(vlistID, nvars, ncvars, ncdims);
4135 if (CDI_Debug) Message("isGeosData %d", isGeosData);
4136 }
4137
4138 if ( CDI_Debug ) cdf_print_vars(ncvars, nvars, "find coordinates vars");
4139
4140 // find coordinates vars
4141 for (int ncdimid = 0; ncdimid < ndims; ncdimid++)
4142 {
4143 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
4144 {
4145 ncvar_t *ncvar = &ncvars[ncvarid];
4146 if (ncvar->ndims == 1)
4147 {
4148 if (timedimid != CDI_UNDEFID && timedimid == ncvar->dimids[0])
4149 {
4150 if (ncvar->varStatus != CoordVar) cdf_set_var(ncvar, DataVar);
4151 }
4152 else
4153 {
4154 // if ( ncvar->varStatus != DataVar ) cdf_set_var(ncvar, CoordVar);
4155 }
4156 // if ( ncvar->varStatus != DataVar ) cdf_set_var(ncvar, CoordVar);
4157
4158 if (ncdimid == ncvar->dimids[0] && ncdims[ncdimid].ncvarid == CDI_UNDEFID)
4159 if (strIsEqual(ncvar->name, ncdims[ncdimid].name))
4160 {
4161 ncdims[ncdimid].ncvarid = ncvarid;
4162 ncvar->varStatus = CoordVar;
4163 }
4164 }
4165 }
4166 }
4167
4168 // find time vars
4169 find_time_vars(nvars, ncvars, ncdims, timedimid, streamptr, &time_has_units, &time_has_bounds, &time_climatology);
4170
4171 leadtime_id = find_leadtime(nvars, ncvars, timedimid);
4172 if ( leadtime_id != CDI_UNDEFID ) ncvars[leadtime_id].varStatus = CoordVar;
4173
4174 // check ncvars
4175 timedimid = cdf_check_variables(streamptr, nvars, ncvars, ntsteps, timedimid);
4176
4177 // verify coordinates vars - first scan (dimname == varname)
4178 bool isHybridCF = false;
4179 verify_coordinates_vars_1(fileID, ndims, ncdims, ncvars, timedimid, &isHybridCF);
4180
4181 // verify coordinates vars - second scan (all other variables)
4182 verify_coordinates_vars_2(streamptr, nvars, ncvars);
4183
4184 if ( CDI_Debug ) cdf_print_vars(ncvars, nvars, "verify_coordinate_vars");
4185
4186 if ( ucla_les ) cdf_set_ucla_dimtype(ndims, ncdims, ncvars);
4187
4188 /*
4189 for (int ncdimid = 0; ncdimid < ndims; ncdimid++)
4190 {
4191 int ncvarid = ncdims[ncdimid].ncvarid;
4192 if ( ncvarid != -1 )
4193 {
4194 printf("coord var %d %s %s\n", ncvarid, ncvar->name, ncvar->units);
4195 if ( ncdims[ncdimid].dimtype == X_AXIS )
4196 printf("coord var %d %s is x dim\n", ncvarid, ncvar->name);
4197 if ( ncdims[ncdimid].dimtype == Y_AXIS )
4198 printf("coord var %d %s is y dim\n", ncvarid, ncvar->name);
4199 if ( ncdims[ncdimid].dimtype == Z_AXIS )
4200 printf("coord var %d %s is z dim\n", ncvarid, ncvar->name);
4201 if ( ncdims[ncdimid].dimtype == T_AXIS )
4202 printf("coord var %d %s is t dim\n", ncvarid, ncvar->name);
4203
4204 if ( ncvar->isLon )
4205 printf("coord var %d %s is lon\n", ncvarid, ncvar->name);
4206 if ( ncvar->isLat )
4207 printf("coord var %d %s is lat\n", ncvarid, ncvar->name);
4208 if ( ncvar->isZaxis )
4209 printf("coord var %d %s is lev\n", ncvarid, ncvar->name);
4210 }
4211 }
4212 */
4213
4214 // Set coordinates varids (att: associate)
4215 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
4216 {
4217 ncvar_t *ncvar = &ncvars[ncvarid];
4218 if (ncvar->varStatus == DataVar && ncvar->ncoordvars)
4219 {
4220 for (int i = 0; i < ncvar->ncoordvars; i++)
4221 {
4222 const int ncvaridc = ncvar->coordvarids[i];
4223 if (ncvaridc != CDI_UNDEFID)
4224 {
4225 ncvar_t *ncvarc = &ncvars[ncvaridc];
4226 if (ncvarc->isLon || ncvarc->isXaxis) ncvar->xvarid = ncvaridc;
4227 else if (ncvarc->isLat || ncvarc->isYaxis) ncvar->yvarid = ncvaridc;
4228 else if (ncvarc->isZaxis) ncvar->zvarid = ncvaridc;
4229 else if (ncvarc->isTaxis) ncvar->tvarid = ncvaridc;
4230 else if (ncvarc->isCharAxis) ncvar->cvarids[i] = ncvaridc;
4231 else if (ncvarc->printWarning)
4232 {
4233 Warning("Coordinates variable %s can't be assigned!", ncvarc->name);
4234 ncvarc->printWarning = false;
4235 }
4236 }
4237 }
4238 }
4239 }
4240
4241 // set dim type
4242 cdf_set_dimtype(nvars, ncvars, ncdims);
4243
4244 // read ECHAM VCT if present
4245 size_t vctsize = 0;
4246 double *vct = NULL;
4247 if (!isHybridCF) read_vct_echam(fileID, nvars, ncvars, ncdims, &vct, &vctsize);
4248
4249
4250 if (CDI_Debug) cdf_print_vars(ncvars, nvars, "cdf_define_all_grids");
4251
4252 // define all grids
4253 gridInfo.timedimid = timedimid;
4254 int status = cdf_define_all_grids(streamptr, streamptr->ncgrid, vlistID, ncdims, nvars, ncvars, &gridInfo);
4255 if (status < 0) return status;
4256
4257 // define all zaxes
4258 status = cdf_define_all_zaxes(streamptr, vlistID, ncdims, nvars, ncvars, vctsize, vct, uuidOfVGrid);
4259 if (vct) Free(vct);
4260 if (status < 0) return status;
4261
4262 // verify vars
4263 cdfVerifyVars(nvars, ncvars, ncdims);
4264
4265 // select vars
4266 int nvars_data = 0;
4267 for (int ncvarid = 0; ncvarid < nvars; ncvarid++)
4268 if (ncvars[ncvarid].varStatus == DataVar) nvars_data++;
4269
4270 if (CDI_Debug) Message("time varid = %d", streamptr->basetime.ncvarid);
4271 if (CDI_Debug) Message("ntsteps = %zu", ntsteps);
4272 if (CDI_Debug) Message("nvars_data = %d", nvars_data);
4273
4274 if ( nvars_data == 0 )
4275 {
4276 streamptr->ntsteps = 0;
4277 Warning("No data arrays found!");
4278 return CDI_EUFSTRUCT;
4279 }
4280
4281 if ( ntsteps == 0 && streamptr->basetime.ncdimid == CDI_UNDEFID && streamptr->basetime.ncvarid != CDI_UNDEFID )
4282 ntsteps = 1;
4283
4284 streamptr->ntsteps = (long)ntsteps;
4285
4286 // define all data variables
4287 cdf_define_all_vars(streamptr, vlistID, instID, modelID, nvars_data, nvars, ncvars, ncdims);
4288
4289
4290 cdiCreateTimesteps(streamptr);
4291
4292 // time varID
4293 int nctimevarid = streamptr->basetime.ncvarid;
4294
4295 if (nctimevarid != CDI_UNDEFID && (!time_has_units || streamptr->basetime.lwrf)) ncvars[nctimevarid].units[0] = 0;
4296 if (nctimevarid != CDI_UNDEFID && time_has_units) streamptr->basetime.lunits = true;
4297
4298 if ( time_has_units )
4299 {
4300 taxis_t *taxis = &streamptr->tsteps[0].taxis;
4301
4302 if ( setBaseTime(ncvars[nctimevarid].units, taxis) == 1 )
4303 {
4304 nctimevarid = CDI_UNDEFID;
4305 streamptr->basetime.ncvarid = CDI_UNDEFID;
4306 streamptr->basetime.lunits = false;
4307 }
4308
4309 if ( leadtime_id != CDI_UNDEFID && taxis->type == TAXIS_RELATIVE )
4310 {
4311 streamptr->basetime.leadtimeid = leadtime_id;
4312 taxis->type = TAXIS_FORECAST;
4313
4314 int timeunit = -1;
4315 if ( ncvars[leadtime_id].units[0] != 0 ) timeunit = scanTimeUnit(ncvars[leadtime_id].units);
4316 if ( timeunit == -1 ) timeunit = taxis->unit;
4317 taxis->fc_unit = timeunit;
4318
4319 setForecastTime(fcreftime, taxis);
4320 }
4321 }
4322
4323 if ( time_has_bounds )
4324 {
4325 streamptr->tsteps[0].taxis.has_bounds = true;
4326 if ( time_climatology ) streamptr->tsteps[0].taxis.climatology = true;
4327 }
4328
4329 if ( nctimevarid != CDI_UNDEFID )
4330 {
4331 taxis_t *taxis = &streamptr->tsteps[0].taxis;
4332 ptaxisDefName(taxis, ncvars[nctimevarid].name);
4333 if (ncvars[nctimevarid].longname[0]) ptaxisDefLongname(taxis, ncvars[nctimevarid].longname);
4334 if (ncvars[nctimevarid].units[0]) ptaxisDefUnits(taxis, ncvars[nctimevarid].units);
4335
4336 int xtype = ncvars[nctimevarid].xtype;
4337 int datatype = (xtype == NC_INT) ? CDI_DATATYPE_INT32 : ((xtype == NC_FLOAT) ? CDI_DATATYPE_FLT32 : CDI_DATATYPE_FLT64);
4338 ptaxisDefDatatype(taxis, datatype);
4339 }
4340
4341 int calendar = CDI_UNDEFID;
4342 if (nctimevarid != CDI_UNDEFID)
4343 if (ncvars[nctimevarid].hasCalendar)
4344 {
4345 char attstring[1024];
4346 cdfGetAttText(fileID, nctimevarid, "calendar", sizeof(attstring), attstring);
4347 strToLower(attstring);
4348 calendar = attribute_to_calendar(attstring);
4349 }
4350
4351 int taxisID;
4352 if ( streamptr->tsteps[0].taxis.type == TAXIS_FORECAST )
4353 {
4354 taxisID = taxisCreate(TAXIS_FORECAST);
4355 }
4356 else if ( streamptr->tsteps[0].taxis.type == TAXIS_RELATIVE )
4357 {
4358 taxisID = taxisCreate(TAXIS_RELATIVE);
4359 }
4360 else
4361 {
4362 taxisID = taxisCreate(TAXIS_ABSOLUTE);
4363 if ( !time_has_units )
4364 {
4365 taxisDefTunit(taxisID, TUNIT_DAY);
4366 streamptr->tsteps[0].taxis.unit = TUNIT_DAY;
4367 }
4368 }
4369
4370
4371 if (calendar == CDI_UNDEFID && streamptr->tsteps[0].taxis.type != TAXIS_ABSOLUTE) calendar = CALENDAR_STANDARD;
4372
4373 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 5)
4374 #pragma GCC diagnostic push
4375 #pragma GCC diagnostic warning "-Wstrict-overflow"
4376 #endif
4377 if ( calendar != CDI_UNDEFID )
4378 {
4379 taxis_t *taxis = &streamptr->tsteps[0].taxis;
4380 taxis->calendar = calendar;
4381 taxisDefCalendar(taxisID, calendar);
4382 }
4383 #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 5)
4384 #pragma GCC diagnostic pop
4385 #endif
4386
4387 vlistDefTaxis(vlistID, taxisID);
4388
4389 streamptr->curTsID = 0;
4390 streamptr->rtsteps = 1;
4391
4392 (void) cdfInqTimestep(streamptr, 0);
4393
4394 cdfCreateRecords(streamptr, 0);
4395
4396 // free ncdims
4397 if (ncdims) Free(ncdims);
4398
4399 // free ncvars
4400 if (ncvars) Free(ncvars);
4401
4402 return 0;
4403 }
4404
4405 static
wrf_read_timestep(int fileID,int nctimevarid,int tsID,taxis_t * taxis)4406 void wrf_read_timestep(int fileID, int nctimevarid, int tsID, taxis_t *taxis)
4407 {
4408 size_t start[2], count[2];
4409 char stvalue[128];
4410 start[0] = (size_t) tsID; start[1] = 0;
4411 count[0] = 1; count[1] = 19;
4412 stvalue[0] = 0;
4413 cdf_get_vara_text(fileID, nctimevarid, start, count, stvalue);
4414 stvalue[19] = 0;
4415 {
4416 int year = 1, month = 1, day = 1 , hour = 0, minute = 0, second = 0;
4417 if ( strlen(stvalue) == 19 )
4418 sscanf(stvalue, "%d-%d-%d_%d:%d:%d", &year, &month, &day, &hour, &minute, &second);
4419 taxis->vdate = cdiEncodeDate(year, month, day);
4420 taxis->vtime = cdiEncodeTime(hour, minute, second);
4421 taxis->type = TAXIS_ABSOLUTE;
4422 }
4423 }
4424
4425 static
get_timevalue(int fileID,int nctimevarid,int tsID,timecache_t * tcache)4426 double get_timevalue(int fileID, int nctimevarid, int tsID, timecache_t *tcache)
4427 {
4428 double timevalue = 0;
4429
4430 if ( tcache )
4431 {
4432 if ( tcache->size == 0 || (tsID < (int)tcache->startid || tsID > (int)(tcache->startid+tcache->size-1)) )
4433 {
4434 size_t maxvals = MAX_TIMECACHE_SIZE;
4435 if (maxvals > tcache->maxvals) maxvals = tcache->maxvals;
4436 tcache->startid = tsID;
4437 //tcache->startid = (tsID/MAX_TIMECACHE_SIZE)*MAX_TIMECACHE_SIZE;
4438 //if ( (tcache->startid + maxvals) > tcache->maxvals ) maxvals = (tcache->maxvals)%MAX_TIMECACHE_SIZE;
4439 tcache->size = maxvals;
4440 //fprintf(stderr, "fill time cache: %d %d %d %d %d\n", tcache->maxvals, tsID, tcache->startid, tcache->startid+maxvals-1, maxvals);
4441 cdf_get_vara_double(fileID, nctimevarid, &tcache->startid, &maxvals, tcache->cache);
4442 }
4443
4444 //timevalue = tcache->cache[tsID%MAX_TIMECACHE_SIZE];
4445 timevalue = tcache->cache[tsID];
4446 }
4447 else
4448 {
4449 size_t index = (size_t) tsID;
4450 cdf_get_var1_double(fileID, nctimevarid, &index, &timevalue);
4451 }
4452
4453 if ( timevalue >= NC_FILL_DOUBLE || timevalue < -NC_FILL_DOUBLE ) timevalue = 0;
4454
4455 return timevalue;
4456 }
4457
4458
cdfInqTimestep(stream_t * streamptr,int tsID)4459 int cdfInqTimestep(stream_t * streamptr, int tsID)
4460 {
4461 if ( CDI_Debug ) Message("streamID = %d tsID = %d", streamptr->self, tsID);
4462
4463 if ( tsID < 0 ) Error("unexpected tsID = %d", tsID);
4464
4465 if ( tsID < streamptr->ntsteps && streamptr->ntsteps > 0 )
4466 {
4467 cdfCreateRecords(streamptr, tsID);
4468
4469 taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
4470 if ( tsID > 0 )
4471 ptaxisCopy(taxis, &streamptr->tsteps[0].taxis);
4472
4473 double timevalue = tsID;
4474
4475 int nctimevarid = streamptr->basetime.ncvarid;
4476 if (nctimevarid != CDI_UNDEFID && streamptr->basetime.lunits)
4477 {
4478 int fileID = streamptr->fileID;
4479 size_t index = (size_t)tsID;
4480
4481 if ( streamptr->basetime.lwrf )
4482 {
4483 wrf_read_timestep(fileID, nctimevarid, tsID, taxis);
4484 }
4485 else
4486 {
4487 #ifdef USE_TIMECACHE
4488 if ( streamptr->basetime.timevar_cache == NULL )
4489 {
4490 streamptr->basetime.timevar_cache = (timecache_t *) Malloc(sizeof(timecache_t));
4491 streamptr->basetime.timevar_cache->size = 0;
4492 streamptr->basetime.timevar_cache->maxvals = streamptr->ntsteps;
4493 }
4494 #endif
4495 timevalue = get_timevalue(fileID, nctimevarid, tsID, streamptr->basetime.timevar_cache);
4496 cdiDecodeTimeval(timevalue, taxis, &taxis->vdate, &taxis->vtime);
4497 }
4498
4499 int nctimeboundsid = streamptr->basetime.ncvarboundsid;
4500 if ( nctimeboundsid != CDI_UNDEFID )
4501 {
4502 size_t start[2], count[2];
4503 start[0] = index; count[0] = 1; start[1] = 0; count[1] = 1;
4504 cdf_get_vara_double(fileID, nctimeboundsid, start, count, &timevalue);
4505 if ( timevalue >= NC_FILL_DOUBLE || timevalue < -NC_FILL_DOUBLE ) timevalue = 0;
4506
4507 cdiDecodeTimeval(timevalue, taxis, &taxis->vdate_lb, &taxis->vtime_lb);
4508
4509 start[0] = index; count[0] = 1; start[1] = 1; count[1] = 1;
4510 cdf_get_vara_double(fileID, nctimeboundsid, start, count, &timevalue);
4511 if ( timevalue >= NC_FILL_DOUBLE || timevalue < -NC_FILL_DOUBLE ) timevalue = 0;
4512
4513 cdiDecodeTimeval(timevalue, taxis, &taxis->vdate_ub, &taxis->vtime_ub);
4514 }
4515
4516 int leadtimeid = streamptr->basetime.leadtimeid;
4517 if ( leadtimeid != CDI_UNDEFID )
4518 {
4519 timevalue = get_timevalue(fileID, leadtimeid, tsID, NULL);
4520 cdiSetForecastPeriod(timevalue, taxis);
4521 }
4522 }
4523 }
4524
4525 streamptr->curTsID = tsID;
4526 long nrecs = streamptr->tsteps[tsID].nrecs;
4527
4528 return (int) nrecs;
4529 }
4530
4531 #endif
4532 /*
4533 * Local Variables:
4534 * c-file-style: "Java"
4535 * c-basic-offset: 2
4536 * indent-tabs-mode: nil
4537 * show-trailing-whitespace: t
4538 * require-trailing-newline: t
4539 * End:
4540 */
4541