1 /* Authored by Don Hooper */
2 /*
3 * Include ./configure's header file
4 */
5 #ifdef HAVE_CONFIG_H
6 #include <config.h>
7 #endif
8
9 #if USESDF == 1
10 #define Success 1
11 #define Failure 0
12
13 /* Functions related to Self-Describing File (SDF) access. -hoop, 95/07/10 */
14 #include <stdlib.h>
15 #include <unistd.h>
16 #include <stdio.h>
17 #include <string.h>
18 #include <math.h>
19 #include "grads.h"
20
21
22 /* expose Mike Fiorino's global struct to these for 365 day calendars */
23 extern struct gamfcmn mfcmn ;
24
25 /* expose descriptor global from gaddes.c here so can share w/routines there */
26 extern FILE *descr ;
27
28 /* kk 020624 ---s */
29 char *gxgnam(char *) ; /* This is also in gx.h */
30 /* kk 020624 ---e */
31 static char pout[512]; /* Build error msgs here */
32
33 /* Read the metadata from a Self-Describing File and fill in the information */
34 /* into a gafile structure. The gafile structure should be allocated and */
35 /* must be initialized by getpfi. If this routine returns an error, release */
36 /* the pfi structure and allocated storage via frepfi. */
37
38 #define XINDEX 0
39 #define YINDEX 1
40 #define ZINDEX 2
41 #define TINDEX 3
42
43
44 int utISinit = 0 ;
45
46 /* For STNDALN, routines included are gadxdf and three others contained therein
47 These are used for parsing the descriptors just to see if GrADS will open them */
48 #ifndef STNDALN
49
50 int
gadsdf(char * name,struct gafile * pfi,char * template,int ntsteps,GASDFPARMS parms)51 gadsdf(char *name, struct gafile *pfi, char *template, int ntsteps, GASDFPARMS parms) {
52 struct gavar *pvar;
53 struct dt tdef, tdefi, dt1, dt2;
54 struct gaindx *pindx;
55 float *vals;
56 int size, rc, len, swpflg, cnt, boole, boole2 ;
57 char *ch, *dn, *pos;
58 char *utname;
59 int flgs[8], cflg, i, j, err, hdrb, trlb, mflflg;
60 int mcnt, maxlv, maxct, inx, xedni;
61 int levs, acum, fpos, recacm, indx, numdvars ;
62 int iyr, imo, idy, ihr, imn, isec, ispress, idim, isDatavar ;
63 int xdimension_id, ydimension_id, zdimension_id, tdimension_id ;
64 short *tempsarry ;
65 float temp, v1, v2, *tempfarry, *tempfarry2 ;
66 FILE *mfile;
67 utUnit timeunit ;
68 long start[1] = {0}, count[1] = {1}, tdimsiz ;
69 double *time1ptr, *time2ptr, time1, time2, dsec, incrfactor, *tempdarry ;
70 VAR_INFO *Xcoord=NULL, *Ycoord=NULL, *Zcoord=NULL, *Tcoord=NULL ;
71 VAR_INFO *lclvar, *DataVar=NULL, *SaveDvar=NULL ;
72 struct attrib_list *positive_attr, *globttl_attr, *missval_attr ;
73 struct attrib_list *deltat_attr, *title_vattr, *timeunits_attr ;
74 struct attrib_list *calendar_attr ;
75 struct attrib_list *find_att(struct attrib_list *first_att, char *attname) ;
76 int read_io_std(char *path, IO_STD *sdf_ptr) ; /* grab metadata */
77 int isdvar(VAR_INFO *var, VAR_INFO *Xcoord, VAR_INFO *Ycoord,
78 VAR_INFO *Zcoord, VAR_INFO *Tcoord, int xdimension_id,
79 int ydimension_id, int zdimension_id, int tdimension_id) ;
80 void free_io_std(IO_STD **sdf_ptr) ;
81 int close_netcdf(int ncid) ;
82 /* these routines try to find the right var struct in the */
83 /* linked list in the IO_STD structure */
84 int findX(IO_STD *sdf_ptr, VAR_INFO **coord) ;
85 int findY(IO_STD *sdf_ptr, VAR_INFO **coord) ;
86 int findZ(IO_STD *sdf_ptr, VAR_INFO **coord, int *ispressptr) ;
87 int findT(IO_STD *sdf_ptr, VAR_INFO **coord) ;
88 int sdfdeflev(struct gafile *pfi, VAR_INFO *coord, int GrADSdimnum, int revflag) ;
89 int trydeflin(struct gafile *pfi, VAR_INFO *coord, int GrADSdimnum, int isX, int revflag) ;
90 int convtype(void *srcptr, nc_type srctype, void *desptr, nc_type destype,
91 int srcindex, int desindex) ;
92 int decode_ud_time(utUnit *unit, double time_val, int *yr, int *mo,
93 int *da, int *hr, int *min, double *sec) ;
94 int decode_standard_time(double time_val, int *year, int *month, int *day,
95 int *hour, int *minn, double *sec) ;
96 int find_dim(IO_STD *std_ptr, char *name) ;
97 int find_dimix(int dimids[], int ndims, int dim) ;
98 VAR_INFO *find_var(IO_STD *std_ptr, char *name) ;
99 int findmatch(char *var_name, char **varnamelist, int namecount, int *inxptr) ;
100 int copy_io_std(IO_STD **std_ptr2, IO_STD *std_ptr1) ;
101
102 /* Enable griping, disable aborting, from within NetCDF library */
103 ncopts = NC_VERBOSE ;
104 if (!utISinit) {
105 utname = gxgnam("udunits.dat") ;
106 if (utname != NULL) {
107 if (utInit(utname)) {
108 gaprnt(0, "UDUNITS package initialization failure.\n") ;
109 return(1) ;
110 }
111 }
112 utISinit = 1 ;
113 }
114
115 /* Grab the metadata */
116 if (!read_io_std(name, pfi->sdf_ptr)) {
117 gaprnt(0, "\nCouldn't ingest SDF metadata.\n") ;
118 #if USEHDF == 0
119 gaprnt(0, "If this was an HDF-SDS file, try gradshdf.\n") ;
120 #endif
121 return(1) ;
122 }
123
124 if (!parms.isxdf) {
125 pfi->calendar = 0 ; /* 365 day kind not available under COARDS */
126 /* But we're going to fake something anyway. -Hoop, 2001/05/31 */
127 /* The code will be in the time coordinate setup section. */
128 }
129 if (parms.xsrch) {
130 if (!findX(pfi->sdf_ptr, &Xcoord)) {
131 gaprnt(0, "SDF file has no discernable X coordinate.\n") ;
132 return(1) ;
133 }
134 } else if (parms.xsetup) {
135 Xcoord = find_var(pfi->sdf_ptr, parms.xdimname) ;
136 if (!Xcoord) {
137 gaprnt(0, "XDF file has no coordinate variable for X dimension.\n") ;
138 return(1) ;
139 }
140 }
141 if (parms.ysrch) {
142 if (!findY(pfi->sdf_ptr, &Ycoord)) {
143 gaprnt(0, "SDF file has no discernable Y coordinate.\n") ;
144 return(1) ;
145 }
146 } else if (parms.ysetup) {
147 Ycoord = find_var(pfi->sdf_ptr, parms.ydimname) ;
148 if (!Ycoord) {
149 gaprnt(0, "XDF file has no coordinate variable for Y dimension.\n") ;
150 return(1) ;
151 }
152 }
153 if (parms.tsrch) {
154 if (!findT(pfi->sdf_ptr, &Tcoord)) {
155 gaprnt(0, "SDF file has no discernable time coordinate.\n") ;
156 }
157 } else if (parms.tsetup) {
158 if (!parms.tdimname) {
159 gaprnt(0, "XDF file has no time dimension, but TDEF is incomplete.\n") ;
160 gaprnt(0, "Failure parsing XDF-style data descriptor file (DDF).\n") ;
161 return(1) ;
162 }
163 Tcoord = find_var(pfi->sdf_ptr, parms.tdimname) ;
164 if (!Tcoord) {
165 gaprnt(0, "XDF file has no coordinate variable for T dimension.\n") ;
166 return(1) ;
167 }
168 }
169 if (parms.zsrch) {
170 (void) findZ(pfi->sdf_ptr, &Zcoord, &ispress) ;
171
172 } else if (parms.zsetup) {
173 Zcoord = find_var(pfi->sdf_ptr, parms.zdimname) ;
174 if (!Zcoord) {
175 gaprnt(0, "XDF file has no coordinate variable for Z dimension.\n") ;
176 return(1) ;
177 }
178 }
179 #ifdef TESTING
180 if (Xcoord) fprintf(stderr, "Using %s as X variable.\n", Xcoord->varnam) ;
181 if (Ycoord) fprintf(stderr, "Using %s as Y variable.\n", Ycoord->varnam) ;
182 if (Tcoord) fprintf(stderr, "Using %s as T variable.\n", Tcoord->varnam) ;
183 if (Zcoord)
184 fprintf(stderr, "Using %s as Z variable.\n", Zcoord->varnam) ;
185 else
186 fprintf(stderr, "Didn't find a Z variable.\n") ;
187
188 #endif
189 if (parms.isxdf && (!parms.tsrch)) {
190 if (!(parms.tdimname)) {
191 tdimension_id = -1 ;
192 } else {
193 tdimension_id = find_dim(pfi->sdf_ptr, parms.tdimname) ;
194 }
195 } else {
196 if (Tcoord) {
197 tdimension_id = find_dim(pfi->sdf_ptr, Tcoord->varnam) ;
198 } else {
199 tdimension_id = -1 ;
200 }
201 }
202 if (parms.isxdf && (!parms.zsrch)) {
203 zdimension_id = find_dim(pfi->sdf_ptr, parms.zdimname) ;
204 } else {
205 if (Zcoord) {
206 zdimension_id = find_dim(pfi->sdf_ptr, Zcoord->varnam) ;
207 } else {
208 zdimension_id = -1 ;
209 }
210 }
211 if (parms.isxdf && (!parms.ysrch)) {
212 ydimension_id = find_dim(pfi->sdf_ptr, parms.ydimname) ;
213 } else {
214 if (Ycoord) {
215 ydimension_id = find_dim(pfi->sdf_ptr, Ycoord->varnam) ;
216 } else {
217 ydimension_id = -1 ;
218 }
219 }
220 if (parms.isxdf && (!parms.xsrch)) {
221 xdimension_id = find_dim(pfi->sdf_ptr, parms.xdimname) ;
222 } else {
223 if (Xcoord) {
224 xdimension_id = find_dim(pfi->sdf_ptr, Xcoord->varnam) ;
225 } else {
226 xdimension_id = -1 ;
227 }
228 }
229 if (parms.dvsrch) {
230 numdvars = 0 ;
231 SaveDvar = pfi->sdf_ptr->var ;
232 for (DataVar = pfi->sdf_ptr->var ; DataVar ; DataVar = DataVar->next) {
233 isDatavar = 0 ;
234 isDatavar = isdvar(DataVar, Xcoord, Ycoord, Zcoord, Tcoord,
235 xdimension_id, ydimension_id, zdimension_id,
236 tdimension_id) ;
237 if (isDatavar) {
238 ++numdvars ;
239 if (numdvars == 1) /* Use first Dvar found as master */ {
240 SaveDvar = DataVar ;
241 }
242 }
243 } /* for each data variable */
244 if (numdvars == 0) /* oh dear... */ {
245 gaprnt(0,"SDF file doesn't appear to have any non-coordinate variables. Yikes!\n") ;
246 return(1) ;
247 } else {
248 parms.dvcount = numdvars ;
249 }
250 } else {
251 numdvars = parms.dvcount ;
252 for (DataVar = pfi->sdf_ptr->var ; DataVar ; DataVar = DataVar->next) {
253 if (!strcmp(parms.dvncnames[0], DataVar->varnam)) { /* Use first Dvar as master */
254 SaveDvar = DataVar ;
255 break ;
256 }
257 }
258 if (!SaveDvar) {
259 sprintf(pout,"Couldn't find \"master\" data variable %s in SDF file.\n", parms.dvncnames[0]);
260 gaprnt(0,pout);
261 return(1);
262 }
263 }
264 if (!parms.isxdf) {
265 getwrd(pfi->dnam, name, 4095) ;
266 }
267 getwrd(pfi->name, name, 4095) ; /* Is this template or first nc name? */
268 /* I dunno about the next 5 lines... -hoop */
269 hdrb = 0;
270 trlb = 0;
271 mflflg = 0; /* no map file to open */
272 pfi->mfile = NULL;
273 mcnt = -1;
274
275 if (parms.xsetup) {
276 pfi->dnum[XINDEX] = pfi->sdf_ptr->dimsiz[Xcoord->vardimid[0]] ;
277 if (!read_one_dimension(pfi->sdf_ptr, Xcoord, Xcoord->vardimid[0])) {
278 gaprnt(0, "Failure checking X coordinate of SDF file.\n") ;
279 return(1) ;
280 }
281 rc = trydeflin(pfi, Xcoord, XINDEX, 1, 0) ;
282 if (rc < 0) {
283 gaprnt(0, "Error evaluating X coordinate of SDF file.\n") ;
284 return(1) ;
285 }
286 if (rc == 0) /* not linear; define levels */ {
287 if (sdfdeflev(pfi, Xcoord, XINDEX, 0)) {
288 gaprnt(0, "Failure defining X levels for SDF file.\n") ;
289 return(1) ;
290 }
291 #ifdef TESTING2
292 } else {
293 fprintf(stderr, "Using Linear X dim. on SDF file.\n") ;
294 #endif
295 }
296 } /* if doing X setup */
297 if (parms.ysetup) {
298 pfi->dnum[YINDEX] = pfi->sdf_ptr->dimsiz[Ycoord->vardimid[0]] ;
299 if ((!read_one_dimension(pfi->sdf_ptr, Ycoord, Ycoord->vardimid[0])) ||
300 (!(Ycoord->data))) {
301 gaprnt(0, "Failure checking Y coordinate of SDF file.\n") ;
302 return(1) ;
303 }
304 /* Next few lines of code only for Y; not analogous for X, tho' should be */
305 tempfarry = (float *) malloc(2 * sizeof(float)) ;
306 if (!convtype(Ycoord->data, Ycoord->vartype, tempfarry, NC_FLOAT, 0, 0)) {
307 gaprnt(0, "Failure interpreting first Y coordinate value of SDF file.\n") ;
308 return(1) ;
309 }
310 if ( (pfi->dnum[YINDEX]) > 1) {
311 if (!convtype(Ycoord->data, Ycoord->vartype, tempfarry, NC_FLOAT, 1, 1)) {
312 gaprnt(0, "Failure interpreting second Y coordinate value of SDF file.\n") ;
313 return(1) ;
314 }
315 if (tempfarry[1] < tempfarry[0]) {
316 pfi->yrflg = 1 ;
317 } else {
318 pfi->yrflg = 0 ;
319 }
320 } else {
321 pfi->yrflg = 0 ;
322 }
323 /* End Y-only code */
324 rc = trydeflin(pfi, Ycoord, YINDEX, 0, pfi->yrflg) ;
325 if (rc < 0) {
326 gaprnt(0, "Error evaluating Y coordinate of SDF file.\n") ;
327 return(1) ;
328 }
329 if (rc == 0) /* not linear; define levels */ {
330 if (sdfdeflev(pfi, Ycoord, YINDEX, pfi->yrflg)) {
331 gaprnt(0, "Failure defining levels for Y coordinate for SDF file.\n") ;
332 return(1) ;
333 }
334 #ifdef TESTING2
335 } else {
336 fprintf(stderr, "Using Linear Y dim. on SDF file.\n") ;
337 #endif
338 }
339 } /* if doing y setup */
340 if (parms.zsetup) {
341 if (Zcoord) {
342 pfi->zrflg = 0 ;
343 }
344 } /* if doing z setup */
345 if (parms.needtitle) {
346 globttl_attr = find_att(pfi->sdf_ptr->first_gattr, "title") ;
347 if (globttl_attr) {
348 if ((globttl_attr->len) > 510) {
349 getstr(pfi->title, (char *)(globttl_attr->data), 511) ;
350 } else {
351 getstr(pfi->title, (char *)(globttl_attr->data), globttl_attr->len) ;
352 }
353 }
354 }
355 if (parms.needundef) {
356 /* the following is a kludge because, by COARDS conventions, missing_value is */
357 /* a per-variable attribute, while GrADS expects it to be global to the data */
358 /* file. So, the first non-coordinate variable is used to set a pseudo- */
359 /* global "undefined" value */
360 missval_attr = find_att(SaveDvar->first_vattr, "missing_value") ;
361 if (!missval_attr) {
362 missval_attr = find_att(SaveDvar->first_vattr, "_FillValue") ;
363 }
364 if (!missval_attr) {
365 pfi->undef = FILL_FLOAT ; /* default _FillValue */
366 } else {
367 if (missval_attr->type == NC_SHORT) /* KLUDGE! */ {
368 pfi->undef = (-1.0) * FILL_FLOAT ; /* CDC default w/packed data */
369 } else /* must be NC_FLOAT or NC_DOUBLE, right? I hope so... */ {
370 if (missval_attr->type == NC_DOUBLE) {
371 tempdarry = (double *) missval_attr->data ;
372 pfi->undef = (float) tempdarry[0] ;
373 } else {
374 tempfarry = (float *) missval_attr->data ;
375 pfi->undef = tempfarry[0] ;
376 }
377 }
378 } /* there was a missing value of some name */
379 } /* if undef needs to be set */ /* next if MAY need to be moved before this line */
380 if (SETMISS) {
381 /* Gads... Will these statements work OK when pfi->undef is -infinity? */
382 pfi->ulow = fabs(pfi->undef/EPSILON);
383 pfi->uhi = pfi->undef + pfi->ulow;
384 pfi->ulow = pfi->undef - pfi->ulow;
385 }
386
387 if (parms.zsetup) {
388 if (!Zcoord) {
389 pfi->dnum[ZINDEX] = 1 ;
390 /* These lines adapted from deflin routine */
391 tempfarry = (float *) malloc(sizeof(float) * 6) ;
392 if (!tempfarry) {
393 gaprnt(0, "Error evaluating fixed Z coordinate in SDF file.\n") ;
394 return(1) ;
395 }
396 tempfarry[0] = 0.0 ;
397 tempfarry[1] = 0.0 ; /* can do better if interp. level_desc attr. */
398 tempfarry[2] = -999.9 ;
399 tempfarry[3] = 1.0 ;
400 tempfarry[4] = 1.0 ;
401 tempfarry[5] = -999.9 ;
402 pfi->grvals[ZINDEX] = tempfarry ;
403 pfi->abvals[ZINDEX] = tempfarry + 3 ;
404 pfi->ab2gr[ZINDEX] = liconv ;
405 pfi->gr2ab[ZINDEX] = liconv ;
406 pfi->linear[ZINDEX] = 1 ;
407 /* end of lines adapted from deflin routine */
408 } else /* not a dummy; an actual Z coordinate variable */ {
409 pfi->dnum[ZINDEX] = pfi->sdf_ptr->dimsiz[Zcoord->vardimid[0]] ;
410 if (!read_one_dimension(pfi->sdf_ptr, Zcoord, Zcoord->vardimid[0])) {
411 gaprnt(0, "Failure checking Z coordinate of SDF file.\n") ;
412 return(1) ;
413 }
414 tempfarry = (float *) malloc(2 * sizeof(float)) ;
415 if (!convtype(Zcoord->data, Zcoord->vartype, tempfarry, NC_FLOAT, 0, 0)) {
416 gaprnt(0, "Failure interpreting first Z coordinate value of SDF file.\n") ;
417 return(1) ;
418 }
419 if (pfi->dnum[ZINDEX] > 1) {
420 if (!convtype(Zcoord->data, Zcoord->vartype, tempfarry,
421 NC_FLOAT, 1, 1)) {
422 gaprnt(0, "Failure interpreting second Z coordinate value of SDF file.\n") ;
423 return(1) ;
424 }
425 positive_attr = find_att(Zcoord->first_vattr, "positive") ;
426 if (positive_attr) {
427 if (!strncmp("down", (char *) (positive_attr->data), 4)) {
428 /* positive:down */
429 if (tempfarry[1] > tempfarry[0]) /* ascending Z values */ {
430 pfi->zrflg = 1 ;
431 } else {
432 pfi->zrflg = 0 ;
433 }
434 } else /* positive:up */ {
435 if (tempfarry[1] > tempfarry[0]) /* ascending Z values */ {
436 pfi->zrflg = 0 ;
437 } else {
438 pfi->zrflg = 1 ;
439 }
440 }
441 }
442 if (sdfdeflev(pfi, Zcoord, ZINDEX, pfi->zrflg)) {
443 #ifdef TESTING
444 fprintf(stderr, "sdfdeflev failed on Zcoord.\n") ;
445 #endif
446 gaprnt(0, "Couldn't define Z-axis usinng LEVELS method.\n") ;
447 return(1) ;
448 } /* sdfdeflev */
449 } else /* 1 Z value, so it linearly */ {
450 /* These lines adapted from deflin routine */
451 tempfarry2 = (float *) malloc(sizeof(float) * 6) ;
452 if (!tempfarry2) {
453 gaprnt(0, "Error evaluating solo Z coordinate in SDF file.\n") ;
454 return(1) ;
455 }
456 tempfarry2[0] = 1.0 ;
457 tempfarry2[1] = tempfarry[0] - 1.0 ;
458 tempfarry2[2] = -999.9 ;
459 tempfarry2[3] = 1.0 ;
460 tempfarry2[4] = 1.0 - tempfarry[0] ;
461 tempfarry2[5] = -999.9 ;
462 pfi->grvals[ZINDEX] = tempfarry2 ;
463 pfi->abvals[ZINDEX] = tempfarry2 + 3 ;
464 pfi->ab2gr[ZINDEX] = liconv ;
465 pfi->gr2ab[ZINDEX] = liconv ;
466 pfi->linear[ZINDEX] = 1 ;
467 /* end of lines adapted from deflin routine */
468 } /* one Z coord. value */
469 } /* setting up actual Z coordinate variable */
470 } /* if parms say setup Z coordinate */
471 if (parms.tsetup) {
472 /* Initial setting only; could be updated if templating */
473 if (!Tcoord) /* fake the time coord */ {
474 pfi->dnum[TINDEX] = 1 ;
475 tempfarry = (float *) malloc(sizeof(float) * 8) ;
476 if (!tempfarry) {
477 gaprnt(0, "Error finding storage to define time coordinate in SDF file.\n") ;
478 return(1) ;
479 }
480 tempfarry[7] = -999.9 ;
481 pfi->grvals[TINDEX] = tempfarry ;
482 pfi->abvals[TINDEX] = tempfarry ;
483 pfi->linear[TINDEX] = 1 ;
484 tempfarry[0] = 1.0 ;
485 tempfarry[1] = 1.0 ;
486 tempfarry[2] = 1.0 ;
487 tempfarry[3] = 0.0 ; /* initial hours */
488 tempfarry[4] = 0.0 ;
489 tempfarry[5] = 0.0 ; /* step in months */
490 tempfarry[6] = 1.0 ; /* step in minutes */
491 } else {
492 calendar_attr = find_att(Tcoord->first_vattr, CALENDAR) ;
493 if (calendar_attr) {
494 if (!strcasecmp(CAL365, (char *) (calendar_attr->data))) {
495 mfcmn.cal365 = pfi->calendar = 1;
496 }
497 if (!strcasecmp(ALTCAL365, (char *) (calendar_attr->data))) {
498 mfcmn.cal365 = pfi->calendar = 1;
499 }
500 }
501 pfi->dnum[TINDEX] = pfi->sdf_ptr->dimsiz[Tcoord->vardimid[0]] ;
502 timeunits_attr = find_att(Tcoord->first_vattr, "units") ;
503 if (!timeunits_attr) {
504 gaprnt(0, "Couldn't find units attribute for Time coordinate of SDF file.\n") ;
505 return(1) ;
506 }
507 #if USEHDF == 1
508 if (!read_variable_data(pfi->sdf_ptr->cdfid, Tcoord, start, count)) {
509 gaprnt(0, "Error reading initial time value in SDF file.\n") ;
510 return(1) ;
511 }
512 time1ptr = &time1 ;
513 if (!convtype(Tcoord->data, Tcoord->vartype, time1ptr, NC_DOUBLE, 0, 0)) {
514 gaprnt(0, "Error interpreting initial time value in SDF file.\n") ;
515 return(1) ;
516 }
517 if (pfi->dnum[TINDEX] > 1) /* deduce increment */ {
518 (start[0])++ ; /* set up to read second time value */
519 if (!read_variable_data(pfi->sdf_ptr->cdfid, Tcoord, start, count)) {
520 gaprnt(0, "Error reading second time value in SDF file.\n") ;
521 return(1) ;
522 }
523 (start[0])-- ; /* reset it in case we need it again */
524 time2ptr = &time2 ;
525 if (!convtype(Tcoord->data, Tcoord->vartype, time2ptr, NC_DOUBLE, 0, 0)) {
526 gaprnt(0, "Error interpreting second time value in SDF file.\n") ;
527 return(1) ;
528 }
529 } else {
530 time2 = time1 + 1 ; /* fake for single timestep files */
531 }
532 #else
533 if (pfi->dnum[TINDEX] > 1) /* deduce increment */ {
534 count[0] = 2 ; /* read both time vals at once for DODS bug */
535 if (!read_variable_data(pfi->sdf_ptr->cdfid, Tcoord, start, count)) {
536 gaprnt(0, "Error reading initial two time values in SDF file.\n") ;
537 return(1) ;
538 }
539 count[0] = 1 ; /* restore for other calls, if any */
540 time1ptr = &time1 ;
541 if (!convtype(Tcoord->data, Tcoord->vartype, time1ptr, NC_DOUBLE, 0, 0)) {
542 gaprnt(0, "Error interpreting initial time value in SDF file.\n") ;
543 return(1) ;
544 }
545 time2ptr = &time2 ;
546 if (!convtype(Tcoord->data, Tcoord->vartype, time2ptr, NC_DOUBLE, 1, 0)) {
547 gaprnt(0, "Error interpreting second time value in SDF file.\n") ;
548 return(1) ;
549 }
550 } else /* fake it when only 1 timestep exists */ {
551 if (!read_variable_data(pfi->sdf_ptr->cdfid, Tcoord, start, count)) {
552 gaprnt(0, "Error reading initial time value in SDF file.\n") ;
553 return(1) ;
554 }
555 time1ptr = &time1 ;
556 if (!convtype(Tcoord->data, Tcoord->vartype, time1ptr, NC_DOUBLE, 0, 0)) {
557 gaprnt(0, "Error interpreting initial time value in SDF file.\n") ;
558 return(1) ;
559 }
560 time2 = time1 + 1 ; /* a kludge, of course */
561 }
562 #endif
563 tempfarry = (float *) malloc(sizeof(float) * 8) ;
564 if (!tempfarry) {
565 gaprnt(0,
566 "Error finding storage to define time coordinate in SDF file.\n") ;
567 return(1) ;
568 }
569 tempfarry[7] = -999.9 ;
570 pfi->grvals[TINDEX] = tempfarry ;
571 pfi->abvals[TINDEX] = tempfarry ;
572 pfi->linear[TINDEX] = 1 ;
573 if ((timeunits_attr->len < 10) &&
574 (!strncmp((char *)timeunits_attr->data, "YYMMDDHH", 8))) /*GSFC time */{
575 tempfarry[0] = (float) (((int) time1) / 1000000) ;
576 tempfarry[1] = (float) ((((int) time1) / 10000) % 100) ;
577 tempfarry[2] = (float) ((((int) time1) / 100) % 100) ;
578 tempfarry[3] = (float) (((int) time1) % 100) ;
579 tempfarry[4] = 0.0 ;
580 if (pfi->dnum[TINDEX] > 1) /* deduce increment */ {
581 dt2.yr = (((int) time2) / 1000000) - (((int) time1) / 1000000) ;
582 dt2.mo = ((((int) time2) / 10000) % 100) -
583 ((((int) time1) / 10000) % 100) ;
584 dt2.dy = ((((int) time2) / 100) % 100) -
585 ((((int) time1) / 100) % 100) ;
586 dt2.hr = (((int) time2) % 100) - (((int) time1) % 100) ;
587 dt2.mn = 0 ;
588 if ((dt2.yr > 0) || (dt2.mo > 0)) {
589 tempfarry[5] = (dt2.yr * 12.0) + dt2.mo ;
590 tempfarry[6] = 0.0 ;
591 } else {
592 tempfarry[5] = 0.0 ;
593 tempfarry[6] = (dt2.dy * 1440.0) + (dt2.hr * 60.0) + dt2.mn ;
594 if (tempfarry[6] < 1.0) {
595 gaprnt(0, "Time unit in SDF file has too small an increment (min. 1 minute).\n") ;
596 return(1) ;
597 }
598 }
599 } else /* one time step exists */ {
600 tempfarry[5] = 0.0 ;
601 tempfarry[6] = 1.0 ; /* KLUDGE for one time-step files */
602 }
603 } else if (pfi->sdf_ptr->time_type == CDC) /* handle YYYYMMDDHHMMSS time */ {
604 if (!decode_standard_time(time1, &iyr, &imo, &idy, &ihr, &imn, &dsec)) {
605 gaprnt(0, "Error deciphering initial time value in SDF file.\n") ;
606 return(1) ;
607 }
608 if (iyr <= 0) iyr = 1 ;
609 if (imo <= 0) imo = 1 ;
610 if (idy <= 0) idy = 1 ;
611 tempfarry[0] = (float) iyr ;
612 tempfarry[1] = (float) imo ;
613 tempfarry[2] = (float) idy ;
614 tempfarry[3] = (float) ihr ;
615 tempfarry[4] = (float) imn ;
616 /* will cop out here and assume delta_t attribute exists, as divining the */
617 /* increment from the first two values is too messy to contemplate */
618 if (pfi->dnum[TINDEX] > 1) {
619 deltat_attr = find_att(Tcoord->first_vattr, "delta_t") ;
620 if (!deltat_attr) {
621 gaprnt(0, "Error in determining time increment in SDF file.\n") ;
622 return(1) ;
623 }
624 ch = (char *) deltat_attr->data ;
625 if (!decode_delta_t(ch, &dt2.yr, &dt2.mo, &dt2.dy, &dt2.hr, &dt2.mn, &isec)) {
626 gaprnt(0, "Error deciphering time increment in SDF file.\n") ;
627 return(1) ;
628 }
629 if ((dt2.yr > 0) || (dt2.mo > 0)) {
630 tempfarry[5] = (dt2.yr * 12.0) + dt2.mo ;
631 tempfarry[6] = 0.0 ;
632 } else {
633 tempfarry[5] = 0.0 ;
634 tempfarry[6] = (dt2.dy * 1440.0) + (dt2.hr * 60.0) + dt2.mn ;
635 if (tempfarry[6] < 1.0) {
636 gaprnt(0, "Time unit in SDF file has too small an increment (min. 1 minute).\n") ;
637 return(1) ;
638 }
639 }
640 /* Yet Another Kludge: any seconds value in SDF file is ignored as minute */
641 /* is the finest resolution in GrADS */
642 } else {
643 tempfarry[5] = 0.0 ;
644 tempfarry[6] = 1.0 ; /* KLUDGE for one time-step files */
645 }
646 } else /* assume udunits time */ {
647 char *time_units = (char *) 0 ;
648 int len_time_units ;
649 #define DFLTORIGIN " since 1-1-1 00:00:0.0"
650
651 if (!strncmp((char *) timeunits_attr->data, "common_year", 11)) {
652 mfcmn.cal365 = pfi->calendar = 1;
653 } /* that's the first step. Now, have add new matches */
654 if (!strstr((char *) timeunits_attr->data, " since ")) /* no origin */ {
655 len_time_units = (int) strlen((char *) timeunits_attr->data)
656 + (int) strlen(DFLTORIGIN) + 1 ;
657 time_units = (char *) malloc((size_t) len_time_units) ;
658 if (!time_units) {
659 gaprnt(0, "Malloc error generating time_units for SDF file.\n") ;
660 return(1) ;
661 }
662 strcpy(time_units, (char *) timeunits_attr->data) ;
663 strcat(time_units, DFLTORIGIN) ;
664 }
665 if (time_units) /* use units with default origin */ {
666 if (utScan(time_units, &timeunit)) {
667 gaprnt(0, "Error parsing time_units for SDF file.\n") ;
668 return(1) ;
669 }
670 } else {
671 if (!get_ud_time_unit(Tcoord, &timeunit)) {
672 gaprnt(0, "Error parsing UDUNITS time units in SDF file.\n") ;
673 return(1) ;
674 }
675 time_units = (char *) timeunits_attr->data ;
676 } /* using file's time units */
677 if (!decode_ud_time(&timeunit, time1, &iyr, &imo, &idy, &ihr, &imn,
678 &dsec)) {
679 gaprnt(0,
680 "Error decoding initial udunits time value in SDF file.\n") ;
681 return(1) ;
682 }
683 if (imo == 0) imo = 1 ;
684 if (idy == 0) idy = 1 ;
685 tempfarry[0] = (float) iyr ;
686 tempfarry[1] = (float) imo ;
687 tempfarry[2] = (float) idy ;
688 tempfarry[3] = (float) ihr ;
689 tempfarry[4] = (float) imn ;
690 if (pfi->dnum[TINDEX] > 1) /* deduce increment */ {
691 int compare_units(char *test_units, char *trunc_units) ;
692 char *trunc_units, *temp_str ;
693 int trunc_point ;
694
695 temp_str = strstr(time_units, " since ") ;
696 /* the above guaranteed to succeed since default origin added if need be */
697 if (!temp_str) /* just in case, though */ {
698 trunc_point = strlen(time_units) ;
699 } else {
700 trunc_point = (int) (temp_str - time_units) ;
701 }
702 trunc_units = (char *) malloc(sizeof(char) * (size_t) (trunc_point+1)) ;
703 strncpy(trunc_units, time_units, trunc_point) ;
704 trunc_units[trunc_point] = '\0' ;
705 start[0] = 1 ;
706 incrfactor = time2 - time1 ;
707 if (compare_units("year", trunc_units)) /* match is 1 */ {
708 tempfarry[5] = 12.0 * ((float) incrfactor) ;
709 if (tempfarry[5] < 1.0) {
710 gaprnt(0, "Time unit in SDF file has too small an increment (min. 1 minute).\n") ;
711 return(1) ;
712 }
713 tempfarry[6] = 0.0 ;
714 } else {
715 /* COARDS conventions say only year, day, hour, and minute are OK, not month */
716 /* But I've accepted Camiel's patch for "months since ..." -Hoop 2K/07/25 */
717 tempfarry[5] = 0.0 ;
718 /* CAS: Added month support */
719 if ((!strncmp(time_units, "month", 5)) ||
720 (!strncmp(time_units, "common_year/12", 14)) ||
721 (!strncmp(time_units, "common_years/12", 15))) {
722 /* if (compare_units("month", trunc_units)) { */
723 /* I doubt udunits will help us with units of "month" */
724 tempfarry[5] = ((float) incrfactor) ;
725 tempfarry[6] = 0.0 ;
726 if (tempfarry[5] < 1.0) {
727 gaprnt(0, "SDF: Fractional months are ill-defined and not supported by GrADS\n") ;
728 return(1) ;
729 }
730 /* CAS: End of month support */
731 } else if (compare_units("day", trunc_units)) {
732 if (incrfactor < 28.0) {
733 tempfarry[6] = ((float) incrfactor) * 24.0 * 60.0 ;
734 if (tempfarry[6] < 1.0) {
735 gaprnt(0, "SDF: Time unit has too small an increment (min. 1 minute).\n") ;
736 return(1) ;
737 }
738 } else if (incrfactor < 360.0) /* assume really months */ {
739 /* This dirty trick should get the right number of months for monthly, */
740 /* bi-monthly, and seasonal data. If there's anything between that and */
741 /* annual data (which should have units of "year(s) since"), we're broken */
742 tempfarry[5] = (float)(((int) (incrfactor + 0.5)) / 28) ;
743 tempfarry[6] = 0.0 ;
744 if (tempfarry[5] < 1.0) {
745 gaprnt(0, "SDF: Time unit has too small an increment (min. 1 minute).\n") ;
746 return(1) ;
747 }
748 } else /* annual or multi-annual w/"days since" */{
749 /* also a dirty trick to figure out how many years & mult. by 12 */
750 tempfarry[5] =
751 12.0 * ((float)(((int) (incrfactor + 0.5))/360)) ;
752 tempfarry[6] = 0.0 ;
753 }
754 } else if (compare_units("hour", trunc_units)) {
755 if (incrfactor < (28.0 * 24.0)) {
756 tempfarry[6] = ((float) incrfactor) * 60.0 ;
757 if (tempfarry[6] < 1.0) {
758 gaprnt(0, "SDF: Time unit has too small an increment (min. 1 minute).\n") ;
759 return(1) ;
760 }
761 } else {
762 if (incrfactor >= (360 * 24)) /* try years? */{
763 tempfarry[5] = 12.0 * ((float) ((int)
764 (((int)(incrfactor + 0.5)) / (360.0 * 24.0)))) ;
765 tempfarry[6] = 0.0 ;
766 } else /* assume really months */ {
767 tempfarry[5] = (float) (((int) (incrfactor + 0.5))/
768 (28 * 24)) ;
769 tempfarry[6] = 0.0 ;
770 }
771 if (tempfarry[5] < 1.0) {
772 gaprnt(0, "SDF: Time unit has too small an increment (min. 1 minute).\n") ;
773 return(1) ;
774 }
775 }
776 } else if (compare_units("minute", trunc_units)) {
777 if (incrfactor < (60.0 * 24.0 * 28.0)) {
778 tempfarry[5] = 0.0 ;
779 tempfarry[6] = (float) incrfactor ;
780 if (tempfarry[6] < 1.0) {
781 gaprnt(0, "SDF: Time unit has too small an increment (min. 1 minute).\n") ;
782 return(1) ;
783 }
784 } else /* monthly or greater */ {
785 if (tempfarry[6] < (60.0 * 24.0 * 360.0)) {
786 tempfarry[5] =
787 (float) (((int) ((incrfactor / (60.0 * 24.0)) + 0.5)) / 28) ;
788 tempfarry[6] = 0.0 ;
789 } else {
790 gaprnt(0,"SDF: Attempt at >monthly increment with 'minutes since' time units attribute value.\n") ;
791 return(1) ;
792 }
793 }
794 } else if (compare_units("seconds", trunc_units)) {
795 if (incrfactor < 60.0) {
796 gaprnt(0,"SDF: Time unit has too small an increment (min. 1 minute).\n") ;
797 return(1) ;
798 } else {
799 if (incrfactor < (60.0 * 60.0 * 24.0 * 28.0)) {
800 /* less than monthly, so use tempfarry[6] */
801 tempfarry[6] = incrfactor / 60.0 ;
802 } else /* monthly or greater */ {
803 if (incrfactor < (60.0 * 60.0 * 24.0 * 360)) {
804 /* assume monthly */
805 tempfarry[5] =
806 (float) (((int) ((incrfactor/(60.0 * 60.0 * 24.0)) + 0.5)) / 28) ;
807 tempfarry[6] = 0.0 ;
808 } else {
809 gaprnt(0,"SDF: Attempt at > monthly interval with 'seconds since' time units attribute.\n") ;
810 return(1) ;
811 }
812 }
813 }
814 } else {
815 gaprnt(0, "Error parsing time units in SDF file.\n") ;
816 return(1) ;
817 }
818 } /* finer than years resolution */
819 if (trunc_units) free(trunc_units) ;
820 } else {
821 tempfarry[5] = 0.0 ;
822 tempfarry[6] = 1.0 ; /* KLUDGE for one time-step files */
823 } /* end deduce increment */
824 #ifdef TESTING
825 fprintf(stderr, "time1=%lf, time2=%lf, incrfactor=%lf.\n", time1, time2,
826 incrfactor) ;
827 fprintf(stderr, "tempfarry={%f, %f, %f, %f, %f, %f, %f}.\n", tempfarry[0],
828 tempfarry[1], tempfarry[2], tempfarry[3], tempfarry[4],
829 tempfarry[5], tempfarry[6]) ;
830 fprintf(stderr, "diminq return=%d, time dimsize=%ld.\n",
831 ncdiminq(pfi->sdf_ptr->cdfid,3,(char *)0,&tdimsiz), tdimsiz) ;
832 #endif
833 } /* end if udunits time */
834 } /* has a T coordinate */
835 } /* doing T setup */
836 if (parms.dvsrch) {
837 pfi->vnum = numdvars ;
838 parms.dvcount = pfi->vnum ;
839 pfi->ivnum = pfi->lvnum = 0 ;
840 size = numdvars * sizeof(struct gavar) ;
841 pvar = (struct gavar *) malloc(size);
842 pfi->pvar1 = pvar;
843 parms.dvcount = numdvars ;
844 parms.dvsetup = (int *) malloc(numdvars * sizeof(int)) ;
845 for (xedni = 0 ; xedni < pfi->vnum ; ++xedni) {
846 parms.dvsetup[xedni] = 1 ;
847 }
848 } else {
849 pvar = pfi->pvar1 ;
850 }
851 if (parms.isxdf) {
852 int xdfopenerr ;
853
854 xdfopenerr = 0 ;
855 if (!parms.tsrch) {
856 if (parms.tdimname) {
857 if (find_dim(pfi->sdf_ptr, parms.tdimname) == -1) {
858 gaprnt(0,"Alleged time dimension (from TDEF) is not an SDF dimension. Cannot continue.\n") ;
859 gaprnt(0,"Time dimension name is:\n") ;
860 gaprnt(0,parms.tdimname) ;
861 gaprnt(0,"\n.\n") ;
862 xdfopenerr = 1 ;
863 }
864 }
865 }
866 if ((!parms.zsrch) && (find_dim(pfi->sdf_ptr, parms.zdimname) == -1)) {
867 gaprnt(0,"Alleged level dimension (from ZDEF) is not an SDF dimension. Cannot continue.\n") ;
868 gaprnt(0,"Level dimension name is:\n") ;
869 gaprnt(0,parms.zdimname) ;
870 gaprnt(0,"\n.\n") ;
871 xdfopenerr = 2 ;
872 }
873 if ((!parms.ysrch) && (find_dim(pfi->sdf_ptr, parms.ydimname) == -1)) {
874 gaprnt(0,"Alleged lon dimension (from YDEF) is not an SDF dimension. Cannot continue.\n") ;
875 gaprnt(0,"Lon dimension name is:\n") ;
876 gaprnt(0,parms.ydimname) ;
877 gaprnt(0,"\n.\n") ;
878 xdfopenerr = 3 ;
879 }
880 if ((!parms.xsrch) && (find_dim(pfi->sdf_ptr, parms.xdimname) == -1)) {
881 gaprnt(0,"Alleged lat dimension (from XDEF) is not an SDF dimension. Cannot continue.\n") ;
882 gaprnt(0,"Lat dimension name is:\n") ;
883 gaprnt(0,parms.xdimname) ;
884 gaprnt(0,"\n.\n") ;
885 xdfopenerr = 4 ;
886 }
887 if (xdfopenerr) {
888 return(1) ;
889 }
890 } /* if it's from xdfopen */
891 if (parms.isxdf && (!parms.tsrch)) {
892 if (!parms.tdimname) {
893 tdimension_id = -1 ;
894 } else {
895 tdimension_id = find_dim(pfi->sdf_ptr, parms.tdimname) ;
896 }
897 } else {
898 if (Tcoord) {
899 tdimension_id = find_dim(pfi->sdf_ptr, Tcoord->varnam) ;
900 } else {
901 tdimension_id = -1 ;
902 }
903 }
904 if (parms.isxdf && (!parms.zsrch)) {
905 zdimension_id = find_dim(pfi->sdf_ptr, parms.zdimname) ;
906 } else {
907 if (Zcoord) {
908 zdimension_id = find_dim(pfi->sdf_ptr, Zcoord->varnam) ;
909 } else {
910 zdimension_id = -1 ;
911 }
912 }
913 if (parms.isxdf && (!parms.ysrch)) {
914 ydimension_id = find_dim(pfi->sdf_ptr, parms.ydimname) ;
915 } else {
916 if (Ycoord) {
917 ydimension_id = find_dim(pfi->sdf_ptr, Ycoord->varnam) ;
918 } else {
919 ydimension_id = -1 ;
920 }
921 }
922 if (parms.isxdf && (!parms.xsrch)) {
923 xdimension_id = find_dim(pfi->sdf_ptr, parms.xdimname) ;
924 } else {
925 if (Xcoord) {
926 xdimension_id = find_dim(pfi->sdf_ptr, Xcoord->varnam) ;
927 } else {
928 xdimension_id = -1 ;
929 }
930 }
931 for (lclvar = pfi->sdf_ptr->var ; lclvar ; lclvar = lclvar->next) {
932 if (parms.isxdf && (!parms.dvsrch)) {
933 boole = findmatch(lclvar->varnam, parms.dvncnames, parms.dvcount,
934 &inx) ;
935 if (boole) {
936 while (strcmp(pvar->abbrv, parms.dvganames[inx])) {
937 pvar++ ;
938 } /* so we mess with right pvar from gadxdf */
939 }
940 } else {
941 boole = isdvar(lclvar, Xcoord, Ycoord, Zcoord, Tcoord,
942 xdimension_id, ydimension_id, zdimension_id,
943 tdimension_id) ;
944 }
945 if (boole) {
946 if (parms.isxdf && (!parms.tsrch)) {
947 if (!(parms.tdimname)) /* tdim-less support */ {
948 lclvar->dimidmap[0] = -1 ;
949 lclvar->dimmap[0] = -1 ;
950 } else {
951 lclvar->dimidmap[0] = tdimension_id ;
952 lclvar->dimmap[0] =
953 find_dimix(lclvar->vardimid, lclvar->nvardims,
954 lclvar->dimidmap[0]) ;
955 }
956 } else {
957 if (Tcoord) {
958 lclvar->dimidmap[0] = tdimension_id ;
959 lclvar->dimmap[0] = find_dimix(lclvar->vardimid,
960 lclvar->nvardims, lclvar->dimidmap[0]) ;
961 } /* leave -1 otherwise and if not in this dvar */
962 }
963 if (parms.isxdf && (!parms.zsrch)) {
964 lclvar->dimidmap[1] = zdimension_id ;
965 lclvar->dimmap[1] =
966 find_dimix(lclvar->vardimid, lclvar->nvardims,
967 lclvar->dimidmap[1]) ;
968 } else {
969 if (Zcoord) {
970 lclvar->dimidmap[1] = zdimension_id ;
971 lclvar->dimmap[1] =
972 find_dimix(lclvar->vardimid, lclvar->nvardims,
973 lclvar->dimidmap[1]) ;
974 } /* leave -1 otherwise and if not in this dvar */
975 }
976 if (parms.isxdf && (!parms.ysrch)) {
977 lclvar->dimidmap[2] = ydimension_id ;
978 } else {
979 lclvar->dimidmap[2] = ydimension_id ;
980 }
981 lclvar->dimmap[2] = find_dimix(lclvar->vardimid, lclvar->nvardims,
982 lclvar->dimidmap[2]) ;
983 if (parms.isxdf && (!parms.xsrch)) {
984 lclvar->dimidmap[3] = xdimension_id ;
985 } else {
986 lclvar->dimidmap[3] = xdimension_id ;
987 }
988 lclvar->dimmap[3] = find_dimix(lclvar->vardimid, lclvar->nvardims,
989 lclvar->dimidmap[3]) ;
990 boole2 = 0 ;
991 if (!(parms.isxdf)) {
992 boole2 = 1 ;
993 } else {
994 if (parms.dvsrch) {
995 boole2 = 1 ;
996 } else {
997 if ((parms.dvcount > inx) && (parms.dvsetup[inx])) {
998 boole2 = 1 ;
999 } else {
1000 boole2 = 0 ;
1001 }
1002 }
1003 }
1004 if (boole2) {
1005 if ((!(parms.isxdf)) || (parms.dvsrch)) {
1006 /* if we got the data variables by searching in this routine */
1007 strncpy(pvar->abbrv, lclvar->varnam, 15) ;
1008 pvar->abbrv[15] = '\0' ;
1009 lowcas(pvar->abbrv) ;
1010 }
1011 strcpy(lclvar->gradsabbr, pvar->abbrv) ;
1012 if (lclvar->dimmap[1] != -1) {
1013 pvar->levels = pfi->dnum[ZINDEX] ; /* by def'n in SDF file */
1014 ++(pfi->lvnum) ;
1015 } else {
1016 pvar->levels = 0 ; /* I hope that's right... */
1017 ++(pfi->ivnum) ;
1018 }
1019 for (j = 0 ; j < 4 ; j++) pvar->units[j] = -999;
1020 if (!(title_vattr = find_att(lclvar->first_vattr, "long_name"))) {
1021 /* gad, I hope it's an old CDC standard file */
1022 if (!(title_vattr = find_att(lclvar->first_vattr, "title"))) {
1023 /* give up */
1024 pvar->varnm[0] = '\0' ;
1025 }
1026 }
1027 if (title_vattr) {
1028 if (title_vattr->len > 79) {
1029 strncpy(pvar->varnm, (char *)(title_vattr->data), 79) ;
1030 pvar->varnm[79] = '\0' ;
1031 } else {
1032 strncpy(pvar->varnm, (char *)(title_vattr->data),
1033 title_vattr->len) ;
1034 pvar->varnm[title_vattr->len] = '\0' ;
1035 }
1036 }
1037 pvar->offset = 0 ;
1038 pvar->recoff = 0 ;
1039 pvar->dfrm = 0 ;
1040 pvar->var_t = 0 ;
1041 pvar->var_z = 0 ;
1042 pvar->y_x = 0 ;
1043 } else {
1044 strcpy(lclvar->gradsabbr, pvar->abbrv) ;
1045 }
1046 if (lclvar->vartype == NC_DOUBLE) {
1047 sprintf(pout,"SDF: The double precision values of the variable %s\n", pvar->abbrv) ;
1048 strcat(pout,"will be automatically converted to single precision.\n") ;
1049 gaprnt(0, pout) ;
1050 }
1051 if (parms.needundef) /* need to divine it now */ {
1052 /* get missing_value attribute value, and stuff it into lclvar->missing_value */
1053 missval_attr = find_att(lclvar->first_vattr, "missing_value") ;
1054 if (!missval_attr) {
1055 missval_attr = find_att(lclvar->first_vattr, "_FillValue") ;
1056 }
1057 /* if the type of the data variable and that of the missing_value attribute */
1058 /* don't match, use the default, as we do if we have no attribute */
1059 if ((!missval_attr) || (missval_attr->type != lclvar->vartype)) {
1060 switch (lclvar->vartype) {
1061 case NC_BYTE:
1062 lclvar->missing_value.bval = FILL_BYTE ;
1063 break ;
1064 case NC_SHORT:
1065 lclvar->missing_value.sval = FILL_SHORT ;
1066 break ;
1067 case NC_LONG:
1068 lclvar->missing_value.lval = FILL_LONG ;
1069 break ;
1070 case NC_FLOAT:
1071 lclvar->missing_value.fval = FILL_FLOAT ;
1072 break ;
1073 case NC_DOUBLE:
1074 lclvar->missing_value.dval = FILL_DOUBLE ;
1075 break ;
1076 } /* end switch */
1077 } else /* have an attribute value and types match */ {
1078 switch (lclvar->vartype) {
1079 case NC_BYTE:
1080 lclvar->missing_value.bval = *((signed char *)missval_attr->data) ;
1081 break ;
1082 case NC_SHORT:
1083 lclvar->missing_value.sval = *((short *)missval_attr->data) ;
1084 break ;
1085 case NC_LONG:
1086 lclvar->missing_value.lval = *((long *)missval_attr->data) ;
1087 break ;
1088 case NC_FLOAT:
1089 lclvar->missing_value.fval = *((float *)missval_attr->data) ;
1090 break ;
1091 case NC_DOUBLE:
1092 lclvar->missing_value.dval = *((double *)missval_attr->data) ;
1093 break ;
1094 } /* end switch */
1095 } /* else is have attribute and types match */
1096 } else /* use value from UNDEF statement in DDF */ {
1097 lclvar->missing_value.fval = pfi->undef ; /* always float */
1098 if (parms.hasDDFundef) {
1099 pfi->sdf_ptr->hasDDFundef = 1 ;
1100 }
1101 } /* end if need to get missing_value attribute value */
1102 if (parms.isxdf && (!parms.dvsrch)) {
1103 pvar = pfi->pvar1 ; /* reset for searching out next one */
1104 } else {
1105 pvar++;
1106 }
1107 }
1108 } /* end for each var in SDF file or descriptor file varlist */
1109 /* Figure out locations of variables within a time group ( N.A. ) */
1110
1111 /* Allocate an I/O buffer the size of one grid */
1112
1113 maxlv = pfi->dnum[ZINDEX] ;
1114 if (pfi->type > 1) {
1115 size = maxlv * sizeof(float);
1116 } else {
1117 size = pfi->dnum[XINDEX] * pfi->dnum[YINDEX] * sizeof(float);
1118 }
1119
1120 /* If the file name is a time series template, figure out
1121 which times go with which files, so we don't waste a lot
1122 of time later opening and closing files unnecessarily. */
1123
1124 if ((!(parms.isxdf)) && (pfi->tmplat)) {
1125 pfi->dnum[TINDEX] = ntsteps ;
1126 strcpy(pfi->name, template) ;
1127 pfi->fnums = (int *) malloc(sizeof(int) * pfi->dnum[TINDEX]) ;
1128 if (pfi->fnums == NULL) {
1129 gaprnt(0, "SDFOPEN error: Memory allocation error.\n") ;
1130 /* what could I do if close_netcdf did bomb? */
1131 (void) close_netcdf(pfi->sdf_ptr->cdfid) ;
1132 if ((pfi->sdf_ptr) == (pfi->first_sdf_ptr)) {
1133 pfi->first_sdf_ptr = (IO_STD *) 0 ;
1134 }
1135 free_io_std(&(pfi->sdf_ptr)) ;
1136 return(1) ;
1137 }
1138 j = 1 ;
1139 gr2t(pfi->grvals[TINDEX], 1.0, &tdefi) ;
1140 ch = gafndt(pfi->name, &tdefi, &tdefi, pfi->abvals[TINDEX], NULL, 1) ;
1141 if (ch == NULL) {
1142 gaprnt(0, "SDFOPEN error: Memory allocation error.\n") ;
1143 (void) close_netcdf(pfi->sdf_ptr->cdfid) ;
1144 if ((pfi->sdf_ptr) == (pfi->first_sdf_ptr)) {
1145 pfi->first_sdf_ptr = (IO_STD *) 0 ;
1146 }
1147 free_io_std(&(pfi->sdf_ptr)) ;
1148 return(1) ;
1149 }
1150 pfi->fnums[0] = j ;
1151 for (i = 2 ; i <= pfi->dnum[TINDEX] ; ++i) {
1152 gr2t(pfi->grvals[TINDEX], (float) i, &tdef) ;
1153 pos = gafndt(pfi->name, &tdef, &tdefi, pfi->abvals[TINDEX], NULL, 1) ;
1154 if (pos == NULL) {
1155 gaprnt(0, "SDFOPEN error: Memory allocation error.\n") ;
1156 (void) close_netcdf(pfi->sdf_ptr->cdfid) ;
1157 free_io_std(&(pfi->sdf_ptr)) ;
1158 return(1) ;
1159 }
1160 if (strcmp(ch, pos) != 0) {
1161 j = i ;
1162 if (ch) {
1163 free(ch) ;
1164 }
1165 ch = pos ;
1166 }
1167 pfi->fnums[i - 1] = j ;
1168 }
1169 if (ch) {
1170 free(ch) ;
1171 }
1172 pfi->fnumc = 0 ;
1173 } /* if templating in use and not XDF */ else if (pfi->tmplat) {
1174 strcpy(pfi->name, template) ;
1175 }
1176 if (copy_io_std(&(pfi->first_sdf_ptr), pfi->sdf_ptr) == Success) {
1177 return(0);
1178 } else {
1179 gaprnt(0, "SDF operational error (copy_io_std failed making current first)\n") ;
1180 return(1) ;
1181 }
1182 }
1183
1184 int
compare_units(char * test_unit,char * trunc_unit)1185 compare_units(char *test_unit, char *trunc_unit) {
1186 utUnit testing_unit, truncated_unit ;
1187 int dequal(double op1, double op2, double threshold) ;
1188 int return_value ;
1189 double slope, intercept ;
1190
1191 return_value = utScan(test_unit, &testing_unit) ;
1192 if (return_value != 0) /* probably shouldn't bitch about an err, but strange */ {
1193 return(0) ;
1194 }
1195 return_value = utScan(trunc_unit, &truncated_unit) ;
1196 if (return_value != 0) /* see above comment */ {
1197 return(0) ;
1198 }
1199 return_value = utConvert(&truncated_unit, &testing_unit, &slope, &intercept) ;
1200 if (return_value != 0) {
1201 return(0) ;
1202 }
1203 if (dequal(slope, 1.0, 0.000001) && dequal(intercept, 0.0, 0.000001)) {
1204 return(1) ; /* hopefully, this will catch all combos from common_year */
1205 } else {
1206 return(0) ;
1207 }
1208 /*NOTREACHED*/
1209 }
1210
1211 int
findmatch(char * var_name,char ** varnamelist,int namecount,int * inxptr)1212 findmatch(char *var_name, char **varnamelist, int namecount, int *inxptr) {
1213 int inx ;
1214
1215 for (inx = 0 ; inx < namecount ; ++inx) {
1216 if (!strcmp(var_name, varnamelist[inx])) {
1217 *inxptr = inx ;
1218 return(1) ;
1219 }
1220 }
1221 *inxptr = -1 ;
1222 return(0) ;
1223 }
1224
1225
1226 void
doparms4sdf(GASDFPARMS * parms)1227 doparms4sdf(GASDFPARMS *parms) {
1228
1229 /* set values in parm struct as appropos for sdfopen (not xdfopen) */
1230
1231 parms->xsrch = parms->ysrch = parms->zsrch = parms->tsrch = 1 ;
1232 parms->dvsrch = 1 ;
1233 parms->isxdf = 0 ;
1234 parms->xsetup = parms->ysetup = parms->zsetup = parms->tsetup = 1 ;
1235 parms->needtitle = parms->needundef = 1 ;
1236 parms->xdimname = parms->ydimname = parms->zdimname = parms->tdimname = (char *) 0 ;
1237 parms->dvcount = -1 ;
1238 parms->dvncnames = parms->dvganames = (char **) 0 ;
1239 parms->dvsetup = (int *) 0 ;
1240 parms->hasDDFundef = 0 ;
1241 return ;
1242 }
1243
1244 int
isdvar(VAR_INFO * var,VAR_INFO * X,VAR_INFO * Y,VAR_INFO * Z,VAR_INFO * T,int xdimid,int ydimid,int zdimid,int tdimid)1245 isdvar(VAR_INFO *var, VAR_INFO *X, VAR_INFO *Y, VAR_INFO *Z, VAR_INFO *T,
1246 int xdimid, int ydimid, int zdimid, int tdimid) {
1247 int idim, hasX, hasY, hasT ;
1248
1249 if ((var->nvardims) < 2) {
1250 return(0) ;
1251 }
1252 hasX = hasY = hasT = 0 ;
1253 for (idim = 0 ; idim < var->nvardims ; ++idim) {
1254 if (var->vardimid[idim] == xdimid) {
1255 hasX = 1 ;
1256 }
1257 if (var->vardimid[idim] == ydimid) {
1258 hasY = 1 ;
1259 }
1260 if (tdimid > -1) {
1261 if (var->vardimid[idim] == tdimid) {
1262 hasT = 1 ;
1263 }
1264 }
1265 }
1266 if ((X && (var == X)) || (Y && (var == Y)) || (Z && (var == Z)) ||
1267 (T && (var == T)) ) {
1268 return(0) ;
1269 }
1270 return(hasX || hasY) ;
1271 }
1272
1273 int
sdfdeflev(struct gafile * pfi,VAR_INFO * coord,int GrADSdimnum,int revflag)1274 sdfdeflev(struct gafile *pfi, VAR_INFO *coord, int GrADSdimnum, int revflag) {
1275 float *tempfarry, *tempfarry2 ;
1276 int indx ;
1277 int convtype(void *srcptr, nc_type srctype, void *desptr, nc_type destype,
1278 int srcindex, int desindex) ;
1279
1280 #ifdef TESTING2
1281 fprintf(stderr, "sdfdeflev(%s, %s, %d, %d);\n", pfi->name, coord->varnam,
1282 GrADSdimnum, revflag);
1283 #endif
1284 /* the following adapted from deflev routine */
1285 tempfarry2 =
1286 (float *) malloc((pfi->dnum[GrADSdimnum] + 5) * sizeof(float)) ;
1287 if (!tempfarry2) {
1288 gaprnt(0, "Failure setting up coordinate levels in SDF file.\n") ;
1289 return(1) ;
1290 }
1291 tempfarry2[0] = pfi->dnum[GrADSdimnum] ;
1292 for (indx = 1 ; indx <= (pfi->dnum[GrADSdimnum]) ; ++indx) {
1293 if (!revflag) {
1294 if (!convtype(coord->data, coord->vartype, tempfarry2, NC_FLOAT,
1295 indx - 1, indx)) {
1296 gaprnt(0, "SDF: Error interpreting coordinate value.\n") ;
1297 return(1) ;
1298 }
1299 } else {
1300 if (!convtype(coord->data, coord->vartype, tempfarry2, NC_FLOAT,
1301 pfi->dnum[GrADSdimnum] - indx, indx)) {
1302 gaprnt(0, "SDF: Error interpreting coordinate value.\n") ;
1303 return(1) ;
1304 }
1305 }
1306 }
1307 tempfarry2[(pfi->dnum[GrADSdimnum]) + 1] = -999.9 ;
1308 pfi->abvals[GrADSdimnum] = tempfarry2 ;
1309 pfi->grvals[GrADSdimnum] = tempfarry2 ;
1310 pfi->ab2gr[GrADSdimnum] = lev2gr ;
1311 pfi->gr2ab[GrADSdimnum] = gr2lev ;
1312 pfi->linear[GrADSdimnum] = 0 ;
1313 /* end of stuff adapted from deflev routine */
1314 return(0) ;
1315 }
1316
1317
1318 int
findX(IO_STD * sdf_ptr,VAR_INFO ** Xcoordptr)1319 findX(IO_STD *sdf_ptr, VAR_INFO **Xcoordptr) {
1320 /* check for coordinate variable whose units are degrees_east, degree_east */
1321 /* degrees_E, or degree_E */
1322 struct attrib_list *units_vattr ;
1323 struct attrib_list *find_att(struct attrib_list *first_attr, char *attname) ;
1324 VAR_INFO *lclvar ;
1325 int iscoordvar, i ;
1326 char *ch ;
1327
1328 *Xcoordptr = NULL ;
1329 for (lclvar = sdf_ptr->var ; lclvar ; lclvar = lclvar->next) {
1330 iscoordvar = 0 ;
1331 if (lclvar->nvardims == 1) {
1332 for (i = 0 ; i < sdf_ptr->ndims ; ++i) {
1333 if (lclvar->vardimid[0] == sdf_ptr->dimids[i]) {
1334 if (!strcmp(sdf_ptr->dimnam[i], lclvar->varnam)) {
1335 iscoordvar = 1 ;
1336 }
1337 }
1338 } /* end for */
1339 }
1340 if (iscoordvar) /* it's a coordinate variable */ {
1341 if (!(units_vattr = find_att(lclvar->first_vattr, "units"))) {
1342 /* return(0) ; */ continue ;
1343 }
1344 ch = (char *) malloc(units_vattr->len + 1) ;
1345 strncpy(ch, (char *) units_vattr->data, units_vattr->len) ;
1346 ch[units_vattr->len] = '\0' ;
1347 /* Use strncmp to make up for ncgen's annoying lack of NULL terminators */
1348 if (strncmp(ch, "degrees_east", 12)) {
1349 if (strncmp(ch, "degree_east", 11)) {
1350 if (strncmp(ch, "degrees_E", 9)) {
1351 if (strncmp(ch, "degree_E", 8)) {
1352 if (ch) {
1353 free(ch) ;
1354 }
1355 continue ;
1356 }
1357 }
1358 }
1359 }
1360 if (ch) {
1361 free(ch) ;
1362 }
1363 *Xcoordptr = lclvar ;
1364 return(1) ; /* got a match on one of them */
1365 }
1366 } /* end for each variable */
1367 return(0) ;
1368 }
1369
1370 int
findY(IO_STD * sdf_ptr,VAR_INFO ** Ycoordptr)1371 findY(IO_STD *sdf_ptr, VAR_INFO **Ycoordptr) {
1372 /* check for coordinate variable whose units are degrees_north, degree_north */
1373 /* degrees_N, or degree_N */
1374 struct attrib_list *units_vattr ;
1375 struct attrib_list *find_att(struct attrib_list *first_attr, char *attname) ;
1376 VAR_INFO *lclvar ;
1377 int iscoordvar, i ;
1378 char *ch ;
1379
1380 *Ycoordptr = NULL ;
1381 for (lclvar = sdf_ptr->var ; lclvar ; lclvar = lclvar->next) {
1382 iscoordvar = 0 ;
1383 if (lclvar->nvardims == 1) {
1384 for (i = 0 ; i < sdf_ptr->ndims ; ++i) {
1385 if (lclvar->vardimid[0] == sdf_ptr->dimids[i]) {
1386 if (!strcmp(sdf_ptr->dimnam[i], lclvar->varnam)) {
1387 iscoordvar = 1 ;
1388 }
1389 }
1390 } /* end for */
1391 }
1392 if (iscoordvar) /* it's a coordinate variable */ {
1393 if (!(units_vattr = find_att(lclvar->first_vattr, "units"))) {
1394 continue ; /* try next variable */
1395 }
1396 ch = (char *) malloc(units_vattr->len + 1) ;
1397 strncpy(ch, (char *) units_vattr->data, units_vattr->len) ;
1398 ch[units_vattr->len] = '\0' ;
1399 /* Use strncmp to make up for ncgen's annoying lack of NULL terminators */
1400 if (strncmp(ch, "degrees_north", 13)) {
1401 if (strncmp(ch, "degree_north", 12)) {
1402 if (strncmp(ch, "degrees_N", 9)) {
1403 if (strncmp(ch, "degree_N", 8)) {
1404 if (ch) {
1405 free(ch) ;
1406 }
1407 continue ;
1408 }
1409 }
1410 }
1411 }
1412 if (ch) {
1413 free(ch) ;
1414 }
1415 *Ycoordptr = lclvar ;
1416 return(1) ; /* got a match on one of them */
1417 }
1418 } /* end for each variable */
1419 return(0) ;
1420 }
1421
1422 int
findZ(IO_STD * sdf_ptr,VAR_INFO ** Zcoordptr,int * ispressptr)1423 findZ(IO_STD *sdf_ptr, VAR_INFO **Zcoordptr, int *ispressptr) {
1424 /* check for coordinate variable whose units are that of pressure, */
1425 /* or another unit approved by COARDS conventions */
1426 /* initially, the pressure units are "millibars" or "pascals" (caseless) */
1427 /* should probably allow for prefixes through udunits package */
1428 /* Will allow exact match on "mb" as a kludge for misbegotten NMC data */
1429 struct attrib_list *units_vattr, *positive_vattr ;
1430 struct attrib_list *find_att(struct attrib_list *first_attr, char *attname) ;
1431 VAR_INFO *lclvar ;
1432 int iscoordvar, i ;
1433 char *ch, *lwrunits ;
1434 struct utUnit feet, thisguy, pascals, kelvins ;
1435 double slope, intcept ;
1436
1437 *Zcoordptr = NULL ;
1438 if (utScan("feet", &feet) != 0) /* bloody hell... */ {
1439 gaprnt(0, "The udunits library doesn't know feet; giving up...\n") ;
1440 return(9999) ;
1441 }
1442 if (utScan("pascals", &pascals) != 0) /* bloody hell... */ {
1443 gaprnt(0, "The udunits library doesn't know pascals; giving up...\n") ;
1444 return(9999) ;
1445 }
1446 if (utScan("kelvins", &kelvins) != 0) /* bloody hell... */ {
1447 gaprnt(0, "The udunits library doesn't know kelvins; giving up...\n") ;
1448 return(9999) ;
1449 }
1450 for (lclvar = sdf_ptr->var ; lclvar ; lclvar = lclvar->next) {
1451 iscoordvar = 0 ;
1452 if (lclvar->nvardims == 1) {
1453 for (i = 0 ; i < sdf_ptr->ndims ; ++i) {
1454 if (lclvar->vardimid[0] == sdf_ptr->dimids[i]) {
1455 if (!strcmp(sdf_ptr->dimnam[i], lclvar->varnam)) {
1456 iscoordvar = 1 ;
1457 }
1458 }
1459 } /* end for */
1460 }
1461 if (iscoordvar) /* it's a coordinate variable */ {
1462 if (!(units_vattr = find_att(lclvar->first_vattr, "units"))) {
1463 continue ;
1464 }
1465 ch = (char *) units_vattr->data ;
1466 lwrunits = (char *) malloc(units_vattr->len + 1) ;
1467 strncpy(lwrunits, ch, units_vattr->len) ;
1468 lwrunits[units_vattr->len] = '\0' ;
1469 lowcas(lwrunits) ;
1470 if (!strcmp(lwrunits, "mb")) {
1471 *Zcoordptr = lclvar ;
1472 *ispressptr = 1 ;
1473 if (lwrunits) {
1474 free(lwrunits) ;
1475 }
1476 return(1) ;
1477 }
1478 /* depth, isothermic, sigma_level, & hybrid_sigma_pressure are the only others */
1479 if ((!strcmp(lwrunits, "sigma_level")) || (!strcmp(lwrunits, "degreesk")) ||
1480 (!strcmp(lwrunits, "degrees_k")) || (!strcmp(lwrunits, "level")) ||
1481 (!strcmp(lwrunits, "layer")) || (!strcmp(lwrunits, "layers"))) {
1482 *Zcoordptr = lclvar ;
1483 *ispressptr = 0 ;
1484 if (lwrunits) {
1485 free(lwrunits) ;
1486 }
1487 return(1) ;
1488 }
1489 if (!strcmp(lwrunits, "hybrid_sigma_pressure")) { /* CSM - NCAR */
1490 *Zcoordptr = lclvar ;
1491 *ispressptr = 1 ;
1492 if (lwrunits) {
1493 free(lwrunits) ;
1494 }
1495 return(1) ;
1496 }
1497 /* if we can convert the units to feet, then it could be depth */
1498 if (utScan(ch, &thisguy) == 0) {
1499 if (utConvert(&thisguy, &feet, &slope, &intcept) == 0) {
1500 *Zcoordptr = lclvar ;
1501 *ispressptr = 0 ;
1502 if (lwrunits) {
1503 free(lwrunits) ;
1504 }
1505 return(1) ;
1506 }
1507 /* if we can convert the units to pascals, then it could be pressure */
1508 if (utConvert(&thisguy, &pascals, &slope, &intcept) == 0) {
1509 *Zcoordptr = lclvar ;
1510 *ispressptr = 1 ;
1511 if (lwrunits) {
1512 free(lwrunits) ;
1513 }
1514 return(1) ;
1515 }
1516 /* if we can convert the units to kelvins, then it could be isothermic */
1517 if (utConvert(&thisguy, &kelvins, &slope, &intcept) == 0) {
1518 *Zcoordptr = lclvar ;
1519 if (lwrunits) {
1520 free(lwrunits) ;
1521 }
1522 *ispressptr = 0 ;
1523 return(1) ;
1524 }
1525 } /* if utScan-able */
1526 if (lwrunits) {
1527 free(lwrunits) ;
1528 }
1529 } /* if a coord. var. */
1530 } /* for each variable in SDF */
1531 return(0) ;
1532 }
1533
1534 int
findT(IO_STD * sdf_ptr,VAR_INFO ** Tcoordptr)1535 findT(IO_STD *sdf_ptr, VAR_INFO **Tcoordptr) {
1536 /* find a coordinate variable whose units mark it as one sort of time or
1537 another. */
1538 struct attrib_list *units_vattr ;
1539 struct attrib_list *find_att(struct attrib_list *first_attr, char *attname) ;
1540 VAR_INFO *lclvar ;
1541 int iscoordvar, i ;
1542 char *ch ;
1543 utUnit timeunit ;
1544
1545 *Tcoordptr = NULL ;
1546 for (lclvar = sdf_ptr->var ; lclvar ; lclvar = lclvar->next) {
1547 iscoordvar = 0 ;
1548 if (lclvar->nvardims == 1) {
1549 for (i = 0 ; i < sdf_ptr->ndims ; ++i) {
1550 if (lclvar->vardimid[0] == sdf_ptr->dimids[i]) {
1551 if (!strcmp(sdf_ptr->dimnam[i], lclvar->varnam)) {
1552 iscoordvar = 1 ;
1553 }
1554 }
1555 } /* end for */
1556 }
1557 if (!iscoordvar) continue ;
1558 units_vattr = find_att(lclvar->first_vattr, "units") ;
1559 if (!units_vattr) continue ;
1560 ch = (char *) malloc(units_vattr->len + 1) ;
1561 strncpy(ch, (char *) units_vattr->data, units_vattr->len) ;
1562 ch[units_vattr->len] = '\0' ;
1563 if (!strcasecmp(ch, "yyyymmddhhmmss")) /* old-style CDC */ {
1564 if (ch) {
1565 free(ch) ;
1566 }
1567 *Tcoordptr = lclvar ;
1568 return(1) ;
1569 }
1570 if (!strcasecmp(ch, "yymmddhh")) /* "DAO format" */ {
1571 if (ch) {
1572 free(ch) ;
1573 }
1574 *Tcoordptr = lclvar ;
1575 return(1) ;
1576 }
1577 if (utScan(ch, &timeunit)) {
1578 if (ch) {
1579 free(ch) ;
1580 }
1581 continue ;
1582 }
1583 if (ch) {
1584 free(ch) ;
1585 }
1586 if (utIsTime(&timeunit)) /* will now supply default unit */ {
1587 *Tcoordptr = lclvar ;
1588 return(1) ;
1589 }
1590 } /* foreach var. */
1591 return(0) ;
1592 }
1593
1594 /* Open a data set by reading the self-describing file for that
1595 data set, and create a gafile structure. Chain the gafile
1596 structure on to the list anchored in the gastat. */
1597
1598 int
gasdfopen(char * args,struct gacmn * pcm)1599 gasdfopen (char *args, struct gacmn *pcm) {
1600 struct gafile *pfi, *pfio;
1601 int size, rc, ntsteps, ichpos, idummy ,i, len;
1602 char pathname[4096], *optargs, template[4096], *cntsteps, *fnpart, *aux , *ch;
1603 GASDFPARMS parms ;
1604 int init_io_std(IO_STD **ptr) ;
1605 int gadsdf(char *name, struct gafile *pfi, char *template, int ntsteps,
1606 GASDFPARMS parms) ;
1607 void doparms4sdf(GASDFPARMS *parms) ;
1608
1609 doparms4sdf(&parms) ; /* set up for sdfopen calling gadsdf */
1610 getwrd(pathname, args, 4095) ;
1611 optargs = nxtwrd(args) ;
1612 gaprnt (2, "Scanning self-describing file: ");
1613 gaprnt (2, pathname);
1614 gaprnt (2, "\n");
1615 #ifdef TESTING3
1616 if (optargs && (*optargs != '\0')) {
1617 gaprnt (2, "Optional args to SDFOPEN command are: ") ;
1618 gaprnt (2, optargs) ;
1619 gaprnt (2, "\n") ;
1620 }
1621 #endif
1622 pfi = getpfi();
1623 if (pfi == NULL) {
1624 gaprnt (0,"Memory Allocation Error: Prior to SDF File Open\n");
1625 return (1);
1626 }
1627 pfi->is_a_SDF = 1 ;
1628 pfi->ncflg = 3; /* set a flag for detecting sdf files that is not #ifdef'd */
1629 if (!init_io_std(&(pfi->sdf_ptr))) {
1630 gaprnt(0, "Memory Allocation Error: On SDF File Open\n") ;
1631 return(1) ;
1632 }
1633 pfi->tmplat = 0 ; /* assume no templating, unless get two more good args */
1634 if (!optargs || (*optargs == '\0')) {
1635 template[0] = '\0' ;
1636 ntsteps = -999 ;
1637 } else {
1638 /* get the template string */
1639 ch = optargs;
1640 len = 0;
1641 while (*(ch+len)!=' ' && *(ch+len)!='\n' && *(ch+len)!='\t') len++;
1642 for (i=0; i<len; i++) *(template+i) = *(ch+i);
1643 *(template+len) = '\0';
1644 /* get the number of time steps */
1645 if ( (ch=nxtwrd(ch))==NULL ) {
1646 template[0] = '\0' ;
1647 ntsteps = -999 ;
1648 gaprnt(1,"Warning: Missing ntsteps arg to SDFOPEN command, ignoring templating info.\n");
1649 } else {
1650 if (!intprs(ch, &ntsteps) || ntsteps < 1) {
1651 template[0] = '\0' ;
1652 ntsteps = -999 ;
1653 gaprnt(1,"Warning: Bad ntsteps arg to SDFOPEN command, ignoring templating info.\n") ;
1654 }
1655 else {
1656 /* everything parsed ok, set templating flag */
1657 pfi->tmplat = 1;
1658 }
1659 }
1660 #ifdef TESTING3
1661 sprintf(pout," template=%s\n ntsteps=%d\n", template,ntsteps);
1662 gaprnt(2,pout);
1663 #endif
1664 } /* if no optargs */
1665
1666 rc = gadsdf(pathname, pfi, template, ntsteps, parms);
1667 #ifdef TESTING
1668 fprintf(stderr, "Return code from gadsdf=%d.\n", rc) ;
1669 fprintf(stderr, "Calling xdumpfile after gadsdf.\n") ;
1670 xdumpfile(pfi) ;
1671 if (pfi->pvar1) {
1672 fprintf(stderr, "Calling xdumpvar on first var after gadsdf.\n") ;
1673 xdumpvar(pfi->pvar1) ;
1674 }
1675 #endif
1676 if (rc) {
1677 frepfi (pfi, 0) ;
1678 return (rc) ;
1679 }
1680 if (pcm->pfi1==NULL) {
1681 pcm->pfi1 = pfi;
1682 } else {
1683 pfio = pcm->pfi1;
1684 while (pfio->pforw!=NULL) pfio = pfio->pforw;
1685 pfio->pforw = pfi;
1686 }
1687 pfi->pforw = NULL;
1688 pcm->fnum++;
1689
1690 if (pcm->fnum==1) {pcm->pfid = pcm->pfi1; pcm->dfnum = 1;}
1691 sprintf (pout,"SDF file %s is open as file %i\n",pfi->name,pcm->fnum);
1692 gaprnt (2,pout);
1693
1694 /* If first file open, set up some default dimension
1695 ranges for the user */
1696
1697 if (pcm->fnum==1) {
1698 if (pfi->type==2 || pfi->wrap) gacmd ("set lon 0 360",pcm,0);
1699 else {
1700 sprintf (pout,"set x 1 %i",pfi->dnum[0]);
1701 gacmd (pout,pcm,0);
1702 }
1703 if (pfi->type==2) {
1704 gacmd ("set lat -90 90",pcm,0);
1705 gacmd ("set lev 500",pcm,0);
1706 } else {
1707 sprintf (pout,"set y 1 %i",pfi->dnum[1]);
1708 gacmd (pout,pcm,0);
1709 gacmd ("set z 1",pcm,0);
1710 }
1711 gacmd ("set t 1",pcm,0);
1712 }
1713
1714 if (pfi->ppflag) {
1715 sprintf (pout,"Notice: Implied interpolation for SDF file %s\n",pathname);
1716 gaprnt (1,pout);
1717 gaprnt (1," Interpolation will be performed on any data ");
1718 gaprnt (1,"displayed from this file\n");
1719 }
1720
1721 return (0);
1722 }
1723
1724 /* Routine to open appropriate file when using file templates */
1725 /* Warning -- changes time value to time with respect to this file */
1726
gaopsdf(int t,struct gagrid * gridptr,VAR_INFO ** data_var_ptr)1727 int gaopsdf(int t, struct gagrid *gridptr, VAR_INFO **data_var_ptr) {
1728 int i, abbrvlen, neednewfile, abbrlen, idim ;
1729 struct dt dtim, dtimi ;
1730 struct gafile *pfi ;
1731 struct attrib_list *missval_attr ;
1732 struct attrib_list *find_att(struct attrib_list *first_att, char *attname) ;
1733 char *fn ;
1734 IO_STD *oldsdf_ptr ;
1735 VAR_INFO *oldlclvar, *newlclvar ;
1736 void free_io_std(IO_STD **sdf_ptr) ;
1737 int close_netcdf(int ncid) ;
1738 int init_io_std(IO_STD **sdf_ptr) ;
1739 int read_io_std(char *path, IO_STD *sdf_ptr) ;
1740 int copy_io_std(IO_STD **sdf_ptr2, IO_STD *std_ptr1) ;
1741
1742 neednewfile = 0 ;
1743 pfi = gridptr->pfile ;
1744 if ((t < 1) || (t > gridptr->pfile->dnum[TINDEX]) ) return(-99999) ;
1745 i = *(gridptr->pfile->fnums + t - 1) ;
1746 if (i != gridptr->pfile->fnumc) /* different file needed */ {
1747 #ifdef TESTING3
1748 fprintf(stderr, "Going to new SDF file for index %d.\n", i) ;
1749 #endif
1750 neednewfile = 1 ;
1751 if ((gridptr->pfile->tempname) && (gridptr->pfile->tempname != NULL)) {
1752 if (gridptr->pfile->tempname) {
1753 free(gridptr->pfile->tempname) ;
1754 }
1755 }
1756 oldsdf_ptr = (IO_STD *) NULL ;
1757 if (copy_io_std(&oldsdf_ptr, gridptr->pfile->sdf_ptr) == Failure) {
1758 gaprnt(1, "SDF templating: Error copying SDF info from default.\n") ;
1759 return(-88888) ;
1760 }
1761 if (oldsdf_ptr == (IO_STD *) 0) {
1762 if (copy_io_std(&oldsdf_ptr, gridptr->pfile->first_sdf_ptr) == Failure) {
1763 gaprnt(1, "SDF templating: Error copying SDF info from first.\n") ;
1764 return(-88888) ;
1765 }
1766 }
1767 gridptr->pfile->sdf_ptr = (IO_STD *) 0 ;
1768 gr2t(gridptr->pfile->grvals[TINDEX], (float) t, &dtim) ;
1769 gr2t(gridptr->pfile->grvals[TINDEX], 1.0, &dtimi) ;
1770 fn = gafndt(gridptr->pfile->name, &dtim, &dtimi, gridptr->pfile->abvals[TINDEX], NULL, 1) ;
1771 #ifdef TESTING3
1772 fprintf(stderr, "dtimi=(%d, %d, %d, %d, %d), dtim=(%d, %d, %d, %d, %d)\n",
1773 dtimi.yr, dtimi.mo, dtimi.dy, dtimi.hr, dtimi.mn, dtim.yr, dtim.mo, dtim.dy,
1774 dtim.hr, dtim.mn) ;
1775 fprintf(stderr, "New SDF filename=%s.\n", fn) ;
1776 #endif
1777 if (fn == NULL) {
1778 /* what could I do if close_netcdf did bomb? */
1779 (void) close_netcdf(oldsdf_ptr->cdfid) ;
1780 if (oldsdf_ptr != gridptr->pfile->first_sdf_ptr) {
1781 free_io_std(&(oldsdf_ptr)) ;
1782 } /* gotta keep that first one for the varlookups later */
1783 /* restore on failure */
1784 if (copy_io_std(&(gridptr->pfile->sdf_ptr),
1785 gridptr->pfile->first_sdf_ptr) == Failure) {
1786 /* SIGH. Two screwups at once! */
1787 gaprnt(1, "SDF templating: Error copying SDF info from first to current.\n") ;
1788 return(-99999) ; /* should be -88888, but... */
1789 }
1790 gaprnt(1, "SDF templating: Error getting new filename.\n") ;
1791 return(-99999) ;
1792 }
1793 if (!init_io_std(&(gridptr->pfile->sdf_ptr))) {
1794 gaprnt(0, "SDF templating: Allocation error opening new SDF file.\n") ;
1795 gridptr->pfile->fnumc = 0 ;
1796 /* what could I do if close_netcdf did bomb? */
1797 (void) close_netcdf(oldsdf_ptr->cdfid) ;
1798 if (oldsdf_ptr != gridptr->pfile->first_sdf_ptr) {
1799 free_io_std(&(oldsdf_ptr)) ;
1800 } /* gotta keep that first one for the varlookups later */
1801 /* restore on failure */
1802 gridptr->pfile->sdf_ptr = (IO_STD *) NULL ;
1803 if (copy_io_std(&(gridptr->pfile->sdf_ptr),
1804 gridptr->pfile->first_sdf_ptr) == Failure) {
1805 gaprnt(1, "SDF templating: Error copying first to current in switch.\n") ;
1806 return(-88888) ;
1807 }
1808 gaprnt(1, "SDF templating: Error moving to new file.\n") ;
1809 return(-88888) ;
1810 }
1811 if (!read_io_std(fn, gridptr->pfile->sdf_ptr)) {
1812 gaprnt(0, "SDF templating: Open of new SDF file failed.\n Filename is ") ;
1813 gaprnt(0, fn) ;
1814 gaprnt(0, ".\n") ;
1815 #if USEHDF == 0
1816 gaprnt(0, "If this was an HDF-SDS file, try gradshdf.\n") ;
1817 #endif
1818 gridptr->pfile->fnumc = 0 ;
1819 /* what could I do if close_netcdf did bomb? */
1820 (void) close_netcdf(oldsdf_ptr->cdfid) ;
1821 /* restore on failure */
1822 gridptr->pfile->sdf_ptr = oldsdf_ptr ;
1823 gaprnt(1, "SDF templating: New open failed.\n") ;
1824 return(-88888) ;
1825 }
1826 /* copy over the gradsabbr values into new file's structs */
1827 for (oldlclvar = oldsdf_ptr->var ; oldlclvar ;
1828 oldlclvar = oldlclvar->next) {
1829 /* find same netCDF var in new file */
1830 for (newlclvar = gridptr->pfile->sdf_ptr->var ; newlclvar ;
1831 newlclvar = newlclvar->next) {
1832 if (!strcmp(oldlclvar->varnam, newlclvar->varnam)) {
1833 /* match */
1834 if (oldlclvar->gradsabbr) {
1835 strcpy(newlclvar->gradsabbr, oldlclvar->gradsabbr) ;
1836 }
1837 for (idim = 0 ; idim < 4 ; ++idim) /* copy dim. maps */ {
1838 newlclvar->dimmap[idim] = oldlclvar->dimmap[idim] ;
1839 newlclvar->dimidmap[idim] = oldlclvar->dimidmap[idim] ;
1840 }
1841 /* Here's the rub: if this is an XDFopen-ed file, and that DDF had an UNDEF */
1842 /* statement, we don't want to look for it in an attribute value. However, */
1843 /* I don't know how to tell this origin info at this point; it was in the */
1844 /* parms structure during data discovery phase, (isxdf and needundef */
1845 /* components), but I don't see any way to access that here. Hoop, 2002/03/13 */
1846 /* Figured out a way to communicate it from DDF parsing code, first via parms */
1847 /* struct to data discovery fcn (first in this file), then via IO_STD struct to */
1848 /* here. hasDDFundef is new structure element in both cases. -Hoop, 2002/03/28 */
1849 /* */
1850 /* get missing_value attribute value, & stuff it into newlclvar->missing_value */
1851 /* Added here from the data_var discovery code in gadsdf(), 2001/04/15 -Hoop */
1852 if (gridptr->pfile->sdf_ptr->hasDDFundef) /* always float */ {
1853 newlclvar->missing_value.fval = gridptr->pfile->undef ;
1854 } else {
1855 missval_attr = find_att(newlclvar->first_vattr, "missing_value") ;
1856 if (!missval_attr) {
1857 missval_attr = find_att(newlclvar->first_vattr, "_FillValue") ;
1858 }
1859 /* if the type of the data variable and that of the missing_value attribute */
1860 /* don't match, use the default, as we do if we have no attribute */
1861 if ((!missval_attr) ||
1862 (missval_attr->type != newlclvar->vartype)) {
1863 switch (newlclvar->vartype) {
1864 case NC_BYTE:
1865 newlclvar->missing_value.bval = FILL_BYTE ;
1866 break ;
1867 case NC_SHORT:
1868 newlclvar->missing_value.sval = FILL_SHORT ;
1869 break ;
1870 case NC_LONG:
1871 newlclvar->missing_value.lval = FILL_LONG ;
1872 break ;
1873 case NC_FLOAT:
1874 newlclvar->missing_value.fval = FILL_FLOAT ;
1875 break ;
1876 case NC_DOUBLE:
1877 newlclvar->missing_value.dval = FILL_DOUBLE ;
1878 break ;
1879 } /* end switch */
1880 } else /* have an attribute value and types match */ {
1881 switch (newlclvar->vartype) {
1882 case NC_BYTE:
1883 newlclvar->missing_value.bval =
1884 *((signed char *)missval_attr->data) ;
1885 break ;
1886 case NC_SHORT:
1887 newlclvar->missing_value.sval =
1888 *((short *)missval_attr->data) ;
1889 break ;
1890 case NC_LONG:
1891 newlclvar->missing_value.lval =
1892 *((long *)missval_attr->data) ;
1893 break ;
1894 case NC_FLOAT:
1895 newlclvar->missing_value.fval =
1896 *((float *)missval_attr->data) ;
1897 break ;
1898 case NC_DOUBLE:
1899 newlclvar->missing_value.dval =
1900 *((double *)missval_attr->data) ;
1901 break ;
1902 } /* end switch */
1903 } /* else is have attribute and types match */
1904 }
1905 break ; /* found this oldlclvar amid new ones, go for next */
1906 } /* end if match found */
1907 } /* end for each new var */
1908 } /* end for each old var */
1909 /* what could I do if close_netcdf did bomb? */
1910 (void) close_netcdf(oldsdf_ptr->cdfid) ;
1911 /* must redetermine data_var */
1912 abbrvlen = (int) strlen(gridptr->pvar->abbrv) ;
1913 for ((*data_var_ptr) = gridptr->pfile->sdf_ptr->var ; (*data_var_ptr) ;
1914 (*data_var_ptr) = (*data_var_ptr)->next) {
1915 if ((*data_var_ptr)->gradsabbr[0] != '\0') {
1916 if (!strcasecmp(gridptr->pvar->abbrv, (*data_var_ptr)->gradsabbr)) {
1917 break ;
1918 }
1919 }
1920 }
1921 if (!(*data_var_ptr)) {
1922 sprintf(pout, "SDF templating: Couldn't find variable %s in SDF file.\n", gridptr->pvar->abbrv) ;
1923 gaprnt(0, pout) ;
1924 (void) close_netcdf(oldsdf_ptr->cdfid) ;
1925 gridptr->pfile->fnumc = 0 ;
1926 /* restore on failure */
1927 if (copy_io_std(&(gridptr->pfile->sdf_ptr),
1928 gridptr->pfile->first_sdf_ptr) == Failure) {
1929 gaprnt(1, "SDF templating: Error copying current to first after failure finding variable.\n") ;
1930 return(-88888) ;
1931 }
1932 return(-88888) ;
1933 }
1934 gridptr->pfile->fnumc = i ;
1935 gridptr->pfile->tempname = fn ;
1936 } /* if i != gridptr->pfile->fnumc */
1937 t = 1 + t - gridptr->pfile->fnumc ;
1938 return(t) ;
1939 }
1940
1941 #endif
1942 int
gadxdf(char * name,struct gafile * pfi,char ** template,int * ntsteps,char ** sdfname,GASDFPARMS * parms)1943 gadxdf(char *name, struct gafile *pfi, char **template, int *ntsteps,
1944 char **sdfname, GASDFPARMS *parms) {
1945 struct gavar *pvar;
1946 struct dt tdef,tdefi,dt1,dt2;
1947 struct gaindx *pindx;
1948 float *vals;
1949 int size,rc,len,swpflg,cnt, idum, diag=0, ichar;
1950 char rec[512], mrec[512], *ch, *dn, *pos, *sname, *misc ; /*mf mf*/
1951 int flgs[1],cflg,i,j,err,hdrb,trlb,mflflg,crayflg=0;
1952 int mcnt,maxlv,maxct, mflag = 0 ;
1953 int levs,acum,acumvz,fpos,recacm;
1954 float temp,v1,v2, fdum;
1955 FILE *mfile = (FILE *) 0 ;
1956 static char *errs[1] = {"DSET"} ;
1957 int nzstride=0,first=1;
1958 int levsua=0,acumstride=0;
1959
1960 #if GRADS_CDUNIF == 1
1961 mfcmn.fullyear=1; /*mf --- define here vice grads.c for cdunif.c mf*/
1962 #endif
1963
1964 char *putncdnm(char *ch, char **dimname) ;
1965 void dupstrval(char *src, char **dest) ;
1966
1967 hdrb = 0;
1968 trlb = 0;
1969 mflflg = 0; /* map file not open */
1970 pfi->mfile = NULL;
1971 mcnt = -1;
1972
1973 /* Initialize variables */
1974 parms->xsrch = parms->ysrch = parms->zsrch = parms->tsrch = 1 ;
1975 parms->dvsrch = 1 ; parms->isxdf = 1 ;
1976 parms->xsetup = parms->ysetup = parms->zsetup = parms->tsetup = 1 ;
1977 parms->needtitle = parms->needundef = 1 ;
1978 parms->dvcount = -1 ;
1979 parms->xdimname = parms->ydimname = parms->zdimname = parms->tdimname = (char *) 0 ;
1980 parms->dvncnames = parms->dvganames = (char **) 0 ;
1981 parms->dvsetup = (int *) 0 ;
1982 parms->hasDDFundef = 0 ;
1983 sname=NULL;
1984
1985 /* Open descriptor file */
1986 descr = fopen (name, "r");
1987 if (descr == NULL) {
1988 /* Add default suffix of .ctl */
1989 sname = (char *)malloc(strlen(name)+5);
1990 if(sname == NULL) {
1991 gaprnt(0,"XDFopen: malloc error in creating data descriptor file name\n");
1992 return(1);
1993 }
1994 for(i=0;i<=strlen(name);i++) *(sname+i)=*(name+i);
1995 strcat(sname,".ctl");
1996 descr = fopen (sname, "r");
1997 }
1998
1999 if (descr == NULL) {
2000 gaprnt (0,"XDFOpen Error: Can't open description file\n");
2001 if (sname) {
2002 free(sname);
2003 }
2004 return(1);
2005 }
2006
2007 /* Copy descriptor file name into gafile structure */
2008 if(sname != NULL) {
2009 getwrd (pfi->dnam,sname,255);
2010 if (sname) {
2011 free(sname);
2012 }
2013 } else {
2014 getwrd (pfi->dnam,name,255);
2015 }
2016
2017 /* initialize error flags */
2018 for (i=0;i<1;i++) flgs[i] = 1;
2019
2020 /*mf initialize the calendar to standard and cray_ieee to 0 mf*/
2021 /* pfi->calendar=0; */ /* Mike stopped doing it here in gaddes.c, so me too */
2022 pfi->cray_ieee=0;
2023
2024 /* Parse the descriptor file */
2025 pfi->vnum = 0 ;
2026 while (fgets(rec,512,descr)!=NULL) {
2027 for (ichar = 0 ; ichar < strlen(rec) ; ++ichar) {
2028 if (rec[ichar] == '\n') {
2029 rec[ichar] = '\0' ;
2030 }
2031 }
2032 strcpy (mrec,rec);
2033 lowcas(rec);
2034
2035 pfi->fhdr = 0;
2036 pfi->xyhdr = 0;
2037 pfi->thdr = 0;
2038
2039 if ( (cmpwrd("*", rec) || cmpwrd("#", rec)) || !isalnum(rec[0]) ) {
2040 cflg = 1;
2041 /* Parse comment if it contains attribute metadata */
2042 if (strncmp("*:attr",mrec,6)==0) {
2043 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
2044 }
2045
2046 } else if (cmpwrd("options",rec)) {
2047 if ( (ch=nxtwrd(rec))==NULL ) {
2048 gaprnt (1,"XDFopen: Description file warning: Missing options keyword\n");
2049 } else {
2050 while (ch != NULL) {
2051 if (cmpwrd("yrev",ch)) pfi->yrflg = 1;
2052 else if (cmpwrd("zrev",ch)) pfi->zrflg = 1;
2053 else if (cmpwrd("template",ch)) pfi->tmplat = 1;
2054 else if (cmpwrd("365_day_calendar",ch)) {
2055 pfi->calendar = 1 ;
2056 mfcmn.cal365 = pfi->calendar ;
2057 } else {
2058 gaprnt (0,"XDFopen Error: Data file option invalid\n");
2059 goto err9;
2060 }
2061 ch = nxtwrd(ch);
2062 }
2063 }
2064
2065 } else if (cmpwrd("title",rec)) {
2066 parms->needtitle = 0 ;
2067 if ( (ch=nxtwrd(mrec))==NULL ) {
2068 gaprnt (1,"XDFopen: Description file warning: Missing title string\n");
2069 pfi->title[0] = '\0' ;
2070 } else {
2071 getstr (pfi->title,ch,511);
2072 }
2073
2074 } else if (cmpwrd("dset",rec)) {
2075 ch = nxtwrd(mrec);
2076 if (ch==NULL) {
2077 gaprnt (0,"XDFOpen Error: Data file name is missing\n");
2078 goto err9;
2079 }
2080 if (*ch=='^' || *ch=='$') {
2081 fnmexp (pfi->name,ch,name);
2082 } else {
2083 getwrd (pfi->name,ch,511);
2084 }
2085 flgs[0] = 0;
2086
2087 } else if (cmpwrd("undef",rec)) {
2088 ch = nxtwrd(rec);
2089 if (ch==NULL) {
2090 gaprnt (0,"XDFOpen Error: Missing undef value\n");
2091 goto err9;
2092 }
2093 pos = valprs(ch,&(pfi->undef));
2094 if (pos==NULL) {
2095 gaprnt (0,"XDFOpen Error: Invalid undef value\n");
2096 goto err9;
2097 }
2098 if (SETMISS) {
2099 pfi->ulow = fabs(pfi->undef/EPSILON);
2100 pfi->uhi = pfi->undef + pfi->ulow;
2101 pfi->ulow = pfi->undef - pfi->ulow;
2102 }
2103 parms->needundef = 0 ;
2104 parms->hasDDFundef = 1 ;
2105
2106 } else if (cmpwrd("xdef",rec)) {
2107 if (pfi->type == 2) continue;
2108 if ( (ch = nxtwrd(mrec)) == NULL) goto err0; /* xdimname must be mixed case version*/
2109 parms->xsrch = 0 ;
2110 ch = putncdnm(ch, &(parms->xdimname)) ; /* mallocs and strcpys */
2111 ch = nxtwrd(rec) ; /* skip over xdef in lowcase version */
2112 if ((ch = nxtwrd(ch)) == NULL) {
2113 parms->xsetup = 1 ;
2114 } else {
2115 if ( (pos = intprs(ch,&(pfi->dnum[0])))==NULL) {
2116 goto err1 ;
2117 }
2118 if (*pos != ' ') goto err1;
2119 if ( (ch = nxtwrd(ch))==NULL) goto err2;
2120 if (cmpwrd("linear",ch)) {
2121 rc = deflin(ch, pfi, 0, 0);
2122 if (rc==-1) goto err8;
2123 if (rc) goto err9;
2124 v2 = *(pfi->grvals[0]);
2125 v1 = *(pfi->grvals[0]+1) + v2;
2126 temp = v1+((float)(pfi->dnum[0]))*v2;
2127 temp=temp-360.0;
2128 if (fabs(temp-v1)<0.01) pfi->wrap = 1;
2129 } else if (cmpwrd("levels",ch)) {
2130 if (pfi->dnum[0]<1) goto err7;
2131 rc = deflev (ch, rec, pfi, 0);
2132 if (rc==-1) goto err8;
2133 if (rc) goto err9;
2134 } else goto err2;
2135 parms->xsetup = 0 ;
2136 }
2137
2138 } else if (cmpwrd("ydef",rec)) {
2139 if (pfi->type == 2) continue;
2140 if ( (ch = nxtwrd(mrec)) == NULL) goto err0; /* ydimname must be mixed case version*/
2141 parms->ysrch = 0 ;
2142 ch = putncdnm(ch, &(parms->ydimname)) ; /* mallocs and strcpys */
2143 ch = nxtwrd(rec) ; /* skip over ydef in lowcase version*/
2144 if ((ch = nxtwrd(ch)) == NULL) { ;
2145 parms->ysetup = 1 ;
2146 } else {
2147 if ( (pos = intprs(ch,&(pfi->dnum[1])))==NULL) {
2148 goto err1 ;
2149 }
2150 if (*pos!=' ') goto err1;
2151 if ( (ch = nxtwrd(ch))==NULL) goto err2;
2152 if (cmpwrd("linear",ch)) {
2153 rc = deflin(ch, pfi, 1, 0);
2154 if (rc==-1) goto err8;
2155 if (rc) goto err9;
2156 } else if (cmpwrd("levels",ch)) {
2157 if (pfi->dnum[1]<1) goto err7;
2158 rc = deflev (ch, rec, pfi, 1);
2159 if (rc==-1) goto err8;
2160 if (rc) goto err9;
2161 } else if (cmpwrd("gausr40",ch)) {
2162 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2163 if ( (pos = intprs(ch,&i))==NULL) goto err3;
2164 pfi->grvals[1] = gagaus(i,pfi->dnum[1]);
2165 if (pfi->grvals[1]==NULL) goto err9;
2166 pfi->abvals[1] = pfi->grvals[1];
2167 pfi->ab2gr[1] = lev2gr;
2168 pfi->gr2ab[1] = gr2lev;
2169 pfi->linear[1] = 0;
2170 } else if (cmpwrd("mom32",ch)) {
2171 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2172 if ( (pos = intprs(ch,&i))==NULL) goto err3;
2173 pfi->grvals[1] = gamo32(i,pfi->dnum[1]);
2174 if (pfi->grvals[1]==NULL) goto err9;
2175 pfi->abvals[1] = pfi->grvals[1];
2176 pfi->ab2gr[1] = lev2gr;
2177 pfi->gr2ab[1] = gr2lev;
2178 pfi->linear[1] = 0;
2179 } else if (cmpwrd("gausr30",ch)) {
2180 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2181 if ( (pos = intprs(ch,&i))==NULL) goto err3;
2182 pfi->grvals[1] = gags30(i,pfi->dnum[1]);
2183 if (pfi->grvals[1]==NULL) goto err9;
2184 pfi->abvals[1] = pfi->grvals[1];
2185 pfi->ab2gr[1] = lev2gr;
2186 pfi->gr2ab[1] = gr2lev;
2187 pfi->linear[1] = 0;
2188 } else if (cmpwrd("gausr20",ch)) {
2189 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2190 if ( (pos = intprs(ch,&i))==NULL) goto err3;
2191 pfi->grvals[1] = gags20(i,pfi->dnum[1]);
2192 if (pfi->grvals[1]==NULL) goto err9;
2193 pfi->abvals[1] = pfi->grvals[1];
2194 pfi->ab2gr[1] = lev2gr;
2195 pfi->gr2ab[1] = gr2lev;
2196 pfi->linear[1] = 0;
2197 } else if (cmpwrd("gausr15",ch)) {
2198 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2199 if ( (pos = intprs(ch,&i))==NULL) goto err3;
2200 pfi->grvals[1] = gags15(i,pfi->dnum[1]);
2201 if (pfi->grvals[1]==NULL) goto err9;
2202 pfi->abvals[1] = pfi->grvals[1];
2203 pfi->ab2gr[1] = lev2gr;
2204 pfi->gr2ab[1] = gr2lev;
2205 pfi->linear[1] = 0;
2206 } else goto err2;
2207 parms->ysetup = 0 ;
2208 }
2209
2210 } else if (cmpwrd("zdef",rec)) {
2211 if (pfi->type == 2) continue;
2212 if ( (ch = nxtwrd(mrec)) == NULL) goto err0; /* get mixed case version */
2213 parms->zsrch = 0 ;
2214 ch = putncdnm(ch, &(parms->zdimname)) ;
2215 ch = nxtwrd(rec) ; /* point past zdef in lowcased version */
2216 if ((ch = nxtwrd(ch)) == NULL) {
2217 parms->zsetup = 1 ;
2218 } else {
2219 if ( (pos = intprs(ch,&(pfi->dnum[2])))==NULL) {
2220 goto err1 ;
2221 }
2222 if (*pos!=' ') goto err1;
2223 if ( (ch = nxtwrd(ch))==NULL) goto err2;
2224 if (cmpwrd("linear",ch)) {
2225 rc = deflin(ch, pfi, 2, 0);
2226 if (rc==-1) goto err8;
2227 if (rc) goto err9;
2228 } else if (cmpwrd("levels",ch)) {
2229 if (pfi->dnum[2]<1) goto err7;
2230 rc = deflev (ch, rec, pfi, 2);
2231 if (rc==-1) goto err8;
2232 if (rc) goto err9;
2233 } else goto err2;
2234 parms->zsetup = 0 ;
2235 }
2236
2237 } else if (cmpwrd("tdef",rec)) {
2238 if ( (ch = nxtwrd(mrec)) == NULL) goto err0; /* get mixed case version */
2239 parms->tsrch = 0 ;
2240 if (!strncasecmp(ch, "%nodim%", 7)) {
2241 parms->tdimname = (char *) 0 ; /* we won't be using any tdimname */
2242 pfi->dnum[TINDEX] = 1 ; /* 1 time step ; be sure not to map any tdim */
2243 } else {
2244 ch = putncdnm(ch, &(parms->tdimname)) ;
2245 }
2246 ch = nxtwrd(rec) ; /* skip over tdef in lowcased version */
2247 if ((ch = nxtwrd(ch)) == NULL) {
2248 if (!(parms->tdimname)) {
2249 vals = (float *)malloc(sizeof(float)*8);
2250 if (vals==NULL) goto err8;
2251 vals[7] = -999.9 ;
2252 pfi->grvals[TINDEX] = vals ;
2253 pfi->abvals[TINDEX] = vals ;
2254 pfi->linear[TINDEX] = 1 ;
2255 vals[0] = 1.0 ;
2256 vals[1] = 1.0 ;
2257 vals[2] = 1.0 ;
2258 vals[3] = 0.0 ; /* initial hours */
2259 vals[4] = 0.0 ;
2260 vals[5] = 0.0 ; /* step in months */
2261 vals[6] = 1.0 ; /* step in minutes */
2262 parms->tsetup = 0 ;
2263 } else {
2264 parms->tsetup = 1 ;
2265 }
2266 } else {
2267 if ( (pos = intprs(ch,&(pfi->dnum[3])))==NULL) {
2268 goto err1 ;
2269 }
2270 if (!(parms->tdimname)) {
2271 if ((pfi->tmplat) == 0) {
2272 /* non-templating, %nodim% case can only have 1 timestep */
2273 if (pfi->dnum[3] != 1) {
2274 gaprnt(0, "TDEF with %nodim% has timestep count != 1; resetting to 1.\n") ;
2275 pfi->dnum[3] = 1 ;
2276 }
2277 }
2278 }
2279 if (*pos!=' ') goto err1;
2280 if ( (ch = nxtwrd(ch))==NULL) goto err2;
2281 if (cmpwrd("linear",ch)) {
2282 if ( (ch = nxtwrd(ch))==NULL) goto err3;
2283 tdef.yr = -1000;
2284 tdef.mo = -1000;
2285 tdef.dy = -1000;
2286 if ( (pos = adtprs(ch,&tdef,&dt1))==NULL) goto err3;
2287 if (*pos!=' ' || dt1.yr == -1000 || dt1.mo == -1000.0 ||
2288 dt1.dy == -1000) goto err3;
2289 if ( (ch = nxtwrd(ch))==NULL) goto err4;
2290 if ( (pos = rdtprs(ch,&dt2))==NULL) goto err4;
2291 v1 = (dt2.yr * 12) + dt2.mo;
2292 v2 = (dt2.dy * 1440) + (dt2.hr * 60) + dt2.mn;
2293 /*mf --- check if 0 dt ---mf*/
2294 if ((v1 == 0) && (v2 == 0)) goto err4a ;
2295 vals = (float *)malloc(sizeof(float)*8);
2296 if (vals==NULL) goto err8;
2297 *(vals) = dt1.yr;
2298 *(vals+1) = dt1.mo;
2299 *(vals+2) = dt1.dy;
2300 *(vals+3) = dt1.hr;
2301 *(vals+4) = dt1.mn;
2302 *(vals+5) = v1;
2303 *(vals+6) = v2;
2304 *(vals+7) = -999.9;
2305 pfi->grvals[3] = vals;
2306 pfi->abvals[3] = vals;
2307 pfi->linear[3] = 1;
2308 } else goto err2;
2309 parms->tsetup = 0 ;
2310 }
2311
2312 } else if (cmpwrd("vars",rec)) {
2313 if ( (ch = nxtwrd(rec)) == NULL) goto err5;
2314 if ( (pos = intprs(ch,&(pfi->vnum)))==NULL) goto err5;
2315 parms->dvsrch = 0 ;
2316 parms->dvcount = pfi->vnum ;
2317 size = pfi->vnum * sizeof(struct gavar) ;
2318 pvar = (struct gavar *)malloc(size);
2319 pfi->pvar1 = pvar;
2320 parms->dvncnames = (char **) malloc(pfi->vnum * sizeof(char *)) ;
2321 parms->dvganames = (char **) malloc(pfi->vnum * sizeof(char *)) ;
2322 parms->dvsetup = (int *) malloc(pfi->vnum * sizeof(int)) ;
2323 i = 0;
2324 while (i<pfi->vnum) {
2325 if (fgets(rec,512,descr)==NULL) {
2326 gaprnt (0,"XDFOpen Error: Unexpected EOF reading variables\n");
2327 sprintf (pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
2328 gaprnt (2,pout);
2329 goto retrn;
2330 }
2331 /* replace newline with null at end of record */
2332 for (ichar = strlen(rec) - 1 ; ichar >= 0 ; --ichar) {
2333 if (rec[ichar] == '\n') {
2334 rec[ichar] = '\0' ;
2335 break ;
2336 }
2337 }
2338 /* Allow comments between VARS and ENDVARS */
2339 strcpy (mrec,rec);
2340 lowcas(rec);
2341 if ( (strncmp(mrec,"*",1)==0) || (strncmp(mrec,"#",1)==0) || !isalnum(*(mrec)) ) {
2342 /* Parse comment if it contains attribute metadata */
2343 if (strncmp("*:attr",mrec,6)==0) {
2344 rc = ddfattr (mrec, pfi);
2345 if (rc == -1) goto retrn;
2346 else continue;
2347 }
2348 else continue;
2349 }
2350 if (cmpwrd("endvars",rec)) {
2351 gaprnt (0,"XDFOpen Error: Unexpected ENDVARS record\n");
2352 sprintf (pout, "Was expecting %i records. Found %i.\n", pfi->vnum, i);
2353 gaprnt (2,pout);
2354 goto err9;
2355 }
2356
2357 /* Get the compound variable name (ncname=>ganame). */
2358 pvar->longnm = NULL;
2359 if ((getvnm(pvar, mrec))!=0) goto err6;
2360
2361 /* copy abbrv to ganame */
2362 parms->dvganames[i] = (char *) malloc( strlen(pvar->abbrv)*sizeof(char)) ;
2363 strcpy(parms->dvganames[i],pvar->abbrv);
2364
2365 /* copy longnm to ncname if it exists, otherwise copy abbrv */
2366 if (pvar->longnm) {
2367 parms->dvncnames[i] = (char *) malloc( strlen(pvar->longnm)*sizeof(char)) ;
2368 strcpy(parms->dvncnames[i],pvar->longnm);
2369 } else {
2370 parms->dvncnames[i] = (char *) malloc( strlen(pvar->abbrv)*sizeof(char)) ;
2371 strcpy(parms->dvncnames[i],pvar->abbrv);
2372 }
2373
2374 if ( (ch=nxtwrd(rec))==NULL) goto err6; /* advance pointer to levels field */
2375
2376 if ((ch == NULL) || (*ch == '\0')) {
2377 parms->dvsetup[i] = 1 ; /* will figure out lvls, etc., in gadsdf */
2378 } else {
2379 parms->dvsetup[i] = 0 ;
2380 /* get the number of vertical levels */
2381 if ( (pos=intprs(ch,&(pvar->levels)))==NULL ) goto err6;
2382 if ( (ch=nxtwrd(ch))==NULL) goto err6;
2383 /* parse the units field */
2384 for (j=0;j<4;j++) pvar->units[j] = -999;
2385 j = 0;
2386 while (1) {
2387 if ( (ch=intprs(ch,&(pvar->units[j])))==NULL ) goto err6;
2388 while ((*ch==' ') || (*ch=='\t')) ch++;
2389 if (*ch=='\0' || *ch=='\n') goto err6;
2390 if (*ch!=',') break;
2391 ch++;
2392 while ((*ch==' ') || (*ch=='\t')) ch++;
2393 if (*ch=='\0' || *ch=='\n') goto err6;
2394 j++;
2395 if (j>3) goto err6;
2396 }
2397 /* points into mixcase version where we left off in lowcased version */
2398 ch = mrec + (ch - rec) ;
2399 getstr (pvar->varnm,ch,127);
2400
2401 /* initialize the var_z counter for NASA GLA format */
2402 pvar->var_z = 1;
2403 if(pvar->units[0] == -1 && pvar->units[1] == 10 ) {
2404 nzstride++;
2405 }
2406
2407 /* var_t is for DRS var-t transforms */
2408 pvar->var_t = 0;
2409 if(pvar->units[0] == -1 && pvar->units[1] == 20 ) pvar->var_t = 1;
2410
2411 /* x-y transpose for lat/lon vice lon/lat data VERY INEFFICIENT!!! */
2412 pvar->y_x = 0;
2413 if(pvar->units[0] == -1 && pvar->units[1] == 30) pvar->y_x = 1;
2414
2415 /* Non-float data types */
2416 pvar->dfrm = 0;
2417 if ((pvar->units[0] == -1 && pvar->units[1] == 40) &&
2418 (pvar->units[2] >= 1 && pvar->units[2] <= 4)) {
2419 pvar->dfrm= pvar->units[2];
2420 if ( pvar->units[2] == 2 ) pvar->dfrm=2;
2421 if ( pvar->units[2] == 2 && pvar->units[3] == -1) pvar->dfrm=-2;
2422 }
2423
2424 #if GRADS_CRAY == 1
2425 /* 32-bit big endian ieee to cray float */
2426 if( (pvar->units[0]==-1 && pvar->units[1]==50) || crayflg ) pvar->dfrm=8;
2427 #endif
2428
2429 }
2430 i++; pvar++;
2431 }
2432 if (fgets(rec,512,descr)==NULL) {
2433 gaprnt (0,"XDFOpen Error: Missing ENDVARS statement.\n");
2434 goto retrn;
2435 }
2436
2437 /* See if record is an attribute comment, otherwise send error message */
2438 strcpy (mrec,rec);
2439 lowcas(rec);
2440 while (!cmpwrd("endvars",rec)) {
2441 if (strncmp("*:attr",mrec,6)==0) {
2442 if ((ddfattr(mrec,pfi)) == -1) goto retrn;
2443 }
2444 else {
2445 sprintf(pout,"Open Error: Looking for \"endvars\", found \"%s\" instead.\n",rec);
2446 gaprnt (0,pout);
2447 goto err9;
2448 }
2449 if (fgets(rec,256,descr)==NULL) {
2450 gaprnt (0,"Open Error: Missing ENDVARS statement.\n");
2451 goto retrn;
2452 }
2453 }
2454
2455 } else {
2456 /* Parse error of descriptor file */
2457 gaprnt (0,"XDFOpen Error: Unknown keyword in description file\n");
2458 goto err9;
2459 }
2460
2461 }
2462
2463 err=0;
2464 for (i=0; i<1; i++) {
2465 if (flgs[i]) {
2466 sprintf (pout,"XDFOpen Error: missing %s record \n", errs[i]);
2467 gaprnt (0,pout);
2468 err=1;
2469 }
2470 }
2471
2472 if (err) goto retrn;
2473
2474 if (pfi->type>1 && mflag) {
2475 if (mcnt==-1) {
2476 gaprnt (0,"XDFOpen Error: missing STNMAP record\n");
2477 err=1;
2478 } else if (mcnt != pfi->dnum[3]) {
2479 gaprnt (0,"XDFOpen Error: Inconsistent time count\n");
2480 sprintf (pout," Count in station map file = %li\n",mcnt);
2481 gaprnt (0,pout);
2482 sprintf (pout," Count in descriptor file = %i\n",pfi->dnum[3]);
2483 gaprnt (0,pout);
2484 err=1;
2485 }
2486 }
2487
2488 if (err) goto retrn;
2489
2490 /* Figure out locations of variables within a time group */
2491 if (pfi->vnum > 0) {
2492 pvar = pfi->pvar1;
2493 }
2494 /* Grid data */
2495 if (pfi->type==1) {
2496 if (((!(parms->xsrch)) && (!(parms->xsetup))) &&
2497 ((!(parms->ysrch)) && (!(parms->ysetup))) ) {
2498 pfi->gsiz = pfi->dnum[0] * pfi->dnum[1];
2499 }
2500
2501 /*mf
2502 add a constant xy header
2503 mf*/
2504 if (pfi->xyhdr) {
2505 pfi->gsiz = pfi->gsiz + pfi->xyhdr;
2506 }
2507
2508 if (pfi->seqflg) {
2509 pfi->gsiz+=2;
2510 if (hdrb>0) hdrb+=2;
2511 pvar->offset = 1+hdrb;
2512 acum = 1+hdrb;
2513 } else {
2514 if (pfi->vnum > 0) {
2515 pvar->offset = hdrb;
2516 }
2517 acum = hdrb;
2518 }
2519 if (pfi->vnum > 0) {
2520 levs = pvar->levels;
2521 if (levs==0) levs=1;
2522 pvar->recoff = 0;
2523 recacm = 0;
2524 pvar++;
2525 }
2526
2527 acumvz=acum;
2528
2529 for (i=1; i<pfi->vnum; i++) {
2530
2531 /* NASA GLA FORMAT CHECKS */
2532 /* upper air fields which var and z are transposed*/
2533
2534 if( (pvar->units[0]==-1) &&
2535 (pvar->units[1]==10) &&
2536 (pvar->units[2]==1) ) {
2537
2538 levsua = pvar->levels;
2539 acum = acum + pfi->gsiz;
2540 pvar->var_z = nzstride;
2541
2542 /* diagnstotic fields AFTER upper air fields which are in GrADS normal order */
2543
2544 } else if( (pvar->units[0]==-1) &&
2545 (pvar->units[1]==10) &&
2546 (pvar->units[2]==2) ) {
2547 if(first) {
2548 acumstride = acumstride + nzstride*levsua*pfi->gsiz;
2549 acum = acumstride + (pvar->levels*pfi->gsiz) ;
2550 first = 0 ;
2551 } else {
2552 acum = acum + (levs*pfi->gsiz);
2553 }
2554
2555 } else if(pvar->var_t) { /* DRS transposition of var and t */
2556
2557 if(pfi->tmplat) { /* time template, read the third unit param
2558 for the size of each chunk in a file */
2559
2560 if(pvar->units[2] != -999) {
2561 acum = acum + levs*(pfi->gsiz)*(pvar->units[1]);
2562 } else {
2563 gaprnt (0,"Using time templat and 4-D variables, # times / file not specified\n");
2564 gaprnt (0,"Defaulting to the number of times in the .ctl file\n");
2565 acum = acum + levs*(pfi->gsiz)*(pfi->dnum[3]);
2566 }
2567 } else {
2568 acum = acum + levs*(pfi->gsiz)*(pfi->dnum[3]);
2569 }
2570
2571 } else { /* simple GrADS */
2572
2573 acum = acum + (levs*pfi->gsiz);
2574 acumstride = acum ;
2575
2576 }
2577
2578 recacm += levs;
2579 pvar->offset = acum;
2580 pvar->recoff = recacm;
2581 levs = pvar->levels;
2582 if (levs==0) levs=1;
2583 pvar++;
2584
2585 }
2586
2587 recacm += levs;
2588
2589 /*mf 960514 correct for case where the last variable is a NASA UA */
2590
2591 if (pfi->vnum > 0) {
2592 pvar--;
2593 if( (pvar->units[0]==-1) &&
2594 (pvar->units[1]==10) &&
2595 (pvar->units[2]==1) ) {
2596 acum = acumvz + recacm*pfi->gsiz;
2597 }else {
2598 acum = acum + (levs*pfi->gsiz);
2599 }
2600 }
2601 /*mf --------------mf*/
2602 /*
2603 time chunk header; the default is 0
2604 */
2605 if(pfi->seqflg && pfi->thdr>0) {
2606 pfi->thdr+=2;
2607 }
2608 pfi->tsiz = acum + pfi->thdr;
2609
2610 /*mf --------------mf*/
2611
2612 pfi->trecs = recacm;
2613 if (pfi->seqflg) pfi->tsiz-=1;
2614 pfi->tsiz += trlb;
2615
2616 /*------------------- non grid data???? -------- */
2617
2618 } else {
2619 if (!(parms->dvsrch)) {
2620 pfi->ivnum = pfi->lvnum = 0 ;
2621 for (i=0; i<pfi->vnum; i++) {
2622 if (!(parms->dvsetup[i])) /* catch rest in gadsdf */ {
2623 if (pvar->levels!=0) {
2624 ++(pfi->ivnum) ;
2625 } else {
2626 ++(pfi->lvnum) ;
2627 }
2628 }
2629 pvar->offset = 0;
2630 pvar++;
2631 } /* for */
2632 } /* if we found dvars here */
2633 }
2634
2635 /*mf---
2636
2637 961204
2638
2639 set the global calendar and check if we are trying to change with a new file...
2640 we do this here to set the calandar for templating
2641
2642 ----mf*/
2643
2644 if(mfcmn.cal365<0) {
2645 mfcmn.cal365=pfi->calendar;
2646 } else {
2647 if(pfi->calendar != mfcmn.cal365) {
2648 gaprnt(0,"Attempt to change the global calendar...\n");
2649 if(mfcmn.cal365) {
2650 gaprnt(0,
2651 "The calendar is NOW 365 DAYS and you attempted to open a standard calendar file\n");
2652 } else {
2653 gaprnt(0,
2654 "The calendar is NOW STANDARD and you attempted to open a 365-day calendar file\n");
2655 }
2656 goto retrn;
2657 }
2658 }
2659
2660
2661 /* if time series templating was specified, & the TDEF line was incomplete, */
2662 /* ERROR! */
2663 if (pfi->tmplat && parms->tsetup) {
2664 goto err1 ;
2665 }
2666 /* If the file name is a time series template, figure out
2667 which times go with which files, so we don't waste a lot
2668 of time later opening and closing files unnecessarily. */
2669
2670 if (pfi->tmplat) {
2671 *ntsteps = pfi->dnum[3] ;
2672 dupstrval(pfi->name, template) ;
2673 pfi->fnums = (int *)malloc(sizeof(int)*pfi->dnum[3]);
2674 if (pfi->fnums==NULL) goto err8;
2675 j = 1;
2676 gr2t(pfi->grvals[3],1.0,&tdefi);
2677 ch = gafndt(pfi->name,&tdefi,&tdefi,pfi->abvals[3],NULL,1);
2678 if (ch==NULL) goto err8;
2679 dupstrval(ch, sdfname) ;
2680 *(pfi->fnums) = j;
2681 for (i=2; i<=pfi->dnum[3]; i++) {
2682 gr2t(pfi->grvals[3],(float)i,&tdef);
2683 pos = gafndt(pfi->name,&tdef,&tdefi,pfi->abvals[3],NULL,1);
2684 if (pos==NULL) goto err8;
2685 if (strcmp(ch,pos)!=0) {
2686 j = i;
2687 if (ch) {
2688 free(ch);
2689 }
2690 ch = pos;
2691 }
2692 *(pfi->fnums+i-1) = j;
2693 }
2694 if (ch) {
2695 free(ch);
2696 }
2697 pfi->fnumc = 0;
2698 } else {
2699 dupstrval(pfi->name, sdfname) ;
2700 }
2701
2702 fclose (descr);
2703 return(0);
2704
2705 errm:
2706 gaprnt(0,"XDFopen Error: Invalid pdef record.\n");
2707 pfi->ppflag = 0;
2708 goto err9;
2709
2710 err0:
2711 gaprnt(0, "XDFopen Error: Missing or invalid dimension name.\n") ;
2712 goto err9;
2713
2714 err1:
2715 gaprnt (0,"XDFopen Error: Missing or invalid dimension size.\n");
2716 goto err9;
2717
2718 err2:
2719 gaprnt (0,"XDFopen Error: Missing or invalid dimension");
2720 gaprnt (0," scaling type\n");
2721 goto err9;
2722
2723 err3:
2724 gaprnt (0,"XDFopen Error: Missing or invalid dimension");
2725 gaprnt (0," starting value\n");
2726 goto err9;
2727
2728 err4:
2729 gaprnt (0,"XDFopen Error: Missing or invalid dimension");
2730 gaprnt (0," increment value\n");
2731 goto err9;
2732
2733 err4a:
2734 gaprnt (0,"XDFopen Error: 0 time increment in tdef\n");
2735 gaprnt (0," use 1 for single time data\n");
2736 goto err9;
2737
2738 err5:
2739 gaprnt (0,"XDFopen Error: Missing or invalid variable");
2740 gaprnt (0," count\n");
2741 goto err9;
2742
2743 err6:
2744 gaprnt (0,"XDFOpen Error: Invalid variable record\n");
2745 goto err9;
2746
2747 err7:
2748 gaprnt (0,"XDFopen Error: Invalid number of levels\n");
2749 goto err9;
2750
2751 err8:
2752 gaprnt (0,"XDFopen Error: Memory allocation Error\n");
2753 goto retrn;
2754
2755 err9:
2756 gaprnt (0," --> The invalid description file record is: \n");
2757 gaprnt (0," --> ");
2758 gaprnt (0,rec);
2759 gaprnt (0,"\n");
2760
2761 retrn:
2762 gaprnt (0," The data file was not opened. \n");
2763 fclose (descr);
2764 if (mflflg) fclose(mfile);
2765 return(1);
2766
2767 }
2768
2769
2770 char *
putncdnm(char * ch,char ** dimname)2771 putncdnm(char *ch, char **dimname) {
2772 int lenstr, lendimname ;
2773
2774 lenstr = strlen(ch) ;
2775 for (lendimname = 0 ; (lendimname < lenstr) && (ch[lendimname] != ' ') ;
2776 ++lendimname) ;
2777 if (!(*dimname = (char *) malloc((lendimname + 1) * sizeof(char)))) {
2778 gaprnt(0, "Malloc err!\n") ;
2779 exit(13) ;
2780 }
2781 strncpy(*dimname, ch, lendimname) ;
2782 (*dimname)[lendimname] = '\0' ;
2783 for (ch += lendimname ; (*ch == ' ') || (*ch == '\t') ; ++ch) ;
2784 return(ch) ;
2785 }
2786
2787 void
dupstrval(char * src,char ** dest)2788 dupstrval(char *src, char **dest) {
2789 int lenstr ;
2790
2791 lenstr = strlen(src) ;
2792 if (!(*dest = (char *) malloc((lenstr + 1) * sizeof(char)))) {
2793 gaprnt(0, "Malloc err!\n") ;
2794 exit(13) ;
2795 }
2796 strcpy(*dest, src) ;
2797 return ;
2798 }
2799
2800 #ifndef STNDALN
2801
2802 /* Open an XDF data set by parsing the descriptor and reading the
2803 self-describing file for that data set and help to fill-in a gafile
2804 structure as need be. Chain the gafile structure on to the list anchored
2805 in the gastat. */
2806
2807 int
gaxdfopen(char * args,struct gacmn * pcm)2808 gaxdfopen (char *args, struct gacmn *pcm) {
2809 struct gafile *pfi, *pfio;
2810 int size, rc, ntsteps, ichpos ;
2811 char pathname[256], *optargs, *template, *cntsteps, *fnpart, *aux, *faux ;
2812 char *sdfpath ;
2813 int i, idummy ;
2814 int gadxdf(char *path, struct gafile *pfi, char **template, int *ntsteps,
2815 char **sdfpath, GASDFPARMS *parmsptr) ;
2816 int gadsdf(char *name, struct gafile *pfi, char *template, int ntsteps,
2817 GASDFPARMS parms) ;
2818 GASDFPARMS parms ;
2819 int init_io_std(IO_STD **ptr) ;
2820
2821 getwrd(pathname, args, 255) ;
2822 optargs = nxtwrd(args) ;
2823 gaprnt (2, "Scanning Descriptor File: ");
2824 gaprnt (2, pathname);
2825 gaprnt (2, "\n");
2826 #ifdef TESTING3
2827 if (optargs && (*optargs != '\0')) {
2828 gaprnt (2, "Optional args to XDFopen command are: ") ;
2829 gaprnt (2, optargs) ;
2830 gaprnt (2, "\n") ;
2831 }
2832 #endif
2833 pfi = getpfi();
2834 if (pfi == NULL) {
2835 gaprnt (0,"Memory Allocation Error: On File Open\n");
2836 return (1);
2837 }
2838 pfi->is_a_SDF = 1 ;
2839 pfi->ncflg = 3; /* set a flag for detecting xdf files that is not #ifdef'd */
2840 if (!init_io_std(&(pfi->sdf_ptr))) {
2841 gaprnt(0, "Memory Allocation Error: On XDF File Open\n") ;
2842 return(1) ;
2843 }
2844 pfi->tmplat = 0 ; /* assume no templating, unless get indication on OPTIONS */
2845 rc = gadxdf(pathname, pfi, &template, &ntsteps, &sdfpath, &parms) ;
2846
2847 #ifdef TESTING
2848 fprintf(stderr, "Return code from gadxdf=%d.\n", rc) ;
2849 fprintf(stderr, "Calling xdumpfile after gadxdf.\n") ;
2850 xdumpfile(pfi) ;
2851 if (pfi->pvar1) {
2852 fprintf(stderr, "Calling xdumpvar on first var after gadxdf.\n") ;
2853 xdumpvar(pfi->pvar1) ;
2854 }
2855 #endif
2856 if (rc) {
2857 frepfi (pfi, 0) ;
2858 return (rc) ;
2859 }
2860 if (pcm->pfi1==NULL) {
2861 pcm->pfi1 = pfi;
2862 } else {
2863 pfio = pcm->pfi1;
2864 while (pfio->pforw!=NULL) pfio = pfio->pforw;
2865 pfio->pforw = pfi;
2866 }
2867 pfi->pforw = NULL;
2868 pcm->fnum++;
2869 if (pcm->fnum==1) {pcm->pfid = pcm->pfi1; pcm->dfnum = 1;}
2870 rc = gadsdf(sdfpath, pfi, template, ntsteps, parms) ;
2871 if (rc) {
2872 sprintf(pout, "SDF file %s was not successfully opened & parsed.\n", sdfpath) ;
2873 gaprnt(0, pout) ;
2874 frepfi (pfi, 0) ;
2875 if (pcm->fnum <= 1) {
2876 pcm->pfid = pcm->pfi1 = (struct gafile *) 0 ;
2877 pcm->fnum = 0 ;
2878 pcm->dfnum = 0 ;
2879 } else {
2880 pcm->fnum-- ;
2881 for (idummy = 1, pfi = pcm->pfi1 ; (idummy < pcm->fnum) && pfi ; ++idummy) {
2882 pfi = pfi->pforw ;
2883 }
2884 if (pfi) {
2885 pfi->pforw = (struct gafile *) 0 ;
2886 }
2887 }
2888 pfi = (struct gafile *) 0 ;
2889 return(rc) ;
2890 } else {
2891 sprintf(pout, "SDF file %s is open as file %i\n", pfi->name, pcm->fnum);
2892 gaprnt(2, pout);
2893 }
2894 /* If first file open, set up some default dimension
2895 ranges for the user */
2896
2897 if (pcm->fnum==1) {
2898 if (pfi->type==2 || pfi->wrap) gacmd ("set lon 0 360",pcm,0);
2899 else {
2900 sprintf (pout,"set x 1 %i",pfi->dnum[0]);
2901 gacmd (pout,pcm,0);
2902 }
2903 if (pfi->type==2) {
2904 gacmd ("set lat -90 90",pcm,0);
2905 gacmd ("set lev 500",pcm,0);
2906 } else {
2907 sprintf (pout,"set y 1 %i",pfi->dnum[1]);
2908 gacmd (pout,pcm,0);
2909 gacmd ("set z 1",pcm,0);
2910 }
2911 gacmd ("set t 1",pcm,0);
2912 }
2913
2914 if (pfi->ppflag) {
2915 sprintf (pout,"Notice: Implied interpolation for file %s\n",pathname);
2916 gaprnt (1,pout);
2917 gaprnt (1," Interpolation will be performed on any data ");
2918 gaprnt (1,"displayed from this file\n");
2919 }
2920
2921 return (0);
2922 }
2923
2924 /* This function takes a grid request struct, fills it from an SDF file, */
2925 /* and returns the fgrid array with floating point values, unpacking the */
2926 /* SDF values if need be */
2927
2928 #define SDFTINDEX 0
2929 #define GXINDEX 0
2930 #define GYINDEX 1
2931 #define GZINDEX 2
2932 #define GTINDEX 3
2933 #define AONAME "add_offset"
2934 #define SFNAME "scale_factor"
2935 #define AONAME2 "intercept"
2936 #define SFNAME2 "slope"
2937 #define MVNAME "missing_value"
2938
2939 int
gagsdf(struct gagrid * gridptr,float fgrid[])2940 gagsdf(struct gagrid *gridptr, float fgrid[]) {
2941 /* This is a wrapper that determines the part of the request that is within */
2942 /* the bounds of the file, and makes a fake request to gafgsdf for that part */
2943 /* It first fills in the entire request with UNDEF values. After gafgsdf */
2944 /* returns, the "good" area is copied into the real grid, and the real grid */
2945 /* request structure is updated for sdf_ptr (might've opened a new file if */
2946 /* templating) and for pfile->fnumc (same reason) */
2947 int gafgsdf(struct gagrid *fakegrdptr, float fakegrd[], int Xinvolved) ;
2948 int rgxstart, rgxstop, rgxsize, rgystart, rgystop, rgysize ;
2949 int rgzstart, rgzstop, rgzsize, rgtstart, rgtstop, rgtsize ;
2950 int gxinx, gyinx, gzinx, gtinx, ginx, indx, retcode, yxsize, iindx, jindx ;
2951 int lxinx, lyinx, lzinx, ltinx, gxsize, gysize, gzsize, gtsize, zyxsize ;
2952 int rzyxsize, ryxsize, gridsize, istart, istop, jstart, jstop, ioff, joff ;
2953 int *iinxptr, *jinxptr, Xinvolved, abbrvlen ;
2954 int dim0size, dim1size, dim2size, dim3size, dim2_3size, dim1_2_3size ;
2955 int *dim0inx, *dim1inx, *dim2inx, *dim3inx, iisx, jisx ;
2956 VAR_INFO *data_var ;
2957 float *fakefgrid, missval ;
2958 struct gagrid *fakegridptr ;
2959 #define MIn(a, b) ((a < b)?a:b)
2960 #define MAx(a, b) ((a > b)?a:b)
2961
2962 for (data_var = gridptr->pfile->sdf_ptr->var ; data_var ;
2963 data_var = data_var->next) {
2964 if (!strcasecmp(gridptr->pvar->abbrv, data_var->gradsabbr)) {
2965 break ;
2966 }
2967 }
2968 if (!data_var) {
2969 if (gridptr->pvar->longnm) {
2970 sprintf(pout, "Couldn't find variable %s (aliased to %s) in SDF file.\n",
2971 gridptr->pvar->longnm, gridptr->pvar->abbrv) ;
2972 gaprnt(0, pout) ;
2973 return(2) ;
2974 } else {
2975 sprintf(pout, "Couldn't find variable %s in SDF file.\n", gridptr->pvar->abbrv) ;
2976 gaprnt(0, pout) ;
2977 return(2) ;
2978 }
2979 }
2980 switch (data_var->vartype) {
2981 case NC_BYTE:
2982 gridptr->undef = (float) data_var->missing_value.bval ;
2983 break ;
2984 case NC_SHORT:
2985 gridptr->undef = (float) data_var->missing_value.sval ;
2986 break ;
2987 case NC_LONG:
2988 gridptr->undef = (float) data_var->missing_value.lval ;
2989 break ;
2990 case NC_FLOAT:
2991 gridptr->undef = (float) data_var->missing_value.fval ;
2992 break ;
2993 case NC_DOUBLE:
2994 gridptr->undef = (float) data_var->missing_value.dval ;
2995 break ;
2996 } /* end switch */
2997 if (gridptr->pfile->sdf_ptr->hasDDFundef) {
2998 gridptr->undef = gridptr->pfile->undef ;
2999 } /* UNDEF stmt in DDF overrides NetCDF metadata, if any */
3000 Xinvolved = (gridptr->idim == GXINDEX) || (gridptr->jdim == GXINDEX) ;
3001 if ((gridptr->idim) >= 0) {
3002 if ((gridptr->jdim) >= 0) {
3003 gridsize = gridptr->isiz * gridptr->jsiz ;
3004 } else {
3005 gridsize = gridptr->isiz ;
3006 }
3007 } else {
3008 if ((gridptr->jdim) >= 0) {
3009 gridsize = gridptr->jsiz ;
3010 } else /* both invalid => 1 grid locs */ {
3011 gridsize = 1 ;
3012 }
3013 }
3014 gxsize = gridptr->dimmax[GXINDEX] - gridptr->dimmin[GXINDEX] + 1 ;
3015 gysize = gridptr->dimmax[GYINDEX] - gridptr->dimmin[GYINDEX] + 1 ;
3016 gzsize = gridptr->dimmax[GZINDEX] - gridptr->dimmin[GZINDEX] + 1 ;
3017 gtsize = gridptr->dimmax[GTINDEX] - gridptr->dimmin[GTINDEX] + 1 ;
3018 yxsize = gysize * gxsize ;
3019 zyxsize = gzsize * yxsize ;
3020 for (indx = 0 ; indx < gridsize ; ++indx) {
3021 fgrid[indx] = gridptr->undef ;
3022 } /* assume all missing initially; this covers out-of-bounds requests */
3023 if ((gridptr->idim != GXINDEX) && (gridptr->jdim != GXINDEX)) {
3024 rgxstart = rgxstop = MAx(gridptr->dimmin[GXINDEX], 1) ;
3025 /* X-dim. not involved this request */
3026 } else {
3027 if (gridptr->pfile->wrap && Xinvolved) {
3028 /* if longitudes wrap, there are no "outside" x locations */
3029 rgxstart = gridptr->dimmin[GXINDEX] ;
3030 rgxstop = gridptr->dimmax[GXINDEX] ;
3031 } else {
3032 if ((gridptr->dimmin[GXINDEX] > gridptr->pfile->dnum[GXINDEX]) ||
3033 (gridptr->dimmax[GXINDEX] < 1)) {
3034 return(0) /* all points are out of bounds */ ;
3035 }
3036 rgxstart = MAx(gridptr->dimmin[GXINDEX], 1) ;
3037 rgxstart = MIn(rgxstart, gridptr->pfile->dnum[GXINDEX]) ;
3038 rgxstop = MIn(gridptr->dimmax[GXINDEX], gridptr->pfile->dnum[GXINDEX]) ;
3039 rgxstop = MAx(rgxstop, 1) ;
3040 }
3041 }
3042 rgxsize = rgxstop - rgxstart + 1 ;
3043 if ((gridptr->idim != GYINDEX) && (gridptr->jdim != GYINDEX)) {
3044 rgystart = rgystop = MAx(gridptr->dimmin[GYINDEX], 1) ;
3045 /* Y-dim. not involved this request */
3046 } else {
3047 if ((gridptr->dimmin[GYINDEX] > gridptr->pfile->dnum[GYINDEX]) ||
3048 (gridptr->dimmax[GYINDEX] < 1)) {
3049 return(0) /* all points are out of bounds */ ;
3050 }
3051 rgystart = MAx(gridptr->dimmin[GYINDEX], 1) ;
3052 rgystart = MIn(rgystart, gridptr->pfile->dnum[GYINDEX]) ;
3053 rgystop = MIn(gridptr->dimmax[GYINDEX], gridptr->pfile->dnum[GYINDEX]) ;
3054 rgystop = MAx(rgystop, 1) ;
3055 }
3056 rgysize = rgystop - rgystart + 1 ;
3057 if ((gridptr->idim != GZINDEX) && (gridptr->jdim != GZINDEX)) {
3058 if (gridptr->pvar->levels > 0) {
3059 rgzstart = rgzstop = MAx(gridptr->dimmin[GZINDEX], 1) ;
3060 } else {
3061 rgzstart = rgzstop = 1 ;
3062 }
3063 /* Z-dim. not involved this request */
3064 } else {
3065 if ((gridptr->dimmin[GZINDEX] > gridptr->pfile->dnum[GZINDEX]) ||
3066 (gridptr->dimmax[GZINDEX] < 1)) {
3067 return(0) /* all points are out of bounds */ ;
3068 }
3069 rgzstart = MAx(gridptr->dimmin[GZINDEX], 1) ;
3070 rgzstart = MIn(rgzstart, gridptr->pfile->dnum[GZINDEX]) ;
3071 rgzstop = MIn(gridptr->dimmax[GZINDEX], gridptr->pfile->dnum[GZINDEX]) ;
3072 rgzstop = MAx(rgzstop, 1) ;
3073 }
3074 rgzsize = rgzstop - rgzstart + 1 ;
3075 if ((gridptr->idim != GTINDEX) && (gridptr->jdim != GTINDEX)) {
3076 rgtstart = rgtstop = MAx(gridptr->dimmin[GTINDEX], 1) ;
3077 /* T-dim. not involved in this request */
3078 } else {
3079 if ((gridptr->dimmin[GTINDEX] > gridptr->pfile->dnum[GTINDEX]) ||
3080 (gridptr->dimmax[GTINDEX] < 1)) {
3081 return(0) /* all points are out of bounds */ ;
3082 }
3083 rgtstart = MAx(gridptr->dimmin[GTINDEX], 1) ;
3084 rgtstart = MIn(rgtstart, gridptr->pfile->dnum[GTINDEX]) ;
3085 rgtstop = MIn(gridptr->dimmax[GTINDEX], gridptr->pfile->dnum[GTINDEX]) ;
3086 rgtstop = MAx(rgtstop, 1) ;
3087 }
3088 rgtsize = rgtstop - rgtstart + 1 ;
3089 switch (gridptr->idim) {
3090 case GXINDEX:
3091 istart = rgxstart ;
3092 istop = rgxstop ;
3093 /* actually, istop should be restrict. by SDF max., but it works, so no change */
3094 iinxptr = &gxinx ;
3095 break ;
3096 case GYINDEX:
3097 istart = rgystart ;
3098 istop = rgystop ;
3099 iinxptr = &gyinx ;
3100 break ;
3101 case GZINDEX:
3102 istart = rgzstart ;
3103 istop = rgzstop ;
3104 iinxptr = &gzinx ;
3105 break ;
3106 case GTINDEX:
3107 istart = rgtstart ;
3108 istop = rgtstop ;
3109 iinxptr = >inx ;
3110 break ;
3111 default:
3112 istart = istop = 0 ; iinxptr = NULL ;
3113 break ;
3114 }
3115 if (gridptr->dimmin[gridptr->idim] < 1) {
3116 ioff = 1 - gridptr->dimmin[gridptr->idim] ;
3117 } else {
3118 ioff = 0 ;
3119 }
3120 if (gridptr->pfile->wrap && (gridptr->idim == GXINDEX)) /* wrap case */ {
3121 ioff = 0 ; /* would otherwise be miss-set, in normal case to 1 */
3122 }
3123 switch (gridptr->jdim) {
3124 case GXINDEX:
3125 jstart = rgxstart ;
3126 jstop = rgxstop ;
3127 jinxptr = &gxinx ;
3128 break ;
3129 case GYINDEX:
3130 jstart = rgystart ;
3131 jstop = rgystop ;
3132 jinxptr = &gyinx ;
3133 break ;
3134 case GZINDEX:
3135 jstart = rgzstart ;
3136 jstop = rgzstop ;
3137 jinxptr = &gzinx ;
3138 break ;
3139 case GTINDEX:
3140 jstart = rgtstart ;
3141 jstop = rgtstop ;
3142 jinxptr = >inx ;
3143 break ;
3144 default:
3145 jstart = jstop = 0 ; jinxptr = NULL ;
3146 break ;
3147 }
3148 if (gridptr->dimmin[gridptr->jdim] < 1) {
3149 joff = 1 - gridptr->dimmin[gridptr->jdim] ;
3150 } else {
3151 joff = 0 ;
3152 }
3153 if (gridptr->pfile->wrap && (gridptr->jdim == GXINDEX)) /* wrap case */ {
3154 joff = 0 ; /* would otherwise be miss-set, in normal case to 1 */
3155 }
3156 ryxsize = rgysize * rgxsize ;
3157 rzyxsize = rgzsize * ryxsize ;
3158 if (gridsize > 0) {
3159 fakefgrid = (float *) malloc(sizeof(float)*gridsize) ;
3160 if (!fakefgrid) {
3161 gaprnt(0, "Memory allocation error handling SDF grid request(1).\n") ;
3162 return(13) ;
3163 }
3164 fakegridptr = (struct gagrid *) malloc(sizeof(struct gagrid)) ;
3165 if (!fakegridptr) {
3166 gaprnt(0, "Memory allocation error handling SDF grid request(2).\n") ;
3167 if (fakefgrid) {
3168 free(fakefgrid) ;
3169 }
3170 return(14) ;
3171 }
3172 } else /* didn't we return on gridsize == 0 already? */ {
3173 return(0) ;
3174 }
3175 /* no doubt, some C compiler on some platform will blow structure assignment, */
3176 /* so this is done the hard way, rather than "*fakegridptr = *gridptr ;" */
3177 fakegridptr->grid = gridptr->grid ;
3178 fakegridptr->pfile = gridptr->pfile ;
3179 fakegridptr->undef = gridptr->undef ;
3180 fakegridptr->rmin = gridptr->rmin ;
3181 fakegridptr->rmax = gridptr->rmax ;
3182 fakegridptr->isiz = gridptr->isiz ;
3183 if ((gridptr->idim) < 0) fakegridptr->isiz = 1 ;
3184 fakegridptr->jsiz = gridptr->jsiz ;
3185 if ((gridptr->jdim) < 0) fakegridptr->jsiz = 1 ;
3186 fakegridptr->idim = gridptr->idim ;
3187 fakegridptr->jdim = gridptr->jdim ;
3188 /* only two (at most!) of these should really vary now */
3189 fakegridptr->dimmin[GXINDEX] = rgxstart ;
3190 fakegridptr->dimmax[GXINDEX] = rgxstop ;
3191 fakegridptr->dimmin[GYINDEX] = rgystart ;
3192 fakegridptr->dimmax[GYINDEX] = rgystop ;
3193 fakegridptr->dimmin[GZINDEX] = rgzstart ;
3194 fakegridptr->dimmax[GZINDEX] = rgzstop ;
3195 fakegridptr->dimmin[GTINDEX] = rgtstart ;
3196 fakegridptr->dimmax[GTINDEX] = rgtstop ;
3197 fakegridptr->pvar = gridptr->pvar ;
3198 fakegridptr->exprsn = gridptr->exprsn ;
3199 fakegridptr->alocf = gridptr->alocf ;
3200 fakegridptr->igrab = gridptr->igrab ;
3201 fakegridptr->jgrab = gridptr->jgrab ;
3202 fakegridptr->iabgr = gridptr->iabgr ;
3203 fakegridptr->jabgr = gridptr->jabgr ;
3204 fakegridptr->ivals = gridptr->ivals ;
3205 fakegridptr->jvals = gridptr->jvals ;
3206 fakegridptr->ilinr = gridptr->ilinr ;
3207 fakegridptr->jlinr = gridptr->jlinr ;
3208 retcode = gafgsdf(fakegridptr, fakefgrid, Xinvolved) ;
3209 if (retcode == 0) {
3210 for (gxinx = rgxstart, lxinx = 0 ; gxinx <= rgxstop ;
3211 ++gxinx, ++lxinx) {
3212 for (gyinx = rgystart, lyinx = 0 ; gyinx <= rgystop ;
3213 ++gyinx, ++lyinx) {
3214 for (gzinx = rgzstart, lzinx = 0 ; gzinx <= rgzstop ;
3215 ++gzinx, ++lzinx) {
3216 for (gtinx = rgtstart, ltinx = 0 ; gtinx <= rgtstop ;
3217 ++gtinx, ++ltinx) {
3218 if (!jinxptr) {
3219 if (!iinxptr) {
3220 ginx = 0 ;
3221 } else {
3222 ginx = (*iinxptr) - istart + ioff ;
3223 }
3224 } else /* jinxptr is good */ {
3225 if (!iinxptr) {
3226 ginx = (*jinxptr) - jstart + joff ;
3227 } else /* both are good */ {
3228 }
3229 }
3230 ginx = (gtinx - gridptr->dimmin[GTINDEX]) * zyxsize +
3231 (gzinx - gridptr->dimmin[GZINDEX]) * yxsize +
3232 (gyinx - gridptr->dimmin[GYINDEX]) * gxsize +
3233 (gxinx - gridptr->dimmin[GXINDEX]) ;
3234 /* doesn't this calculation make the t-z-y-x assumption? */
3235 indx = ltinx * rzyxsize + lzinx * ryxsize +
3236 lyinx * rgxsize + lxinx ;
3237 fgrid[ginx] = fakefgrid[indx] ;
3238 } /* time */
3239 } /* z */
3240 } /* y */
3241 } /* x */
3242 } /* if there's any point in copying over the results */
3243 if (fakefgrid) {
3244 free(fakefgrid) ;
3245 }
3246 if (fakegridptr) {
3247 free(fakegridptr) ;
3248 }
3249 return(retcode) ;
3250 }
3251
3252 /* This is the fake grid request handler for SDF files. It is wrapped by */
3253 /* gagsdf, which passes a fake grid request to this that is within the */
3254 /* file's bounds */
3255 int
gafgsdf(struct gagrid * gridptr,float fgrid[],int Xinvolved)3256 gafgsdf(struct gagrid *gridptr, float fgrid[], int Xinvolved) {
3257 long *start, *count, gridcnt ;
3258 int is4d, ispacked, abbrvlen, /* whichdim[4], */ gridinx ;
3259 int xstart, ystart, zstart, tstart, xstop, ystop, zstop, tstop ;
3260 int xincr, yincr, zincr, tincr, xinx, yinx, zinx, tinx ;
3261 int gxinx, gyinx, gzinx, gtinx, ztsize, yztsize ;
3262 int gxstart, gxstop, gxincr, gtstart, gtstop, gtincr, grinx, inx ;
3263 int gystart, gystop, gyincr, gzstart, gzstop, gzincr ;
3264 int xindex, yindex, zindex, tindex, yxsize, zyxsize, mininx, maxinx ;
3265 int timeindex, tcnt, zcnt, ycnt, xcnt ;
3266 int xdimindx, ydimindx, zdimindx, tdimindx, tcount, starttime, newtime ;
3267 int willwrap, littlex, bigx, sdfxcount, usingFV=0, wasmiss ;
3268 int fillwithmiss, alltimesinthisfile, alltimesin1file ;
3269 float *fdata, fmiss, myao, mysf, fFillValue ;
3270 short *sdata, smiss, sFillValue ;
3271 double *ddata, dFillValue, dmiss ;
3272
3273 struct attrib_list *aoattr, *sfattr, *missvalattr, *FVattr ;
3274 VAR_INFO *data_var ;
3275 union att_val ao, sf, junk ;
3276 nc_type unpktype, pktype ;
3277 struct attrib_list *find_att(struct attrib_list *first_att, char *attname) ;
3278 int are_data_packed(VAR_INFO *in_var, union att_val *scale,
3279 union att_val *offset, nc_type *unpktype,
3280 nc_type *pktype) ;
3281 int sdfwrap(VAR_INFO **dvar_ptr, struct attrib_list *missvalattr, int ispacked,
3282 int is4d, float ao, float sf, int tdimindx, int tindex,
3283 int zdimindx, int zindex, int ydimindx, int yindex, int xdimindx,
3284 int xindex, struct gagrid *gridptr, float *fgrid,
3285 nc_type pktype, nc_type unpktype, int sdfxcount, int usingFV) ;
3286 int gaopsdf(int timeindex, struct gagrid *gridptr, VAR_INFO **data_var_ptr) ;
3287 void init_start_count(long **start, long **count, int ndims) ;
3288 int fequal(float op1, float op2, float delta) ;
3289 int dequal(double op1, double op2, double delta) ;
3290
3291 /* Search for data_var based on max. 15 char match on gridptr->pvar->abbrv */
3292 /* Fill start and count based on gridptr->dimmin and gridptr->dimmax and is4d */
3293 /* Call read_variable_data(gridptr->pfile->sdf_ptr->cdfid, data_var, start, */
3294 /* count) */
3295 /* Check data_var->vartype: NC_SHORT => ispacked=1, NC_FLOAT => ispacked=0 */
3296 /* If ispacked, unpack into grid; else copy into grid. */
3297 /* Addenda to above: if (ispacked || (pktype != FLOAT)), treat as packed */
3298
3299 /* abbrvlen = (int) strlen(gridptr->pvar->abbrv) ; */
3300 /* netCDF name may be longer */
3301 fillwithmiss = 0 ;
3302 for (data_var = gridptr->pfile->sdf_ptr->var ; data_var ;
3303 data_var = data_var->next) {
3304 if (!strcasecmp(gridptr->pvar->abbrv, data_var->gradsabbr)) {
3305 break ;
3306 }
3307 }
3308 if (!data_var) {
3309 sprintf(pout, "Couldn't find variable %s in SDF file.\n", gridptr->pvar->abbrv) ;
3310 gaprnt(0, pout) ;
3311 return(2) ;
3312 }
3313 init_start_count(&start, &count, data_var->nvardims) ;
3314 /* Here's the rub: if this is an XDFopen-ed file, and that DDF had an UNDEF */
3315 /* statement, we don't want to look for it in an attribute value. However, */
3316 /* I don't know how to tell this origin info at this point; it was in the */
3317 /* parms structure during data discovery phase, (isxdf and needundef */
3318 /* components), but I don't see any way to access that here. Hoop, 2002/03/13 */
3319 missvalattr = (struct attrib_list *) 0 ;
3320 if (gridptr->pfile->sdf_ptr->hasDDFundef) {
3321 fmiss = gridptr->pfile->undef ;
3322 } else {
3323 missvalattr = find_att(data_var->first_vattr, MVNAME) ; /* ok if not found? */
3324 if (!missvalattr) {
3325 missvalattr = find_att(data_var->first_vattr, "_FillValue") ;
3326 usingFV = 1 ; /* don't need to make two checks if using this one */
3327 }
3328 }
3329 ispacked = are_data_packed(data_var, &sf, &ao, &unpktype, &pktype) ;
3330 if (pktype == NC_FLOAT) {
3331 if (!(gridptr->pfile->sdf_ptr->hasDDFundef)) {
3332 if (missvalattr) {
3333 fdata = (float *) missvalattr->data ;
3334 fmiss = *fdata ;
3335 } else {
3336 fmiss = FILL_FLOAT ;
3337 }
3338 }
3339 } else if (pktype == NC_DOUBLE) {
3340 if (!(gridptr->pfile->sdf_ptr->hasDDFundef)) {
3341 if (missvalattr) {
3342 ddata = (double *) missvalattr->data ;
3343 dmiss = *ddata ;
3344 } else {
3345 dmiss = FILL_DOUBLE ;
3346 }
3347 }
3348 }
3349 if ((pktype != NC_FLOAT) &&(pktype != NC_DOUBLE)) ispacked = 1 ;
3350 if (ispacked && (pktype != NC_SHORT)) /* bail for now */ {
3351 gaprnt(0, "SDF data variable has unsupported packed data type.\n") ;
3352 return(6) ;
3353 }
3354 if (ispacked && (pktype == unpktype)) {
3355 /* trouble is, this call doesn't happen in are_data_packed if pktype=unpktype */
3356 /* Another trouble is that we don't want to do this if sf or ao are present */
3357 if (!(sfattr = find_att(data_var->first_vattr, SFNAME))) {
3358 if (!(sfattr = find_att(data_var->first_vattr, SFNAME2))) {
3359 set_default_scale_offset(unpktype, &sf, &junk) ;
3360 }
3361 }
3362 if (!(aoattr = find_att(data_var->first_vattr, AONAME))) {
3363 if (!(aoattr = find_att(data_var->first_vattr, AONAME2))) {
3364 set_default_scale_offset(unpktype, &junk, &ao) ;
3365 }
3366 }
3367 }
3368 if (ispacked) /* figure myao, mysf */ {
3369 switch (unpktype) {
3370 case NC_SHORT:
3371 myao = (float) ao.sval ;
3372 mysf = (float) sf.sval ;
3373 break ;
3374 case NC_LONG:
3375 myao = (float) ao.lval ;
3376 mysf = (float) sf.lval ;
3377 break ;
3378 case NC_FLOAT:
3379 myao = ao.fval ;
3380 mysf = sf.fval ;
3381 break ;
3382 case NC_DOUBLE:
3383 myao = (float) ao.dval ;
3384 mysf = (float) sf.dval ;
3385 break ;
3386 default:
3387 gaprnt(0, "SDF data variable has unsupported unpacked data type.\n") ;
3388 break ;
3389 } /* end switch */
3390 }
3391 tdimindx = data_var->dimidmap[0] ;
3392 tindex = data_var->dimmap[0] ;
3393 zdimindx = data_var->dimidmap[1] ;
3394 zindex = data_var->dimmap[1] ;
3395 ydimindx = data_var->dimidmap[2] ;
3396 yindex = data_var->dimmap[2] ;
3397 xdimindx = data_var->dimidmap[3] ;
3398 xindex = data_var->dimmap[3] ;
3399 is4d = (zindex != -1) ? 1 : 0 ;
3400 if (Xinvolved) {
3401 willwrap = 1 ;
3402 littlex = gridptr->dimmin[GXINDEX] ;
3403 bigx = gridptr->dimmax[GXINDEX] ;
3404 sdfxcount = gridptr->pfile->sdf_ptr->dimsiz[xdimindx] ;
3405 if ((littlex <= sdfxcount) && (bigx <= sdfxcount)) {
3406 if ((littlex > 0) && (bigx > 0)) {
3407 willwrap = 0 ;
3408 }
3409 }
3410 } else {
3411 willwrap = 0 ;
3412 }
3413 if ((gridptr->pfile->wrap && Xinvolved) && (willwrap)) {
3414 return(sdfwrap(&data_var, missvalattr, ispacked, is4d, myao, mysf,
3415 tdimindx, tindex, zdimindx, zindex, ydimindx, yindex,
3416 xdimindx, xindex, gridptr, fgrid, pktype, unpktype,
3417 sdfxcount, usingFV)) ;
3418 } else {
3419 if (!ispacked) /* just copy */ {
3420 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3421 /* need to check for _FillValue as well as missing_value */
3422 FVattr = find_att(data_var->first_vattr, "_FillValue") ;
3423 if (pktype == NC_FLOAT) {
3424 if (FVattr) {
3425 fdata = (float *) FVattr->data ;
3426 fFillValue = *fdata ;
3427 } else {
3428 fFillValue = FILL_FLOAT ;
3429 }
3430 } else /* must be double, right? */ {
3431 if (FVattr) {
3432 ddata = (double *) FVattr-> data ;
3433 dFillValue = *ddata ;
3434 } else {
3435 dFillValue = FILL_DOUBLE ;
3436 }
3437 } /* is double, I guess */
3438 } /* !usingFV */
3439 } else /* must unpack */ {
3440 if ((!missvalattr) && (pktype != unpktype)) {
3441 smiss = (int) (((gridptr->undef - myao) / mysf) + 0.5) ;
3442 } else {
3443 if ((pktype == unpktype) && (!missvalattr)) {
3444 smiss = FILL_SHORT ;
3445 } else {
3446 sdata = (short *) missvalattr->data ;
3447 smiss = *sdata ;
3448 }
3449 }
3450 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3451 FVattr = find_att(data_var->first_vattr, "_FillValue") ;
3452 if (FVattr) {
3453 sdata = (short *) FVattr->data ;
3454 sFillValue = *sdata ;
3455 } else {
3456 sFillValue = FILL_SHORT ;
3457 }
3458 } /* !using FV */
3459 } /* is it packed? */
3460 start[xindex] = gridptr->dimmin[GXINDEX] - 1 ; /* dimmin be 1-based */
3461 if (gridptr->pfile->yrflg) {
3462 maxinx = gridptr->pfile->sdf_ptr->dimsiz[ydimindx]
3463 - gridptr->dimmin[GYINDEX] ;
3464 mininx = gridptr->pfile->sdf_ptr->dimsiz[ydimindx]
3465 - gridptr->dimmax[GYINDEX] ;
3466 start[yindex] = mininx ;
3467 } else {
3468 start[yindex] = gridptr->dimmin[GYINDEX] - 1 ;
3469 }
3470 count[yindex] = gridptr->dimmax[GYINDEX] - gridptr->dimmin[GYINDEX] + 1 ;
3471 if (count[yindex] > gridptr->pfile->sdf_ptr->dimsiz[ydimindx]) {
3472 count[yindex] = gridptr->pfile->sdf_ptr->dimsiz[ydimindx] ;
3473 }
3474 count[xindex] = gridptr->dimmax[GXINDEX] - gridptr->dimmin[GXINDEX] + 1 ;
3475 if (count[xindex] > gridptr->pfile->sdf_ptr->dimsiz[xdimindx]) {
3476 count[xindex] = gridptr->pfile->sdf_ptr->dimsiz[xdimindx] ;
3477 }
3478 gxstart = xstart = 0 ;
3479 gxstop = xstop = count[xindex] ;
3480 gxincr = xincr = 1 ;
3481 gystart = 0 ;
3482 gystop = count[yindex] ;
3483 gyincr = 1 ;
3484 if (gridptr->pfile->yrflg) {
3485 ystart = gystop - 1 ;
3486 ystop = gystart - 1 ;
3487 yincr = -1 ;
3488 } else {
3489 ystart = gystart ;
3490 ystop = gystop ;
3491 yincr = 1 ;
3492 }
3493 gridcnt = count[xindex] * count[yindex] ;
3494 if (is4d) {
3495 if (gridptr->pfile->zrflg) {
3496 /* Should determine z dim index empirically */
3497 maxinx = gridptr->pfile->sdf_ptr->dimsiz[zdimindx]
3498 - gridptr->dimmin[GZINDEX] ;
3499 mininx = gridptr->pfile->sdf_ptr->dimsiz[zdimindx]
3500 - gridptr->dimmax[GZINDEX] ;
3501 start[zindex] = mininx ;
3502 } else {
3503 start[zindex] = gridptr->dimmin[GZINDEX] - 1 ;
3504 }
3505 count[zindex] = gridptr->dimmax[GZINDEX] - gridptr->dimmin[GZINDEX]
3506 + 1 ;
3507 if (count[zindex] > gridptr->pfile->sdf_ptr->dimsiz[zdimindx]) {
3508 count[zindex] = gridptr->pfile->sdf_ptr->dimsiz[zdimindx] ;
3509 }
3510 gridcnt *= count[zindex] ;
3511 gzstart = 0 ;
3512 gzstop = count[zindex] ;
3513 gzincr = 1 ;
3514 if (gridptr->pfile->zrflg) {
3515 zstart = gzstop - 1 ;
3516 zstop = gzstart - 1 ;
3517 zincr = -1 ;
3518 } else {
3519 zstart = gzstart ;
3520 zstop = gzstop ;
3521 zincr = 1 ;
3522 }
3523 } else {
3524 gzstart = zstart = 0 ;
3525 gzstop = zstop = 1 ;
3526 gzincr = zincr = 1 ;
3527 }
3528 starttime = gridptr->dimmin[GTINDEX] ;
3529 tcount = gridptr->dimmax[GTINDEX] - starttime + 1 ;
3530 if (tindex != -1) {
3531 count[tindex] = 1 ; /* read each time step separately for templating */
3532 }
3533 gtstart = tstart = 0 ;
3534 gtstop = tstop = tcount ;
3535 gtincr = tincr = 1 ;
3536 yxsize = count[yindex] * count[xindex] ;
3537 zyxsize = (gzstop - gzstart) * yxsize ;
3538 gridcnt *= tcount ; /*???*/
3539 #ifdef TESTING
3540 fprintf(stderr, "Xloop runs from %d to %d step %d (%d to %d step %d).\n",
3541 gxstart, gxstop, gxincr, xstart, xstop, xincr) ;
3542 fprintf(stderr, "Yloop runs from %d to %d step %d (%d to %d step %d).\n",
3543 gystart, gystop, gyincr, ystart, ystop, yincr) ;
3544 fprintf(stderr, "Zloop runs from %d to %d step %d (%d to %d step %d).\n",
3545 gzstart, gzstop, gzincr, zstart, zstop, zincr) ;
3546 fprintf(stderr, "Tloop runs from %d to %d step %d (%d to %d step %d).\n",
3547 gtstart, gtstop, gtincr, tstart, tstop, tincr) ;
3548 #endif
3549 /* experiment to increase performance when not templating Hoop 2K/11/22 */
3550 alltimesinthisfile = 0 ;
3551 alltimesin1file = 0 ;
3552 if (gridptr->pfile->tmplat) {
3553 if ((*(gridptr->pfile->fnums + starttime - 1)) == gridptr->pfile->fnumc) {
3554 /* this time slice starts in the current file, so let's see if */
3555 /* it ends in the same file */
3556 if ((*(gridptr->pfile->fnums + starttime - 1 + tstop - tstart - 1)) ==
3557 gridptr->pfile->fnumc) /* ends in same file */ {
3558 alltimesinthisfile = 1 ;
3559 alltimesin1file = 1 ;
3560 }
3561 } else /* first time not in current file */ {
3562 if ((*(gridptr->pfile->fnums + starttime - 1)) ==
3563 (*(gridptr->pfile->fnums + starttime - 1 + tstop - tstart - 1))) {
3564 alltimesin1file = 1 ;
3565 }
3566 } /* is first time in current file */
3567 } /* if templating */
3568 if ((!(gridptr->pfile->tmplat)) || alltimesinthisfile) {
3569 if (tindex != -1) {
3570 start[tindex] = starttime - 1 ;
3571 count[tindex] = tstop - tstart ;
3572 } /* if we have a valid tindex => we have a time dim. */
3573 if (gridptr->pfile->tmplat) {
3574 if ((gridptr->pfile->fnums[0]) != (gridptr->pfile->fnumc)) {
3575 /* if we are not on the first file, we still need to call gaopsdf to */
3576 /* adulterate the starttime to be one relative to the current file */
3577 /* -Hoop, 2001/04/14 */
3578 newtime = gaopsdf(starttime, gridptr, &data_var) ;
3579 if (newtime == -99999) /* bad surprise */ {
3580 gaprnt(0, "Bad time index in SDF file.\n") ;
3581 return(1) ; /* bad time index */
3582 } /* bad surprise */
3583 if (newtime == -88888) /* bad file, bad surprise */ {
3584 for (grinx = 0 ; grinx < gridcnt ; ++grinx) {
3585 fgrid[grinx] = gridptr->undef ;
3586 }
3587 return(0) ;
3588 } /* bad file, bad surprise */
3589 if (tindex != -1) {
3590 start[tindex] = newtime - 1 ;
3591 }
3592 } /* if not 1st file, thus need new starttime */
3593 } /* if templating */
3594 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid,
3595 data_var, start, count)) {
3596 gaprnt(0, "Failure reading SDF data. #1\n") ;
3597 return(3) ;
3598 }
3599 } /* end if not templating, read all times at once */ else {
3600 if (alltimesin1file) {
3601 if (tindex != -1) {
3602 start[tindex] = starttime - 1 ;
3603 count[tindex] = tstop - tstart ;
3604 } /* if tindex is valid => we have a time dim. */
3605 newtime = gaopsdf(starttime, gridptr, &data_var) ;
3606 if (newtime == -99999) {
3607 gaprnt(0, "Bad time index in SDF file.\n") ;
3608 return(1) ; /* bad time index */
3609 }
3610 if (newtime == -88888) /* bad file */ {
3611 for (grinx = 0 ; grinx < gridcnt ; ++grinx) {
3612 fgrid[grinx] = gridptr->undef ;
3613 }
3614 return(0) ;
3615 } /* if bad file */
3616 if (tindex != -1) {
3617 start[tindex] = newtime - 1 ;
3618 }
3619 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid,
3620 data_var, start, count)) {
3621 gaprnt(0, "Failure reading SDF data. #2\n") ;
3622 return(3) ;
3623 }
3624 } /* still do all in one read, but after gaopsdf */
3625 } /* templating, not all in this file */
3626 for (timeindex = starttime, gtinx = gtstart, tinx = tstart ;
3627 (gtinx < gtstop) && (tinx < tstop) ;
3628 gtinx += gtincr, tinx += tincr, timeindex += tincr) {
3629 newtime = timeindex ;
3630 if (gridptr->pfile->tmplat) {
3631 if ((!alltimesinthisfile) && (!alltimesin1file)) {
3632 newtime = gaopsdf(timeindex, gridptr, &data_var) ;
3633 if (newtime == -99999) {
3634 gaprnt(0, "Bad time index in SDF file.\n") ;
3635 return(1) ; /* bad time index */
3636 }
3637 if (newtime == -88888) /* bad file */ {
3638 for (grinx = 0 ; grinx < gridcnt ; ++grinx) {
3639 fgrid[grinx] = gridptr->undef ;
3640 }
3641 return(0) ;
3642 } /* if bad file */
3643 } /* if not alltimesinthisfile */
3644 if (tindex != -1) {
3645 start[tindex] = newtime - 1 ;
3646 }
3647 } /* if templating in use */
3648 #ifdef TESTING
3649 fprintf(stderr,
3650 "read_variable_data(%d, %s, [%d, %d, %d, %d], [%d, %d, %d, %d])\n",
3651 gridptr->pfile->sdf_ptr->cdfid, data_var->varnam, start[0], start[1],
3652 start[2], start[3], count[0], count[1], count[2], count[3]) ;
3653 #endif
3654 /* experiment to increase performance when not templating Hoop 2K/11/22 */
3655 if (gridptr->pfile->tmplat && (!alltimesinthisfile) &&
3656 (!alltimesin1file)) {
3657 /* no need to read per timestep if not templating (or all times in one file); */
3658 /* all read above */
3659 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid, data_var,
3660 start, count)) {
3661 gaprnt(0, "Failure reading SDF data. #3\n") ;
3662 return(3) ;
3663 }
3664 } /* end if templating and not all in one file */
3665 if (!ispacked) /* just copy */ {
3666 if (pktype == NC_DOUBLE) {
3667 ddata = (double *) data_var->data ;
3668 } else {
3669 fdata = (float *) data_var->data ;
3670 }
3671 } else /* must unpack */ {
3672 sdata = (short *) data_var->data ;
3673 } /* if unpacking needed */
3674 for (gzinx = gzstart, zinx = zstart ;
3675 (gzinx < gzstop) && ((zincr > 0)?(zinx < zstop):(zinx > zstop)) ;
3676 gzinx += gzincr, zinx += zincr) {
3677 for (gyinx = gystart, yinx = ystart ;
3678 (gyinx < gystop) && ((yincr > 0)?(yinx < ystop):(yinx > ystop)) ;
3679 gyinx += gyincr, yinx += yincr) {
3680 for (gxinx = gxstart, xinx = xstart ;
3681 (gxinx < gxstop) && (xinx < xstop) ;
3682 gxinx += gxincr, xinx += xincr) {
3683 grinx = (gtinx * zyxsize) + (gzinx * yxsize) +
3684 (gyinx * count[xindex]) + gxinx ;
3685 xcnt = xinx - xstart ;
3686 if (gridptr->pfile->yrflg) {
3687 ycnt = yinx - ystop - 1 ;
3688 } else {
3689 ycnt = yinx - ystart ;
3690 }
3691 if (is4d) {
3692 if (gridptr->pfile->zrflg) {
3693 zcnt = zinx - zstop - 1 ;
3694 } else {
3695 zcnt = zinx - zstart ;
3696 }
3697 } else {
3698 zcnt = 1 ;
3699 }
3700 tcnt = tinx - tstart ;
3701 if (is4d) {
3702 inx = (zcnt * yxsize) + (ycnt * count[xindex]) + xcnt ;
3703 } else {
3704 inx = (ycnt * count[xindex]) + xcnt ;
3705 }
3706 /* experiment to increase performance when not templating Hoop 2K/11/22 */
3707 if ((!(gridptr->pfile->tmplat)) || alltimesinthisfile
3708 || alltimesin1file) {
3709 inx += (tinx * zyxsize) ;
3710 }
3711 if (ispacked) {
3712 /* handle missing */
3713 wasmiss = 0 ;
3714 if (sdata[inx] == smiss) {
3715 wasmiss = 1 ;
3716 } else if ((!usingFV)&&(!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3717 /* check FV too */
3718 if (sdata[inx] == sFillValue) {
3719 wasmiss = 1 ;
3720 }
3721 }
3722 if (!wasmiss) {
3723 fgrid[grinx] = (mysf * ((float) sdata[inx]))
3724 + myao ;
3725 } else {
3726 fgrid[grinx] = gridptr->undef ;
3727 }
3728 } else {
3729 /* May want to swap missing values */
3730 /* Do want to. -Hoop, 2K/07/25 */
3731 if (pktype == NC_DOUBLE) {
3732 wasmiss = 0 ;
3733 if (dequal(ddata[inx], dmiss, (double)1.0e-5)) {
3734 wasmiss = 1 ;
3735 }
3736 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3737 if (dequal(ddata[inx], dFillValue, (double)1.0e-5)) {
3738 wasmiss = 1 ;
3739 }
3740 }
3741 if (!wasmiss) {
3742 fgrid[grinx] = (float) ddata[inx] ;
3743 } else {
3744 fgrid[grinx] = gridptr->undef ;
3745 }
3746 } else {
3747 wasmiss = 0 ;
3748 if (fequal(fdata[inx], fmiss, 1.0e-5)) {
3749 wasmiss = 1 ;
3750 }
3751 if ((!usingFV)&&(!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3752 if (fequal(fdata[inx], fFillValue, 1.0e-5)) {
3753 wasmiss = 1 ;
3754 }
3755 }
3756 if (!wasmiss) {
3757 fgrid[grinx] = fdata[inx] ;
3758 } else {
3759 fgrid[grinx] = gridptr->undef ;
3760 }
3761 }
3762 } /* end if packed */
3763 } /* end for x indexing */
3764 } /* end for y indexing */
3765 } /* end for z indexing */
3766 } /* end for time indexing */
3767 return(0) ;
3768 } /* end if not wrapped data */
3769 }
3770
3771 /* below changed 1998/04/20 to buffer in all longitudes of a grid and pick */
3772 /* through the buffering to achieve wrapping effect -Hoop */
3773 int
sdfwrap(VAR_INFO ** data_var_ptr,struct attrib_list * missvalattr,int ispacked,int is4d,float ao,float sf,int tdimindx,int tindex,int zdimindx,int zindex,int ydimindx,int yindex,int xdimindx,int xindex,struct gagrid * gridptr,float * fgrid,nc_type pktype,nc_type unpktype,int sdfxcount,int usingFV)3774 sdfwrap(VAR_INFO **data_var_ptr, struct attrib_list *missvalattr, int ispacked,
3775 int is4d, float ao, float sf, int tdimindx, int tindex,
3776 int zdimindx, int zindex, int ydimindx, int yindex, int xdimindx,
3777 int xindex, struct gagrid *gridptr, float *fgrid, nc_type pktype,
3778 nc_type unpktype, int sdfxcount, int usingFV) {
3779 long start[4] = {0, 0, 0, 0}, count[4] = {1, 1, 1, 1}, gridcnt, yxsize ;
3780 long zyxsize, xsize ;
3781 float *fdata, fmiss, *ftemp, fFillValue ;
3782 short *sdata, smiss, *stemp, sFillValue ;
3783 double *ddata, dmiss, *dtemp, dFillValue ;
3784 struct attrib_list *FVattr ;
3785 VAR_INFO *data_var = *data_var_ptr ;
3786 int gridinx, xstart, ystart, zstart, tstart, xstop, ystop, zstop, tstop ;
3787 int xincr, yincr, zincr, tincr, xinx, yinx, zinx, tinx, gxinx, gyinx ;
3788 int gzinx, gtinx, ztsize, yztsize, gxstart, gxstop, gxincr, gystart ;
3789 int gystop, gyincr, gzstart, gzstop, gzincr, mininx, maxinx, xact ;
3790 int xreqmin, xreqmax, xreq, xcount, xpos, gtstart, gtstop, gtincr ;
3791 int inx, grinx, starttime, newtime, tcount, timeindex, gyxsize, gzyxsize ;
3792 int tcnt, zcnt, ycnt, xcnt, wasmiss ;
3793 int fillwithmiss ; /* under orders from gaopsdf */
3794 int alltimesinthisfile, alltimesin1file ;
3795 int gaopsdf(int timeindex, struct gagrid *gridptr, VAR_INFO **data_var_ptr) ;
3796 struct attrib_list *find_att(struct attrib_list *first_att, char *attname) ;
3797 int fequal(float op1, float op2, float delta) ;
3798 int dequal(double op1, double op2, double delta) ;
3799
3800 fillwithmiss = 0 ;
3801 xcount = gridptr->pfile->dnum[GXINDEX] ;
3802 xreqmin = gridptr->dimmin[GXINDEX] ;
3803 xreqmax = gridptr->dimmax[GXINDEX] ;
3804 if (gridptr->pfile->yrflg) {
3805 maxinx = gridptr->pfile->sdf_ptr->dimsiz[ydimindx]
3806 - gridptr->dimmin[GYINDEX] ;
3807 mininx = gridptr->pfile->sdf_ptr->dimsiz[ydimindx]
3808 - gridptr->dimmax[GYINDEX] ;
3809 start[yindex] = mininx ;
3810 } else {
3811 start[yindex] = gridptr->dimmin[GYINDEX] - 1 ;
3812 }
3813 count[yindex] = gridptr->dimmax[GYINDEX] - gridptr->dimmin[GYINDEX] + 1 ;
3814 if (count[yindex] > gridptr->pfile->sdf_ptr->dimsiz[ydimindx]) {
3815 count[yindex] = gridptr->pfile->sdf_ptr->dimsiz[ydimindx] ;
3816 }
3817 start[xindex] = 0 ;
3818 count[xindex] = sdfxcount ;
3819 xsize = gridptr->dimmax[GXINDEX] - gridptr->dimmin[GXINDEX] + 1 ;
3820 gystart = 0 ;
3821 gystop = count[yindex] ;
3822 gyincr = 1 ;
3823 if (gridptr->pfile->yrflg) {
3824 ystart = gystop - 1 ;
3825 ystop = gystart - 1 ;
3826 yincr = -1 ;
3827 } else {
3828 ystart = gystart ;
3829 ystop = gystop ;
3830 yincr = 1 ;
3831 }
3832 gridcnt = count[xindex] * count[yindex] ;
3833 if (is4d) {
3834 if (gridptr->pfile->zrflg) {
3835 /* Should determine z dim index empirically */
3836 maxinx = gridptr->pfile->sdf_ptr->dimsiz[zdimindx]
3837 - gridptr->dimmin[GZINDEX] ;
3838 mininx = gridptr->pfile->sdf_ptr->dimsiz[zdimindx]
3839 - gridptr->dimmax[GZINDEX] ;
3840 start[zindex] = mininx ;
3841 } else {
3842 start[zindex] = gridptr->dimmin[GZINDEX] - 1 ;
3843 }
3844 count[zindex] = gridptr->dimmax[GZINDEX] - gridptr->dimmin[GZINDEX] + 1 ;
3845 if (count[zindex] > gridptr->pfile->sdf_ptr->dimsiz[zdimindx]) {
3846 count[zindex] = gridptr->pfile->sdf_ptr->dimsiz[zdimindx] ;
3847 }
3848 gridcnt *= count[zindex] ;
3849 gzstart = 0 ;
3850 gzstop = count[zindex] ;
3851 gzincr = 1 ;
3852 if (gridptr->pfile->zrflg) {
3853 zstart = gzstop - 1 ;
3854 zstop = gzstart - 1 ;
3855 zincr = -1 ;
3856 } else {
3857 zstart = gzstart ;
3858 zstop = gzstop ;
3859 zincr = 1 ;
3860 }
3861 } else {
3862 gzstart = zstart = 0 ;
3863 gzstop = zstop = 1 ;
3864 gzincr = zincr = 1 ;
3865 }
3866 starttime = gridptr->dimmin[GTINDEX] ;
3867 tcount = gridptr->dimmax[GTINDEX] - starttime + 1 ;
3868 if (tindex != -1) {
3869 count[tindex] = 1 ; /* read each time step separately for templating */
3870 }
3871 gtstart = tstart = 0 ;
3872 gtstop = tstop = tcount ;
3873 gtincr = tincr = 1 ;
3874 yxsize = count[yindex] * count[xindex] ;
3875 gyxsize = count[yindex] * xsize ; /*output array has full xdim */
3876 zyxsize = (gzstop - gzstart) * yxsize ;
3877 gzyxsize = (gzstop - gzstart) * gyxsize ;
3878 gridcnt *= tcount ; /*???*/
3879 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3880 /* need to check for _FillValue as well as missing_value */
3881 FVattr = find_att((*data_var_ptr)->first_vattr, "_FillValue") ;
3882 }
3883 if (!ispacked) {
3884 if (pktype == NC_FLOAT) {
3885 if (missvalattr) {
3886 fdata = (float *) missvalattr->data ;
3887 fmiss = *fdata ;
3888 } else {
3889 fmiss = FILL_FLOAT ;
3890 }
3891 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3892 if (FVattr) {
3893 fdata = (float *) FVattr->data ;
3894 fFillValue = *fdata ;
3895 } else {
3896 fFillValue = FILL_FLOAT ;
3897 }
3898 } /* !usingFV */
3899 } else if (pktype == NC_DOUBLE) {
3900 if (missvalattr) {
3901 ddata = (double *) missvalattr->data ;
3902 dmiss = *ddata ;
3903 } else {
3904 dmiss = FILL_DOUBLE ;
3905 }
3906 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3907 if (FVattr) {
3908 ddata = (double *) FVattr->data ;
3909 dFillValue = *ddata ;
3910 } else {
3911 dFillValue = FILL_DOUBLE ;
3912 }
3913 } /* !usingFV */
3914 /* else not yet a supported unpacked type */
3915 }
3916 } else {
3917 if ((!missvalattr) && (pktype != unpktype)) {
3918 smiss = (int) (((gridptr->undef - ao) / sf) + 0.5) ;
3919 } else {
3920 if ((!missvalattr) && (pktype == unpktype)) {
3921 smiss = FILL_SHORT ;
3922 } else {
3923 sdata = (short *) missvalattr->data ;
3924 smiss = *sdata ;
3925 }
3926 }
3927 if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))) {
3928 if (FVattr) {
3929 sdata = (short *) FVattr->data ;
3930 sFillValue = *sdata ;
3931 } else {
3932 sFillValue = FILL_SHORT ;
3933 }
3934 } /* !usingFV */
3935 } /* ispacked */
3936 alltimesinthisfile = 0 ;
3937 alltimesin1file = 0 ;
3938 if (gridptr->pfile->tmplat) {
3939 if ((*(gridptr->pfile->fnums + starttime - 1)) == gridptr->pfile->fnumc) {
3940 /* this time slice starts in the current file, so let's see if */
3941 /* it ends in the same file */
3942 if ((*(gridptr->pfile->fnums + starttime - 1 + tstop - tstart - 1)) ==
3943 gridptr->pfile->fnumc) /* ends in same file */ {
3944 alltimesinthisfile = 1 ;
3945 alltimesin1file = 1 ;
3946 }
3947 } else /* first time not in current file */ {
3948 if ((*(gridptr->pfile->fnums + starttime - 1)) ==
3949 (*(gridptr->pfile->fnums + starttime - 1 + tstop - tstart - 1))) {
3950 alltimesin1file = 1 ;
3951 }
3952 } /* is first time in current file */
3953 } /* if templating */
3954 if ((!(gridptr->pfile->tmplat)) || alltimesinthisfile) {
3955 if (tindex != -1) {
3956 start[tindex] = starttime - 1 ;
3957 count[tindex] = tstop - tstart ;
3958 } /* if tindex is valid => we have a time dim. */
3959 if (gridptr->pfile->tmplat) {
3960 if ((gridptr->pfile->fnums[0]) != (gridptr->pfile->fnumc)) {
3961 /* if we are not on the first file, we still need to call gaopsdf to */
3962 /* adulterate the starttime to be one relative to the current file */
3963 /* -Hoop, 2001/04/14 */
3964 newtime = gaopsdf(starttime, gridptr, data_var_ptr) ;
3965 if (newtime == -99999) /* bad surprise */ {
3966 gaprnt(0, "Bad time index in SDF file.\n") ;
3967 return(1) ; /* bad time index */
3968 } /* bad surprise */
3969 if (newtime == -88888) /* bad file, bad surprise */ {
3970 for (grinx = 0 ; grinx < gridcnt ; ++grinx) {
3971 fgrid[grinx] = gridptr->undef ;
3972 }
3973 return(0) ;
3974 } /* bad file, bad surprise */
3975 if (tindex != -1) {
3976 start[tindex] = newtime - 1 ;
3977 }
3978 } /* if not 1st file, thus need new starttime */
3979 } /* if templating */
3980 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid,
3981 (*data_var_ptr), start, count)) {
3982 gaprnt(0, "Failure reading SDF data. #4\n") ;
3983 return(3) ;
3984 }
3985 } /* end if not templating or allin1, read all times at once */ else {
3986 if (alltimesin1file) {
3987 if (tindex != -1) {
3988 start[tindex] = starttime - 1 ;
3989 count[tindex] = tstop - tstart ;
3990 } /* if tindex is valid => we have a time dim. */
3991 newtime = gaopsdf(starttime, gridptr, data_var_ptr) ;
3992 if (newtime == -99999) {
3993 gaprnt(0, "Bad time index in SDF file.\n") ;
3994 return(1) ; /* bad time index */
3995 }
3996 if (newtime == -88888) /* bad file */ {
3997 for (grinx = 0 ; grinx < gridcnt ; ++grinx) {
3998 fgrid[grinx] = gridptr->undef ;
3999 }
4000 return(0) ;
4001 } /* if bad file */
4002 data_var = *data_var_ptr ;
4003 if (tindex != -1) {
4004 start[tindex] = newtime - 1 ;
4005 }
4006 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid,
4007 (*data_var_ptr), start, count)) {
4008 gaprnt(0, "Failure reading SDF data. #5\n") ;
4009 return(3) ;
4010 }
4011 } /* still do all in one read, but after gaopsdf */
4012 } /* templating, not all in this file */
4013 for (timeindex = starttime, gtinx = gtstart, tinx = tstart ;
4014 (gtinx < gtstop) && (tinx < tstop) ;
4015 gtinx += gtincr, tinx += tincr, timeindex += tincr) {
4016 newtime = timeindex ;
4017 if (gridptr->pfile->tmplat) {
4018 if ((!alltimesinthisfile) && (!alltimesin1file)) {
4019 newtime = gaopsdf(timeindex, gridptr, data_var_ptr) ;
4020 if (newtime == -99999) {
4021 gaprnt(0, "Bad time index in SDF file.\n") ;
4022 return(1) ; /* bad time index */
4023 }
4024 if (newtime == -88888) /* bad file */ {
4025 fillwithmiss = 1 ;
4026 } /* if bad file */
4027 data_var = *data_var_ptr ;
4028 }
4029 } /* if templating in use */
4030 if (!fillwithmiss) {
4031 if (tindex != -1) {
4032 start[tindex] = newtime - 1 ;
4033 }
4034 }
4035 #ifdef TESTING
4036 fprintf(stderr,
4037 "read_variable_data(%d, %s, [%ld, %ld, %ld, %ld], [%ld, %ld, %ld, %ld])\n",
4038 gridptr->pfile->sdf_ptr->cdfid, (*data_var_ptr)->varnam, start[0], start[1],
4039 start[2], start[3], count[0], count[1], count[2], count[3]) ;
4040 #endif
4041 if (!fillwithmiss) {
4042 /* experiment to increase performance when not templating Hoop 2K/11/22 */
4043 if (gridptr->pfile->tmplat && (!alltimesinthisfile) &&
4044 (!alltimesin1file)) {
4045 /* no need to read per time step if not templating or all times in one file; */
4046 /* all read above */
4047 if (!read_variable_data(gridptr->pfile->sdf_ptr->cdfid,
4048 (*data_var_ptr), start, count)) {
4049 gaprnt(0, "Failure reading SDF data. #6\n") ;
4050 return(3) ;
4051 }
4052 } /* end if templating and not all in one file */
4053 if (!ispacked) {
4054 if (pktype == NC_DOUBLE) {
4055 ddata = (double *) (*data_var_ptr)->data ;
4056 } else {
4057 fdata = (float *) (*data_var_ptr)->data ;
4058 }
4059 } else /* must unpack */ {
4060 sdata = (short *) (*data_var_ptr)->data ;
4061 } /* is it packed? */
4062 } /* only do above if not fillwithmiss */
4063 for (gzinx = gzstart, zinx = zstart ;
4064 (gzinx < gzstop) && ((zincr > 0)?(zinx < zstop):(zinx > zstop)) ;
4065 gzinx += gzincr, zinx += zincr) {
4066 for (gyinx = gystart, yinx = ystart ;
4067 (gyinx < gystop) && ((yincr > 0)?(yinx < ystop):(yinx > ystop)) ;
4068 gyinx += gyincr, yinx += yincr) {
4069 for (xreq = xreqmin, xpos = 0 ; xreq <= xreqmax ; ++xreq, ++xpos) {
4070 xact = xreq ;
4071 while (xact < 1) { xact += xcount ; }
4072 while (xact > xcount) { xact -= xcount ; }
4073 grinx = (gtinx * gzyxsize) + (gzinx * gyxsize) +
4074 (gyinx * xsize) + xpos ;
4075 if (gridptr->pfile->yrflg) {
4076 ycnt = yinx - ystop - 1 ;
4077 } else {
4078 ycnt = yinx - ystart ;
4079 }
4080 if (is4d) {
4081 if (gridptr->pfile->zrflg) {
4082 zcnt = zinx - zstop - 1 ;
4083 } else {
4084 zcnt = zinx - zstart ;
4085 }
4086 } else {
4087 zcnt = 1 ;
4088 }
4089 tcnt = tinx - tstart ;
4090 if (is4d) {
4091 inx = (zcnt * yxsize) + (ycnt * sdfxcount) + xact - 1 ;
4092 } else {
4093 inx = (ycnt * sdfxcount) + xact - 1 ;
4094 }
4095 /* experiment to increase performance when not templating Hoop 2K/11/22 */
4096 if ((!(gridptr->pfile->tmplat)) || alltimesinthisfile
4097 || alltimesin1file) {
4098 inx += (tinx * zyxsize) ;
4099 }
4100 #ifdef TESTING2
4101 fprintf(stderr, "grinx=%d, inx=%d.\n", grinx, inx) ;
4102 #endif
4103 if (fillwithmiss) {
4104 fgrid[grinx] = gridptr->undef ;
4105 } else {
4106 if (ispacked) {
4107 /* handle missing */
4108 wasmiss = 0 ;
4109 if (sdata[inx] == smiss) {
4110 wasmiss = 1 ;
4111 } else if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))){
4112 if (sdata[inx] == sFillValue) {
4113 wasmiss = 1 ;
4114 }
4115 }
4116 if (!wasmiss) {
4117 fgrid[grinx] = (sf * ((float) sdata[inx])) + ao ;
4118 } else {
4119 fgrid[grinx] = gridptr->undef ;
4120 }
4121 } else {
4122 /* may wish to swap missing values */
4123 /* Do want to. -Hoop, 2K/07/25 */
4124 if (pktype == NC_DOUBLE) {
4125 wasmiss = 0 ;
4126 if (dequal(ddata[inx], dmiss, (double)1.0e-5)) {
4127 wasmiss = 1 ;
4128 } else if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))){
4129 if (dequal(ddata[inx], dFillValue, (double)1.0e-5)) {
4130 wasmiss = 1 ;
4131 }
4132 }
4133 if (!wasmiss) {
4134 fgrid[grinx] = (float) ddata[inx] ;
4135 } else {
4136 fgrid[grinx] = gridptr->undef ;
4137 }
4138 } else if (pktype == NC_FLOAT) {
4139 wasmiss = 0 ;
4140 if (fequal(fdata[inx], fmiss, 1.0e-5)) {
4141 wasmiss = 1 ;
4142 } else if ((!usingFV) && (!(gridptr->pfile->sdf_ptr->hasDDFundef))){
4143 if(fequal(fdata[inx], fFillValue, 1.0e-5)) {
4144 wasmiss = 1 ;
4145 }
4146 }
4147 if (!wasmiss) {
4148 fgrid[grinx] = fdata[inx] ;
4149 } else {
4150 fgrid[grinx] = gridptr->undef ;
4151 }
4152 /* else using unsupported (not)packed type */
4153 }
4154 } /* end if packed */
4155 } /* end if fillwithmiss */
4156 } /* for x if wrapped */
4157 } /* end for y indexing */
4158 } /* end for z indexing */
4159 } /* end for time indexing */
4160 return(0) ;
4161 }
4162
4163 int
convtype(void * srcptr,nc_type srctype,void * desptr,nc_type destype,int srcindex,int desindex)4164 convtype(void *srcptr, nc_type srctype, void *desptr, nc_type destype,
4165 int srcindex, int desindex) {
4166 short *ssrcptr, *sdesptr ;
4167 long *lsrcptr, *ldesptr ;
4168 float *fsrcptr, *fdesptr ;
4169 double *dsrcptr, *ddesptr ;
4170
4171 switch (srctype) {
4172 case NC_SHORT:
4173 ssrcptr = (short *) srcptr ;
4174 switch (destype) {
4175 case NC_SHORT: /* stupid, but... */
4176 sdesptr = (short *) desptr ;
4177 sdesptr[desindex] = (short) (ssrcptr[srcindex]) ;
4178 break ;
4179 case NC_LONG:
4180 ldesptr = (long *) desptr ;
4181 ldesptr[desindex] = (long) (ssrcptr[srcindex]) ;
4182 break ;
4183 case NC_FLOAT:
4184 fdesptr = (float *) desptr ;
4185 fdesptr[desindex] = (float) (ssrcptr[srcindex]) ;
4186 break ;
4187 case NC_DOUBLE:
4188 ddesptr = (double *) desptr ;
4189 ddesptr[desindex] = (double) (ssrcptr[srcindex]) ;
4190 break ;
4191 default:
4192 return(0) ; /* inappropriate type for conversion */
4193 break ;
4194 } /* end destype switch */
4195 break ;
4196 case NC_LONG:
4197 lsrcptr = (long *) srcptr ;
4198 switch (destype) {
4199 case NC_SHORT:
4200 sdesptr = (short *) desptr ;
4201 sdesptr[desindex] = (short) (lsrcptr[srcindex]) ;
4202 break ;
4203 case NC_LONG: /* stupid, but... */
4204 ldesptr = (long *) desptr ;
4205 ldesptr[desindex] = (long) (lsrcptr[srcindex]) ;
4206 break ;
4207 case NC_FLOAT:
4208 fdesptr = (float *) desptr ;
4209 fdesptr[desindex] = (float) (lsrcptr[srcindex]) ;
4210 break ;
4211 case NC_DOUBLE:
4212 ddesptr = (double *) desptr ;
4213 ddesptr[desindex] = (double) (lsrcptr[srcindex]) ;
4214 break ;
4215 default:
4216 return(0) ; /* inappropriate type for conversion */
4217 break ;
4218 } /* end destype switch */
4219 break ;
4220 case NC_FLOAT:
4221 fsrcptr = (float *) srcptr ;
4222 switch (destype) {
4223 case NC_SHORT:
4224 sdesptr = (short *) desptr ;
4225 sdesptr[desindex] = (short) (fsrcptr[srcindex]) ;
4226 break ;
4227 case NC_LONG:
4228 ldesptr = (long *) desptr ;
4229 ldesptr[desindex] = (long) (fsrcptr[srcindex]) ;
4230 break ;
4231 case NC_FLOAT: /* stupid, but... */
4232 fdesptr = (float *) desptr ;
4233 fdesptr[desindex] = (float) (fsrcptr[srcindex]) ;
4234 break ;
4235 case NC_DOUBLE:
4236 ddesptr = (double *) desptr ;
4237 ddesptr[desindex] = (double) (fsrcptr[srcindex]) ;
4238 break ;
4239 default:
4240 return(0) ; /* inappropriate type for conversion */
4241 break ;
4242 } /* end destype switch */
4243 break ;
4244 case NC_DOUBLE:
4245 dsrcptr = (double *) srcptr ;
4246 switch (destype) {
4247 case NC_SHORT:
4248 sdesptr = (short *) desptr ;
4249 sdesptr[desindex] = (short) (dsrcptr[srcindex]) ;
4250 break ;
4251 case NC_LONG:
4252 ldesptr = (long *) desptr ;
4253 ldesptr[desindex] = (long) (dsrcptr[srcindex]) ;
4254 break ;
4255 case NC_FLOAT:
4256 fdesptr = (float *) desptr ;
4257 fdesptr[desindex] = (float) (dsrcptr[srcindex]) ;
4258 break ;
4259 case NC_DOUBLE: /* stupid, but... */
4260 ddesptr = (double *) desptr ;
4261 ddesptr[desindex] = (double) (dsrcptr[srcindex]) ;
4262 break ;
4263 default:
4264 return(0) ; /* inappropriate type for conversion */
4265 break ;
4266 } /* end destype switch */
4267 break ;
4268 default:
4269 return(0) ; /* inappropriate type for conversion */
4270 break ;
4271 } /* end srctype switch */
4272 return(1) ; /* conversion successful */
4273 /* NOTREACHED */
4274 }
4275
4276 #include "gasdf_std_time.h"
4277 /* Old RCS ID: ReadNetCDF.c,v 1.1 1995/05/02 20:23:12 jac Exp jac $ */
4278 /* Revision 1.1 1995/05/02 20:23:12 jac (Julia Collins) */
4279 /* Initial revision */
4280
4281 /* open and read a netcdf file */
4282 /* print variable and dimension information to standard output */
4283
4284 #define True 1
4285 #define False 0
4286
4287 int
print_dim_info(std_ptr,dim,all)4288 print_dim_info (std_ptr, dim, all)
4289 IO_STD *std_ptr;
4290 int dim;
4291 int all; /* whether to print all values. If False, */
4292 /* only the first and last values are shown */
4293 {
4294 VAR_INFO *var = NULL;
4295 union att_val dim_val;
4296 int i;
4297 int is_time = False;
4298 char tstring[128];
4299 utUnit unit;
4300
4301 int decode_time ();
4302 void set_data_to_type ();
4303 void adjust_data_type ();
4304 VAR_INFO *find_var ();
4305
4306 /* print dim info if it exists */
4307 if (dim >= 0)
4308 {
4309 if ((var = find_var (std_ptr, std_ptr->dimnam[dim])) == NULL)
4310 return Failure;
4311 if (read_one_dimension (std_ptr, var, dim) != Success)
4312 return Failure;
4313
4314 /* print the dimension name */
4315 printf ("%s\n", var->varnam);
4316
4317 if (!strcmp (var->varnam, dims[TIME_IX]))
4318 {
4319 is_time = True;
4320 if (std_ptr->time_type == COOP)
4321 get_ud_time_unit (var, &unit);
4322 }
4323
4324 /* print the data values */
4325 if (all)
4326 {
4327 for (i = 0; i < std_ptr->dimsiz[dim]; i++)
4328 {
4329 set_data_to_type (var->vartype, i, var->data, &dim_val);
4330 if (is_time)
4331 {
4332 decode_time (dim_val.dval, std_ptr->time_type, var, tstring, &unit);
4333 printf ("%s\n", tstring);
4334 }
4335 else
4336 {
4337 adjust_data_type (var->vartype, NC_FLOAT, &dim_val);
4338 printf ("%.2f\n", dim_val.fval);
4339 }
4340 }
4341 }
4342 else
4343 {
4344 /* print first value */
4345 set_data_to_type (var->vartype, 0, var->data, &dim_val);
4346 if (is_time)
4347 {
4348 decode_time (dim_val.dval, std_ptr->time_type, var, tstring, &unit);
4349 printf ("%s\n", tstring);
4350 }
4351 else
4352 {
4353 adjust_data_type (var->vartype, NC_FLOAT, &dim_val);
4354 printf ("%.2f\n", dim_val.fval);
4355 }
4356
4357 /* print last value */
4358 set_data_to_type (var->vartype, std_ptr->dimsiz[dim] - 1,
4359 var->data, &dim_val);
4360 if (is_time)
4361 {
4362 decode_time (dim_val.dval, std_ptr->time_type, var, tstring, &unit);
4363 printf ("%s\n", tstring);
4364 }
4365 else
4366 {
4367 adjust_data_type (var->vartype, NC_FLOAT, &dim_val);
4368 printf ("%.2f\n", dim_val.fval);
4369 }
4370 }
4371 printf ("END\n");
4372 }
4373
4374 return Success;
4375 }
4376
4377 int
decode_time(tval,time_type,var,tstring,unit)4378 decode_time (tval, time_type, var, tstring, unit)
4379 double tval; /* value to decode */
4380 int time_type; /* whether time is formatted according to "old" CDC */
4381 /* standards or according to udunits */
4382 VAR_INFO *var;
4383 char *tstring;
4384 utUnit *unit;
4385 {
4386 int yr, mo, da, hr, min;
4387 double sec;
4388 struct attrib_list *delta;
4389 int pix;
4390 char delta_t_str[128];
4391
4392 struct attrib_list *find_att ();
4393 int decode_ud_time() ;
4394 int decode_standard_time() ;
4395 void set_default_time_components ();
4396
4397 /* get delta_t */
4398 delta = find_att (var->first_vattr, time_atts[DELTA_T_IX]);
4399 ((char *) delta->data)[delta->len] = '\0';
4400
4401 /* get delta_t values */
4402 if (decode_delta_t ((char *) delta->data, &yr, &mo, &da, &hr,
4403 &min, &sec) != Success)
4404 pix = -1;
4405 else if (get_period (yr, mo, da, hr, &pix) != Success)
4406 pix = -1;
4407
4408 if (time_type == CDC)
4409 {
4410 if (decode_standard_time (tval, &yr, &mo, &da,
4411 &hr, &min, &sec) != Success)
4412 return Failure;
4413 }
4414 else
4415 {
4416 if (decode_ud_time (unit, tval, &yr, &mo, &da, &hr, &min, &sec) != Success)
4417 return Failure;
4418 }
4419
4420 /* now format the string */
4421 set_default_time_components (pix, &yr, &mo, &da, &hr, &min, &sec);
4422 encode_time (yr, mo, da, hr, min, (int) sec, tstring);
4423
4424 return Success;
4425 }
4426
4427 int
get_period(year,month,day,hour,pix)4428 get_period (year, month, day, hour, pix)
4429 int year;
4430 int month;
4431 int day;
4432 int hour;
4433 int *pix;
4434 {
4435 if ((year != MISSING) && (year > 0))
4436 *pix = PYEARS;
4437 else if ((month != MISSING) && (month > 0))
4438 {
4439 /* this can be months or seasons */
4440 if (month % MONTHS_PER_SEASON)
4441 /* month is multiple of size of season, so assume seasonal file */
4442 *pix = PMONTHS;
4443 else
4444 *pix = PSEASONS;
4445 }
4446 else if ((day != MISSING) && (day > 0))
4447 {
4448 /* This can be days or pentads */
4449 if (day % DAYS_PER_PENTAD)
4450 *pix = PDAYS;
4451 else
4452 *pix = PPENTADS;
4453 }
4454 else if ((hour != MISSING) && (hour > 0))
4455 *pix = PHOURS;
4456 else
4457 {
4458 *pix = MISSING;
4459 return Failure;
4460 }
4461
4462 return Success;
4463 }
4464
4465 /* set any missing components to a default beginning value */
4466 void
set_default_time_components(pix,yr,mo,day,hour,min,sec)4467 set_default_time_components (pix, yr, mo, day, hour, min, sec)
4468 int pix;
4469 int *yr;
4470 int *mo;
4471 int *day;
4472 int *hour;
4473 int *min;
4474 double *sec;
4475 {
4476 if (pix < 0)
4477 return;
4478
4479 switch (pix)
4480 {
4481 case PMONTHS:
4482 case PSEASONS:
4483 if (*mo < 0)
4484 *mo = 1;
4485 *day = MISSING;
4486 *hour = MISSING;
4487 *min = MISSING;
4488 *sec = (double) MISSING;
4489 break;
4490
4491 case PDAYS:
4492 case PPENTADS:
4493 if (*day < 0)
4494 *day = 1;
4495 *hour = MISSING;
4496 *min = MISSING;
4497 *sec = (double) MISSING;
4498 break;
4499
4500 case PHOURS:
4501 if (*hour < 0)
4502 *hour = 0;
4503 *min = MISSING;
4504 *sec = (double) MISSING;
4505 break;
4506
4507 default:
4508 break;
4509 }
4510 return;
4511 }
4512
4513 int
encode_time(yr,mo,da,hr,min,sec,tstring)4514 encode_time (yr, mo, da, hr, min, sec, tstring)
4515 int yr;
4516 int mo;
4517 int da;
4518 int hr;
4519 int min;
4520 int sec;
4521 char *tstring;
4522 {
4523 static char *month[13] = {" ", "Jan", "Feb", "Mar", "Apr",
4524 "May", "June", "July", "Aug",
4525 "Sep", "Oct", "Nov", "Dec" };
4526 char tstr[32];
4527
4528 strcpy (tstring, " ");
4529 if (mo != MISSING)
4530 {
4531 sprintf (tstring, "%s ", month[mo]);
4532 }
4533 if (da != MISSING)
4534 {
4535 sprintf (tstr, "%d ", da);
4536 strcat (tstring, tstr);
4537 }
4538 if (hr != MISSING)
4539 {
4540 sprintf (tstr, "%.2d", hr);
4541 strcat (tstring, tstr);
4542 if (min != MISSING)
4543 {
4544 sprintf (tstr, ":%.2d", min);
4545 strcat (tstring, tstr);
4546 if (sec != MISSING)
4547 {
4548 sprintf (tstr, "%.2d ", sec);
4549 strcat (tstring, tstr);
4550 }
4551 }
4552 }
4553 if (yr != MISSING)
4554 {
4555 sprintf (tstr, "%d", yr);
4556 strcat (tstring, tstr);
4557 }
4558
4559 return Success;
4560 }
4561
4562 /* get the udunit structure for time */
get_ud_time_unit(time,unit)4563 get_ud_time_unit (time, unit)
4564 VAR_INFO *time; /* time variable information */
4565 utUnit *unit;
4566 {
4567 struct attrib_list *units_str;
4568 struct attrib_list *find_att ();
4569
4570 /* get units string */
4571 if ((units_str = find_att (time->first_vattr,
4572 time_atts[T_UNITS_IX])) == NULL)
4573 return Failure;
4574
4575 /* scan the units */
4576 if (utScan ((char *) units_str->data, unit))
4577 return Failure;
4578
4579 return Success;
4580 }
4581
4582 /* use the udunits library to decode a time value */
4583 int
decode_ud_time(unit,time_val,yr,mo,da,hr,min,sec)4584 decode_ud_time (unit, time_val, yr, mo, da, hr, min, sec)
4585 utUnit *unit; /* time unit information */
4586 double time_val; /* input time value */
4587 int *yr; /* (returned) decoded year */
4588 int *mo; /* (returned) decoded month */
4589 int *da; /* (returned) decoded day */
4590 int *hr; /* (returned) decoded hour */
4591 int *min; /* (returned) decoded minute */
4592 double *sec; /* (returned) decoded second */
4593 {
4594
4595 /* convert the time to integer components */
4596 if (utCalendar (time_val, unit, yr, mo, da, hr, min, (float *) sec))
4597 return Failure;
4598
4599 return Success;
4600 }
4601
4602 /* Strip the given number of characters from the string and return the
4603 new length. If the number of characters to strip is less than
4604 or equal to zero, or if the strip length is greater than the
4605 string length, return Failure. Otherwise, return Success. */
4606 int
strip_char(strip_num,str1,int_len)4607 strip_char (strip_num, str1, int_len)
4608 int strip_num; /* Number of characters to strip. */
4609 char *str1; /* (input and output) String */
4610 int *int_len; /* (input and output) String length */
4611 {
4612 int slen;
4613
4614 slen = strlen (str1);
4615 if (strip_num <= 0)
4616 return Failure;
4617 slen -= strip_num;
4618 if (slen < 0)
4619 return Failure;
4620 *int_len = slen;
4621 str1[slen] = '\0';
4622
4623 return Success;
4624 }
4625
4626 /* Decode standard time. Return Success if OK, Failure if error. */
4627 int
decode_standard_time(time_val,year,month,day,hour,minn,sec)4628 decode_standard_time (time_val, year, month, day, hour, minn, sec)
4629 double time_val; /* Time value to be decoded. */
4630 int *year; /* (returned) year */
4631 int *month; /* (returned) month */
4632 int *day; /* (returned) day */
4633 int *hour; /* (returned) hour */
4634 int *minn; /* (returned) minute */
4635 double *sec; /* (returned) second */
4636 {
4637 char str1[20];
4638 int i,
4639 slen,
4640 int_len;
4641
4642 /* Make time value into character string. */
4643 sprintf (str1, "%f", time_val);
4644
4645 /* Find decimal point. */
4646 slen = strlen (str1);
4647 for (i = 0; i < slen; i++)
4648 {
4649 if (str1[i] == '.')
4650 {
4651 int_len = i;
4652 break;
4653 }
4654 }
4655 if (int_len == 0)
4656 {
4657 return Failure;
4658 }
4659
4660 /* Get second. */
4661 if (int_len <= 2)
4662 {
4663 sscanf (str1, "%lf", sec);
4664 int_len = 0;
4665 }
4666 else
4667 {
4668 sscanf (&str1[int_len - 2], "%lf", sec);
4669 str1[int_len] = '\0';
4670 (void) strip_char (2, str1, &int_len);
4671 }
4672
4673 /* Get minute. */
4674 *minn = MISSING;
4675 if (int_len > 0)
4676 {
4677 if (int_len <= 2)
4678 {
4679 sscanf (str1, "%d", minn);
4680 int_len = 0;
4681 }
4682 else
4683 {
4684 sscanf (&str1[int_len - 2], "%d", minn);
4685 str1[int_len] = '\0';
4686 (void) strip_char (2, str1, &int_len);
4687 }
4688 }
4689
4690 /* Get hour. */
4691 *hour = MISSING;
4692 if (int_len > 0)
4693 {
4694 if (int_len <= 2)
4695 {
4696 sscanf (str1, "%d", hour);
4697 int_len = 0;
4698 }
4699 else
4700 {
4701 sscanf (&str1[int_len - 2], "%d", hour);
4702 str1[int_len] = '\0';
4703 (void) strip_char (2, str1, &int_len);
4704 }
4705 }
4706
4707 /* Get day. */
4708 *day = MISSING;
4709 if (int_len > 0)
4710 {
4711 if (int_len <= 2)
4712 {
4713 sscanf (str1, "%d", day);
4714 int_len = 0;
4715 }
4716 else
4717 {
4718 sscanf (&str1[int_len - 2], "%d", day);
4719 str1[int_len] = '\0';
4720 (void) strip_char (2, str1, &int_len);
4721 }
4722 }
4723
4724 /* Get month. */
4725 *month = MISSING;
4726 if (int_len > 0)
4727 {
4728 if (int_len <= 2)
4729 {
4730 sscanf (str1, "%d", month);
4731 int_len = 0;
4732 }
4733 else
4734 {
4735 sscanf (&str1[int_len - 2], "%d", month);
4736 str1[int_len] = '\0';
4737 (void) strip_char (2, str1, &int_len);
4738 }
4739 }
4740
4741 /* Get year. A year of 0000 or 9999 defaults to missing. */
4742 *year = MISSING;
4743 if (int_len > 0)
4744 sscanf (str1, "%d", year);
4745 if ((*year == 0) || (*year == 9999))
4746 *year = MISSING;
4747
4748 /* All OK. */
4749 return Success;
4750 }
4751
4752
4753 void
adjust_data_type(old_type,new_type,val)4754 adjust_data_type (old_type, new_type, val)
4755 nc_type old_type; /* current data type */
4756 nc_type new_type; /* type to cast data TO */
4757 union att_val *val; /* current values, updated and returned */
4758
4759 {
4760 union att_val save_val;
4761
4762 switch (new_type)
4763 {
4764 case NC_BYTE:
4765 switch (old_type)
4766 {
4767 case NC_SHORT:
4768 save_val.sval = val->sval;
4769 if (val->sval == SMISS)
4770 val->bval = BMISS;
4771 else if (val->sval == SFILL)
4772 val->bval = BFILL;
4773 else
4774 val->bval = (char) save_val.sval;
4775 break;
4776 case NC_LONG:
4777 save_val.lval = val->lval;
4778 if (val->lval == LMISS)
4779 val->bval = BMISS;
4780 else if (val->lval == LFILL)
4781 val->bval = BFILL;
4782 else
4783 val->bval = (char) save_val.lval;
4784 break;
4785 case NC_FLOAT:
4786 save_val.fval = val->fval;
4787 if (dbl_equal (val->fval, FMISS))
4788 val->bval = BMISS;
4789 else if (dbl_equal (val->fval, FFILL))
4790 val->bval = BFILL;
4791 else
4792 val->bval = (char) save_val.fval;
4793 break;
4794 case NC_DOUBLE:
4795 save_val.dval = val->dval;
4796 if (dbl_equal (val->dval, DMISS))
4797 val->bval = BMISS;
4798 else if (dbl_equal (val->dval, DFILL))
4799 val->bval = BFILL;
4800 else
4801 val->bval = (char) save_val.dval;
4802 break;
4803 default:
4804 break;
4805 }
4806 break;
4807 case NC_SHORT:
4808 switch (old_type)
4809 {
4810 case NC_BYTE:
4811 save_val.bval = val->bval;
4812 if (val->bval == BMISS)
4813 val->sval = SMISS;
4814 else if (val->bval == BFILL)
4815 val->sval = SFILL;
4816 else
4817 val->sval = (short) save_val.bval;
4818 break;
4819 case NC_LONG:
4820 save_val.lval = val->lval;
4821 if (val->lval == LMISS)
4822 val->sval = SMISS;
4823 else if (val->lval == LFILL)
4824 val->sval = SFILL;
4825 else
4826 val->sval = (short) save_val.lval;
4827 break;
4828 case NC_FLOAT:
4829 save_val.fval = val->fval;
4830 if (dbl_equal (val->fval, FMISS))
4831 val->sval = SMISS;
4832 else if (dbl_equal (val->fval, FFILL))
4833 val->sval = SFILL;
4834 else
4835 val->sval = (short) save_val.fval;
4836 break;
4837 case NC_DOUBLE:
4838 save_val.dval = val->dval;
4839 if (dbl_equal (val->dval, DMISS))
4840 val->sval = SMISS;
4841 else if (dbl_equal (val->dval, DFILL))
4842 val->sval = SFILL;
4843 else
4844 val->sval = (short) save_val.dval;
4845 break;
4846 default:
4847 break;
4848 }
4849 break;
4850 case NC_LONG:
4851 switch (old_type)
4852 {
4853 case NC_SHORT:
4854 save_val.sval = val->sval;
4855 if (val->sval == SMISS)
4856 val->lval = LMISS;
4857 else if (val->sval == SFILL)
4858 val->lval = LFILL;
4859 else
4860 val->lval = (long) save_val.sval;
4861 break;
4862 case NC_BYTE:
4863 save_val.bval = val->bval;
4864 if (val->bval == BMISS)
4865 val->lval = LMISS;
4866 else if (val->bval == BFILL)
4867 val->lval = LFILL;
4868 else
4869 val->lval = (long) save_val.bval;
4870 break;
4871 case NC_FLOAT:
4872 save_val.fval = val->fval;
4873 if (dbl_equal (val->fval, FMISS))
4874 val->lval = LMISS;
4875 else if (dbl_equal (val->fval, FFILL))
4876 val->lval = LFILL;
4877 else
4878 val->lval = (long) save_val.fval;
4879 break;
4880 case NC_DOUBLE:
4881 save_val.dval = val->dval;
4882 if (dbl_equal (val->dval, DMISS))
4883 val->lval = LMISS;
4884 else if (dbl_equal (val->dval, DFILL))
4885 val->lval = LFILL;
4886 else
4887 val->lval = (long) save_val.dval;
4888 break;
4889 default:
4890 break;
4891 }
4892 break;
4893 case NC_FLOAT:
4894 switch (old_type)
4895 {
4896 case NC_SHORT:
4897 save_val.sval = val->sval;
4898 if (val->sval == SMISS)
4899 val->fval = FMISS;
4900 else if (val->sval == SFILL)
4901 val->fval = FFILL;
4902 else
4903 val->fval = (float) save_val.sval;
4904 break;
4905 case NC_LONG:
4906 save_val.lval = val->lval;
4907 if (val->lval == LMISS)
4908 val->fval = FMISS;
4909 else if (val->lval == LFILL)
4910 val->fval = FFILL;
4911 else
4912 val->fval = (float) save_val.lval;
4913 break;
4914 case NC_BYTE:
4915 save_val.bval = val->bval;
4916 if (val->bval == BMISS)
4917 val->fval = FMISS;
4918 else if (val->bval == BFILL)
4919 val->fval = FFILL;
4920 else
4921 val->fval = (float) save_val.bval;
4922 break;
4923 case NC_DOUBLE:
4924 save_val.dval = val->dval;
4925 if (dbl_equal (val->dval, DMISS))
4926 val->fval = FMISS;
4927 else if (dbl_equal (val->dval, DFILL))
4928 val->fval = FFILL;
4929 else
4930 val->fval = (float) save_val.dval;
4931 break;
4932 default:
4933 break;
4934 }
4935 break;
4936 case NC_DOUBLE:
4937 switch (old_type)
4938 {
4939 case NC_SHORT:
4940 save_val.sval = val->sval;
4941 if (val->sval == SMISS)
4942 val->dval = DMISS;
4943 else if (val->sval == SFILL)
4944 val->dval = DFILL;
4945 else
4946 val->dval = (double) save_val.sval;
4947 break;
4948 case NC_LONG:
4949 save_val.lval = val->lval;
4950 if (val->lval == LMISS)
4951 val->dval = DMISS;
4952 else if (val->lval == LFILL)
4953 val->dval = DFILL;
4954 else
4955 val->dval = (double) save_val.lval;
4956 break;
4957 case NC_FLOAT:
4958 save_val.fval = val->fval;
4959 if (dbl_equal (val->fval, FMISS))
4960 val->dval = DMISS;
4961 else if (dbl_equal (val->fval, FFILL))
4962 val->dval = DFILL;
4963 else
4964 val->dval = (double) save_val.fval;
4965 break;
4966 case NC_BYTE:
4967 save_val.bval = val->bval;
4968 if (val->bval == BMISS)
4969 val->dval = DMISS;
4970 else if (val->bval == BFILL)
4971 val->dval = DFILL;
4972 else
4973 val->dval = (double) save_val.bval;
4974 break;
4975 default:
4976 break;
4977 }
4978 break;
4979 default:
4980 break;
4981 }
4982
4983 return;
4984 }
4985
4986 /* cast void data to correct type */
4987 void
set_data_to_type(type,ix,in_data,out_data)4988 set_data_to_type (type, ix, in_data, out_data)
4989 nc_type type; /* data type */
4990 int ix; /* if data is greater than one item long,
4991 * index into data array */
4992 void *in_data; /* data to be cast */
4993 union att_val *out_data; /* output data (returned) */
4994
4995 {
4996
4997 switch (type)
4998 {
4999 case NC_BYTE:
5000 out_data->bval = ((char *) in_data)[ix];
5001 break;
5002 case NC_SHORT:
5003 out_data->sval = ((short *) in_data)[ix];
5004 break;
5005 case NC_LONG:
5006 out_data->lval = ((long *) in_data)[ix];
5007 break;
5008 case NC_FLOAT:
5009 out_data->fval = ((float *) in_data)[ix];
5010 break;
5011 case NC_DOUBLE:
5012 out_data->dval = ((double *) in_data)[ix];
5013 break;
5014 default:
5015 break;
5016 }
5017
5018 return;
5019 }
5020
5021 /* initialize a netCDF standard file structure */
5022 int
init_io_std(std_ptr)5023 init_io_std (std_ptr)
5024 IO_STD **std_ptr; /* pointer to netCDF data structure */
5025 {
5026 int init_dim_info ();
5027 void free_io_std ();
5028
5029 if ((*std_ptr) != NULL)
5030 free_io_std (std_ptr);
5031
5032 if (((*std_ptr) = (IO_STD *) malloc (sizeof (IO_STD))) == NULL) {
5033 #if USEHDF == 1
5034 printf ("Could not allocate memory for netCDF/HDF-SDS structure.\n");
5035 #else
5036 printf ("Could not allocate memory for netCDF structure.\n") ;
5037 #endif
5038 return Failure;
5039 }
5040
5041 /* initialize contents of structure */
5042 (*std_ptr)->cdfid = -1;
5043 (*std_ptr)->ndims = 0;
5044 (*std_ptr)->nvars = 0;
5045 (*std_ptr)->ngatts = 0;
5046 (*std_ptr)->recdim = -1;
5047 (*std_ptr)->time_type = CDC;
5048 (*std_ptr)->hasDDFundef = 0 ;
5049
5050 (*std_ptr)->first_gattr = NULL;
5051
5052 /* don't create variable list quite yet */
5053 (*std_ptr)->var = NULL;
5054
5055 /* intialize dimension information */
5056 if (init_dim_info ((*std_ptr)->dimids, (*std_ptr)->dimnam,
5057 (*std_ptr)->dimsiz) != Success)
5058 return Failure;
5059
5060 return Success;
5061 }
5062
5063 /* copy a netCDF standard file structure */
5064 int
copy_io_std(IO_STD ** std_ptr2,IO_STD * std_ptr1)5065 copy_io_std (IO_STD **std_ptr2, IO_STD *std_ptr1)
5066 /* pointers to netCDF data structure */
5067 {
5068 int indx ;
5069 int init_dim_info ();
5070 int copy_attr_list() ;
5071 int copy_var_list() ;
5072 void free_io_std ();
5073
5074 if ((*std_ptr2) != NULL)
5075 free_io_std (std_ptr2);
5076
5077 if (((*std_ptr2) = (IO_STD *) malloc (sizeof (IO_STD))) == NULL) {
5078 #if USEHDF == 1
5079 printf ("Could not allocate memory for netCDF/HDF-SDS structure.\n");
5080 #else
5081 printf ("Could not allocate memory for netCDF structure.\n") ;
5082 #endif
5083 return Failure;
5084 }
5085
5086 /* initialize contents of structure */
5087 (*std_ptr2)->cdfid = std_ptr1->cdfid ;
5088 (*std_ptr2)->ndims = std_ptr1->ndims ;
5089 (*std_ptr2)->nvars = std_ptr1->nvars ;
5090 (*std_ptr2)->ngatts = std_ptr1->ngatts ;
5091 (*std_ptr2)->recdim = std_ptr1->recdim ;
5092 (*std_ptr2)->time_type = std_ptr1->time_type ;
5093 if (copy_attr_list(&((*std_ptr2)->first_gattr), std_ptr1->first_gattr) ==
5094 Failure) {
5095 return Failure ;
5096 }
5097 /* don't create variable list quite yet */
5098 (*std_ptr2)->var = NULL; /* copy the linked list.... */
5099
5100 /* initialize dimension information */
5101 for (indx = 0 ; indx < MAX_NC_DIMS ; ++indx) {
5102 (*std_ptr2)->dimids[indx] = std_ptr1->dimids[indx] ; /* cp the list.... */
5103 strcpy((*std_ptr2)->dimnam[indx], std_ptr1->dimnam[indx]) ;
5104 (*std_ptr2)->dimsiz[indx] = std_ptr1->dimsiz[indx] ;
5105 } /* end for all dims */
5106
5107 if (copy_var_list(&((*std_ptr2)->var), std_ptr1->var) == Failure) {
5108 return Failure ;
5109 }
5110
5111 return Success;
5112 }
5113
5114 int
copy_var_list(VAR_INFO ** newvarlistptr,VAR_INFO * oldvarlist)5115 copy_var_list(VAR_INFO **newvarlistptr, VAR_INFO *oldvarlist) {
5116 int indx ;
5117 char *charptr1, *charptr2 ;
5118 short *shortptr1, *shortptr2 ;
5119 long *longptr1, *longptr2 ;
5120 float *floatptr1, *floatptr2 ;
5121 double *dblptr1, *dblptr2 ;
5122 VAR_INFO **lclnewvarlistptr ;
5123 int copy_attr_list() ;
5124
5125 if (oldvarlist == (VAR_INFO *) NULL) {
5126 return Failure ;
5127 }
5128 lclnewvarlistptr = newvarlistptr ;
5129 while (oldvarlist != (VAR_INFO *) NULL) {
5130 if (((*lclnewvarlistptr) = (VAR_INFO *)malloc(sizeof(VAR_INFO))) ==
5131 (VAR_INFO *) NULL) {
5132 return Failure ;
5133 }
5134 (*lclnewvarlistptr)->varid = oldvarlist->varid ;
5135 strcpy((*lclnewvarlistptr)->varnam, oldvarlist->varnam) ;
5136 strcpy((*lclnewvarlistptr)->gradsabbr, oldvarlist->gradsabbr) ;
5137 (*lclnewvarlistptr)->vartype = oldvarlist->vartype ;
5138 (*lclnewvarlistptr)->nvardims = oldvarlist->nvardims ;
5139 for (indx = 0 ; indx < oldvarlist->nvardims ; ++indx) {
5140 (*lclnewvarlistptr)->vardimid[indx] = oldvarlist->vardimid[indx] ;
5141 }
5142 for (indx = 0 ; indx < 4 ; ++indx) {
5143 (*lclnewvarlistptr)->dimmap[indx] = oldvarlist->dimmap[indx] ;
5144 (*lclnewvarlistptr)->dimidmap[indx] = oldvarlist->dimidmap[indx] ;
5145 }
5146 (*lclnewvarlistptr)->nvatts = oldvarlist->nvatts ;
5147 if (oldvarlist->nvatts > 0) {
5148 if (copy_attr_list(&((*lclnewvarlistptr)->first_vattr),
5149 oldvarlist->first_vattr) == Failure) {
5150 return Failure ;
5151 }
5152 } else {
5153 (*lclnewvarlistptr)->first_vattr = (struct attrib_list *) NULL ;
5154 }
5155 /* forget about messing with data member of struct */
5156 oldvarlist = oldvarlist->next ;
5157 if (oldvarlist != (VAR_INFO *) NULL) {
5158 lclnewvarlistptr = &((*lclnewvarlistptr)->next) ;
5159 }
5160 /* following added by jma to fix redhat linux problem */
5161 else sprintf(pout, " ") ;
5162 } /* end while there's still a list cell to copy */
5163 return Success ;
5164 }
5165
5166 int
copy_attr_list(struct attrib_list ** first_list1,struct attrib_list * first_list2)5167 copy_attr_list(struct attrib_list **first_list1,
5168 struct attrib_list *first_list2) {
5169 struct attrib_list *this_struct_ptr, *new_struct_ptr ;
5170 int indx ;
5171 char *charptr1, *charptr2 ;
5172 short *shortptr1, *shortptr2 ;
5173 long *longptr1, *longptr2 ;
5174 float *floatptr1, *floatptr2 ;
5175 double *dblptr1, *dblptr2 ;
5176
5177 if (first_list2 == (struct attrib_list*)NULL) {
5178 /* return Failure ; */
5179 *first_list1 == (struct attrib_list*)NULL;
5180 return Success ;
5181 }
5182 if ((*first_list1 =
5183 (struct attrib_list *) malloc(sizeof(struct attrib_list))) ==
5184 (struct attrib_list *) NULL) {
5185 return Failure ;
5186 }
5187 for(this_struct_ptr = first_list2, new_struct_ptr = *first_list1;
5188 this_struct_ptr != (struct attrib_list *)NULL ;
5189 this_struct_ptr = this_struct_ptr->next) {
5190 new_struct_ptr->next = this_struct_ptr->next ;
5191 new_struct_ptr->name = strdup(this_struct_ptr->name) ;
5192 new_struct_ptr->type = this_struct_ptr->type ;
5193 new_struct_ptr->len = this_struct_ptr->len ;
5194 if ((new_struct_ptr->data =
5195 (void *) malloc((size_t)(new_struct_ptr->len *
5196 nctypelen(new_struct_ptr->type))))
5197 == (void *) NULL) {
5198 return Failure ;
5199 }
5200 switch (new_struct_ptr->type) {
5201 case NC_BYTE:
5202 case NC_CHAR:
5203 charptr1 = (char *) new_struct_ptr->data ;
5204 charptr2 = (char *) this_struct_ptr->data ;
5205 break ;
5206 case NC_SHORT:
5207 shortptr1 = (short *) new_struct_ptr->data ;
5208 shortptr2 = (short *) this_struct_ptr->data ;
5209 break ;
5210 case NC_LONG:
5211 longptr1 = (long *) new_struct_ptr->data ;
5212 longptr2 = (long *) this_struct_ptr->data ;
5213 break ;
5214 case NC_FLOAT:
5215 floatptr1 = (float *) new_struct_ptr->data ;
5216 floatptr2 = (float *) this_struct_ptr->data ;
5217 break ;
5218 case NC_DOUBLE:
5219 dblptr1 = (double *) new_struct_ptr->data ;
5220 dblptr2 = (double *) this_struct_ptr->data ;
5221 break ;
5222 } /* esac */
5223 for (indx = 0 ; indx < new_struct_ptr->len ; ++indx) {
5224 switch (new_struct_ptr->type) {
5225 case NC_BYTE:
5226 case NC_CHAR:
5227 charptr1[indx] = charptr2[indx] ;
5228 break ;
5229 case NC_SHORT:
5230 shortptr1[indx] = shortptr2[indx] ;
5231 break ;
5232 case NC_LONG:
5233 longptr1[indx] = longptr2[indx] ;
5234 break ;
5235 case NC_FLOAT:
5236 floatptr1[indx] = floatptr2[indx] ;
5237 break ;
5238 case NC_DOUBLE:
5239 dblptr1[indx] = dblptr2[indx] ;
5240 break ;
5241 } /* esac */
5242 } /* end for */
5243 if (this_struct_ptr->next != (struct attrib_list *)NULL) {
5244 if ((new_struct_ptr->next =
5245 (struct attrib_list *) malloc(sizeof(struct attrib_list)))
5246 == (struct attrib_list *)NULL) {
5247 return Failure ;
5248 } /* if malloc failure */
5249 } /* if there's more list to do, malloc next cell */
5250 } /* while linked list isn't over */
5251 return Success ;
5252 }
5253
5254
5255 /* initialize dimension information */
5256 int
init_dim_info(dimids,dimnam,dimsiz)5257 init_dim_info (dimids, dimnam, dimsiz)
5258 int *dimids; /* array of dimension IDs */
5259 char dimnam[MAX_NC_DIMS][MAX_NC_NAME + 1]; /* array of dimension
5260 * names */
5261 long *dimsiz; /* array of dimension sizes */
5262 {
5263 int i;
5264
5265 /* set all dimension information to bogus values */
5266 for (i = 0; i < MAX_NC_DIMS; i++)
5267 {
5268 dimids[i] = -1;
5269 dimnam[i][0] = '\0';
5270 dimsiz[i] = -1;
5271 }
5272
5273 return Success;
5274 }
5275
5276 /* initialize a variable structure */
5277 int
init_var_info(var)5278 init_var_info (var)
5279 VAR_INFO **var; /* variable data structure */
5280 {
5281 int i;
5282
5283 void free_var_info ();
5284
5285 if (*var != NULL)
5286 free_var_info (var);
5287
5288 /* get memory */
5289 if ((*var = (VAR_INFO *) malloc (sizeof (VAR_INFO))) == NULL)
5290 {
5291 printf ("Could not allocate memory for variable information.\n");
5292 return Failure;
5293 }
5294
5295 /* initialize contents to bogus values */
5296 (*var)->next = NULL;
5297 (*var)->varid = -1;
5298 (*var)->varnam[0] = '\0';
5299 (*var)->gradsabbr[0] = '\0' ;
5300 (*var)->vartype = -1;
5301 (*var)->nvardims = 0;
5302 (*var)->nvatts = 0;
5303
5304 for (i = 0; i < MAX_VAR_DIMS; i++)
5305 (*var)->vardimid[i] = -1;
5306 for (i = 0 ; i < 4 ; ++i) {
5307 (*var)->dimmap[i] = -1 ;
5308 (*var)->dimidmap[i] = -1 ;
5309 }
5310
5311 (*var)->first_vattr = NULL;
5312 (*var)->data = NULL;
5313 (*var)->dimmap[0] = (*var)->dimmap[1] = (*var)->dimmap[2] =
5314 (*var)->dimmap[3] = -1 ;
5315 (*var)->dimidmap[0] = (*var)->dimidmap[1] = (*var)->dimidmap[2] =
5316 (*var)->dimidmap[3] = -1 ;
5317 return Success;
5318 }
5319
5320 /* free a netCDF standard file structure */
5321 void
free_io_std(std_ptr)5322 free_io_std (std_ptr)
5323 IO_STD **std_ptr;
5324 {
5325 void free_var_info ();
5326 int free_netcdf_att_list ();
5327
5328 if (*std_ptr != NULL)
5329 {
5330 /* if the variable list exists, free it */
5331 if ((*std_ptr)->var != NULL)
5332 free_var_info (&((*std_ptr)->var));
5333
5334 /* free the global attributes if they have been allocated */
5335 if (((*std_ptr)->first_gattr != NULL) && ((*std_ptr)->ngatts > 0)) {
5336 free_netcdf_att_list (&((*std_ptr)->first_gattr));
5337 }
5338
5339 /* now free the IO_STD structure itself */
5340 if (*std_ptr) {
5341 free (*std_ptr);
5342 }
5343 *std_ptr = NULL;
5344 }
5345
5346 return;
5347 }
5348
5349 /* free a variable structure */
5350 void
free_var_info(var)5351 free_var_info (var)
5352 VAR_INFO **var;
5353 {
5354 void free_var_info ();
5355 int free_netcdf_att_list ();
5356
5357 /* go through list to the end and free all variable structures */
5358 if ((*var)->next != NULL)
5359 free_var_info (&((*var)->next));
5360
5361 if (((*var)->first_vattr != NULL) && ((*var)->nvatts > 0)) {
5362 free_netcdf_att_list (&((*var)->first_vattr));
5363 }
5364
5365 /* if data space has been allocated, free it */
5366 if ((*var)->data != NULL) {
5367 free ((*var)->data);
5368 }
5369
5370 /* now free the variable structure itself */
5371 if (*var) {
5372 free (*var);
5373 }
5374 *var = NULL;
5375
5376 return;
5377 }
5378
5379 /* Free a netCDF attribute list. */
5380 int
free_netcdf_att_list(first_attr)5381 free_netcdf_att_list (first_attr)
5382 struct attrib_list **first_attr; /* First attribute in list. */
5383 {
5384 struct attrib_list *temp_attr,
5385 *save_attr;
5386
5387 temp_attr = *first_attr;
5388 while (temp_attr != NULL) {
5389 save_attr = temp_attr->next;
5390 if (temp_attr->data != NULL) {
5391 free (temp_attr->data);
5392 }
5393 free (temp_attr);
5394 temp_attr = save_attr;
5395 }
5396 *first_attr = NULL;
5397 return Success;
5398 }
5399
5400 /* open and read the contents of a netCDF file */
5401 /* std_ptr should already be initialized by init_io_std */
5402 int
read_io_std(path,std_ptr)5403 read_io_std (path, std_ptr)
5404 char *path; /* file pathname */
5405 IO_STD *std_ptr; /* data structure to be filled */
5406 {
5407 int inq_netcdf (),
5408 open_netcdf (),
5409 close_netcdf (),
5410 read_atts (),
5411 read_dims (),
5412 read_vars ();
5413 void handle_error ();
5414 int set_time_type ();
5415
5416 if ((path == NULL) || (strlen (path) == 0))
5417 return Failure;
5418
5419 if (open_netcdf (path, 0, &(std_ptr->cdfid)) != Success)
5420 {
5421 #if USEHDF == 1
5422 printf ("%s does not exist or is not a netCDF/HDF-SDS file.\n", path);
5423 #else
5424 printf ("%s does not exist or is not a netCDF file.\n", path) ;
5425 #endif
5426 return Failure;
5427 }
5428
5429 /* Get general information. If error, return. */
5430 if (inq_netcdf (std_ptr->cdfid, &(std_ptr->ndims), &(std_ptr->nvars),
5431 &(std_ptr->ngatts), &(std_ptr->recdim)) != Success)
5432 {
5433 (void) close_netcdf (std_ptr->cdfid);
5434 return Failure;
5435 }
5436 if (std_ptr->ndims > MAX_NC_DIMS) {
5437 std_ptr->ndims = MAX_NC_DIMS ;
5438 fprintf(stderr, "Only first %d dimensions processed; library limit.\n",
5439 MAX_NC_DIMS) ;
5440 }
5441 /* Create the global attribute list. */
5442 if (read_atts (std_ptr->cdfid, NC_GLOBAL, std_ptr->ngatts,
5443 &(std_ptr->first_gattr)) != Success)
5444 {
5445 (void) close_netcdf (std_ptr->cdfid);
5446 return Failure;
5447 }
5448
5449 /* get dimension information */
5450 if (read_dims (std_ptr) != Success)
5451 {
5452 (void) close_netcdf (std_ptr->cdfid);
5453 return Failure;
5454 }
5455
5456 /* inquire about variable(s) */
5457 if (read_vars (std_ptr->cdfid, std_ptr->nvars, &(std_ptr->var)) != Success)
5458 {
5459 (void) close_netcdf (std_ptr->cdfid);
5460 return Failure;
5461 }
5462
5463 /* determine if new or old units are being used */
5464 if (set_time_type (std_ptr) != Success)
5465 ;
5466 /* return Failure; */
5467
5468 /* set up standard tables according to time unit being used */
5469 if (init_standard_arrays (std_ptr->time_type) != Success)
5470 ;
5471 /* return Failure; */
5472
5473 return Success;
5474 }
5475
5476 int
init_standard_arrays(time_type)5477 init_standard_arrays (time_type)
5478 int time_type;
5479 {
5480 if (time_type == CDC)
5481 {
5482 dims = cdc_dims;
5483 vars = cdc_vars;
5484 var_type = cdc_var_type;
5485 var_atts = cdc_var_atts;
5486 var_atts_type = cdc_var_atts_type;
5487 var_atts_val = cdc_var_atts_val;
5488 obs_atts_val = cdc_obs_atts_val;
5489 vatts_abbrev = cdc_vatts_abbrev;
5490 time_atts = cdc_time_atts;
5491 time_atts_val = cdc_time_atts_val;
5492 latlon_atts = cdc_latlon_atts;
5493 num_reqd_vatts = NUM_REQD_VATTS;
5494 num_reqd_vars = NUM_REQD_VARS;
5495 num_reqd_dims = NUM_REQD_DIMS;
5496 }
5497 else
5498 {
5499 dims = coop_dims;
5500 vars = coop_vars;
5501 var_type = coop_var_type;
5502 var_atts = coop_var_atts;
5503 var_atts_type = coop_var_atts_type;
5504 var_atts_val = coop_var_atts_val;
5505 obs_atts_val = coop_obs_atts_val;
5506 vatts_abbrev = coop_vatts_abbrev;
5507 time_atts = coop_time_atts;
5508 time_atts_val = coop_time_atts_val;
5509 latlon_atts = coop_latlon_atts;
5510 num_reqd_vatts = NUM_REQD_COOP_VATTS;
5511 num_reqd_vars = NUM_REQD_COOP_VARS;
5512 num_reqd_dims = NUM_REQD_COOP_DIMS;
5513 }
5514
5515 return Success;
5516 }
5517
5518 int
set_time_type(in_ptr)5519 set_time_type (in_ptr)
5520 IO_STD *in_ptr;
5521 {
5522 VAR_INFO *time = NULL, *lclvar;
5523 char *ch ;
5524 utUnit timeunit ;
5525 struct attrib_list *units;
5526 struct attrib_list *find_att ();
5527 VAR_INFO *find_var ();
5528
5529 /* the name of the time variable and its units attribute is the same */
5530 /*for both the cdc and cooperative standards. Since the "vars" and */
5531 /* "time_atts" arrays are not yet initialized, use the cdc variable */
5532 /* and attributes definitions */
5533 time = find_var (in_ptr, cdc_vars[TIME_IX]);
5534 if (!time) {
5535 /* Apply COARDSian heuristics on units to find a time-type coordvar */
5536 /* If that fails, complain and return failure */
5537 for (lclvar = in_ptr->var ; lclvar ; lclvar = lclvar->next) {
5538 if ((lclvar->nvardims) != 1) continue ;
5539 units = find_att(lclvar->first_vattr, cdc_time_atts[T_UNITS_IX]) ;
5540 if (!units) continue ;
5541 ch = (char *) malloc(units->len + 1) ;
5542 strncpy(ch, (char *) units->data, units->len) ;
5543 ch[units->len] = '\0' ;
5544 if (utScan(ch, &timeunit)) {
5545 if (ch) {
5546 free(ch) ;
5547 }
5548 continue ;
5549 }
5550 if (ch) {
5551 free(ch) ;
5552 }
5553 if (utIsTime(&timeunit)) {
5554 time = lclvar ;
5555 break ;
5556 }
5557 } /* end for */
5558 } /* end if !time */
5559 if (!time) /* still no time; give up */ {
5560 return Failure ;
5561 }
5562 if ((units = find_att (time->first_vattr, cdc_time_atts[T_UNITS_IX]))
5563 == NULL)
5564 return Failure;
5565
5566 if (!strncasecmp ("yyyy", (char *) units->data, 4))
5567 /* it's the old style */
5568 in_ptr->time_type = CDC;
5569 else
5570 in_ptr->time_type = COOP;
5571
5572 return Success;
5573 }
5574
5575 /* return an attribute structure, given the name of the attribute */
5576 struct attrib_list *
find_att(first_att,attname)5577 find_att (first_att, attname)
5578 struct attrib_list *first_att; /* attribute list to search */
5579 char *attname; /* attribute to find */
5580 {
5581 static struct attrib_list *temp = NULL;
5582
5583 temp = first_att;
5584 while (temp != NULL) {
5585 if (!strcmp (attname, temp->name)) break;
5586 temp = temp->next;
5587 }
5588 return (struct attrib_list *) temp;
5589 }
5590
5591 /* Open the netCDF pathname for input or output. Return <0 if error;
5592 else return 0. */
5593 int
open_netcdf(pathname,cflag,cdfid)5594 open_netcdf (pathname, cflag, cdfid)
5595 char *pathname; /* file to open */
5596 int cflag; /* for output, clobber flag (NC_CLOBBER,NC_NOCLOBBER) */
5597 int *cdfid; /* (returned) netCDF id */
5598 {
5599 #if USEHDF ==1
5600 int oldncopts ; /* to save and restore */
5601 oldncopts = ncopts ;
5602 ncopts = 0; /* Turn off automatic error handling. */
5603 /* Open input file */
5604 *cdfid = ncopen (pathname, NC_NOWRITE);
5605 if (*cdfid == -1) {
5606 ncopts = oldncopts ;
5607 return Failure;
5608 }
5609 ncopts = oldncopts ;
5610 #else
5611 int error_code ; /* version 3 NetCDF API error handling */
5612 error_code = NC_NOERR ;
5613 /* Open input file */
5614 error_code = nc_open(pathname, NC_NOWRITE, cdfid) ;
5615 if (error_code != NC_NOERR) {
5616 fprintf(stderr, "nc_open failure:\n%s\n", nc_strerror(error_code)) ;
5617 return(Failure) ;
5618 }
5619 #endif
5620
5621 return Success;
5622 }
5623
5624
5625 /* Close the netCDF file. Return <0 if error; else return 0. */
5626 int
close_netcdf(cdfid)5627 close_netcdf (cdfid)
5628 int cdfid; /* netCDF id */
5629 {
5630 int oldncopts ; /* to save and restore */
5631
5632 oldncopts = ncopts ;
5633 /* Turn off automatic error handling. */
5634 ncopts = 0;
5635
5636 if (ncclose (cdfid) == -1) {
5637 ncopts = oldncopts ;
5638 return Failure;
5639 }
5640
5641 ncopts = oldncopts ;
5642 return Success;
5643 }
5644
5645
5646 /* Inquire about the netCDF file. Return <0 if error; else return 0. */
5647 int
inq_netcdf(cdfid,ndims,nvars,ngatts,recdim)5648 inq_netcdf (cdfid, ndims, nvars, ngatts, recdim)
5649 int cdfid; /* netCDF id */
5650 int *ndims; /* (returned) number of dimensions */
5651 int *nvars; /* (returned) number of variables */
5652 int *ngatts; /* (returned) number of global attributes */
5653 int *recdim; /* (returned) record dimension, if any */
5654 {
5655 int oldncopts ; /* to save and restore */
5656
5657 oldncopts = ncopts ;
5658 /* Turn off automatic error handling. */
5659 ncopts = 0;
5660
5661 if (ncinquire (cdfid, ndims, nvars, ngatts, recdim) == -1) {
5662 /* Possible errors include NC_EBADID (invalid netCDF id), which could
5663 * mean the netCDF id is bad or the file is closed. */
5664 ncopts = oldncopts ;
5665 return Failure;
5666 }
5667
5668 ncopts = oldncopts ;
5669 return Success;
5670 }
5671
5672 /* Inquire about a netCDF variable. Return <0 if error; else return 0. */
5673 int
inq_netcdf_var(cdfid,varid,varname,vartype,nvardims,vardims,natts)5674 inq_netcdf_var (cdfid, varid, varname, vartype, nvardims, vardims, natts)
5675 int cdfid; /* netCDF id */
5676 int varid; /* variable id */
5677 char *varname; /* (returned) variable name */
5678 nc_type *vartype;
5679
5680 /* (returned) variable data type */
5681 int *nvardims; /* (returned) number of dimensions */
5682 int *vardims; /* (returned) array of dimension numbers */
5683 int *natts; /* (returned) number of attributes */
5684 {
5685 int oldncopts ; /* to save and restore */
5686
5687 oldncopts = ncopts ;
5688 /* Turn off automatic error handling. */
5689 ncopts = 0;
5690
5691 /* Inquire about the variable. */
5692 if (ncvarinq (cdfid, varid, varname, vartype, nvardims,
5693 vardims, natts) == -1) {
5694 ncopts = oldncopts ;
5695 return Failure;
5696 }
5697
5698 ncopts = oldncopts ;
5699 return Success;
5700 }
5701
5702 void
init_att(new_att,msg)5703 init_att (new_att, msg)
5704 struct attrib_list **new_att; /* attribute to create and initialize */
5705 char *msg; /* error message string */
5706 {
5707 if ((*new_att = (struct attrib_list *) malloc (sizeof (struct attrib_list)))
5708 == NULL)
5709 return;
5710
5711 /* initialize all attribute components */
5712 (*new_att)->next = NULL;
5713 (*new_att)->name = NULL;
5714 (*new_att)->type = NC_UNSPECIFIED;
5715 (*new_att)->len = 0;
5716 (*new_att)->data = NULL;
5717
5718 return;
5719 }
5720
5721 /* read attribute information for a file or variable */
5722 int
read_atts(cdfid,varid,natts,first_attr)5723 read_atts (cdfid, varid, natts, first_attr)
5724 int cdfid; /* netCDF file ID */
5725 int varid; /* variable ID or NC_GLOBAL for global
5726 * attribute read */
5727 int natts; /* number of attributes to retrieve */
5728 struct attrib_list **first_attr;/* first attribute in linked list */
5729 {
5730 int i;
5731 struct attrib_list *curr,
5732 *prev;
5733 char name[MAX_NC_NAME + 1];
5734
5735 int get_netcdf_attname (),
5736 get_netcdf_att_info ();
5737 int free_netcdf_att_list ();
5738 void init_att ();
5739
5740 curr = prev = NULL;
5741
5742 /* first make sure nothing is already in the list */
5743 if (*first_attr != NULL)
5744 free_netcdf_att_list (first_attr);
5745
5746 /* loop on all attributes */
5747 for (i = 0; i < natts; i++) {
5748 /* allocate and initialize attribute space */
5749 init_att (&curr, "read_atts");
5750
5751 /* update pointers */
5752 if (*first_attr == NULL)
5753 *first_attr = curr;
5754 else
5755 prev->next = curr;
5756
5757 curr->next = NULL;
5758 prev = curr;
5759
5760 /* get attribute name */
5761 if (get_netcdf_attname (cdfid, varid, i, name) != Success)
5762 return Failure;
5763 copy_char_str (&(curr->name), name);
5764
5765 /* now get attribute information */
5766 if (get_netcdf_att_info (cdfid, varid, curr->name,
5767 &(curr->type), &(curr->len), &(curr->data))
5768 != Success)
5769 return Failure;
5770 }
5771 return Success;
5772 }
5773
5774 /* Get attribute name. Return <0 if error; else return 0. */
5775 int
get_netcdf_attname(cdfid,varid,attnum,name)5776 get_netcdf_attname (cdfid, varid, attnum, name)
5777 int cdfid; /* netcdf id */
5778 int varid; /* variable id */
5779 int attnum; /* attribute number */
5780 char *name; /* (returned) attribute name */
5781 {
5782 int oldncopts ; /* to save and restore */
5783
5784 oldncopts = ncopts ;
5785 /* Turn off automatic error handling. */
5786 ncopts = 0;
5787
5788 if (ncattname (cdfid, varid, attnum, name) == -1) {
5789 ncopts = oldncopts ;
5790 return Failure;
5791 }
5792
5793 /* All OK. */
5794 ncopts = oldncopts ;
5795 return Success;
5796 }
5797
5798
5799 /* Get attribute information. If error, return <0; else, return 0. */
5800 int
get_netcdf_att_info(cdfid,varid,aname,atype,alen,adata)5801 get_netcdf_att_info (cdfid, varid, aname, atype, alen, adata)
5802 int cdfid; /* netcdf id */
5803 int varid; /* variable id */
5804 char *aname; /* attribute name */
5805 nc_type *atype; /* (returned) attribute type */
5806 int *alen; /* (returned) attribute length */
5807 void **adata; /* (returned) attribute data */
5808 {
5809 int oldncopts ; /* to save and restore */
5810 char *char_ptr;
5811
5812 oldncopts = ncopts ;
5813 /* Turn off automatic error handling. */
5814 ncopts = 0;
5815
5816 if (ncattinq (cdfid, varid, aname, atype, alen) == -1) {
5817 ncopts = oldncopts ;
5818 return Failure;
5819 }
5820
5821 /* Allocate space for values. */
5822 if ((*atype == NC_BYTE) || (*atype == NC_CHAR)) {
5823 *adata = (char *) malloc (sizeof (char) * (*alen + 1));
5824 } else if (*atype == NC_SHORT) {
5825 *adata = (short *) malloc (sizeof (short) * *alen);
5826 } else if (*atype == NC_LONG) {
5827 *adata = (long *) malloc (sizeof (long) * *alen);
5828 } else if (*atype == NC_FLOAT) {
5829 *adata = (float *) malloc (sizeof (float) * *alen);
5830 } else if (*atype == NC_DOUBLE) {
5831 *adata = (double *) malloc (sizeof (double) * *alen);
5832 } else {
5833 ncopts = oldncopts ;
5834 return Failure;
5835 }
5836
5837 if (*adata == NULL) {
5838 ncopts = oldncopts ;
5839 return Failure;
5840 }
5841
5842 /* Retrieve values. */
5843 if (ncattget (cdfid, varid, aname, (void *) (*adata)) == -1) {
5844 ncopts = oldncopts ;
5845 return Failure;
5846 }
5847
5848 /* If character, set last byte to null. */
5849 if (*atype == NC_CHAR) {
5850 char_ptr = (char *) *adata;
5851 char_ptr[*alen] = '\0';
5852 }
5853
5854 /* All OK. */
5855 ncopts = oldncopts ;
5856 return Success;
5857 }
5858
5859 /* read all dimensions in a file */
5860 int
read_dims(std_ptr)5861 read_dims (std_ptr)
5862 IO_STD *std_ptr; /* pointer to netCDF data structure */
5863 {
5864 int i;
5865 int inq_netcdf_dim ();
5866 int oldncopts ; /* to save and restore */
5867
5868 oldncopts = ncopts ;
5869 for (i = 0; i < std_ptr->ndims; i++) {
5870 std_ptr->dimids[i] = i;
5871 if (inq_netcdf_dim (std_ptr->cdfid, std_ptr->dimids[i],
5872 std_ptr->dimnam[i], &(std_ptr->dimsiz[i])) != Success) {
5873 ncopts = oldncopts ;
5874 return Failure;
5875 }
5876 }
5877
5878 ncopts = oldncopts ;
5879 return Success;
5880 }
5881
5882 /* Inquire about a netCDF variable. Return <0 if error; else return 0. */
5883 int
inq_netcdf_dim(cdfid,dimid,dimname,dimsize)5884 inq_netcdf_dim (cdfid, dimid, dimname, dimsize)
5885 int cdfid; /* netCDF id */
5886 int dimid; /* dimension id */
5887 char *dimname; /* (returned) dimension name */
5888 long *dimsize; /* (returned) dimension size */
5889 {
5890 int oldncopts ; /* to save and restore */
5891
5892 oldncopts = ncopts ;
5893 /* Turn off automatic error handling. */
5894 ncopts = 0;
5895
5896 /* Inquire about the dimension. */
5897 if (ncdiminq (cdfid, dimid, dimname, (long *) dimsize) == -1) {
5898 /* Possible errors include NC_EBADID (invalid netCDF id), which could
5899 * mean the netCDF id is bad or the file is closed. */
5900 ncopts = oldncopts ;
5901 return Failure;
5902 }
5903
5904 ncopts = oldncopts ;
5905 return Success;
5906 }
5907
5908 /* inquire about variable(s). Note that the actual data aren't */
5909 /* read here, but is read via a call to read_variable_data */
5910 int
read_vars(cdfid,nvars,var)5911 read_vars (cdfid, nvars, var)
5912 int cdfid; /* netCDF file ID */
5913 int nvars; /* number of variables in the file */
5914 VAR_INFO **var; /* pointer to first variable (not
5915 * initialized, will be filled and returned) */
5916 {
5917 int i;
5918 int oldncopts ; /* to save and restore */
5919 VAR_INFO *curr_var = NULL,
5920 *prev_var = NULL;
5921
5922 int init_var_info (),
5923 inq_netcdf_var (),
5924 read_atts ();
5925 void handle_error ();
5926
5927 prev_var = NULL;
5928 oldncopts = ncopts ;
5929
5930 /* loop through all variables and get the information for each */
5931 for (i = 0; i < nvars; i++) {
5932
5933 /* initialize the variable pointer */
5934 curr_var = NULL;
5935 if (init_var_info (&curr_var) != Success) {
5936 ncopts = oldncopts ;
5937 return Failure;
5938 }
5939
5940 /* update pointers */
5941 if (*var == NULL) {
5942 *var = curr_var;
5943 prev_var = curr_var;
5944 } else {
5945 prev_var->next = curr_var;
5946 prev_var = prev_var->next;
5947 }
5948
5949 curr_var->next = NULL;
5950 curr_var->varid = i;
5951
5952 /* inquire about the variable information */
5953 if (inq_netcdf_var (cdfid, i, curr_var->varnam, &(curr_var->vartype),
5954 &(curr_var->nvardims), curr_var->vardimid, &(curr_var->nvatts))
5955 != Success) {
5956 ncopts = oldncopts ;
5957 return Failure;
5958 }
5959
5960 /* get variable attribute information */
5961 if (read_atts (cdfid, curr_var->varid, curr_var->nvatts,
5962 &(curr_var->first_vattr)) != Success) {
5963 if (close_netcdf (cdfid) != Success) {
5964 ncopts = oldncopts ;
5965 return Failure;
5966 }
5967 }
5968
5969 }
5970
5971 ncopts = oldncopts ;
5972 return Success;
5973
5974 }
5975
5976 int
get_additional_dimension_ix(ndims,dimids,timedim,latdim,londim)5977 get_additional_dimension_ix (ndims, dimids, timedim, latdim, londim)
5978 int ndims;
5979 int *dimids;
5980 int timedim;
5981 int latdim;
5982 int londim;
5983 {
5984 int dim_ix = -1;
5985 int i;
5986
5987 for (i = 0; i < ndims; i++) {
5988 if ((dimids[i] != timedim) &&
5989 (dimids[i] != latdim) &&
5990 (dimids[i] != londim)) {
5991 dim_ix = i;
5992 break;
5993 }
5994 }
5995
5996 return dim_ix;
5997 }
5998
5999 /* find a dimension id, given the dimension name */
6000 int
find_dim(std_ptr,name)6001 find_dim (std_ptr, name)
6002 IO_STD *std_ptr; /* pointer to netCDF file information */
6003 char *name; /* name of dimension to find. */
6004 {
6005 int dim = -1;
6006 int i;
6007
6008 for (i = 0; i < std_ptr->ndims; i++)
6009 {
6010 if (!strcmp (std_ptr->dimnam[i], name))
6011 {
6012 dim = std_ptr->dimids[i];
6013 break;
6014 }
6015 }
6016
6017 return dim;
6018 }
6019
6020 /* find the index in the dimension id array, given the dimension id */
6021 int
find_dimix(dimids,ndims,dim)6022 find_dimix (dimids, ndims, dim)
6023 int dimids[]; /* array of dimension ids to search */
6024 int ndims; /* number of dimensions in list */
6025 int dim; /* dimension id to find */
6026 {
6027 int dimix = -1;
6028 int i;
6029
6030 for (i = 0; i < ndims; i++)
6031 {
6032 if (dimids[i] == dim)
6033 {
6034 dimix = i;
6035 break;
6036 }
6037 }
6038
6039 return dimix;
6040 }
6041
6042 int
read_one_dimension(in_ptr,var,dimid)6043 read_one_dimension (in_ptr, var, dimid)
6044 IO_STD *in_ptr;
6045 VAR_INFO *var;
6046 int dimid; /* dimension to read */
6047 {
6048 long *start ;
6049 long *count ;
6050
6051 void init_start_count ();
6052
6053 /* set start and count to read coordinate variable values */
6054 init_start_count (&start, &count, var->nvardims) ;
6055 count[0] = in_ptr->dimsiz[dimid] ; /* traditional unidata coord. var. defn. */
6056
6057 if (read_variable_data (in_ptr->cdfid, var, start, count) != Success)
6058 return Failure;
6059
6060 return Success;
6061 }
6062
6063 /* initialize start and count arrays */
6064 void
init_start_count(start,count,ndims)6065 init_start_count (start, count, ndims)
6066 long **start;
6067 long **count;
6068 int ndims ;
6069
6070 {
6071 int i ;
6072
6073 if (!((*start) = (long *) malloc(sizeof(long) * ndims))) {
6074 fprintf(stderr, "Memory allocation error!\n") ;
6075 exit(13) ;
6076 }
6077 if (!((*count) = (long *) malloc(sizeof(long) * ndims))) {
6078 fprintf(stderr, "Memory allocation error!\n") ;
6079 exit(13) ;
6080 }
6081 for (i = 0 ; i < ndims ; i++) {
6082 (*start)[i] = 0 ;
6083 (*count)[i] = 1 ;
6084 }
6085
6086 return ;
6087 }
6088
6089 /* read variable data */
6090 int
read_variable_data(cdfid,var,start,count)6091 read_variable_data (cdfid, var, start, count)
6092 int cdfid; /* netCDF ID */
6093 VAR_INFO *var; /* pointer to variable structure */
6094 long *start; /* array of start indices */
6095 long *count; /* array of frame counts for each variable
6096 * dimension */
6097 {
6098
6099 int read_netcdf_data ();
6100 int alloc_variable_data ();
6101
6102 if (alloc_variable_data (var, count) != Success)
6103 return Failure;
6104
6105 if (read_netcdf_data (cdfid, var->varid, start, count, var->data)
6106 != Success)
6107 return Failure;
6108
6109 return Success;
6110 }
6111
6112 /* return size of hyperslab */
6113 int
get_hyperslab_size(count,ndims)6114 get_hyperslab_size (count, ndims)
6115 long *count; /* array of hyperslab edge counts */
6116 int ndims;
6117
6118 {
6119 int size;
6120 int i;
6121
6122 size = 1;
6123 for (i = 0; i < ndims; i++)
6124 size *= count[i];
6125
6126 return size;
6127 }
6128
6129 int
alloc_variable_data(var,count)6130 alloc_variable_data (var, count)
6131 VAR_INFO *var; /* variable structure which needs memory */
6132 long *count; /* array of frame counts for each var
6133 * dimension */
6134 {
6135 int malloc_len;
6136 int alloc_data_array ();
6137
6138 /* get rid of any old memory */
6139 if (var->data != NULL)
6140 free (var->data);
6141
6142 /* figure out how much space we need */
6143 malloc_len = get_hyperslab_size (count, var->nvardims);
6144
6145 /* now allocate the data array */
6146 if (alloc_data_array (malloc_len, var->vartype, &(var->data)) != Success)
6147 return Failure;
6148
6149 return Success;
6150 }
6151
6152 /* Read data from a variable in a netCDF input file. Return <0 if error;
6153 else return 0. */
6154 int
read_netcdf_data(cdfid,varid,start,count,data_array)6155 read_netcdf_data (cdfid, varid, start, count, data_array)
6156 int cdfid; /* netCDF id */
6157 int varid; /* variable id */
6158 long *start; /* hyperslab corner */
6159 long *count; /* hyperslab edges */
6160 void *data_array; /* (returned) hyperslab of data values */
6161 {
6162 int rcode;
6163 int oldncopts ; /* to save and restore */
6164
6165 oldncopts = ncopts ;
6166 /* Turn off automatic error handling. */
6167 ncopts = 0;
6168
6169 /* Read the data. */
6170
6171 if (ncvarget (cdfid, varid, (long *) start, (long *) count, data_array) == -1) {
6172 ncopts = oldncopts ;
6173 return Failure;
6174 }
6175
6176 ncopts = oldncopts ;
6177 return Success;
6178 }
6179
6180 /* Allocate one data array. If OK, return Success. If error, return Failure.*/
6181 int
alloc_data_array(malloc_len,datatype,data_array)6182 alloc_data_array (malloc_len, datatype, data_array)
6183 int malloc_len; /* data array length */
6184 nc_type datatype; /* data type */
6185 void **data_array; /* (returned) data array pointer */
6186
6187 {
6188 char str[256];
6189
6190 /* Allocate data array. */
6191 switch (datatype)
6192 {
6193 case (NC_BYTE):
6194 *data_array = (void *) malloc (sizeof (char) * malloc_len);
6195 break;
6196 case (NC_SHORT):
6197 *data_array = (void *) malloc (sizeof (short) * malloc_len);
6198 break;
6199 case (NC_LONG):
6200 *data_array = (void *) malloc (sizeof (long) * malloc_len);
6201 break;
6202 case (NC_FLOAT):
6203 *data_array = (void *) malloc (sizeof (float) * malloc_len);
6204 break;
6205 case (NC_DOUBLE):
6206 *data_array = (void *) malloc (sizeof (double) * malloc_len);
6207 break;
6208 default:
6209 return Failure;
6210 };
6211
6212 if (*data_array == NULL)
6213 return Failure;
6214
6215 /* All OK. */
6216 return Success;
6217 }
6218
6219 /* find a variable pointer, given the variable name */
6220 VAR_INFO *
find_var(std_ptr,var_name)6221 find_var (std_ptr, var_name)
6222 IO_STD *std_ptr;
6223 char *var_name;
6224 {
6225 int i;
6226 static VAR_INFO *found = NULL;
6227 VAR_INFO *curr_var;
6228
6229 found = NULL;
6230 curr_var = std_ptr->var;
6231
6232 for (i = 0; i < std_ptr->nvars; i++)
6233 {
6234 if (curr_var == NULL)
6235 break;
6236
6237 if (!strcmp (var_name, curr_var->varnam))
6238 {
6239 found = curr_var;
6240 break;
6241 }
6242 curr_var = curr_var->next;
6243 }
6244
6245 return (VAR_INFO *) found;
6246 }
6247
6248 /* Check if two single precision values are equal. If they are within */
6249 /* 1/10000000 of each other they are considered equal and 1 is */
6250 /* returned. Else 0 is returned. */
6251 int
flt_equal(f1,f2)6252 flt_equal(f1, f2)
6253 float f1, f2;
6254 {
6255 float ftemp, ftarg1 = -0.0000001, ftarg2 = 0.0000001, flarge, af1, af2;
6256 double dtemp;
6257 int iscale, i;
6258
6259 /* First check if the two numbers are exactly equal. */
6260 if (f1 == f2) {
6261 return 1;
6262 }
6263 /* Next, check if their significant digits are within the specified range. */
6264 af1 = f1;
6265 af2 = f2;
6266 if (af1 < 0) {
6267 af1 *= -1.0 ;
6268 }
6269 if (af2 < 0) {
6270 af2 *= -1.0 ;
6271 }
6272 if (af1 >= af2) {
6273 flarge = af1;
6274 } else {
6275 flarge = af2;
6276 }
6277 /* If the larger number is less than zero, return. */
6278 if (flarge < 1) {
6279 return 0;
6280 }
6281 /* Otherwise, check for significant digits. */
6282 iscale = 1;
6283 for (i = 0; i <= 6; i++) {
6284 dtemp = pow(10.0, (double) i);
6285 if (flarge < (float) dtemp) {
6286 break;
6287 } else {
6288 iscale *= 10;
6289 }
6290 }
6291 ftarg1 *= (float) iscale;
6292 ftarg2 *= (float) iscale;
6293 ftemp = f1 - f2;
6294 if ((ftemp >= ftarg1) && (ftemp <= ftarg2)) {
6295 return 1;
6296 } else {
6297 return 0;
6298 }
6299 }
6300
6301 /* Check if two double precision values are equal. If they are within
6302 1/10000000 of each other they are considered equal and
6303 True is returned. Otherwise False is returned. */
6304 int
dbl_equal(d1,d2)6305 dbl_equal (d1, d2)
6306 double d1,
6307 d2;
6308 {
6309 double dtemp,
6310 dtarg1 = -0.0000001,
6311 dtarg2 = 0.0000001;
6312
6313 /* First check if the two numbers are exactly equal. */
6314 if (d1 == d2)
6315 return Failure;
6316
6317 /* Next, check if they are within the specified range. */
6318 dtemp = d1 - d2;
6319 if ((dtemp >= dtarg1) && (dtemp <= dtarg2))
6320 return Failure;
6321 else
6322 return Success;
6323 }
6324
6325 int
decode_delta_t(delta_t_str,year,month,day,hour,minn,sec)6326 decode_delta_t (delta_t_str, year, month, day, hour, minn, sec)
6327 char *delta_t_str;
6328 int *year;
6329 int *month;
6330 int *day;
6331 int *hour;
6332 int *minn;
6333 int *sec;
6334 {
6335 int year_mark = 4,
6336 month_mark = 7,
6337 day_mark = 10,
6338 hour_mark = 13,
6339 minute_mark = 16,
6340 second_mark = 19;
6341 int delta_t_len;
6342 char temp_str[100];
6343
6344 *year = *month = *day = *hour = *minn = *sec = MISSING;
6345
6346 delta_t_len = strlen (delta_t_str);
6347 if ((delta_t_len > day_mark) && (delta_t_str[day_mark] != ' '))
6348 return Failure;
6349 if ((delta_t_len > hour_mark) && (delta_t_str[hour_mark] != ':'))
6350 return Failure;
6351 if ((delta_t_len > minute_mark) && (delta_t_str[minute_mark] != ':'))
6352 {
6353 return Failure;
6354 }
6355
6356 /* Get year. */
6357 strcpy (temp_str, delta_t_str);
6358 temp_str[year_mark] = '\0';
6359 sscanf (temp_str, "%d", year);
6360
6361 /* Get month. */
6362 strcpy (temp_str, &delta_t_str[year_mark + 1]);
6363 temp_str[month_mark - year_mark - 1] = '\0';
6364 sscanf (temp_str, "%d", month);
6365
6366 /* Get day. */
6367 strcpy (temp_str, &delta_t_str[month_mark + 1]);
6368 temp_str[day_mark - month_mark - 1] = '\0';
6369 sscanf (temp_str, "%d", day);
6370
6371 /* Get other fields if present. */
6372 if (delta_t_len > day_mark)
6373 {
6374
6375 /* Get hour. */
6376 strcpy (temp_str, &delta_t_str[day_mark + 1]);
6377 temp_str[hour_mark - day_mark - 1] = '\0';
6378 sscanf (temp_str, "%d", hour);
6379
6380 /* Get minute. */
6381 strcpy (temp_str, &delta_t_str[hour_mark + 1]);
6382 temp_str[minute_mark - hour_mark - 1] = '\0';
6383 sscanf (temp_str, "%d", minn);
6384
6385 /* Get second. */
6386 strcpy (temp_str, &delta_t_str[minute_mark + 1]);
6387 temp_str[second_mark - minute_mark - 1] = '\0';
6388 sscanf (temp_str, "%d", sec);
6389
6390 }
6391 else
6392 *hour = *minn = *sec = MISSING;
6393
6394 return Success;
6395 }
6396
6397 /* Copy character strings in read structure. */
copy_char_str(dest,source)6398 copy_char_str (dest, source)
6399 char **dest;
6400 char *source;
6401 {
6402 int malloc_len;
6403
6404 /* Free existing string. */
6405 if (*dest != NULL)
6406 free (*dest);
6407 *dest = NULL;
6408
6409 /* If source string exists, allocate corresponding dest string and copy
6410 * source. */
6411 if (source != NULL) {
6412 malloc_len = strlen (source) + 1;
6413 *dest = (char *) malloc (sizeof (char) * malloc_len);
6414 if (*dest == NULL)
6415 return Failure;
6416 strcpy (*dest, source);
6417 }
6418 return Success;
6419 }
6420
6421 /* see if variable data are packed. */
6422 /* Hoop, 2000/02/01 short-circuit to "not-packed" if sf&ao are default valued */
6423 int
are_data_packed(in_var,scale,offset,converted_type,input_type)6424 are_data_packed (in_var, scale, offset, converted_type, input_type)
6425 VAR_INFO *in_var; /* input variable structure */
6426 union att_val *scale; /* scale value, cast to type */
6427 union att_val *offset; /* offset value, cast to type */
6428 nc_type *converted_type; /* unpacked data type (returned) */
6429 nc_type *input_type; /* packed data type (returned) */
6430
6431 {
6432 int convert = 0, dfltsf_ao = 0 ;
6433 struct attrib_list *scale_att, *offset_att;
6434 union att_val junk;
6435 void set_data_to_type();
6436
6437 /* initialize types */
6438 *input_type = in_var->vartype;
6439 *converted_type = NC_UNSPECIFIED;
6440
6441 /* get scale and offset attributes */
6442 scale_att = find_att(in_var->first_vattr, SFNAME) ;
6443 if (!scale_att) {
6444 scale_att = find_att(in_var->first_vattr, SFNAME2) ;
6445 }
6446 offset_att = find_att(in_var->first_vattr, AONAME) ;
6447 if (!offset_att) {
6448 offset_att = find_att(in_var->first_vattr, AONAME2) ;
6449 }
6450
6451 if (scale_att == NULL && offset_att == NULL) {
6452 /* scale and offset do not exist; assume data are unpacked */
6453 *converted_type = *input_type;
6454 convert = 0;
6455 } else {
6456 if (scale_att != NULL) {
6457 *converted_type = scale_att->type;
6458 set_data_to_type(scale_att->type, 0, scale_att->data, scale);
6459 }
6460 if (offset_att != NULL) {
6461 *converted_type = offset_att->type;
6462 set_data_to_type(offset_att->type, 0, offset_att->data, offset);
6463 }
6464
6465 if (scale_att == NULL) {
6466 set_default_scale_offset(*converted_type, scale, &junk);
6467 }
6468 if (offset_att == NULL) {
6469 set_default_scale_offset(*converted_type, &junk, offset);
6470 }
6471 switch (*converted_type) {
6472 case NC_BYTE:
6473 if (scale->bval != B_SCALE || offset->bval != B_OFFSET) {
6474 convert = 1;
6475 } else {
6476 dfltsf_ao = 1 ;
6477 }
6478 break;
6479 case NC_SHORT:
6480 if (scale->sval != S_SCALE || offset->sval != S_OFFSET) {
6481 convert = 1;
6482 } else {
6483 dfltsf_ao = 1 ;
6484 }
6485 break;
6486 case NC_LONG:
6487 if (scale->lval != L_SCALE || offset->lval != L_OFFSET) {
6488 convert = 1;
6489 } else {
6490 dfltsf_ao = 1 ;
6491 }
6492 break;
6493 case NC_FLOAT:
6494 if (!flt_equal(scale->fval, F_SCALE) ||
6495 !flt_equal(offset->fval, F_OFFSET)) {
6496 convert = 1;
6497 } else {
6498 dfltsf_ao = 1 ;
6499 }
6500 break;
6501 case NC_DOUBLE:
6502 if (!dbl_equal(scale->dval, D_SCALE) ||
6503 !dbl_equal(offset->dval, D_OFFSET)) {
6504 convert = 1;
6505 } else {
6506 dfltsf_ao = 1 ;
6507 }
6508 break;
6509 default:
6510 break;
6511 } /* end switch */
6512
6513 /* also check to see if type conversion is necessary */
6514 if (*converted_type != *input_type) {
6515 if (dfltsf_ao == 1) {
6516 convert = 0 ;
6517 } else {
6518 convert = 1;
6519 }
6520 }
6521 } /* end else (either scale or offset or both exist in else ) */
6522
6523 return(convert) ;
6524 }
6525
6526 int
set_default_scale_offset(type,scale,offset)6527 set_default_scale_offset (type, scale, offset)
6528 nc_type type; /* scale and offset data type */
6529 union att_val *scale; /* pointer to scale data */
6530 union att_val *offset; /* pointer to offset data */
6531
6532 {
6533 switch (type)
6534 {
6535 case NC_BYTE:
6536 scale->bval = B_SCALE;
6537 offset->bval = B_OFFSET;
6538 break;
6539 case NC_SHORT:
6540 scale->sval = S_SCALE;
6541 offset->sval = S_OFFSET;
6542 break;
6543 case NC_LONG:
6544 scale->lval = L_SCALE;
6545 offset->lval = L_OFFSET;
6546 break;
6547 case NC_FLOAT:
6548 scale->fval = F_SCALE;
6549 offset->fval = F_OFFSET;
6550 break;
6551 case NC_DOUBLE:
6552 scale->dval = D_SCALE;
6553 offset->dval = D_OFFSET;
6554 break;
6555 default:
6556 break;
6557 }
6558 return Success;
6559 }
6560
6561 #ifdef TESTING
6562 void
xdumpreq(struct gagrid * gridptr)6563 xdumpreq(struct gagrid *gridptr) {
6564 void xdumpfile(struct gafile *fileptr) ;
6565 void xdumpvar(struct gavar *varptr) ;
6566
6567 fprintf(stderr, "grid request address=0x%X;\n", gridptr) ;
6568 fprintf(stderr, "grid=0x%X; undef=%E; rmin=%E; rmax=%E;\n",
6569 gridptr->grid, gridptr->undef, gridptr->rmin, gridptr->rmax) ;
6570 fprintf(stderr, "isiz=%d; jsiz=%d; idim=%d; jdim=%d;\n",
6571 gridptr->isiz, gridptr->jsiz, gridptr->idim, gridptr->jdim) ;
6572 fprintf(stderr, "dimmin: %d; %d; %d; %d;\n", gridptr->dimmin[0],
6573 gridptr->dimmin[1], gridptr->dimmin[2], gridptr->dimmin[3]) ;
6574 fprintf(stderr, "dimmax: %d; %d; %d; %d;\n", gridptr->dimmax[0],
6575 gridptr->dimmax[1], gridptr->dimmax[2], gridptr->dimmax[3]) ;
6576 if (gridptr->exprsn) {
6577 fprintf(stderr, "exprsn=%s;\n", gridptr->exprsn) ;
6578 } else {
6579 fprintf(stderr, "exprsn is null;\n") ;
6580 }
6581 fprintf(stderr, "alocf=%d;\n", gridptr->alocf) ;
6582 fprintf(stderr, "igrab=0x%X; jgrab=0x%X; iabgr=0x%X; jabgr=0x%X;\n",
6583 gridptr->igrab, gridptr->jgrab, gridptr->iabgr, gridptr->jabgr) ;
6584 fprintf(stderr, "ivals=0x%X; jvals=0x%X; iavals=0x%X; javals=0x%X;\n",
6585 gridptr->ivals, gridptr->jvals, gridptr->iavals, gridptr->javals) ;
6586 fprintf(stderr, "ilinr=%d; jlinr=%d;\n", gridptr->ilinr, gridptr->jlinr) ;
6587 fprintf(stderr, "Calling xdumpfile from xdumpreq.\n") ;
6588 xdumpfile(gridptr->pfile) ;
6589 fprintf(stderr, "Calling xdumpvar from xdumpreq.\n") ;
6590 xdumpvar(gridptr->pvar) ;
6591 }
6592
6593 void
xdumpfile(struct gafile * fileptr)6594 xdumpfile(struct gafile *fileptr) {
6595 int varinx ;
6596
6597 fprintf(stderr, "======== start of gafile structure info dump ============\n");
6598 fprintf(stderr, "file address=0x%X;\n", fileptr) ;
6599 fprintf(stderr, "pforw=0x%X; infile=0x%X; type=%d; undef=%E;\n",
6600 fileptr->pforw, fileptr->infile, fileptr->type, fileptr->undef) ;
6601 if (fileptr->name) {
6602 fprintf(stderr, "name=%s;\n", fileptr->name) ;
6603 } else {
6604 fprintf(stderr, "name is null.\n") ;
6605 }
6606 if (fileptr->dnam) {
6607 fprintf(stderr, "dnam=%s;\n", fileptr->dnam) ;
6608 } else {
6609 fprintf(stderr, "dnam is null.\n") ;
6610 }
6611 if (fileptr->mnam) {
6612 fprintf(stderr, "mnam=%s;\n", fileptr->mnam) ;
6613 } else {
6614 fprintf(stderr, "mnam is null.\n") ;
6615 }
6616 if (fileptr->title) {
6617 fprintf(stderr, "title=%s;\n", fileptr->title) ;
6618 } else {
6619 fprintf(stderr, "title is null.\n") ;
6620 }
6621 fprintf(stderr, "ulow=%E; uhi=%E;\n", fileptr->ulow, fileptr->uhi) ;
6622 fprintf(stderr, "rbuf=0x%X; pbuf=0x%X; bbuf=0x%X;\n", fileptr->rbuf,
6623 fileptr->pbuf, fileptr->bbuf) ;
6624 fprintf(stderr, "bswap=%d; mtype=%d; tstrt=0x%X; tcnt=0x%X;\n",
6625 fileptr->bswap, fileptr->mtype, fileptr->tstrt, fileptr->tcnt) ;
6626 fprintf(stderr, "stcnt=%d; stpos=%d; mfile=0x%X;\n", fileptr->stcnt,
6627 fileptr->stpos, fileptr->mfile) ;
6628 fprintf(stderr, "dnum: %d; %d; %d; %d;\n", fileptr->dnum[0],
6629 fileptr->dnum[1], fileptr->dnum[2], fileptr->dnum[3]) ;
6630 fprintf(stderr, "tlpflg=%d; tlpst=%d; vnum=%d; ivnum=%d; lvnum=%d;\n",
6631 fileptr->tlpflg, fileptr->vnum, fileptr->ivnum, fileptr->lvnum) ;
6632 for (varinx = 0 ; varinx < fileptr->vnum ; ++varinx) {
6633 fprintf(stderr, "var #%d=0x%X;\n", varinx, &(fileptr->pvar1[varinx])) ;
6634 }
6635 fprintf(stderr, "gsiz=%d; tsiz=%d; trecs=%d; fhdr=%d;\n",
6636 fileptr->gsiz, fileptr->tsiz, fileptr->trecs, fileptr->fhdr) ;
6637 fprintf(stderr, "wrap=%d; seqflg=%d; yrflg=%d; zrflg=%d;\n",
6638 fileptr->wrap, fileptr->seqflg, fileptr->yrflg, fileptr->zrflg) ;
6639 fprintf(stderr, "ppflag=%d; ppwrot=%d; ppisiz=%d; ppjsiz=%d;\n",
6640 fileptr->ppflag, fileptr->ppwrot, fileptr->ppisiz, fileptr->ppjsiz) ;
6641 fprintf(stderr, "ppi=0x%X; ppf[0]=0x%X; ppf[1]=0x%X;\n", fileptr->ppi,
6642 fileptr->ppf[0], fileptr->ppf[1]) ;
6643 fprintf(stderr, "gr2ab[0]=0x%X; ab2gr[0]=0x%X; grvals[0]=0x%X; abvals[0]=0x%X;\n",
6644 fileptr->gr2ab[0], fileptr->ab2gr[0], fileptr->grvals[0], fileptr->abvals[0]) ;
6645 fprintf(stderr, "linear: %d; %d; %d; %d;\n", fileptr->linear[0],
6646 fileptr->linear[1], fileptr->linear[2], fileptr->linear[3]) ;
6647 fprintf(stderr, "dimoff: %d; %d; %d; %d;\n", fileptr->dimoff[0],
6648 fileptr->dimoff[1], fileptr->dimoff[2], fileptr->dimoff[3]) ;
6649 fprintf(stderr, "climo=%d; cysiz=%d; idxflg=%d, grbgrd=%d;\n",
6650 fileptr->climo, fileptr->cysiz, fileptr->idxflg, fileptr->grbgrd) ;
6651 fprintf(stderr, "pindx=0x%X; tmplat=%d; fnums=0x%X; fnumc=%d;\n",
6652 fileptr->pindx, fileptr->tmplat, fileptr->fnums, fileptr->fnumc) ;
6653 fprintf(stderr, "errcnt=%d; errflg=%d;\n", fileptr->errcnt, fileptr->errflg) ;
6654 fprintf(stderr, "======== end of gafile structure info dump ============\n");
6655 }
6656
6657 void
xdumpvar(struct gavar * varptr)6658 xdumpvar(struct gavar *varptr) {
6659
6660 fprintf(stderr, "var. address=0x%X;\n", varptr) ;
6661 if (varptr->varnm) {
6662 fprintf(stderr, "varnm=%s;\n", varptr->varnm) ;
6663 } else {
6664 fprintf(stderr, "varnm is null.\n") ;
6665 }
6666 if (varptr->abbrv) {
6667 fprintf(stderr, "abbrv=%s;\n", varptr->abbrv) ;
6668 } else {
6669 fprintf(stderr, "abbrv is null.\n") ;
6670 }
6671 fprintf(stderr, "units: %d; %d; %d; %d;\n", varptr->units[0],
6672 varptr->units[1], varptr->units[2], varptr->units[3]) ;
6673 fprintf(stderr, "offset=%d; recoff=%d; levels=%d;\n", varptr->offset,
6674 varptr->recoff, varptr->levels) ;
6675 }
6676 #endif
6677
6678 /* function to attempt a linear definition of a dimension. Returns 1 */
6679 /* if successful, 0 if not linearizable, and <0 if an error occurs */
6680
6681 int
trydeflin(struct gafile * pfi,VAR_INFO * coord,int GrADSdimnum,int isX,int revflag)6682 trydeflin(struct gafile *pfi, VAR_INFO *coord, int GrADSdimnum, int isX,
6683 int revflag) {
6684 float *deltas, coordvals[2], *vals, firstcoordval ;
6685 int islin, deltaidx ;
6686 int fequal(float op1, float op2, float tolerance) ;
6687 int convtype(void *srcptr, nc_type srctype, void *desptr, nc_type destype,
6688 int srcindex, int desindex) ;
6689
6690 if (pfi->dnum[GrADSdimnum] < 2) {
6691 /* WRONG! return(0) ; *//* Use levels */ /* gaddes code doesn't do that */
6692 if (!convtype(coord->data, coord->vartype, (void *) &(coordvals[0]),
6693 NC_FLOAT, 0, 0)) {
6694 gaprnt(0, "Couldn't convert coordinate value to float.\n") ;
6695 return(-1) ;
6696 }
6697 firstcoordval = coordvals[0] ;
6698 /* Ok, now we emulate deflin code */
6699 if (!(vals = (float *) malloc(sizeof(float) * 6))) {
6700 gaprnt(0, "Memory allocation error eval. lin. coord. for SDF.\n") ;
6701 return(-1) ;
6702 }
6703 vals[0] = 1.0 ;
6704 vals[1] = 0.0 ;
6705 vals[2] = -999.9 ;
6706 vals[3] = 1.0 ;
6707 vals[4] = 0.0 ;
6708 vals[5] = -999.9 ;
6709 pfi->grvals[GrADSdimnum] = vals ;
6710 pfi->abvals[GrADSdimnum] = vals + 3 ;
6711 pfi->ab2gr[GrADSdimnum] = liconv ;
6712 pfi->gr2ab[GrADSdimnum] = liconv ;
6713 pfi->linear[GrADSdimnum] = 1 ;
6714 if (isX) /* a coordinate of extent one CANNOT wrap the globe */ {
6715 pfi->wrap = 0 ;
6716 }
6717 return(1) ;
6718 }
6719 deltas = (float *) malloc(sizeof(float)*(pfi->dnum[GrADSdimnum] - 1)) ;
6720 if (!deltas) {
6721 gaprnt(0, "Memory allocation error in trydeflin.\n") ;
6722 return(-1) ;
6723 }
6724 if (!revflag) {
6725 if (!convtype(coord->data, coord->vartype, (void *) &(coordvals[0]),
6726 NC_FLOAT, 0, 0)) {
6727 gaprnt(0, "Couldn't convert coordinate value to float.\n") ;
6728 return(-1) ;
6729 }
6730 } else {
6731 if (!convtype(coord->data, coord->vartype, (void *) &(coordvals[0]),
6732 NC_FLOAT, pfi->dnum[GrADSdimnum] - 1, 0)) {
6733 gaprnt(0, "Couldn't convert coordinate value to float.\n") ;
6734 return(-1) ;
6735 }
6736 }
6737 firstcoordval = coordvals[0] ;
6738 for (deltaidx = 1 ; deltaidx < pfi->dnum[GrADSdimnum] ; ++deltaidx) {
6739 if (!revflag) {
6740 if (!convtype(coord->data, coord->vartype, (void *) &(coordvals[0]),
6741 NC_FLOAT, deltaidx, 1)) {
6742 gaprnt(0, "Couldn't convert coordinate value to float.\n") ;
6743 return(-1) ;
6744 }
6745 } else {
6746 if (!convtype(coord->data, coord->vartype, (void *) &(coordvals[0]),
6747 NC_FLOAT, (pfi->dnum[GrADSdimnum] - 1) - deltaidx,
6748 1)) {
6749 gaprnt(0, "Couldn't convert coordinate value to float.\n") ;
6750 return(-1) ;
6751 }
6752 }
6753 deltas[deltaidx - 1] = coordvals[1] - coordvals[0] ;
6754 if (isX) {
6755 if (coordvals[1] < coordvals[0]) {
6756 deltas[deltaidx - 1] = (360.0 + coordvals[1]) - coordvals[0] ;
6757 }
6758 } /* longitudes can wrap around Greenwich and yet be linear */
6759 coordvals[0] = coordvals[1] ;
6760 } /* end for */
6761 islin = 1 ; /* assume coordinate is linear for now */
6762 for (deltaidx = 1 ; deltaidx < (pfi->dnum[GrADSdimnum] - 1) ; ++deltaidx) {
6763 if (!fequal(deltas[deltaidx - 1], deltas[deltaidx],
6764 ((float)fabs(deltas[deltaidx])) / 1000.0)) {
6765 islin = 0 ;
6766 break ;
6767 }
6768 } /* end for */
6769 if (!islin) {
6770 return(0) ; /* use levels */
6771 }
6772 /* Ok, now we emulate deflin code */
6773 if (!(vals = (float *) malloc(sizeof(float) * 6))) {
6774 gaprnt(0, "Memory allocation error eval. lin. coord. for SDF.\n") ;
6775 return(-1) ;
6776 }
6777 vals[0] = deltas[0] ;
6778 vals[1] = firstcoordval - deltas[0] ;
6779 vals[2] = -999.9 ;
6780 vals[3] = 1.0 / deltas[0] ;
6781 vals[4] = -1.0 * ((firstcoordval - deltas[0]) / deltas[0]) ;
6782 vals[5] = -999.9 ;
6783 pfi->grvals[GrADSdimnum] = vals ;
6784 pfi->abvals[GrADSdimnum] = vals + 3 ;
6785 pfi->ab2gr[GrADSdimnum] = liconv ;
6786 pfi->gr2ab[GrADSdimnum] = liconv ;
6787 pfi->linear[GrADSdimnum] = 1 ;
6788 if (isX) /* does it wrap the globe? */ {
6789 /* fabs & fmod take and return doubles, hence the casts */
6790 if (fequal((float) fabs(fmod((double)(coordvals[1] + deltas[0] - firstcoordval ),
6791 (double) 360.0)),
6792 (double) 0.0, ((float) fabs((double) deltas[0])) / 2.0)) {
6793 pfi->wrap = 1 ;
6794 } else {
6795 pfi->wrap = 0 ;
6796 }
6797 }
6798 return(1) ;
6799 }
6800
6801 int
fequal(float op1,float op2,float tolerance)6802 fequal(float op1, float op2, float tolerance) {
6803 /*mf ---- added casts to be consistent with ANSI --- */
6804 if (fabs((double)op1 - (double)op2) <= (double)tolerance) {
6805 return(1) ;
6806 } else {
6807 return(0) ;
6808 }
6809 }
6810
6811 int
dequal(double op1,double op2,double tolerance)6812 dequal(double op1, double op2, double tolerance) {
6813 if (fabs(op1 - op2) <= tolerance) {
6814 return(1) ;
6815 } else {
6816 return(0) ;
6817 }
6818 }
6819
6820 #endif
6821
6822 #endif
6823